!#define _DEBUG_
module hybrid_ensemble_isotropic
!$$$   module documentation block
!                .      .    .                                       .
! module:    hybrid_ensemble_isotropic
!   prgmmr: parrish          org: np22                date: 2009-09-28
!
! abstract: contains routines for localization of the hybrid ensemble
!            control variable a_en.  this application is for 
!            localization with a regional model with homogeneous scale in horizontal
!            and vertical scale a function of vertical only.
!
! program history log:
!   2009-09-28  parrish, initial documentation.
!   2010-02-26  parrish, remove redundant special spectral and sub2grid/grid2sub code.  replaced
!                 by general purpose code as part of adding dual resolution option.
!   2010-03-17  zhu  - use vlevs from gridmod
!   2010-04-06  parrish - fix dimension error in ensemble_forward_model_ad_dual_res and
!                           add array deallocation in ensemble_forward_model_ad_dual_res and
!                            ensemble_forward_model_dual_res
!   2010-05-14  parrish - repair code to get dual resolution option working again.  remove
!                           error stops when dual resolution invoked.
!   2010-05-20  todling - renamed all cstate to bundle to avoid confusion; the
!                         bundles here are not necessarily idendical to the control vector,
!                         but rather what this module take the ensemble to be composed of
!   2010-07-29  kleist  - combine hybrid_ensemble_isotropic global/regional modules 
!   2011-02-28  parrish - introduce more complete use of gsi_bundlemod to eliminate hard-wired variables
!   2011-06-28  parrish - add code to allow sqrt of localization correlation so hybrid ensemble option
!                            can be extended to sqrt(B) minimization option.
!   2011-08-31  todling - revisit en_perts (single-prec) in light of extended bundle
!   2011-10-26  kleist - add time dimension to en_perts for 4d extensions
!   2011-11-01  kleist - 4d-ensemble-var available as part of ensemble forward models(s) and adjoints,
!                        as described in Liu et al. (2008), MWR, vol.36 and Buehner et al. (2010), 
!                        MWR, vol.138.  This specific implementation uses the linearized observation 
!                        operators as in Buehner et al., but also has the capability to allow for a 
!                        contribution from a static, time-invariant (3DVAR) B.  Capability also exists
!                        for actual hybrid Ens-4DVAR (with TL/AD).
!   2012-01-17  wu      - option "pwgtflg": psfc with vertically integrated contribution 
!                          in forward and adjoint routines
!                       - option "betaflg": height dependent beta in regional
!   2012-02-08  parrish - add changes to allow regional dual res 
!   2012-02-08  parrish - cleanup
!   2012-10-11  wu      - dual resolution for regional hybens options; 
!                         use ensemble dimensions on control variable: alpha
!   2013-04-17  wu      - bug fix in normalizing the recursive filter
!
! subroutines included:
!   sub init_rf_z                         - initialize localization recursive filter (z direction)
!   sub init_rf_x                         - initialize localization recursive filter (x direction)
!   sub init_rf_y                         - initialize localization recursive filter (y direction)
!   sub new_factorization_rf_z            - localization recursive filter (z direction)
!   sub new_factorization_rf_x            - localization recursive filter (x direction)
!   sub new_factorization_rf_y            - localization recursive filter (y direction)
!   sub normal_new_factorization_rf_z     - normalize localization recursive filter (z direction)
!   sub normal_new_factorization_rf_x     - normalize localization recursive filter (x direction)
!   sub normal_new_factorization_rf_y     - normalize localization recursive filter (y direction)
!   sub create_ensemble                   - allocate space for ensemble perturbations
!   sub load_ensemble                     - read/generate ensemble perturbations
!   sub rescale_ensemble_rh_perturbations - for internal generated perturbations only, normalize by ges q
!   sub ensemble_forward_model            - add ensemble contribution to analysis increment
!   sub ensemble_forward_model_ad         - adjoint of ensemble_forward_model
!   sub ensemble_forward_model_dual_res   - dual resolution version of ensemble_forward_model
!   sub ensemble_forward_model_ad_dual_res- adjoint of ensemble_forward_model_dual_res
!
! Variable Definitions:
!   def yyyy      - what yyyy is
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind,r_single
  use mpimod, only: mype
  use gsi_bundlemod, only: gsi_bundlecreate
  use gsi_bundlemod, only: gsi_bundleset
  use gsi_bundlemod, only: gsi_bundleunset
  use gsi_bundlemod, only: gsi_grid
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use gsi_bundlemod, only: gsi_bundledestroy
  use gsi_bundlemod, only: gsi_gridcreate

  use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d
  use hybrid_ensemble_parameters, only: nval_lenz_en
  implicit none

! set default to private
  private
! set subroutines to public
  public :: init_rf_z
  public :: init_rf_x
  public :: init_rf_y
  public :: new_factorization_rf_z
  public :: new_factorization_rf_x
  public :: new_factorization_rf_y
  public :: normal_new_factorization_rf_z
  public :: normal_new_factorization_rf_x
  public :: normal_new_factorization_rf_y
  public :: create_ensemble
  public :: load_ensemble
  public :: rescale_ensemble_rh_perturbations
  public :: ensemble_forward_model
  public :: ensemble_forward_model_dual_res
  public :: ensemble_forward_model_ad
  public :: ensemble_forward_model_ad_dual_res
  public :: beta12mult
  public :: init_sf_xy
  public :: sf_xy
  public :: sqrt_sf_xy
  public :: sqrt_sf_xy_ad

! set passed variables to public
  public :: nelen
  public :: en_perts,ps_bar
  public :: region_lat_ens,region_lon_ens
  public :: region_dx_ens,region_dy_ens

  character(len=*),parameter::myname='hybrid_ensemble_isotropic'
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!      following variables are filter parameters for isotropic
!       homogeneous localization of hybrid control variable a_en

  real(r_kind),allocatable:: fmatz(:,:,:,:)
  real(r_kind),allocatable:: fmat0z(:,:,:)
  real(r_kind),allocatable:: fmatx(:,:,:,:)
  real(r_kind),allocatable:: fmat0x(:,:,:)
  real(r_kind),allocatable:: fmaty(:,:,:)
  real(r_kind),allocatable:: fmat0y(:,:)
  real(r_kind),allocatable:: znorm_new(:,:)
  real(r_kind),allocatable:: xnorm_new(:,:)
  real(r_kind),allocatable:: ynorm_new(:)
  real(r_kind),allocatable:: psbar(:)

! Other local variables for horizontal/spectral localization
  real(r_kind),allocatable,dimension(:,:)  :: spectral_filter,sqrt_spectral_filter
  integer(i_kind),allocatable,dimension(:) :: k_index

! following is for storage of ensemble perturbations:

!   def en_perts            - array of ensemble perturbations
!   def nelen               - length of one ensemble perturbation vector

  integer(i_kind) nelen
  type(gsi_bundle),save,allocatable :: en_perts(:,:)
  real(r_single),dimension(:,:,:),allocatable:: ps_bar

!    following is for interpolation of global ensemble to regional ensemble grid

  real(r_kind),allocatable:: region_lat_ens(:,:),region_lon_ens(:,:)
  real(r_kind),allocatable:: region_dx_ens(:,:),region_dy_ens(:,:)

!    following is for special subdomain to slab variables used when internally generating ensemble members

  integer(i_kind) nval2f,nscl
  integer(i_kind) nh_0,nh_1,nv_0,nv_1
  integer(i_kind),allocatable,dimension(:):: nsend_sd2h,ndsend_sd2h,nrecv_sd2h,ndrecv_sd2h
  integer(i_kind),allocatable,dimension(:):: i_recv,k_recv

  logical,parameter:: debug=.false.


contains

subroutine init_rf_z(z_len)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_rf_z    initialize vertical recursive filter
!   prgmmr: parrish          org: np22                date: 2009-12-16
!
! abstract: initialize vertical recursive filter for hybrid ensemble control variable a_en
!            call this one first, then init_rf_x, init_rf_y
!
! program history log:
!   2009-12-16  parrish
!   2010-04-06  parrish - add 2nd option for units of vertical localization:
!                             if z_len < 0, then abs(z_len) is vertical localization length scale in
!                             units of ln(p).  otherwise, when z_len > 0, localization is in vertical
!                             grid units.  The ln(p) distance measurement is approximate, based on a
!                             fixed surface pressure of 1000mb.  This is because at the point where this
!                             is currently called, the background 3d pressure field is not yet available.
!                             A later version will correct this.
!                             For the current s_ens_v > 0, the measure is vertical grid units.  
!                             s_ens_v = 20 and s_ens_v = -0.44 are roughly comparable, and
!                             connection of .44 is .44 = (sqrt(.15)/sqrt(2))*1.6, where 1.6 is the value used
!                             by Jeff Whitaker for his distance in which the Gaspari-Cohn function 1st = 0.
!   2011-07-19  tong    - add the calculation of pressure vertical profile for regional model,
!                             when vertical localization length scale is in units of ln(p)
!
!   input argument list:
!     z_len    - filter length scale in grid units
!
!   output argument list:
!
! remarks:  see modules used
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use gridmod, only: nsig,ak5,bk5,aeta1_ll,aeta2_ll,pt_ll,pdtop_ll,twodvar_regional, &
                     wrf_nmm_regional,nems_nmmb_regional,wrf_mass_regional,cmaq_regional, &
                     regional
  use constants, only: half,one,rd_over_cp,zero,one_tenth,ten,two
  use hybrid_ensemble_parameters, only: grd_ens,s_ens_v

  real(r_kind)   ,intent(in) :: z_len(grd_ens%nsig)

  integer(i_kind) k,km,kp,nxy,i,ii,jj
  real(r_kind) aspect(nsig),p_interface(nsig+1),ln_p_int(nsig+1)
  real(r_kind) p_layer,dlnp,kap1,kapr,d1,d2,rnsig

  kap1=rd_over_cp+one
  kapr=one/rd_over_cp
  nxy=grd_ens%latlon11
  rnsig=float(nsig)

!    use new factorization:
  allocate(fmatz(2,nsig,2,nxy),fmat0z(nsig,2,nxy))

!   for z_len < zero, use abs val z_len and assume localization scale is in units of ln(p)
  if(s_ens_v > zero) then

!  z_len is in grid units
     do k=1,nsig
        aspect(k)=z_len(k)**2
     end do

     do i=1,nxy
       call get_new_alpha_beta(aspect,nsig,fmatz(1,1,1,i),fmat0z(1,1,i))
     end do

  else

!  abs(z_len) is in units of ln(p)
!              put in approximate vertical scale which depends on ln(p)
!              use ensemble mean ps (ps_bar) to construct pressure
!              profiles using GFS sigma-p coordinate [ add mods to be able
!              to use fully generalized coordinate later** ]
!
     i=0
     do jj=1,grd_ens%lon2
        do ii=1,grd_ens%lat2
           i=i+1

           if(regional)then

              do k=1,nsig+1
                 if (wrf_nmm_regional.or.nems_nmmb_regional.or.cmaq_regional) then
   	            p_interface(k)= one_tenth* &
                          (aeta1_ll(k)*pdtop_ll + &
                           aeta2_ll(k)*(ten*ps_bar(ii,jj,1)-pdtop_ll-pt_ll) + &
                           pt_ll)
                 endif
                 if (wrf_mass_regional .or. twodvar_regional) then
		    p_interface(k)=one_tenth*(aeta1_ll(k)*(ten*ps_bar(ii,jj,1)-pt_ll)+pt_ll)
                 endif
	         ln_p_int(k)=log(max(p_interface(k),0.0001_r_kind))
              end do

           else
              do k=1,nsig+1
                 p_interface(k)=ak5(k)+(bk5(k)*ps_bar(ii,jj,1))
                 ln_p_int(k)=log(max(p_interface(k),0.0001_r_kind))
              end do
           endif

           do k=1,nsig
              dlnp=abs(ln_p_int(k)-ln_p_int(k+1))
 	      d1=abs(z_len(k))/dlnp
              d1=min(rnsig,d1)
              aspect(k)=d1**2
!!            if(mype == 0) write(400,'(" k, vertical localization in grid units for ln(p) scaling =",i4,f10.2,f10.2,f10.2)') &
!!                                         k,sqrt(aspect(k))
           enddo
   
           call get_new_alpha_beta(aspect,nsig,fmatz(1,1,1,i),fmat0z(1,1,i))
        end do
     end do
  end if

end subroutine init_rf_z

subroutine init_rf_x(x_len)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_rf_x    initialize x direction recursive filter
!   prgmmr: parrish          org: np22                date: 2009-12-16
!
! abstract: initialize longitude recursive filters for hybrid ensemble control variable a_en
!
! program history log:
!   2009-12-16  parrish
!
!   input argument list:
!     x_len -- filter length scale in grid units
!
!   output argument list:
!
! remarks:  see modules used
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use hybrid_ensemble_parameters, only: grd_ens
  use constants, only: half

  real(r_kind),intent(in   ) :: x_len

  integer(i_kind) i,j,k,l
  real(r_kind) aspect(grd_ens%nlon)
  real(r_kind) fmatc(2,grd_ens%nlon,2),fmat0c(grd_ens%nlon,2)

!    use new factorization:
  if(allocated(fmatx)) deallocate(fmatx)
  if(allocated(fmat0x)) deallocate(fmat0x)
  allocate(fmatx(grd_ens%nlat,2,grd_ens%nlon,2),fmat0x(grd_ens%nlat,grd_ens%nlon,2))
  do i=1,grd_ens%nlat
     do j=1,grd_ens%nlon
        aspect(j)=(x_len*region_dy_ens(grd_ens%nlat/2,grd_ens%nlon/2)/region_dx_ens(i,j))**2 ! only works for rotated lat-lon grids
     end do
     call get_new_alpha_beta(aspect,grd_ens%nlon,fmatc,fmat0c)
     do k=1,2
        do j=1,grd_ens%nlon
           do l=1,2
              fmatx(i,l,j,k)=fmatc(l,j,k)
           end do
           fmat0x(i,j,k)=fmat0c(j,k)
        end do
     end do
  end do

end subroutine init_rf_x

subroutine init_rf_y(y_len)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_rf_y    initialize y direction recursive filter
!   prgmmr: parrish          org: np22                date: 2009-12-16
!
! abstract: initialize latitude recursive filters for hybrid ensemble control variable a_en
!
! program history log:
!   2009-12-16  parrish
!
!   input argument list:
!     y_len -- filter length scale in grid units
!
!   output argument list:
!
! remarks:  see modules used
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use hybrid_ensemble_parameters, only: grd_ens
  use constants, only: half

  real(r_kind),intent(in   ) :: y_len

  real(r_kind) aspect(grd_ens%nlat)
  integer(i_kind) i

!    use new factorization:
  if(allocated(fmaty)) deallocate(fmaty)
  if(allocated(fmat0y)) deallocate(fmat0y)
  allocate(fmaty(2,grd_ens%nlat,2),fmat0y(grd_ens%nlat,2))
  do i=1,grd_ens%nlat
     aspect(i)=y_len**2
  end do
  call get_new_alpha_beta(aspect,grd_ens%nlat,fmaty,fmat0y)

end subroutine init_rf_y

subroutine new_factorization_rf_z(f,iadvance,iback)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    new_factorization_rf_z
!
!   prgrmmr:  parrish        org: np22                              date: 2009-09-28
!
! abstract:  apply new factorization Purser 1-d high-order filter in z (vertical) dimension.
!
! program history log:
!   2009-09-28  parrish  initial documentation
!   2010-02-20  parrish  modifications for dual resolution
!   2010-04-14  kleist   modifications for ln(ps) localization option
!
!   input argument list:
!     f        - input field to be filtered
!     iadvance - =1  for forward operator, =2 for adjoint operator
!     iback    - =2  for forward operator, =1 for adjoint operator
!
!   output argument list:
!     f        - filtered output
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

  use hybrid_ensemble_parameters, only: grd_ens
  implicit none

  integer(i_kind),intent(in   ) :: iadvance,iback
  real(r_kind)   ,intent(inout) :: f(grd_ens%latlon11,grd_ens%nsig)

  integer(i_kind) i,j,k,l,nxy,nz

  nxy=grd_ens%latlon11 ; nz=grd_ens%nsig
  if(iadvance == 1) then
     do k=1,nz
        do i=1,nxy
           f(i,k)=znorm_new(i,k)*f(i,k)
        end do
     end do
  end if
  do k=1,nz
     do l=1,min(2,k-1)
        do i=1,nxy
           f(i,k)=f(i,k)-fmatz(l,k,iadvance,i)*f(i,k-l)
        end do
     end do
     do i=1,nxy
        f(i,k)=fmat0z(k,iadvance,i)*f(i,k)
     end do
  end do
  do k=nz,1,-1
     do l=1,min(2,nz-k)
        do i=1,nxy
           f(i,k)=f(i,k)-fmatz(l,k+l,iback,i)*f(i,k+l)
        end do
     end do
     do i=1,nxy
        f(i,k)=fmat0z(k,iback,i)*f(i,k)
     end do
  end do
  if(iadvance == 2) then
     do k=1,nz
        do i=1,nxy
           f(i,k)=znorm_new(i,k)*f(i,k)
        end do
     end do
  end if

end subroutine new_factorization_rf_z

subroutine new_factorization_rf_x(f,iadvance,iback,nlevs)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    new_factorization_rf_x
!
!   prgrmmr:  parrish        org: np22                              date: 2009-09-28
!
! abstract:  apply new factorization Purser 1-d high-order filter in x (longitude) direction.
!
! program history log:
!   2009-09-28  parrish  initial documentation
!   2010-02-20  parrish  modifications for dual resolution
!   2010-03-11  parrish  adjust for possibility that nlevs=0
!
!   input argument list:
!     f        - input field to be filtered
!     iadvance - =1  for forward operator, =2 for adjoint operator
!     iback    - =2  for forward operator, =1 for adjoint operator
!
!   output argument list:
!     f        - filtered output
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block
  use hybrid_ensemble_parameters, only: grd_ens
  implicit none

  integer(i_kind),intent(in   ) :: iadvance,iback,nlevs
  real(r_kind)   ,intent(inout) :: f(grd_ens%nlat,grd_ens%nlon,max(nlevs,1))

  integer(i_kind) i,j,k,l,ny,nx,nz

  ny=grd_ens%nlat ; nx=grd_ens%nlon ; nz=nlevs
  do k=1,nz

     if(iadvance == 1) then
        do j=1,nx
           do i=1,ny
              f(i,j,k)=xnorm_new(i,j)*f(i,j,k)
           end do
        end do
     end if

     do j=1,nx
        do l=1,min(2,j-1)
           do i=1,ny
              f(i,j,k)=f(i,j,k)-fmatx(i,l,j,iadvance)*f(i,j-l,k)
           end do
        end do
        do i=1,ny
           f(i,j,k)=fmat0x(i,j,iadvance)*f(i,j,k)
        end do
     end do

     do j=nx,1,-1
        do l=1,min(2,nx-j)
           do i=1,ny
              f(i,j,k)=f(i,j,k)-fmatx(i,l,j+l,iback)*f(i,j+l,k)
           end do
        end do
        do i=1,ny
           f(i,j,k)=fmat0x(i,j,iback)*f(i,j,k)
        end do
     end do

     if(iadvance == 2) then
        do j=1,nx
           do i=1,ny
              f(i,j,k)=xnorm_new(i,j)*f(i,j,k)
           end do
        end do
     end if

  end do

end subroutine new_factorization_rf_x

subroutine new_factorization_rf_y(f,iadvance,iback,nlevs)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    new_factorization_rf_y
!
!   prgrmmr:  parrish        org: np22                              date: 2009-09-28
!
! abstract:  apply new factorization Purser 1-d high-order filter in y (latitude) direction.
!
! program history log:
!   2009-09-28  parrish  initial documentation
!   2010-02-20  parrish  modifications for dual resolution
!   2010-03-11  parrish  adjust for possibility that nlevs=0
!
!   input argument list:
!     f        - input field to be filtered
!     iadvance - =1  for forward operator, =2 for adjoint operator
!     iback    - =2  for forward operator, =1 for adjoint operator
!
!   output argument list:
!     f        - filtered output
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block
  use hybrid_ensemble_parameters, only: grd_ens
                      !                  use mpimod, only: mype
  implicit none

  integer(i_kind),intent(in   ) :: iadvance,iback,nlevs
  real(r_kind)   ,intent(inout) :: f(grd_ens%nlat,grd_ens%nlon*max(nlevs,1))

  integer(i_kind) i,k,l,nx,ny,nz

  nx=grd_ens%nlon ; ny=grd_ens%nlat ; nz=nlevs
  do k=1,nx*nz

     if(iadvance == 1) then
        do i=1,ny
           f(i,k)=ynorm_new(i)*f(i,k)
        end do
     end if

     do i=1,ny
        do l=1,min(2,i-1)
           f(i,k)=f(i,k)-fmaty(l,i,iadvance)*f(i-l,k)
        end do
        f(i,k)=fmat0y(i,iadvance)*f(i,k)
     end do

     do i=ny,1,-1
        do l=1,min(2,ny-i)
           f(i,k)=f(i,k)-fmaty(l,i+l,iback)*f(i+l,k)
        end do
        f(i,k)=fmat0y(i,iback)*f(i,k)
     end do

     if(iadvance == 2) then
        do i=1,ny
           f(i,k)=ynorm_new(i)*f(i,k)
        end do
     end if

  end do

end subroutine new_factorization_rf_y

subroutine normal_new_factorization_rf_z
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    normal_new_factorization_rf_z get normalization factor in z direction
!
!   prgrmmr:  parrish        org: np22                              date: 2009-09-28
!
! abstract:  compute normalization factor in z (vertical dimension)
!
! program history log:
!   2009-09-28  parrish  initial documentation
!   2010-02-20  parrish  modifications for dual resolution
!   2010-04-14  kleist   modifications for ln(p) localization
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  use hybrid_ensemble_parameters, only: grd_ens
  use constants, only: zero,one
  implicit none

  integer(i_kind) k,iadvance,iback,nxy
  real(r_kind) f(grd_ens%latlon11,grd_ens%nsig),diag(grd_ens%latlon11,grd_ens%nsig)

  if(allocated(znorm_new)) deallocate(znorm_new)
  allocate(znorm_new(grd_ens%latlon11,grd_ens%nsig))

  nxy=grd_ens%latlon11

  znorm_new=one

   do k=1,grd_ens%nsig
      f=zero
      f(:,k)=one

      iadvance=1 ; iback=2
      call new_factorization_rf_z(f,iadvance,iback)
      iadvance=2 ; iback=1
      call new_factorization_rf_z(f,iadvance,iback)

      diag(:,k)=sqrt(one/f(:,k))
   end do

   do k=1,grd_ens%nsig
      znorm_new(:,k)=diag(:,k)
   end do

! Check result:
  do k=1,grd_ens%nsig
     f=zero
     f(:,k)=one

     iadvance=1 ; iback=2
     call new_factorization_rf_z(f,iadvance,iback)
     iadvance=2 ; iback=1
     call new_factorization_rf_z(f,iadvance,iback)

     diag(:,k)=sqrt(one/f(:,k))
  end do

  if(debug) write(6,*)'in normal_new_factorization_rf_z, min,max(diag)=',minval(diag),maxval(diag)

end subroutine normal_new_factorization_rf_z

subroutine normal_new_factorization_rf_x
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    normal_new_factorization_rf_x get normalization factor in x direction
!
!   prgrmmr:  parrish        org: np22                              date: 2009-09-28
!
! abstract:  compute normalization factor in longitude direction
!
! program history log:
!   2009-09-28  parrish  initial documentation
!   2010-02-20  parrish  modifications for dual resolution
!   2010-03-11  parrish  correct error that can lead to infinite loop, and introduce grd_ens%kend_alloc
!                         in dimension statements
!   2013-04-17  wu       use grd_loc instead of grd_ens when defining normalization factor
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  use hybrid_ensemble_parameters, only: grd_loc
  use constants, only: zero,one

  integer(i_kind) i,j,k,kcount,lcount,iadvance,iback
  real(r_kind) f(grd_loc%nlat,grd_loc%nlon,grd_loc%kend_alloc+1-grd_loc%kbegin_loc)
  real(r_kind) diag(grd_loc%nlat,grd_loc%nlon)

!                       possible to have kend_loc - kbegin_loc-1 for processors not involved
!                          which results in infinite loops

  if(grd_loc%kend_loc < grd_loc%kbegin_loc) return

  if(allocated(xnorm_new)) deallocate(xnorm_new)
  allocate(xnorm_new(grd_loc%nlat,grd_loc%nlon))
  xnorm_new=one

  kcount=0
  lcount=0
  do
     f=zero
     do k=1,min(grd_loc%nlon,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
        kcount=kcount+1
        do i=1,grd_loc%nlat
           f(i,kcount,k)=one
        end do
        if(kcount == grd_loc%nlon) exit
     end do
     iadvance=1 ; iback=2
     call new_factorization_rf_x(f,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
     iadvance=2 ; iback=1
     call new_factorization_rf_x(f,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
     do k=1,min(grd_loc%nlon,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
        lcount=lcount+1
        do i=1,grd_loc%nlat
           diag(i,lcount)=sqrt(one/f(i,lcount,k))
        end do
        if(lcount == grd_loc%nlon) exit
     end do
     if(lcount == grd_loc%nlon) exit
  end do
  do j=1,grd_loc%nlon
     do i=1,grd_loc%nlat
        xnorm_new(i,j)=diag(i,j)
     end do
  end do

!           check accuracy of xnorm
  kcount=0
  lcount=0
  do
     f=zero
     do k=1,min(grd_loc%nlon,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
        kcount=kcount+1
        do i=1,grd_loc%nlat
           f(i,kcount,k)=one
        end do
        if(kcount == grd_loc%nlon) exit
     end do
     iadvance=1 ; iback=2
     call new_factorization_rf_x(f,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
     iadvance=2 ; iback=1
     call new_factorization_rf_x(f,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
     do k=1,min(grd_loc%nlon,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
        lcount=lcount+1
        do i=1,grd_loc%nlat
           diag(i,lcount)=f(i,lcount,k)
        end do
        if(lcount == grd_loc%nlon) exit
     end do
     if(lcount == grd_loc%nlon) exit
  end do
  if(debug) write(6,*)' in normal_new_factorization_rf_x, min,max(diag)=',minval(diag),maxval(diag)

end subroutine normal_new_factorization_rf_x

subroutine normal_new_factorization_rf_y
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    normal_new_factorization_rf_y  get normalization factor in y direction
!
!   prgrmmr:  parrish        org: np22                              date: 2009-09-28
!
! abstract:  compute normalization factor in latitude direction
!
! program history log:
!   2009-09-28  parrish  initial documentation
!   2010-02-20  parrish  modifications for dual resolution
!   2013-04-17  wu       use grd_loc instead of grd_ens when defining normalization factor
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  use hybrid_ensemble_parameters, only: grd_loc
  use constants, only: zero,one
  use mpimod, only: mype
  implicit none

  integer(i_kind) k,kcount,lcount,iadvance,iback
  real(r_kind) f(grd_loc%nlat,grd_loc%nlon*(grd_loc%kend_alloc+1-grd_loc%kbegin_loc)),diag(grd_loc%nlat)

!                       possible to have kend_loc - kbegin_loc-1 for processors not involved
!                          which results in infinite loops

  if(grd_loc%kend_loc < grd_loc%kbegin_loc) return

  if(allocated(ynorm_new)) deallocate(ynorm_new)
  allocate(ynorm_new(grd_loc%nlat))
  ynorm_new=one

  kcount=0
  lcount=0
  do
     f=zero
     do k=1,min(grd_loc%nlon*(grd_loc%kend_loc+1-grd_loc%kbegin_loc),grd_loc%nlat)
        kcount=kcount+1
        f(kcount,k)=one
        if(kcount == grd_loc%nlat) exit
     end do
     iadvance=1 ; iback=2
     call new_factorization_rf_y(f,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
     iadvance=2 ; iback=1
     call new_factorization_rf_y(f,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
     do k=1,min(grd_loc%nlon*(grd_loc%kend_loc+1-grd_loc%kbegin_loc),grd_loc%nlat)
        lcount=lcount+1
        diag(lcount)=sqrt(one/f(lcount,k))
        if(lcount == grd_loc%nlat) exit
     end do
     if(lcount == grd_loc%nlat) exit
  end do
  do k=1,grd_loc%nlat
     ynorm_new(k)=diag(k)
  end do

!               check that ynorm is corect
  kcount=0
  lcount=0
  do
     f=zero
     do k=1,min(grd_loc%nlon*(grd_loc%kend_loc+1-grd_loc%kbegin_loc),grd_loc%nlat)
        kcount=kcount+1
        f(kcount,k)=one
        if(kcount == grd_loc%nlat) exit
     end do
     iadvance=1 ; iback=2
     call new_factorization_rf_y(f,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
     iadvance=2 ; iback=1
     call new_factorization_rf_y(f,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
     do k=1,min(grd_loc%nlon*(grd_loc%kend_loc+1-grd_loc%kbegin_loc),grd_loc%nlat)
        lcount=lcount+1
        diag(lcount)=f(lcount,k)
        if(lcount == grd_loc%nlat) exit
     end do
     if(lcount == grd_loc%nlat) exit
  end do
  if(debug) write(6,*)' in normal_new_factorization_rf_y, min,max(diag)=',minval(diag),maxval(diag)

end subroutine normal_new_factorization_rf_y

  subroutine create_ensemble
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    create_ensemble        allocate space for ensembles
!   prgmmr: parrish          org: np22                date: 2009-06-16
!
! abstract: allocate space for ensemble perturbations used with the 
!             hybrid ensemble option.
!
! program history log:
!   2009-06-16  parrish
!   2010-02-20  parrish  modifications for dual resolution
!   2011-02-28  parrish - introduce more complete use of gsi_bundlemod to eliminate hard-wired variables
!   2011-08-31  todling - revisit en_perts (single-prec) in light of extended bundle
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use hybrid_ensemble_parameters, only: n_ens,grd_ens,ntlevs_ens

    implicit none

    type(gsi_grid)  :: grid_ens

    integer(i_kind) n,istatus,m
    character(len=*),parameter::myname_=trim(myname)//'*create_ensemble'

    nelen=grd_ens%latlon11*(max(0,nc3d)*grd_ens%nsig+max(0,nc2d))
!   create ensemble perturbations bundles (using newly added r_single capability

    allocate(en_perts(n_ens,ntlevs_ens))
    call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
 
   do m=1,ntlevs_ens
       do n=1,n_ens
          call gsi_bundlecreate(en_perts(n,m),grid_ens,'ensemble perts',istatus, &
                                names2d=cvars2d,names3d=cvars3d,bundle_kind=r_single)
          if(istatus/=0) then
             write(6,*)trim(myname_),': trouble creating en_perts bundle'
             call stop2(999)
          endif
       end do
    end do


    allocate(ps_bar(grd_ens%lat2,grd_ens%lon2,ntlevs_ens) )
    if(debug) write(6,*)' in create_ensemble, grd_ens%latlon11,grd_ens%latlon1n,n_ens,ntlevs_ens=', &
                                    grd_ens%latlon11,grd_ens%latlon1n,n_ens,ntlevs_ens
    if(debug) write(6,*)' in create_ensemble, total bytes allocated=',4*nelen*n_ens*ntlevs_ens

  end subroutine create_ensemble

  subroutine load_ensemble
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    load_ensemble    read/generate ensemble perturbations
!   prgmmr: parrish          org: np22                date: 2009-09-11
!
! abstract: read or generate (if generate_ens=.true.) ensemble
!             perturbations used for hybrid ensemble option.
!
! program history log:
!   2009-09-11  parrish
!   2010-02-20  parrish  modifications for dual resolution
!   2010-03-15  zhu      make changes using cstate
!   2010-03-28  todling  update to use gsi_bundle
!   2010-06-03  parrish  multiple fixes for dual resolution
!   2011-02-28  parrish - introduce more complete use of gsi_bundlemod to eliminate hard-wired variables
!   2011-07-05  carley  - add call to get_nmmb_ensperts
!   2011-12-07  tong    - add the option to read wrf_nmm ensemble
!   2012-01-30  parrish - remove wrf_nmm_regional,wrf_mass_regional,netcdf,nems_nmmb_regional
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use gridmod, only: regional
    use constants, only: zero,one
    use hybrid_ensemble_parameters, only: n_ens,generate_ens,grd_ens,grd_anl,ntlevs_ens, &
                                          pseudo_hybens,regional_ensemble_option
    use mpimod, only: mype,ierror
    implicit none

    type(gsi_bundle),allocatable:: en_bar(:)
    type(gsi_bundle):: bundle_anl,bundle_ens
    type(gsi_grid)  :: grid_anl,grid_ens
    integer(i_kind) i,j,k,n,ii,ic2,m
    integer(i_kind) istatus,ier
    real(r_kind),allocatable:: seed(:,:)
    real(r_kind),pointer,dimension(:,:)   :: cv_ps
    real(r_kind) sig_norm,bar_norm
    character(len=*),parameter::myname_=trim(myname)//'*load_ensemble'
    character(50) title

!      create simple regular grid

    if(generate_ens) then


       call gsi_gridcreate(grid_anl,grd_anl%lat2,grd_anl%lon2,grd_anl%nsig)
       call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)

       allocate(en_bar(ntlevs_ens))

       do m=1,ntlevs_ens
          call gsi_bundlecreate(en_bar(m),grid_ens,'ensemble',istatus, &
                                names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind)
       end do
       if(istatus/=0) then
          write(6,*)trim(myname_),': trouble creating en_bar bundle'
          call stop2(999)
       endif
       sig_norm=sqrt(one/max(one,n_ens-one))
       bar_norm=one/n_ens
       if(n_ens == 1) bar_norm=zero

!                        initialize subdomain to slab routine special_sd2h
       call special_sd2h0
       allocate(seed(nval2f,nscl))
       seed=-one
 
       do m=1,ntlevs_ens
         en_bar(m)%values=zero
       end do

!      create two internal bundles, one on analysis grid and one on ensemble grid

       call gsi_bundlecreate (bundle_anl,grid_anl,'ensemble work',istatus, &
                                 names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind)
       if(istatus/=0) then
          write(6,*)trim(myname_),': trouble creating bundle_anl bundle'
          call stop2(999)
       endif
       call gsi_bundlecreate (bundle_ens,grid_ens,'ensemble work ens',istatus, &
                                 names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind)
       if(istatus/=0) then
          write(6,*)trim(myname_),': trouble creating bundle_ens bundle'
          call stop2(999)
       endif

       do m=1,ntlevs_ens
          do n=1,n_ens
             call generate_one_ensemble_perturbation(bundle_anl,bundle_ens,seed)
             do ii=1,nelen
                en_perts(n,m)%valuesr4(ii)=bundle_ens%values(ii)
                en_bar(m)%values(ii)=en_bar(m)%values(ii)+bundle_ens%values(ii)
             end do
          end do
       
! Load ps_bar for use with vertical localization later
          call gsi_bundlegetpointer (en_bar(m),'ps' ,cv_ps ,istatus)
          do j=1,grd_ens%lon2
             do i=1,grd_ens%lat2
                ps_bar(i,j,m)=cv_ps(i,j)*bar_norm
             end do
          end do
       end do

! do some cleanning
       call gsi_bundledestroy(bundle_anl,istatus)
       if(istatus/=0) then
          write(6,*)trim(myname_),': trouble destroying bundle_anl bundle'
          call stop2(999)
       endif
       call gsi_bundledestroy(bundle_ens,istatus)
       if(istatus/=0) then
          write(6,*)trim(myname_),': trouble destroying bundle_ens bundle'
          call stop2(999)
       endif
!                          remove mean, which is locally significantly non-zero, due to sample size.
!                           with real ensembles, the mean of the actual sample will be removed.

       do m=1,ntlevs_ens
          do n=1,n_ens
             do ii=1,nelen
                en_perts(n,m)%valuesr4(ii)=(en_perts(n,m)%valuesr4(ii)-en_bar(m)%values(ii)*bar_norm)*sig_norm
             end do
          end do

          call gsi_bundledestroy(en_bar(m),istatus)
          if(istatus/=0) then
          write(6,*)trim(myname_),': trouble destroying en_bar bundle'
          call stop2(999)
	  end if
       end do

       deallocate(en_bar)

    else

!            read in ensembles

       if (.not.regional) then

          call get_gefs_ensperts_dualres

       else

          if(regional_ensemble_option < 1 .or. regional_ensemble_option > 4) then
             if(mype==0) then
                write(6,'(" IMPROPER CHOICE FOR ENSEMBLE INPUT IN SUBROUTINE LOAD_ENSEMBLE")')
                write(6,'(" regional_ensemble_option = ",i5)') regional_ensemble_option
                write(6,'(" allowed values of regional_ensemble_option:")')
                write(6,'("   regional_ensemble_option=1:")')
                write(6,'("                    use GEFS internally interpolated to ensemble grid.")')
                write(6,'("                    can add pseudo ensemble hybrid option for hwrf.")')
                write(6,'("   regional_ensemble_option=2:")')
                write(6,'("                    ensembles are WRF NMM (HWRF) format.")')
                write(6,'("   regional_ensemble_option=3:")')
                write(6,'("                    ensembles are ARW netcdf format.")')
                write(6,'("   regional_ensemble_option=4:")')
                write(6,'("                    ensembles are NEMS NMMB format.")')
             end if
             call stop2(999)
          end if

          select case(regional_ensemble_option)

             case(1)

!     regional_ensemble_option = 1: use GEFS internally interpolated to ensemble grid.

                call get_gefs_for_regional

!     pseudo_hybens = .true.: pseudo ensemble hybrid option for hwrf
!                             GEFS ensemble perturbations in TC vortex area
!                             are replaced with TC vortex library perturbations
                if (pseudo_hybens) then
                   call get_pseudo_ensperts
                end if

             case(2)

!     regional_ensemble_option = 2: ensembles are WRF NMM (HWRF) format

                call get_wrf_nmm_ensperts

             case(3)

!     regional_ensemble_option = 3: ensembles are ARW netcdf format.

                call get_wrf_mass_ensperts_netcdf

             case(4)

!     regional_ensemble_option = 4: ensembles are NEMS NMMB format.

                call get_nmmb_ensperts

          end select

       end if

    end if

  end subroutine load_ensemble

  subroutine generate_one_ensemble_perturbation(bundle_anl,bundle_ens,seed)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    generate_one_ensemble_perturbation
!
!   prgrmmr:  parrish        org: np22                              date: 2009-09-28
!
! abstract:  compute normalization factor in latitude direction
!
! program history log:
!   2009-09-28  parrish  initial documentation
!   2010-02-20  parrish  modifications for dual resolution
!   2010-03-14  zhu     - make changes using cstate
!   2010-03-15  derber  - fix qvar3d for ensemble
!   2010-04-28  todling - update to use gsi_bundle
!   2010-05-18  parrish  reactivate dual resolution
!   2010-06-04  parrish  multiple fixes for dual resolution
!   2011-02-28  parrish - introduce more complete use of gsi_bundlemod to eliminate hard-wired variables
!   2011-09-27  parrish - ckgcov is called here, which uses sqrt_smootrf.  The variable nval_lenz is passed
!                         to sqrt_smoothrf through module jfunc, but if lsqrtb=F, then nval_lenz=0.
!                         To fix this, nval_lenz is passed in to this subroutine through module jfunc,
!                         then is temporarily replaced with the proper non-zero value, and restored back
!                         to its original value after exiting from the call to ckgcov.  Also, this is
!                         now an argument in ckgcov, which was missing in this call to ckgcov.
!   2012-06-12  parrish - remove variable nvar_pe (not used)
!   2013-01-12  parrish - remove extra argument nnnn1o from call ckgcov--no longer used
!
!   input argument list:
!     seed     - old random number seeds (used for bit reproducibility of
!                 generated random ensemble perturbations on different 
!                 numbers of processors)
!
!   output argument list:
!     seed     - new random number seeds
!     st       - stream function part of generated ensemble perturbation
!     vp       - velocity potential part of generated ensemble perturbation
!     t        - virtual temperature part of generated ensemble perturbation
!     rh       - relative humidity part of generated ensemble perturbation
!     oz       - ozone part of generated ensemble perturbation
!     cw       - cloud water part of generated ensemble perturbation
!     p        - surface pressure part of generated ensemble perturbation
!     sst      - skin temperature part of generated ensemble perturbation
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

    use kinds, only: r_kind,i_kind,i_llong
    use gridmod, only: vlevs,nnnn1o,regional
    use mpimod, only: mype,mpi_rtype,mpi_comm_world,ierror
    use hybrid_ensemble_parameters, only: uv_hyb_ens,grd_ens,grd_anl,p_e2a
    use general_sub2grid_mod, only: general_suba2sube
    use constants, only: zero,one
    use jfunc, only: nval_lenz
    implicit none

    character(len=*),parameter::myname_=myname//'*generate_one_ensemble_perturbation'
    real(r_kind)    ,intent(inout) :: seed(nval2f,nscl)
    type(gsi_bundle),intent(inout) :: bundle_anl,bundle_ens

    real(r_kind),dimension(nval2f,nnnn1o,nscl):: z
    real(r_kind) vert1(vlevs)
    integer(i_llong) iseed
    integer(i_kind) nvert,i,j,ii,is,naux,k,ic3
    integer(i_kind) istat_st,istat_vp
    integer(i_kind) nval_lenz_save
    real(r_kind) aux
    real(r_kind),dimension(nh_0:nh_1,vlevs,nscl):: zsub
    real(r_kind),dimension(:,:,:),allocatable:: ua,va
    real(r_kind),pointer,dimension(:,:,:):: st,vp

    naux=0
    nvert=vlevs
    if(maxval(seed) <  zero) then

!       create initial seed for random numbers for each horizontal location.
  
       if(mype == 0) then
          call random_number(seed)
          do is=1,nscl
             do i=1,nval2f
                iseed=1+nint(seed(i,is)*1234567._r_kind)
                seed(i,is)=iseed
             end do
          end do
       end if
       call mpi_bcast(seed,nval2f*nscl,mpi_rtype,0,mpi_comm_world,ierror)

    end if

    do is=1,nscl
       do i=nh_0,nh_1
#ifdef ibm_sp
          call dnrand(seed(i,is),nvert,vert1,aux,naux)
#else
! Generate a Vector of Normally Distributed Random Numbers (function from Lapack lib)
          call dlarnv(3,seed(i,is),nvert,vert1)
#endif
          do k=1,nvert
             zsub(i,k,is)=vert1(k)
          end do
       end do
    end do
    call special_sd2h(zsub,z)

!     if this is a global run, then need to fix tropical belt part of z so periodic overlap is correct
    if(.not.regional) call fix_belt(z)

!     temporarily redefine nval_lenz
    nval_lenz_save=nval_lenz
    nval_lenz=nval2f*nnnn1o*nscl
    call ckgcov(z,bundle_anl,nval_lenz)
!     restore nval_lenz
    nval_lenz=nval_lenz_save

!     if uv_hyb_ens=.true., then convert st,vp to u,v
    if(uv_hyb_ens) then
       allocate(ua(grd_anl%lat2,grd_anl%lon2,grd_anl%nsig))
       allocate(va(grd_anl%lat2,grd_anl%lon2,grd_anl%nsig))
       istat_st=-999
       istat_vp=-999
       do ic3=1,nc3d
          if(trim(cvars3d(ic3))=='sf') call gsi_bundlegetpointer (bundle_anl, cvars3d(ic3),st, istat_st)
          if(trim(cvars3d(ic3))=='vp') call gsi_bundlegetpointer (bundle_anl, cvars3d(ic3),vp, istat_vp)
       end do
       if(istat_st/=0.or.istat_vp/=0) then
          write(6,*) myname_,': error getting sf/vp pointers, aborting ...'
          call stop2(999)
       endif
       call getuv(ua,va,st,vp,0)
       st=ua
       vp=va
       deallocate(va)
       deallocate(ua)
    end if

    if(p_e2a%identity) then
       bundle_ens%values=bundle_anl%values
    else
       call general_suba2sube(grd_anl,grd_ens,p_e2a,bundle_anl%values,bundle_ens%values,regional)
    end if


  end subroutine generate_one_ensemble_perturbation

  subroutine fix_belt(z)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    fix_belt
!
!   prgrmmr:  parrish        org: np22                              date: 2009-09-28
!
! abstract: when generating random vector representation of control variable z for
!            global case, need to make adjustment to tropical belt part to
!            properly account for periodicity in overlap zone
!
! program history log:
!   2009-09-28  parrish  initial documentation
!   2010-02-20  parrish  modifications for dual resolution
!
!   input argument list:
!     z        - field to be adjusted
!
!   output argument list:
!     z        - adjusted field
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block


    use kinds, only: r_kind,i_kind
    use gridmod, only: nnnn1o
    use berror, only: nx,ny,nf
    use hybrid_ensemble_parameters, only: grd_ens
    implicit none

    real(r_kind),intent(inout) :: z(nval2f,nnnn1o,nscl)

    real(r_kind) zloc1(ny,nx)
    integer(i_kind) i,ii,j,jj,k

    do j=1,nscl
       do k=1,nnnn1o
          i=0
          do jj=1,nx
             do ii=1,ny
                i=i+1
                zloc1(ii,jj)=z(i,k,j)
             end do
          end do
          do jj=grd_ens%nlon+1,nx
             do ii=1,ny
                zloc1(ii,jj)=zloc1(ii,jj-grd_ens%nlon)
             end do
          end do
          i=0
          do jj=1,nx
             do ii=1,ny
                i=i+1
                z(i,k,j)=zloc1(ii,jj)
             end do
          end do
       end do
    end do
    
  end subroutine fix_belt

  subroutine rescale_ensemble_rh_perturbations
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    rescale_ensemble_rh_perturbations
!   prgmmr: parrish          org: np22                date: 2009-10-15
!
! abstract: rescale internally generated ensemble rh perturbations
!
! program history log:
!   2009-10-15  parrish
!   2010-02-20  parrish  modifications for dual resolution
!   2011-02-28  parrish  bundle changes
!   2011-07-08  todling  adapt to fit upd interface to general_sub2grid
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use kinds, only: r_kind,i_kind
    use gridmod, only: regional
    use hybrid_ensemble_parameters, only: n_ens,grd_ens,grd_anl,grd_a1,grd_e1,p_e2a,ntlevs_ens
    use general_sub2grid_mod, only: general_suba2sube
    use berror, only: qvar3d
    implicit none

    integer(i_kind) i,ii,j,k,n,istatus,m
    real(r_kind) qvar3d_ens(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig,1)
    real(r_single),pointer,dimension(:,:,:):: w3

    call gsi_bundlegetpointer(en_perts(1,1),'q',w3,istatus)
    if(istatus/=0) then
       write(6,*)' rh variable not available, skip subroutine rescale_ensemble_rh_perturbations'
       return
    end if

    if(grd_anl%latlon11 == grd_ens%latlon11) then
       qvar3d_ens=reshape(qvar3d,(/size(qvar3d,1),size(qvar3d,2),size(qvar3d,3),1/))
    else
       call general_suba2sube(grd_a1,grd_e1,p_e2a, &
            reshape(qvar3d,(/size(qvar3d,1),size(qvar3d,2),size(qvar3d,3),1/)),qvar3d_ens,regional)
    end if
!$omp parallel do schedule(dynamic,1) private(n,ii,i,j,k,w3,istatus)
    do m=1,ntlevs_ens
       do n=1,n_ens
          call gsi_bundlegetpointer(en_perts(n,m),'q',w3,istatus)
          if(istatus/=0) then
             write(6,*)' error retrieving pointer to rh variable for ensemble number ',n
             call stop2(999)
          end if
          do k=1,grd_ens%nsig
             do j=1,grd_ens%lon2
                do i=1,grd_ens%lat2
                   w3(i,j,k)=qvar3d_ens(i,j,k,1)*w3(i,j,k)
                end do
             end do
          end do
       end do
    end do
 
  end subroutine rescale_ensemble_rh_perturbations

  subroutine destroy_ensemble
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    destroy_ensemble       deallocate space for ensembles
!   prgmmr: parrish          org: np22                date: 2009-06-16
!
! abstract: deallocate space for ensemble perturbations used with the 
!             hybrid ensemble option.
!
! program history log:
!   2009-06-16  parrish
!   2011-02-28  parrish, replace specific ensemble perturbation arrays with pseudo-bundle en_perts array
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use hybrid_ensemble_parameters, only: l_hyb_ens,n_ens,ntlevs_ens
    implicit none

    integer(i_kind) istatus,n,m

    if(l_hyb_ens) then
       do m=ntlevs_ens,1,-1
          do n=n_ens,1,-1
             call gsi_bundleunset(en_perts(n,m),istatus)
             if(istatus/=0) then
                write(6,*)'in destroy_ensemble: trouble destroying en_perts bundle'
                call stop2(999)
             endif
          end do
       end do
       deallocate(ps_bar)
       deallocate(en_perts)
    end if

  end subroutine destroy_ensemble

  subroutine ensemble_forward_model(cvec,a_en,ibin)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ensemble_forward_model  add ensemble part to anl vars
!   prgmmr: parrish          org: np22                date: 2009-06-16
!
! abstract: For the hybrid ensemble method, add ensemble contribution
!             to standard analysis control variables.  (This follows,
!             method outlined in Wang et al, MWR, 2008).

! program history log:
!   2009-09-11  parrish
!   2010-02-20  parrish  modifications for dual resolution
!   2010-03-23  zhu - use cstate
!   2010-05-07  parrish - remove error stop for dual resolution
!   2010-04-28  todling - update to use gsi_bundle
!   2011-02-28  parrish - bundle changes
!   2011-10-03  wu - add option to weight ensemble contribution to surface pressure with vertical profile
!   2011-10-26  kleist - 4d capability for ensemble/hybrid
!
!   input argument list:
!     cvec      - 
!     a_en      - 
!     ibin      - integer bin number for ensemble perturbations  
! 
!   output argument list:
!     cvec      - 
!
! remarks:  
!    need to reconcile grid in gsi_bundle w/ grid_ens/grid_anl
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use hybrid_ensemble_parameters, only: n_ens,grd_ens,grd_anl
    use hybrid_ensemble_parameters, only: pwgt,pwgtflg

    implicit none
    type(gsi_bundle),intent(inout) :: cvec
    type(gsi_bundle),intent(in)    :: a_en(:)
    integer,intent(in)		   :: ibin

    character(len=*),parameter::myname_=trim(myname)//'*ensemble_forward_model'
    logical nogood
    integer(i_kind) i,j,k,n,im,jm,km,ic2,ic3
    integer(i_kind) ipc3d(nc3d),ipc2d(nc2d),ipe(1),istatus

!   Check resolution consistency between static and ensemble components
    nogood=.not.(cvec%grid%im==a_en(1)%grid%im.and.&
                 cvec%grid%jm==a_en(1)%grid%jm.and.&
                 cvec%grid%km==a_en(1)%grid%km)
    if (nogood) then
       write(6,*) myname_,': static&ensemble vectors have inconsistent dims'
       call stop2(999)
    endif

!   Request ensemble-corresponding fields from control vector
!    NOTE:  because ensemble perturbation bundle structure is same as control vector, use same ipc3d and
!             ipc2d indices for cvec and en_perts bundles.
    call gsi_bundlegetpointer (cvec,cvars3d,ipc3d,istatus)
    if(istatus/=0) then
      write(6,*) myname_,': cannot find 3d pointers'
      call stop2(999)
    endif
    call gsi_bundlegetpointer (cvec,cvars2d,ipc2d,istatus)
    if(istatus/=0) then
      write(6,*) myname_,': cannot find 2d pointers'
      call stop2(999)
    endif
 
    im=cvec%grid%im
    jm=cvec%grid%jm
    km=cvec%grid%km
    ipe(1)=1
       do n=1,n_ens
          do ic3=1,nc3d
             do k=1,km
                do j=1,jm
                   do i=1,im
                      cvec%r3(ipc3d(ic3))%q(i,j,k)=cvec%r3(ipc3d(ic3))%q(i,j,k) &
                            +a_en(n)%r3(ipe(1))%q(i,j,k)*en_perts(n,ibin)%r3(ipc3d(ic3))%qr4(i,j,k)
                   end do
                end do
             end do
          end do

          do ic2=1,nc2d

             select case (trim(cvars2d(ic2)))
 
                case('ps','PS')
 
        if(pwgtflg)then
                   do k=1,km
                      do j=1,jm
                         do i=1,im
                         cvec%r2(ipc2d(ic2))%q(i,j)=cvec%r2(ipc2d(ic2))%q(i,j) &
                            +a_en(n)%r3(ipe(1))%q(i,j,k)*pwgt(i,j,k)*en_perts(n,ibin)%r2(ipc2d(ic2))%qr4(i,j)
                         end do
                      end do
                   end do
        else
                   do j=1,jm
                      do i=1,im
                         cvec%r2(ipc2d(ic2))%q(i,j)=cvec%r2(ipc2d(ic2))%q(i,j) &
                            +a_en(n)%r3(ipe(1))%q(i,j,1)*en_perts(n,ibin)%r2(ipc2d(ic2))%qr4(i,j)
                      end do
                   end do
        endif
 
                case('sst','SST')
 
                   do j=1,jm
                      do i=1,im
                         cvec%r2(ipc2d(ic2))%q(i,j)=cvec%r2(ipc2d(ic2))%q(i,j) &
                            +a_en(n)%r3(ipe(1))%q(i,j,1)*en_perts(n,ibin)%r2(ipc2d(ic2))%qr4(i,j)
                      end do
                   end do
 
             end select

          end do
       end do !end do n_ens

  end subroutine ensemble_forward_model

  subroutine ensemble_forward_model_dual_res(cvec,a_en,ibin)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ensemble_forward_model_dual_res  use for dualres option
!   prgmmr: parrish          org: np22                date: 2010-02-20
!
! abstract: Copy of ensemble_forward_model for use with dual resolution.

! program history log:
!   2010-02-20  parrish
!   2010-03-23  zhu - use cstate
!   2010-04-06  parrish - add deallocate(sube_vars)
!   2010-04-28  todling - update to use gsi_bundle
!   2010-05-07  parrish - remove error stop for dual resolution, and add cstate again
!                         to analysis variables after interpolation from
!                         ensemble grid.
!                         the ensemble part is still not updated for
!                         generalized control variable
!   2010-05-18  todling - revisited bundle usage in light of Dave's change (2010-05-07)
!   2011-02-28  parrish - bundle changes
!   2011-10-03  wu - add option to weight ensemble contribution to surface pressure with vertical profile
!   2011-10-26  kleist  - 4d capability for ensemble/hybrid
!
!   input argument list:
!     cvec      -
!     a_en      -
!     ibin      - integer bin number for ensemble perturbations
!
!   output argument list:
!     cvec      - 
!
! remarks:  
!    need to reconcile grid in gsi_bundle w/ grid_ens/grid_anl
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use hybrid_ensemble_parameters, only: n_ens,grd_ens,grd_anl,p_e2a,uv_hyb_ens
    use hybrid_ensemble_parameters, only: pwgt,pwgtflg
    use general_sub2grid_mod, only: general_sube2suba
    use gridmod,only: regional
    use constants, only: zero
    implicit none

    type(gsi_bundle),intent(inout) :: cvec
    type(gsi_bundle),intent(in)    :: a_en(:)
    integer,intent(in)             :: ibin

    character(len=*),parameter::myname_=trim(myname)//'*ensemble_forward_model_dual_res'
    integer(i_kind) i,j,k,n,im,jm,km,ic2,ic3
    type(gsi_grid)  :: grid_ens,grid_anl
    type(gsi_bundle)  :: work_ens,work_anl
    integer(i_kind) ipc2d(nc2d),ipc3d(nc3d),ipe(1),istatus

!   Request ensemble-corresponding fields from control vector
!    NOTE:  because ensemble perturbation bundle structure is same as control vector, use same ipc3d and
!             ipc2d indices for cvec and en_perts bundles.
    call gsi_bundlegetpointer (cvec,cvars3d,ipc3d,istatus)
    if(istatus/=0) then
      write(6,*) myname_,': cannot find 3d pointers'
      call stop2(999)
    endif
    call gsi_bundlegetpointer (cvec,cvars2d,ipc2d,istatus)
    if(istatus/=0) then
      write(6,*) myname_,': cannot find 2d pointers'
      call stop2(999)
    endif

    call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
    call gsi_bundlecreate (work_ens,grid_ens,'ensemble work',istatus, &
                              names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind)
    if(istatus/=0) then
       write(6,*)trim(myname_),': trouble creating work_ens bundle'
       call stop2(999)
    endif
    call gsi_gridcreate(grid_anl,grd_anl%lat2,grd_anl%lon2,grd_anl%nsig)
    call gsi_bundlecreate (work_anl,grid_anl,'analysis work',istatus, &
                              names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind)
    if(istatus/=0) then
       write(6,*)trim(myname_),': trouble creating work_anl bundle'
       call stop2(999)
    endif


    ipe(1)=1
    work_ens%values(:)=zero
    im=work_ens%grid%im
    jm=work_ens%grid%jm
    km=work_ens%grid%km
       do n=1,n_ens
          do ic3=1,nc3d
             do k=1,km
                do j=1,jm
                   do i=1,im
                      work_ens%r3(ipc3d(ic3))%q(i,j,k)=work_ens%r3(ipc3d(ic3))%q(i,j,k) &
                         +a_en(n)%r3(ipe(1))%q(i,j,k)*en_perts(n,ibin)%r3(ipc3d(ic3))%qr4(i,j,k)
                   end do
                end do
             end do
          end do
          do ic2=1,nc2d

             select case (trim(cvars2d(ic2)))

                case('ps','PS')

        if(pwgtflg)then
                   do k=1,km
                      do j=1,jm
                         do i=1,im
                         work_ens%r2(ipc2d(ic2))%q(i,j)=work_ens%r2(ipc2d(ic2))%q(i,j) &
                            +a_en(n)%r3(ipe(1))%q(i,j,k)*pwgt(i,j,k)*en_perts(n,ibin)%r2(ipc2d(ic2))%qr4(i,j)
                         end do
                      end do
                   end do
        else
                   do j=1,jm
                      do i=1,im
                         work_ens%r2(ipc2d(ic2))%q(i,j)=work_ens%r2(ipc2d(ic2))%q(i,j) &
                            +a_en(n)%r3(ipe(1))%q(i,j,1)*en_perts(n,ibin)%r2(ipc2d(ic2))%qr4(i,j)
                      end do
                   end do
        endif

                case('sst','SST')

                   do j=1,jm
                      do i=1,im
                         work_ens%r2(ipc2d(ic2))%q(i,j)=work_ens%r2(ipc2d(ic2))%q(i,j) &
                            +a_en(n)%r3(ipe(1))%q(i,j,1)*en_perts(n,ibin)%r2(ipc2d(ic2))%qr4(i,j)
                      end do
                   end do

             end select

          end do
       end do !n_ens

    call general_sube2suba(grd_ens,grd_anl,p_e2a,work_ens%values,work_anl%values,regional)
    call gsi_bundledestroy(work_ens,istatus)
    if(istatus/=0) then
       write(6,*)trim(myname_),': trouble destroying work ens bundle'
       call stop2(999)
    endif
    do ic3=1,nc3d
       cvec%r3(ipc3d(ic3))%q=cvec%r3(ipc3d(ic3))%q+work_anl%r3(ipc3d(ic3))%q
    end do
    do ic2=1,nc2d
       cvec%r2(ipc2d(ic2))%q=cvec%r2(ipc2d(ic2))%q+work_anl%r2(ipc2d(ic2))%q
    end do
    call gsi_bundledestroy(work_anl,istatus)
    if(istatus/=0) then
       write(6,*)trim(myname_),': trouble destroying work anl bundle'
       call stop2(999)
    endif

  end subroutine ensemble_forward_model_dual_res

  subroutine ensemble_forward_model_ad(cvec,a_en,ibin)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ensemble_forward_model  add ensemble part to anl vars
!   prgmmr: parrish          org: np22                date: 2009-06-16
!
! abstract: For the hybrid ensemble method, add ensemble contribution
!             to standard analysis control variables.  (This follows,
!             method outlined in Wang et al, MWR, 2008).

! program history log:
!   2009-09-11  parrish
!   2010-02-20  parrish - adapt for dual resolution
!   2010-03-23  zhu - use cstate
!   2010-04-28  todling - update to use gsi_bundle
!   2011-10-03  wu - add option to weight ensemble contribution to surface pressure with vertical profile
!   2011-11-01  kleist  - 4d capability for ensemble/hybrid
!
!   input argument list:
!     cvec      -
!     a_en      -
!     ibin      - integer bin number for ensemble perturbations
!
!   output argument list:
!     cvec      -
!     a_en      -
!
! remarks:  
!    need to reconcile grid in gsi_bundle w/ grid_ens/grid_anl
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

    use hybrid_ensemble_parameters, only: n_ens,grd_ens,grd_anl
    use hybrid_ensemble_parameters, only: pwgt,pwgtflg
    implicit none

    type(gsi_bundle),intent(inout) :: cvec
    type(gsi_bundle),intent(inout) :: a_en(:)
    integer,intent(in)             :: ibin

    character(len=*),parameter::myname_=trim(myname)//'*ensemble_forward_model_ad'
    logical nogood
    integer(i_kind) i,j,k,n,im,jm,km,ic2,ic3
    integer(i_kind) ipc3d(nc3d),ipc2d(nc2d),ipe(1),istatus

!   Check resolution consistency between static and ensemble components
    nogood=.not.(cvec%grid%im==a_en(1)%grid%im.and.&
                 cvec%grid%jm==a_en(1)%grid%jm.and.&
                 cvec%grid%km==a_en(1)%grid%km)
    if (nogood) then
       write(6,*) myname_,': static/ensemble vectors have inconsistent dims'
       call stop2(999)
    endif

!   Request ensemble-corresponding fields from control vector
!    NOTE:  because ensemble perturbation bundle structure is same as control vector, use same ipc3d and
!             ipc2d indices for cvec and en_perts bundles.
    call gsi_bundlegetpointer (cvec,cvars3d,ipc3d,istatus)
    if(istatus/=0) then
      write(6,*) myname_,': cannot find 3d pointers'
      call stop2(999)
    endif
    call gsi_bundlegetpointer (cvec,cvars2d,ipc2d,istatus)
    if(istatus/=0) then
      write(6,*) myname_,': cannot find 2d pointers'
      call stop2(999)
    endif

    im=cvec%grid%im
    jm=cvec%grid%jm
    km=cvec%grid%km
    ipe(1)=1
       do n=1,n_ens
          do ic3=1,nc3d
             do k=1,km
                do j=1,jm
                   do i=1,im
                         a_en(n)%r3(ipe(1))%q(i,j,k)=a_en(n)%r3(ipe(1))%q(i,j,k) &
                               +cvec%r3(ipc3d(ic3))%q(i,j,k)*en_perts(n,ibin)%r3(ipc3d(ic3))%qr4(i,j,k)
                   end do
                end do
             end do
          end do
          do ic2=1,nc2d

             select case (trim(cvars2d(ic2)))
 
                case('ps','PS')
 
                 if(pwgtflg)then
                   do k=1,km
                      do j=1,jm
                         do i=1,im
                            a_en(n)%r3(ipe(1))%q(i,j,k)=a_en(n)%r3(ipe(1))%q(i,j,k) &
                               +cvec%r2(ipc2d(ic2))%q(i,j)*en_perts(n,ibin)%r2(ipc2d(ic2))%qr4(i,j)*pwgt(i,j,k)
                         end do
                      end do
                   end do
                 else
                   do j=1,jm
                      do i=1,im
                         a_en(n)%r3(ipe(1))%q(i,j,1)=a_en(n)%r3(ipe(1))%q(i,j,1) &
                            +cvec%r2(ipc2d(ic2))%q(i,j)*en_perts(n,ibin)%r2(ipc2d(ic2))%qr4(i,j)
                      end do
                   end do
                 endif
  
                case('sst','SST')
  
                   do j=1,jm
                      do i=1,im
                         a_en(n)%r3(ipe(1))%q(i,j,1)=a_en(n)%r3(ipe(1))%q(i,j,1) &
                            +cvec%r2(ipc2d(ic2))%q(i,j)*en_perts(n,ibin)%r2(ipc2d(ic2))%qr4(i,j)
                      end do
                   end do
 
             end select
          end do
       end do !n_ens
  end subroutine ensemble_forward_model_ad

  subroutine ensemble_forward_model_ad_dual_res(cvec,a_en,ibin)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ensemble_forward_model_ad_dual_res  use for dualres option
!   prgmmr: parrish          org: np22                date: 2010-02-20
!
! abstract: Copy of ensemble_forward_model_ad for use with dual resolution.

! program history log:
!   2010-02-20  parrish
!   2010-03-23  zhu - use cstate
!   2010-04-06  parrish - correct dimensions of st,vp,t,rh,oz,cw,p,sst. add deallocate(suba_vars)
!   2010-04-28  todling - update to use gsi_bundle
!   2010-05-07  parrish - remove error stop for dual resolution, and add cstate again
!                           to analysis variables before adjoint interpolation
!                           to ensemble grid.
!                           the ensemble part is still not updated for
!                           generalized control variable
!   2010-05-18  todling - revisited bundle usage in light of Dave's change (2010-05-07)
!   2011-02-28  parrish - bundle changes
!   2011-10-03  wu - add option to weight ensemble contribution to surface pressure with vertical profile
!   2011-11-01  kleist  - 4d capability for ensemble/hybrid
!
!   input argument list:
!     cvec      -
!     a_en      -
!     ibin      - integer bin number for ensemble perturbations
!
!   output argument list:
!     cvec      -
!     a_en      -
!
! remarks:  
!    need to reconcile grid in gsi_bundle w/ grid_ens/grid_anl
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use constants, only: zero
    use hybrid_ensemble_parameters, only: n_ens,grd_ens,grd_anl,p_e2a,uv_hyb_ens
    use hybrid_ensemble_parameters, only: pwgt,pwgtflg
    use general_sub2grid_mod, only: general_sube2suba_ad
    use gridmod,only: regional
    implicit none

    type(gsi_bundle),intent(inout) :: cvec
    type(gsi_bundle),intent(inout) :: a_en(:)
    integer,intent(in)             :: ibin

    character(len=*),parameter::myname_=trim(myname)//'*ensemble_forward_model_ad_dual_res'
    integer(i_kind) i,j,k,n,im,jm,km,ic2,ic3
    type(gsi_grid)  :: grid_ens,grid_anl
    type(gsi_bundle)  :: work_ens,work_anl
    integer(i_kind) ipc2d(nc2d),ipc3d(nc3d),ipe(1),istatus

!   Request ensemble-corresponding fields from control vector
!    NOTE:  because ensemble perturbation bundle structure is same as control vector, use same ipc3d and
!             ipc2d indices for cvec and en_perts bundles.
    call gsi_bundlegetpointer (cvec,cvars3d,ipc3d,istatus)
    if(istatus/=0) then
      write(6,*) myname_,': cannot find 3d pointers'
      call stop2(999)
    endif
    call gsi_bundlegetpointer (cvec,cvars2d,ipc2d,istatus)
    if(istatus/=0) then
      write(6,*) myname_,': cannot find 2d pointers'
      call stop2(999)
    endif

    call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
    call gsi_bundlecreate (work_ens,grid_ens,'ensemble work',istatus, &
                              names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind)
    if(istatus/=0) then
       write(6,*)trim(myname_),': trouble creating work_ens bundle'
       call stop2(999)
    endif
    call gsi_gridcreate(grid_anl,grd_anl%lat2,grd_anl%lon2,grd_anl%nsig)
    call gsi_bundlecreate (work_anl,grid_anl,'analysis work',istatus, &
                              names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind)
    if(istatus/=0) then
       write(6,*)trim(myname_),': trouble creating work_anl bundle'
       call stop2(999)
    endif

    do ic3=1,nc3d
       work_anl%r3(ipc3d(ic3))%q=cvec%r3(ipc3d(ic3))%q
    end do
    do ic2=1,nc2d
       work_anl%r2(ipc2d(ic2))%q=cvec%r2(ipc2d(ic2))%q
    end do
    work_ens%values=zero
    call general_sube2suba_ad(grd_ens,grd_anl,p_e2a,work_ens%values,work_anl%values,regional)
    call gsi_bundledestroy(work_anl,istatus)
    if(istatus/=0) then
       write(6,*)trim(myname_),': trouble destroying work anl bundle'
       call stop2(999)
    endif

    ipe(1)=1
    im=a_en(1)%grid%im
    jm=a_en(1)%grid%jm
    km=a_en(1)%grid%km
       do n=1,n_ens
          do ic3=1,nc3d
             do k=1,km
                do j=1,jm
                   do i=1,im
                      a_en(n)%r3(ipe(1))%q(i,j,k)=a_en(n)%r3(ipe(1))%q(i,j,k) &
                               +work_ens%r3(ipc3d(ic3))%q(i,j,k)*en_perts(n,ibin)%r3(ipc3d(ic3))%qr4(i,j,k)
                   end do
                end do
             end do
          end do
          do ic2=1,nc2d

             select case (trim(cvars2d(ic2)))

                case('ps','PS')

                 if(pwgtflg)then
                   do k=1,km
                      do j=1,jm
                         do i=1,im
                            a_en(n)%r3(ipe(1))%q(i,j,k)=a_en(n)%r3(ipe(1))%q(i,j,k) &
                               +work_ens%r2(ipc2d(ic2))%q(i,j)*en_perts(n,ibin)%r2(ipc2d(ic2))%qr4(i,j)*pwgt(i,j,k)
                         end do
                      end do
                   end do
                 else
                   do j=1,jm
                      do i=1,im
                         a_en(n)%r3(ipe(1))%q(i,j,1)=a_en(n)%r3(ipe(1))%q(i,j,1) &
                            +work_ens%r2(ipc2d(ic2))%q(i,j)*en_perts(n,ibin)%r2(ipc2d(ic2))%qr4(i,j)
                      end do
                   end do
                 endif

                case('sst','SST')

                   do j=1,jm
                      do i=1,im
                         a_en(n)%r3(ipe(1))%q(i,j,1)=a_en(n)%r3(ipe(1))%q(i,j,1) &
                            +work_ens%r2(ipc2d(ic2))%q(i,j)*en_perts(n,ibin)%r2(ipc2d(ic2))%qr4(i,j)
                      end do
                   end do

             end select
          end do
       end do !n_ens
    call gsi_bundledestroy(work_ens,istatus)
    if(istatus/=0) then
       write(6,*)trim(myname_),': trouble destroying work ens bundle'
       call stop2(999)
    endif

  end subroutine ensemble_forward_model_ad_dual_res

  subroutine special_sd2h0
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    special_sd2h0  initialize subroutine special_sd2h
!   prgmmr: parrish          org: np22                date: 2009-06-16
!
! abstract: initialize subroutine special_sd2h (subdomain to slab for  
!             variable a_en).
!
! program history log:
!   2009-06-16  parrish
!   2010-02-10  parrish, correct allocate error on ndrecv_sd2h, found by Arthur Mizzi.
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

    use kinds, only: r_kind,i_kind
    use mpimod, only: npe,mype,mpi_comm_world,ierror,mpi_rtype
    use gridmod, only: nlat,nlon,nsig,nnnn1o,regional,vlevs
    use berror, only: nx,ny,nf
    implicit none

    integer(i_kind),dimension(0:npe-1):: nh_0_all,nh_1_all,nv_0_all,nv_1_all
    integer(i_kind) nvert,nh_tot,nh_this,nn,nv_tot,nv_this,kchk,n,kk,i,k
    real(r_kind),allocatable:: zsub(:,:),z(:)

!    set nval2f (loosely paraphrased from jfunc.f90)

    nscl=3          !  hard-wired here, later generalize when generalizing control variables
    if(regional) then
       nval2f=nlat*nlon
    else
       nval2f=ny*nx + 2*(2*nf+1)*(2*nf+1)
    end if


    allocate(nsend_sd2h(0:npe-1),ndsend_sd2h(0:npe),nrecv_sd2h(0:npe-1),ndrecv_sd2h(0:npe))
    allocate(i_recv(nval2f*nnnn1o),k_recv(nval2f*nnnn1o))
    nvert=vlevs

!  compute nv_0,nv_1

    nv_tot=nvert
    nv_this=nv_tot/npe
    if(mod(nv_tot,npe)/=0) nv_this=nv_this+1
    if(mod(nv_tot,npe)==0) then
       kchk=npe
    else
       kchk=mod(nv_tot,npe)
    end if

    nv_0_all=-1
    nv_1_all=-2
    nn=0
    do n=1,npe
       if(n<=kchk) then
          kk=nv_this
       else
          kk=nv_this-1
       end if
       if(kk>0) then
          nv_0_all(n-1)=nn+1
          nv_1_all(n-1)=nn+kk
       end if
       nn=nn+kk
    end do
    nv_0=nv_0_all(mype)
    nv_1=nv_1_all(mype)

!     compute nh_0, nh_1

    nh_tot=nval2f
    nh_this=nh_tot/npe
    if(mod(nh_tot,npe)/=0) nh_this=nh_this+1
    if(mod(nh_tot,npe)==0) then
       kchk=npe
    else
       kchk=mod(nh_tot,npe)
    end if

    nh_0_all=-1
    nh_1_all=-2
    nn=0
    do n=1,npe
       if(n<=kchk) then
          kk=nh_this
       else
          kk=nh_this-1
       end if
       if(kk>0) then
          nh_0_all(n-1)=nn+1
          nh_1_all(n-1)=nn+kk
       end if
       nn=nn+kk
    end do
    nh_0=nh_0_all(mype)
    nh_1=nh_1_all(mype)

!   compute nsend_sd2h,ndsend_sd2h,nrecv_sd2h,ndrecv_sd2h

    ndsend_sd2h(0)=0
    ndrecv_sd2h(0)=0
    do n=0,npe-1
       nsend_sd2h(n)=max(0,(nv_1_all(n)-nv_0_all(n)+1)*(nh_1-nh_0+1))
       ndsend_sd2h(n+1)=ndsend_sd2h(n)+nsend_sd2h(n)
       nrecv_sd2h(n)=max(0,(nv_1-nv_0+1)*(nh_1_all(n)-nh_0_all(n)+1))
       ndrecv_sd2h(n+1)=ndrecv_sd2h(n)+nrecv_sd2h(n)
    end do
    allocate(zsub(nh_0:nh_1,nvert),z(nval2f*(nv_1-nv_0+1)))
    do k=1,nvert
       do i=nh_0,nh_1
          zsub(i,k)=i
       end do
    end do
    call mpi_alltoallv(zsub,nsend_sd2h,ndsend_sd2h,mpi_rtype,&
                       z,nrecv_sd2h,ndrecv_sd2h,mpi_rtype,mpi_comm_world,ierror)
    do i=1,nval2f*(nv_1-nv_0+1)
       i_recv(i)=nint(z(i))
    end do

    do k=1,nvert
       do i=nh_0,nh_1
          zsub(i,k)=k
       end do
    end do
    call mpi_alltoallv(zsub,nsend_sd2h,ndsend_sd2h,mpi_rtype,&
                       z,nrecv_sd2h,ndrecv_sd2h,mpi_rtype,mpi_comm_world,ierror)
    do i=1,nval2f*(nv_1-nv_0+1)
       k_recv(i)=nint(z(i))
    end do

    deallocate(zsub,z)

  end subroutine special_sd2h0

  subroutine special_sd2h(zsub,z)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    special_sd2h  subdomain to slab for variable a_en
!   prgmmr: parrish          org: np22                date: 2009-06-16
!
! abstract: subdomain to slab for variable a_en.
!
! program history log:
!   2009-06-16  parrish
!
!   input argument list:
!     zsub     - input array on "subdomains"
!
!   output argument list:
!     z        - output array on slabs (form expected for input argument to ckgcov)
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind,i_kind
  use gridmod, only: nnnn1o,nsig,vlevs
  use constants, only: zero
  use mpimod, only: mype,mpi_rtype,ierror,mpi_comm_world
  implicit none

  real(r_kind),dimension(nh_0:nh_1,vlevs,nscl),intent(in   ) :: zsub
  real(r_kind),dimension(nval2f,nv_0:nv_1,nscl)         ,intent(  out) :: z

  real(r_kind) zsub1(nh_0:nh_1,vlevs),work(nval2f*(nv_1-nv_0+1))
  integer(i_kind) i,ii,is,k
! integer(i_kind) ibadp,ibadm,kbadp,kbadm
! logical good

!      1 <= nh_0 <= nh_1 <= nval2f

!      1 <= nv_0 <= nv_1 <= vlevs

  z=zero
  do is=1,nscl
     do k=1,vlevs
        do i=nh_0,nh_1
           zsub1(i,k)=zsub(i,k,is)
        end do
     end do
     call mpi_alltoallv(zsub1,nsend_sd2h,ndsend_sd2h,mpi_rtype,&
                        work,nrecv_sd2h,ndrecv_sd2h,mpi_rtype,mpi_comm_world,ierror)
     do ii=1,nval2f*(nv_1-nv_0+1)
        i=i_recv(ii) ; k=k_recv(ii)
        z(i,k,is)=work(ii)
     end do
  end do
  
end subroutine special_sd2h

subroutine beta12mult(grady)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    beta12mult  multiply grady by beta1inv and beta2inv        
!   prgmmr: parrish          org: np22                date: 2009-09-17
!
! abstract: when the hybrid ensemble option is turned on (l_hyb_ens=.true.)
!            the gradient vector grady contains two parts: the first is
!            the gradient with respect to the control variable associated
!            with the static background error covariance; the second is the
!            gradient with respect to the new ensemble control vector a_en.
!            the first is multiplied by beta1_inv, and the second by beta2_inv
!            where beta1_inv+beta2_inv=1.
!            adjusting beta1_inv between 0 and 1 allows tuning for optimal
!            blend between information provided by static background B and
!            ensemble based background.  beta1_inv=1 gives full weight to B
!            and betainv=0 gives full weight to ensemble.
!
! program history log:
!   2009-10-12  parrish  initial documentation
!   2010-03-29  kleist   comment out beta1_inv for SST
!   2010-04-28  todling  update to use gsi_bundle
!   2011-06-13  wu       used height dependent beta for regional
!
!   input argument list:
!     grady    - input field  grady_x1 : grady_a_en
!
!   output
!     grady    - beta1_inv*grady_x1 : beta2_inv*grady_a_en
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
  use kinds, only: r_kind,i_kind
  use gsi_4dvar, only: nsubwin
  use hybrid_ensemble_parameters, only: beta1_inv,n_ens,oz_univ_static
  use hybrid_ensemble_parameters, only: beta1wgt,beta2wgt,betaflg,grd_ens
  use constants, only:  one
  use control_vectors
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use timermod, only: timer_ini,timer_fnl

  use gridmod, only: nsig,regional,lat2,lon2

  implicit none

! Declare passed variables
  type(control_vector),intent(inout) :: grady

! Declare local variables
  character(len=*),parameter::myname_=myname//'*beta12mult'
  integer(i_kind) ii,nn,ic2,ic3
  integer(i_kind) i,j,k
  real(r_kind) beta2_inv
  integer(i_kind) ipc3d(nc3d),ipc2d(nc2d),istatus

! Initialize timer
  call timer_ini('beta12mult')
if(regional .and. betaflg) then
! Request CV pointers to vars pertinent to ensemble
  call gsi_bundlegetpointer ( grady%step(1), cvars3d, ipc3d, istatus )
  if(istatus/=0) then
     write(6,*) myname_,': cannot proceed, CV does not contain ens-required 3d fields'
     call stop2(999)
  endif
  call gsi_bundlegetpointer ( grady%step(1), cvars2d, ipc2d, istatus )
  if(istatus/=0) then
     write(6,*) myname_,': cannot proceed, CV does not contain ens-required 2d fields'
     call stop2(999)
  endif
  do ii=1,nsubwin

!    multiply by beta1_inv first:
     do ic3=1,nc3d
!  check for ozone and skip if oz_univ_static = true
        if((trim(cvars3d(ic3))=='oz'.or.trim(cvars3d(ic3))=='OZ').and.oz_univ_static) cycle
        do k=1,nsig
           do j=1,lon2
              do i=1,lat2
                 grady%step(ii)%r3(ipc3d(ic3))%q(i,j,k) =beta1wgt(k)*grady%step(ii)%r3(ipc3d(ic3))%q(i,j,k)
              enddo
           enddo
        enddo
     enddo
     do ic2=1,nc2d
! Default to static B estimate for SST
        if(trim(cvars2d(ic2))=='sst'.or.trim(cvars2d(ic2))=='SST') cycle
        do j=1,lon2
           do i=1,lat2
              grady%step(ii)%r2(ipc2d(ic2))%q(i,j) =beta1wgt(1)*grady%step(ii)%r2(ipc2d(ic2))%q(i,j)
           enddo
        enddo
     enddo

!    next multiply by beta2inv:
     do nn=1,n_ens
        do k=1,grd_ens%nsig
           do j=1,grd_ens%lon2
              do i=1,grd_ens%lat2
                 grady%aens(ii,nn)%r3(1)%q(i,j,k) =beta2wgt(k)*grady%aens(ii,nn)%r3(1)%q(i,j,k)
              enddo
           enddo
        enddo
     enddo


  end do
else
  beta2_inv=one-beta1_inv

! Request CV pointers to vars pertinent to ensemble
  call gsi_bundlegetpointer ( grady%step(1), cvars3d, ipc3d, istatus )
  if(istatus/=0) then
     write(6,*) myname_,': cannot proceed, CV does not contain ens-required 3d fields'
     call stop2(999)
  endif
  call gsi_bundlegetpointer ( grady%step(1), cvars2d, ipc2d, istatus )
  if(istatus/=0) then
     write(6,*) myname_,': cannot proceed, CV does not contain ens-required 2d fields'
     call stop2(999)
  endif
  do ii=1,nsubwin

!    multiply by beta1_inv first:
     do ic3=1,nc3d
!  check for ozone and skip if oz_univ_static = true
        if((trim(cvars3d(ic3))=='oz'.or.trim(cvars3d(ic3))=='OZ').and.oz_univ_static) cycle
        grady%step(ii)%r3(ipc3d(ic3))%q =beta1_inv*grady%step(ii)%r3(ipc3d(ic3))%q
     enddo
     do ic2=1,nc2d
! Default to static B estimate for SST
        if(trim(cvars2d(ic2))=='sst'.or.trim(cvars2d(ic2))=='SST') cycle 
        grady%step(ii)%r2(ipc2d(ic2))%q =beta1_inv*grady%step(ii)%r2(ipc2d(ic2))%q
     enddo

!    next multiply by beta2inv:
!$omp parallel do schedule(dynamic,1) private(nn)
     do nn=1,n_ens
        grady%aens(ii,nn)%r3(1)%q =beta2_inv*grady%aens(ii,nn)%r3(1)%q
     enddo
     
 
  end do
endif ! regional
  call timer_fnl('beta12mult')

  return
end subroutine beta12mult

subroutine init_sf_xy(jcap_in)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_sf_xy
!   prgmmr: parrish          org: np23                date: 2009-12-17
!
! abstract: initialize horizontal spectral localization sf_xy
!
! program history log:
!   2009-12-17  parrish
!   2010-06-29  parrish, modify so localization length can be different for each vertical level.
!                 to do this, add new variable array s_ens_hv(nsig), which is read in if readin_localization=.true.
!                 Otherwise, s_ens_hv is set equal to s_ens_h.
!
!   input argument list:
!     jcap_in - maximum spectral truncation allowed
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$

  use kinds, only: r_kind,i_kind,r_single
  use hybrid_ensemble_parameters,only: s_ens_hv,sp_loc,grd_ens,grd_loc
  use hybrid_ensemble_parameters,only: generate_ens,eqspace_ensgrid
  use general_specmod, only: general_init_spec_vars
  use constants, only: zero,half,one,two,three,rearth,pi
  use mpimod, only: mype
  implicit none

  integer(i_kind),intent(in   ) :: jcap_in

  integer(i_kind) i,ii,j,k,l,n,jcap
  real(r_kind),allocatable::g(:),gsave(:),errmax(:)
  real(r_kind) factor
  real(r_kind) rkm(grd_ens%nlat),f(grd_ens%nlat,grd_ens%nlon),f0(grd_ens%nlat,grd_ens%nlon)
  real(r_single) out1(grd_ens%nlon,grd_ens%nlat)
  real(r_single),allocatable::pn0_npole(:)
  real(r_kind) s_ens_h_min
  real(r_single) s_ens_hv4(grd_ens%nsig)
  real(r_kind) hmin,hmax
  character(5) mapname
  logical make_test_maps

  make_test_maps=.false.

!    make sure s_ens_hv is within allowable range  ( pi*rearth*.001/jcap_in <= s_ens_hv <= 5500 )

  s_ens_h_min=pi*rearth*.001_r_kind/jcap_in
  do k=1,grd_ens%nsig
     if(s_ens_hv(k) <  s_ens_h_min) then
        if(mype == 0) write(6,*)' s_ens_hv(',k,') = ',s_ens_hv(k),' km--too small, min value = ', &
                                        s_ens_h_min,' km.'
        if(mype == 0) write(6,*)' s_ens_hv(',k,') reset to min value'
        s_ens_hv(k)=s_ens_h_min
     else if(s_ens_hv(k) >  5500._r_kind) then
        if(mype == 0) write(6,*)' s_ens_hv(',k,') = ',s_ens_hv(k),' km--too large, max value = 5500 km.'
        if(mype == 0) write(6,*)' s_ens_hv(',k,') reset to max value'
        s_ens_hv(k)=5500._r_kind
     end if
  end do


  jcap=nint(1.2_r_kind*pi*rearth*.001_r_kind/minval(s_ens_hv))
  jcap=min(jcap,jcap_in)

!  set up spectral variables for jcap

  call general_init_spec_vars(sp_loc,jcap,jcap,grd_ens%nlat,grd_ens%nlon,eqspace_ensgrid)

!    the following code is used to compute the desired spectrum to get a
!     gaussian localization of desired length-scale.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   how spectrum is obtained:
!
!     Correlation matrix:  C = Y*D*Ytrans
!
!       where Y is matrix of spherical harmonics, evaluated on gaussian grid, and D is a diagonal matrix
!
!     To obtain D, exploit fact that for D a function only of total wave-number n, then C is homogeneous
!       and isotropic on the sphere.
!
!     So look at the special case of a test point centered on the north pole.  The correlation function
!       is then only a function of latitude, call it c(phi), where c(pi/2) = 1.
!
!     Now we have C = P*D*Ptrans, where we have reduced the problem to 1 dimension, latitude, and in
!       spectral space, total wave number n.  P is the zonal component only of Y.
!
!     Next, form the product
!                             C*e1 =P*D*Ptrans*e1,
!
!            where e1 is a vector of all 0, except for 1 at the north pole.
!
!     Then have P*D*Ptrans*e1 = sum(n) p(n,j)*d(n)*p(n,1) = c(j,1)
!
!        where j=1 corresponds to north pole point in this formulation.
!
!     Now if we have available C(j,1), a gaussian of desired length-scale evaluated (note, doesn't have to
!            be gaussian!) on the gaussian grid, then applying the inverse transform subroutine g2s0 to
!     C yields the product
!
!             Chat(n) = d(n)*p(n,1)
!
!     So finally the desired spectrum is
!
!               d(n) = chat(n)/p(n,1)
!
!     To create the full spectral transform of the desired correlation, d(n) is copied to all non-zero
!      zonal wave numbers on lines of constant total wave number n, multiplied by 1/2 to account for
!      two degrees of freedom for non-zero zonal wave numbers.
!
!     Note that while creating this routine, a bug was discovered in s2g0 routine for evaluation of pole
!       value.  There was a missing factor of 1/sqrt(2) apparently.
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!                 create reference gaussian centered on north pole with desired length-scale.

!     compute latitudes in km from north pole
  rkm(grd_ens%nlat)=zero
  rkm(1)=two*asin(one)*rearth*.001_r_kind
  do i=1,(grd_ens%nlat-2)/2
     rkm(grd_ens%nlat-i)=(asin(one)-asin(sp_loc%slat(i)))*rearth*.001_r_kind
     rkm(1+i)=(asin(one)+asin(sp_loc%slat(i)))*rearth*.001_r_kind
  end do
  if(mype == 0) write(6,*)' in init_sf_xy, lat,max(dlat)=', &
           rkm(1+(grd_ens%nlat-2)/2), &
          -rkm(grd_ens%nlat-(grd_ens%nlat-2)/2)+rkm(1+(grd_ens%nlat-2)/2),' km'

  allocate(spectral_filter(sp_loc%nc,grd_ens%nsig))
  allocate(sqrt_spectral_filter(sp_loc%nc,grd_ens%nsig))
  allocate(g(sp_loc%nc),gsave(sp_loc%nc))
  allocate(pn0_npole(0:sp_loc%jcap))
  do k=1,grd_ens%nsig
     do i=1,grd_ens%nlat
        f0(i,1)=exp(-half*(rkm(i)/s_ens_hv(k))**2)
     end do

     do j=2,grd_ens%nlon
        do i=1,grd_ens%nlat
           f0(i,j)=f0(i,1)
        end do
     end do

     call general_g2s0(grd_ens,sp_loc,g,f0)

     call general_s2g0(grd_ens,sp_loc,g,f)

!    adjust so value at np = 1
     f=f/f(grd_ens%nlat,1)
     call general_g2s0(grd_ens,sp_loc,g,f)
     call general_s2g0(grd_ens,sp_loc,g,f)
     if(mype == 0) write(6,*)' in init_sf_xy, jcap,s_ens_hv(',k,'), max diff(f0-f)=', &
                                        sp_loc%jcap,s_ens_hv(k),maxval(abs(f0-f))

!            correct spectrum by dividing by pn0_npole
     gsave=g

!    obtain pn0_npole
     do n=0,sp_loc%jcap
        g=zero
        g(2*n+1)=one
        call general_s2g0(grd_ens,sp_loc,g,f)
        pn0_npole(n)=f(grd_ens%nlat,1)
     end do

     g=zero
     do n=0,sp_loc%jcap
        g(2*n+1)=gsave(2*n+1)/pn0_npole(n)
     end do

!    obtain spectral_filter

     ii=0
     do l=0,sp_loc%jcap
        factor=one
        if(l >  0) factor=half
        do n=l,sp_loc%jcap
           ii=ii+1
           spectral_filter(ii,k)=factor*g(2*n+1)
           ii=ii+1
           if(l == 0) then
              spectral_filter(ii,k)=zero
           else
              spectral_filter(ii,k)=factor*g(2*n+1)
           end if
        end do
     end do
  end do
  sqrt_spectral_filter=sqrt(spectral_filter)
  deallocate(g,gsave,pn0_npole)

!  assign array k_index for each processor, based on grd_loc%kbegin_loc,grd_loc%kend_loc

  allocate(k_index(grd_loc%kbegin_loc:grd_loc%kend_alloc))
  k_index=0
  do k=grd_loc%kbegin_loc,grd_loc%kend_loc
     k_index(k)=1+mod(k-1,grd_loc%nsig)
!!     write(6,*) 'k_index(',k,')=',k_index(k)
  end do

  if(make_test_maps) then
   do k=1,grd_ens%nsig
      if(k>=grd_loc%kbegin_loc.and.k<=grd_loc%kend_loc) then
         f=zero
         f(grd_ens%nlat/2,grd_ens%nlon/2)=one
         call sf_xy(f,k,k)
         do j=1,grd_ens%nlon
            do i=1,grd_ens%nlat
               out1(j,i)=f(i,j)
            end do
         end do
         write(mapname,'("out",i2.2)')1+mod(k-1,grd_ens%nsig)
         call outgrads1(out1,grd_ens%nlon,grd_ens%nlat,mapname)
      end if
   end do
  end if

end subroutine init_sf_xy

subroutine sf_xy(f,k_start,k_end)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    sf_xy       spectral isotropic localization for global domain
!
!   prgrmmr:  parrish        org: np22                              date: 2009-09-30
!
! abstract:  use spherical harmonic transform to do horizontal isotropic localization of hybrid
!             ensemble control variable a_en.
!
! program history log:
!   2009-10-16  parrish  initial documentation
!   2010-03-11  parrish - adjust dimensions for f to allow for nlevs=0
!
!   input argument list:
!     f        - input field to be filtered
!     k_start  - starting horizontal slab index
!     k_end    - ending horizontal slab index    (k_end can be less than k_start, meaning there is
!                                                   no work on this processor)
!                 NOTE: above args allow horizontal localization length to vary in vertical
!
!   output argument list:
!     f        - filtered output
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  use hybrid_ensemble_parameters, only: grd_ens,sp_loc
  implicit none

  integer(i_kind),intent(in   ) :: k_start,k_end
  real(r_kind)   ,intent(inout) :: f(grd_ens%nlat*grd_ens%nlon,k_start:max(k_start,k_end))

  real(r_kind) g(sp_loc%nc)
  integer(i_kind) k

!  Cannot get this loop to thread properly on Zeus - remove threading for now.
!!$omp parallel do schedule(dynamic,1) private(k,g)
  do k=k_start,k_end
     call general_s2g0_ad(grd_ens,sp_loc,g,f(:,k))
     g(:)=g(:)*spectral_filter(:,k_index(k))
!    g(:)=g(:)*sqrt_spectral_filter(:,k_index(k))
!    g(:)=g(:)*sqrt_spectral_filter(:,k_index(k))
     call general_s2g0(grd_ens,sp_loc,g,f(:,k))
  end do

end subroutine sf_xy

subroutine sqrt_sf_xy(z,f,k_start,k_end)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    sqrt_sf_xy  sqrt(sf_xy)
!
!   prgrmmr:  parrish        org: np22                              date: 2011-06-28
!
! abstract:  sqrt of spectral localization operator sf_xy.
!
! program history log:
!   2011-06-28  parrish  initial documentation
!
!   input argument list:
!     z        - input spectral space variable
!     k_start  - starting horizontal slab index
!     k_end    - ending horizontal slab index    (k_end can be less than k_start, meaning there is
!                                                   no work on this processor)
!                 NOTE: above args allow horizontal localization length to vary in vertical
!
!   output argument list:
!     f        - output grid space variable
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  use hybrid_ensemble_parameters, only: grd_ens,sp_loc
  implicit none

  integer(i_kind),intent(in   ) :: k_start,k_end
  real(r_kind)   ,intent(in   ) :: z(sp_loc%nc,k_start:max(k_start,k_end))
  real(r_kind)   ,intent(  out) :: f(grd_ens%nlat*grd_ens%nlon,k_start:max(k_start,k_end))

  real(r_kind) g(sp_loc%nc)
  integer(i_kind) k

  do k=k_start,k_end
     g(:)=z(:,k)*sqrt_spectral_filter(:,k_index(k))
     call general_s2g0(grd_ens,sp_loc,g,f(:,k))
  end do

end subroutine sqrt_sf_xy

subroutine sqrt_sf_xy_ad(z,f,k_start,k_end)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    sqrt_sf_xy_ad  adjoint of sqrt_sf_xy
!
!   prgrmmr:  parrish        org: np22                              date: 2011-06-28
!
! abstract:  adjoint of sqrt_sf_xy.
!
! program history log:
!   2011-06-28  parrish  initial documentation
!
!   input argument list:
!     f        - grid space variable
!     k_start  - starting horizontal slab index
!     k_end    - ending horizontal slab index    (k_end can be less than k_start, meaning there is
!                                                   no work on this processor)
!                 NOTE: above args allow horizontal localization length to vary in vertical
!
!   output argument list:
!     z        - spectral space variable
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  use hybrid_ensemble_parameters, only: grd_ens,sp_loc
  implicit none

  integer(i_kind),intent(in   ) :: k_start,k_end
  real(r_kind)   ,intent(inout) :: z(sp_loc%nc,k_start:max(k_start,k_end))
  real(r_kind)   ,intent(inout) :: f(grd_ens%nlat*grd_ens%nlon,k_start:max(k_start,k_end))

  real(r_kind) g(sp_loc%nc)
  integer(i_kind) k

  do k=k_start,k_end
     call general_s2g0_ad(grd_ens,sp_loc,g,f(:,k))
     z(:,k)=g(:)*sqrt_spectral_filter(:,k_index(k))
  end do

end subroutine sqrt_sf_xy_ad

end module hybrid_ensemble_isotropic

subroutine get_new_alpha_beta(aspect,ng,fmat_out,fmat0_out)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    get_new_alpha_beta
!
!   prgrmmr:
!
! abstract:  compute various constants for new factorization Purser 1-d high-order filter.
!            adapted as simplification from new_alpha_betaa4 in raflib.f90 for use with
!            simple homogeneous isotropic localization filter.
!
! program history log:
!   2009-09-28  parrish  initial documentation
!
!   input argument list:
!     aspect   - squared correlation scale, in grid units squared
!     ng       - length of string
!     m        - filter order
!
!   output argument list:
!     fmat_out,fmat0_out - filter parameters for special factorization
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind,r_double
  use constants, only: one
  use raflib, only: stringop
  implicit none

  integer(i_kind)            , intent(in   ) :: ng
  real(r_kind), dimension(ng), intent(in   ) :: aspect
  real(r_kind)               , intent(  out) :: fmat_out(2,ng,2),fmat0_out(ng,2)

  integer(i_kind) i,j
  real(r_double) sig(0:ng-1),fmat(0:ng-1,-2:0,2)

  do i=1,ng
     sig(i-1)=sqrt(aspect(i))
  end do
  call stringop(ng-1,sig,fmat)

  do i=1,ng
     fmat_out(2,i,1)=fmat(i-1,-2,1)
     fmat_out(1,i,1)=fmat(i-1,-1,1)
     fmat0_out(i,1)=one/fmat(i-1,0,1)
     fmat_out(2,i,2)=fmat(i-1,-2,2)
     fmat_out(1,i,2)=fmat(i-1,-1,2)
     fmat0_out(i,2)=one/fmat(i-1,0,2)
  end do

end subroutine get_new_alpha_beta

subroutine bkerror_a_en(gradx,grady)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    bkerror_a_en  copy of bkerror for hybrid ensemble          
!   prgmmr: parrish          org: np22                date: 2009-09-17
!
! abstract: copy of bkerror for applying localization recursive filter
!            to hybrid ensemble control variable a_en.
!
! program history log:
!   2009-09-17  parrish  initial creation of code from a copy of bkerror
!   2010-05-20  todling  update to use bundle
!
!   input argument list:
!     gradx    - input field  
!
!   output
!     grady    - background structure * gradx 
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
  use kinds, only: r_kind,i_kind
  use gsi_4dvar, only: nsubwin, lsqrtb
  use control_vectors
  use timermod, only: timer_ini,timer_fnl
  use hybrid_ensemble_parameters, only: n_ens
  use hybrid_ensemble_parameters, only: nval_lenz_en
  use gsi_bundlemod,only: gsi_bundlegetpointer
  implicit none

! Declare passed variables
  type(control_vector),intent(inout) :: gradx
  type(control_vector),intent(inout) :: grady

! Declare local variables
  integer(i_kind) ii,nn,ip,istatus
  real(r_kind),allocatable::z(:)

  if (lsqrtb) then
     write(6,*)'bkerror_a_en: not for use with lsqrtb'
     call stop2(317)
  end if

! Initialize timer
  call timer_ini('bkerror_a_en')

! Put things in grady first since operations change input variables
  call gsi_bundlegetpointer ( grady%aens(1,1),'a_en',ip,istatus)
  if(istatus/=0) then
     write(6,*)'bkerror_a_en: trouble getting pointer to ensemble CV'
     call stop2(317)
  endif
  do ii=1,nsubwin
!$omp parallel do schedule(dynamic,1) private(nn)
     do nn=1,n_ens
        grady%aens(ii,nn)%r3(ip)%q=gradx%aens(ii,nn)%r3(ip)%q
     enddo
  end do

! Apply variances, as well as vertical & horizontal parts of background error
  do ii=1,nsubwin
    !if(test_sqrt_localization) then
    !        write(6,*)' using ckgcov_a_en_new_factorization'
    !   allocate(z(nval_lenz_en))
    !   call ckgcov_a_en_new_factorization_ad(z,grady%aens(ii,1:n_ens))
    !   call ckgcov_a_en_new_factorization   (z,grady%aens(ii,1:n_ens))
    !   deallocate(z)
    !else
    !        write(6,*)' using bkgcov_a_en_new_factorization'
        call bkgcov_a_en_new_factorization(grady%aens(ii,1:n_ens))
    !end if
  end do

! Finalize timer
  call timer_fnl('bkerror_a_en')

  return
end subroutine bkerror_a_en

subroutine bkgcov_a_en_new_factorization(a_en)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    bkgcov_a_en copy of bkgcov for hybrid ens var a_en 
!   prgmmr: parrish        org: np22                date: 2009-09-17
!
! abstract: copy of bkgcov to apply localization with recursive filters
!            to hybrid ensemble control variable a_en.
!
! program history log:
!   2009-09-17  parrish
!   2010-02-20  parrish, adapt for dual resolution
!   2010-05-21  todling, update to use bundle
!   2010-07-02  parrish, modify arguments to call sf_xy to allow for vertical variation
!                 of horizontal localization length
!
!   input argument list:
!     a_en     - control variable for ensemble contribution to background error
!     nlevs    - number of vertical levels for smoothing
!
!   output argument list:
!                 all after smoothing, combining scales
!     a_en     - control variable for ensemble contribution to background error
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: regional
  use hybrid_ensemble_parameters, only: n_ens,grd_loc
  use hybrid_ensemble_isotropic, only: &
         new_factorization_rf_z,new_factorization_rf_x,new_factorization_rf_y,sf_xy
  use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer

  implicit none

! Passed Variables
! real(r_kind),dimension(grd_loc%latlon1n,n_ens),intent(inout) :: a_en
  type(gsi_bundle),intent(inout) :: a_en(n_ens)

! Local Variables
  integer(i_kind) i,ii,j,k,n,iflg,iadvance,iback,is,ie,ipnt,istatus
  real(r_kind) hwork(grd_loc%inner_vars,grd_loc%nlat,grd_loc%nlon,grd_loc%kbegin_loc:grd_loc%kend_alloc)
  real(r_kind),allocatable,dimension(:):: a_en_work

  iflg=1

  call gsi_bundlegetpointer(a_en(1),'a_en',ipnt,istatus)
  if(istatus/=0) then
     write(6,*)'bkgcov_a_en_new_factorization: trouble getting pointer to ensemble CV'
     call stop2(999)
  endif

! Apply vertical smoother on each ensemble member
! To avoid my having to touch the general sub2grid and grid2sub,
! get copy for ensemble components to work array
  allocate(a_en_work(n_ens*a_en(1)%ndim),stat=istatus)
  if(istatus/=0) then
     write(6,*)'bkgcov_a_en_new_factorization: trouble in alloc(a_en_work)'
     call stop2(999)
  endif
  iadvance=1 ; iback=2
!$omp parallel do schedule(dynamic,1) private(k,ii,is,ie)
  do k=1,n_ens
     call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback)
     ii=(k-1)*a_en(1)%ndim
     is=ii+1
     ie=ii+a_en(1)%ndim
     a_en_work(is:ie)=a_en(k)%values(1:a_en(k)%ndim)
  enddo

! Convert from subdomain to full horizontal field distributed among processors
  call general_sub2grid(grd_loc,a_en_work,hwork)

! Apply horizontal smoother for number of horizontal scales
  if(regional) then
     iadvance=1 ; iback=2
     call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
     call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
     iadvance=2 ; iback=1
     call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
     call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
  else
     call sf_xy(hwork,grd_loc%kbegin_loc,grd_loc%kend_loc)
  end if

! Put back onto subdomains
  call general_grid2sub(grd_loc,hwork,a_en_work)

! Retrieve ensemble components from long vector
! Apply vertical smoother on each ensemble member
  iadvance=2 ; iback=1
!$omp parallel do schedule(dynamic,1) private(k,ii,is,ie)
  do k=1,n_ens
     ii=(k-1)*a_en(1)%ndim
     is=ii+1
     ie=ii+a_en(1)%ndim
     a_en(k)%values(1:a_en(k)%ndim)=a_en_work(is:ie)
     call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback)
  enddo
  deallocate(a_en_work)

  return
end subroutine bkgcov_a_en_new_factorization

subroutine ckgcov_a_en_new_factorization(z,a_en)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ckgcov_a_en_new_factorization sqrt(bkgcov_a_en_new_factorization)
!   prgmmr: parrish        org: np22                date: 2011-06-27
!
! abstract: make a copy of bkgcov_a_en_new_factorization and form sqrt.
!
! program history log:
!   2011-06-27  parrish, initial documentation
!
!   input argument list:
!     z        - long vector containing sqrt control vector for ensemble extended control variable
!
!   output argument list:
!     a_en     - bundle containing intermediate control variable after multiplication by sqrt(S), the
!                    ensemble localization correlation.
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only: zero
  use gridmod, only: regional
  use hybrid_ensemble_parameters, only: n_ens,grd_loc
  use hybrid_ensemble_isotropic, only: &
         new_factorization_rf_z,new_factorization_rf_x,new_factorization_rf_y
#ifdef LATER
  use hybrid_ensemble_isotropic, only: sqrt_sf_xy
#endif /*LATER*/
  use hybrid_ensemble_parameters, only: nval_lenz_en
  use general_sub2grid_mod, only: general_grid2sub
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer

  implicit none

! Passed Variables
  type(gsi_bundle),intent(inout) :: a_en(n_ens)
  real(r_kind),dimension(nval_lenz_en),intent(in   ) :: z

! Local Variables
  integer(i_kind) i,ii,k,iadvance,iback,is,ie,ipnt,istatus
  real(r_kind) hwork(grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1))
!NOTE:   nval_lenz_en = nhoriz*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1)
!      and nhoriz = grd_loc%nlat*grd_loc%nlon for regional,
!          nhoriz = (sp_loc%jcap+1)*(sp_loc%jcap+2) for global
!   but internal array hwork always has
!      dimension grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1)
!    which just happens to match up with nval_lenz_en for regional case, but not global.
  real(r_kind),allocatable,dimension(:):: a_en_work

  call gsi_bundlegetpointer(a_en(1),'a_en',ipnt,istatus)
  if(istatus/=0) then
     write(6,*)'bkgcov_a_en_new_factorization: trouble getting pointer to ensemble CV'
     call stop2(999)
  endif


  if(grd_loc%kend_loc+1-grd_loc%kbegin_loc==0) then
!     no work to be done on this processor, but hwork still has allocated space, since
!                     grd_loc%kend_alloc = grd_loc%kbegin_loc in this case, so set to zero.
     hwork=zero
  else
! Apply horizontal smoother for number of horizontal scales
     if(regional) then
! Make a copy of input variable z to hwork
        hwork=z
        iadvance=2 ; iback=1
        call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
        call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
     else
#ifdef LATER
        call sqrt_sf_xy(z,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc)
#else
        write(6,*) ' problem with ibm compiler with "use hybrid_ensemble_isotropic, only: sqrt_sf_xy"'
#endif /*LATER*/
     end if
  end if

! Put back onto subdomains
  allocate(a_en_work(n_ens*a_en(1)%ndim),stat=istatus)
  if(istatus/=0) then
     write(6,*)'bkgcov_a_en_new_factorization: trouble in alloc(a_en_work)'
     call stop2(999)
  endif
  call general_grid2sub(grd_loc,hwork,a_en_work)

! Retrieve ensemble components from long vector
  ii=0
  do k=1,n_ens
     is=ii+1
     ie=ii+a_en(1)%ndim
     a_en(k)%values(1:a_en(k)%ndim)=a_en_work(is:ie)
     ii=ii+a_en(1)%ndim
  enddo
  deallocate(a_en_work)

! Apply vertical smoother on each ensemble member
  do k=1,n_ens

     iadvance=2 ; iback=1
     call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback)

  end do

  return
end subroutine ckgcov_a_en_new_factorization

subroutine ckgcov_a_en_new_factorization_ad(z,a_en)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ckgcov_a_en_new_factorization_ad adjoint of ckgcov_a_en_new_factorization
!   prgmmr: parrish        org: np22                date: 2011-06-27
!
! abstract: adjoint of ckgcov_a_en_new_factorization.  Calling ckgcov_a_en_new_factorization_ad,
!            followed by ckgcov_a_en_new_factorization is the equivalent of one call to
!            subroutine bkgcov_a_en_new_factorization.
!
! program history log:
!   2011-06-27  parrish, initial documentation
!
!   input argument list:
!     z        - long vector containing sqrt control vector for ensemble extended control variable
!     a_en     - bundle containing intermediate control variable after multiplication by sqrt(S), the
!                    ensemble localization correlation.
!
!   output argument list:
!     z        - long vector containing sqrt control vector for ensemble extended control variable
!     a_en     - bundle containing intermediate control variable after multiplication by sqrt(S), the
!                    ensemble localization correlation.
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only: zero
  use gridmod, only: regional
  use hybrid_ensemble_parameters, only: n_ens,grd_loc
  use hybrid_ensemble_isotropic, only: &
         new_factorization_rf_z,new_factorization_rf_x,new_factorization_rf_y,sqrt_sf_xy_ad
  use hybrid_ensemble_parameters, only: nval_lenz_en
  use general_sub2grid_mod, only: general_sub2grid
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer

  implicit none

! Passed Variables
  type(gsi_bundle),intent(inout) :: a_en(n_ens)
  real(r_kind),dimension(nval_lenz_en),intent(inout) :: z

! Local Variables
  integer(i_kind) i,ii,k,iadvance,iback,is,ie,ipnt,istatus
  real(r_kind) hwork(grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1))
!NOTE:   nval_lenz_en = nhoriz*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1)
!      and nhoriz = grd_loc%nlat*grd_loc%nlon for regional,
!          nhoriz = (sp_loc%jcap+1)*(sp_loc%jcap+2) for global
!   but internal array hwork always has
!      dimension grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1)
!    which just happens to match up with nval_lenz_en for regional case, but not global.
  real(r_kind),allocatable,dimension(:):: a_en_work

  call gsi_bundlegetpointer(a_en(1),'a_en',ipnt,istatus)
  if(istatus/=0) then
     write(6,*)'bkgcov_a_en_new_factorization: trouble getting pointer to ensemble CV'
     call stop2(999)
  endif

! Apply vertical smoother on each ensemble member
  do k=1,n_ens

     iadvance=1 ; iback=2
     call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback)
 
  end do

! To avoid my having to touch the general sub2grid and grid2sub,
! get copy for ensemble components to work array
  allocate(a_en_work(n_ens*a_en(1)%ndim),stat=istatus)
  if(istatus/=0) then
     write(6,*)'bkgcov_a_en_new_factorization: trouble in alloc(a_en_work)'
     call stop2(999)
  endif
  ii=0
  do k=1,n_ens
     is=ii+1
     ie=ii+a_en(1)%ndim
     a_en_work(is:ie)=a_en(k)%values(1:a_en(k)%ndim)
     ii=ii+a_en(1)%ndim
  enddo

! Convert from subdomain to full horizontal field distributed among processors
  call general_sub2grid(grd_loc,a_en_work,hwork)
  deallocate(a_en_work)

  if(grd_loc%kend_loc+1-grd_loc%kbegin_loc==0) then
!     no work to be done on this processor, but z still has allocated space, since
!                     grd_loc%kend_alloc = grd_loc%kbegin_loc in this case, so set to zero.
     z=zero
  else
! Apply horizontal smoother for number of horizontal scales
     if(regional) then
        iadvance=1 ; iback=2
        call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
        call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
        z=hwork
     else
        call sqrt_sf_xy_ad(z,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc)
     end if
  end if

  return
end subroutine ckgcov_a_en_new_factorization_ad

! ------------------------------------------------------------------------------
! ------------------------------------------------------------------------------

subroutine hybens_grid_setup
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    hybrid_ensemble_setup  initialize everything for hybrid ensemble
!   prgmmr: parrish          org: np22                date: 2009-09-17
!
! abstract: put everything for initializing hybrid ensemble in one subroutine.
!
! program history log:
!   2009-09-17  parrish
!   2010-02-20  parrish, adapt for dual resolution
!   2011-01-30  parrish, fix so regional application depends only on parameters regional
!                  and dual_res.  Rename subroutine get_regional_gefs_grid to get_regional_dual_res_grid.
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind
  use hybrid_ensemble_parameters, only: aniso_a_en,generate_ens,n_ens,&
                      s_ens_h,nlon_ens,nlat_ens,jcap_ens,jcap_ens_test,&
                      grd_ens,grd_loc,grd_a1,grd_e1,grd_anl,sp_ens,p_e2a,&
                      dual_res,uv_hyb_ens,eqspace_ensgrid,grid_ratio_ens
  use gridmod,only: regional,nsig,nlon,nlat,rlats,rlons
  use general_sub2grid_mod, only: general_sub2grid_create_info
  use general_specmod, only: general_init_spec_vars
  use egrid2agrid_mod,only: g_create_egrid2agrid,create_egrid2agrid
  use mpimod, only: mype,ierror,npe
  use constants, only: zero,one
  use hybrid_ensemble_isotropic, only: region_lat_ens,region_lon_ens
  use hybrid_ensemble_isotropic, only: region_dx_ens,region_dy_ens
  use control_vectors, only: cvars3d,nc2d,nc3d
  use gridmod, only: region_lat,region_lon,region_dx,region_dy

  implicit none

  integer(i_kind) inner_vars,num_fields,ic3,k
  integer(i_kind) nord_e2a
  logical,allocatable::vector(:)
  real(r_kind) eps,r_e
  real(r_kind) rlon_a(nlon),rlat_a(nlat),rlon_e(nlon),rlat_e(nlat)

             nord_e2a=4       !   soon, move this to hybrid_ensemble_parameters

  if(aniso_a_en) then
     if(mype == 0) write(6,*)' anisotropic option not available yet for hybrid ensemble localization'
     if(mype >  -10) then
        call stop2(999)
        stop
     end if
  end if

  eps=zero
  r_e=grid_ratio_ens

! Initialize dual_res logical
  if (.not.regional) then
    dual_res = (nlon /= nlon_ens .or. nlat /= nlat_ens)
  else
    dual_res=.false.
  end if
  if(mype==0) write(6,*)' before compute nlat_ens,nlon_ens, nlat,nlon,nlat_ens,nlon_ens,r_e,eps=',&
                          nlat,nlon,nlat_ens,nlon_ens,r_e,eps

!  if regional set up for possible dual-res application:

  if(regional) then
     call get_regional_dual_res_grid(eps,r_e,nlon,nlon_ens,rlon_a,rlon_e)
     call get_regional_dual_res_grid(eps,r_e,nlat,nlat_ens,rlat_a,rlat_e)
     call create_egrid2agrid(nlat,rlat_a,nlon,rlon_a,nlat_ens,rlat_e,nlon_ens,rlon_e,nord_e2a,p_e2a)
     dual_res=.not.p_e2a%identity
!  NOTE:  key dual-res on egrid2agrid parameter p_e2a%identity, not nlat_ens==nlat .or. nlon_ens==nlon
     allocate(region_lat_ens(nlat_ens,nlon_ens))
     allocate(region_lon_ens(nlat_ens,nlon_ens))
     allocate(region_dx_ens(nlat_ens,nlon_ens))
     allocate(region_dy_ens(nlat_ens,nlon_ens))
     if(dual_res) then
        call get_region_lat_lon_ens(region_lat_ens,region_lon_ens, &
                                       rlat_e,rlon_e,nlat_ens,nlon_ens)
     else
        region_lon_ens=region_lon
        region_lat_ens=region_lat
     end if
  end if
                 if(mype==0) write(6,*)' dual_res,nlat,nlon,nlat_ens,nlon_ens,r_e,eps=',&
                                                     dual_res,nlat,nlon,nlat_ens,nlon_ens,r_e,eps

  if(nlon_ens<=0 .or. nlat_ens<=0) then
     nlon_ens=nlon ; nlat_ens=nlat
  end if

!     2.  create grid info for ensemble, including stuff for ensemble general_sub2grid, general_grid2sub
  num_fields=nsig*n_ens
  inner_vars=1
  allocate(vector(num_fields))
  vector=.false.
  call general_sub2grid_create_info(grd_loc,inner_vars,nlat_ens,nlon_ens,nsig,num_fields,regional,vector)
  num_fields=max(0,nc3d)*nsig+max(0,nc2d)
  deallocate(vector)
  allocate(vector(num_fields))
  vector=.false.
  do ic3=1,nc3d
     if(trim(cvars3d(ic3))=='sf'.or.trim(cvars3d(ic3))=='vp') then
        do k=1,nsig
           vector((ic3-1)*nsig+k)=uv_hyb_ens
        end do
     end if
  end do
  call general_sub2grid_create_info(grd_ens,inner_vars,nlat_ens,nlon_ens,nsig,num_fields,regional,vector)
  call general_sub2grid_create_info(grd_anl,inner_vars,nlat,nlon,nsig,num_fields,regional,vector)
  deallocate(vector)
  num_fields=nsig
  allocate(vector(num_fields))
  vector=.false.
  call general_sub2grid_create_info(grd_e1,inner_vars,nlat_ens,nlon_ens,nsig,num_fields,regional,vector)
  call general_sub2grid_create_info(grd_a1,inner_vars,nlat,nlon,nsig,num_fields,regional,vector)
  deallocate(vector)

!    setup dual-resolution feature, if needed

  if(.not.regional) then
     call general_init_spec_vars(sp_ens,jcap_ens,jcap_ens_test,grd_ens%nlat,grd_ens%nlon,eqspace_ensgrid)
     call g_create_egrid2agrid(nlat,rlats,nlon,rlons,grd_ens%nlat,sp_ens%rlats,grd_ens%nlon,sp_ens%rlons, &
                               nord_e2a,p_e2a)
  else
     if(dual_res) then
        call get_region_dx_dy_ens(region_dx_ens,region_dy_ens)
     else
        region_dx_ens=region_dx
        region_dy_ens=region_dy
     end if
  end if

  return
end subroutine hybens_grid_setup

subroutine hybens_localization_setup
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    hybens_localization_setup   setup for hybrid localization
!   prgmmr: kleist          org: np22                date: 2010-07-30
!
! abstract: setup arrays used for control variable localization
!
! program history log:
!   2010-07-30  kleist
!   2011-10-03  wu - add call to setup_ens_pwgt, which computes vertical weighting for ensemble contribution
!                     to psfc.
!   2012-10-16  wu - only call setup_ens_pwgt if necessary
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind
  use hybrid_ensemble_parameters, only: grd_ens,s_ens_v,jcap_ens,s_ens_vv,&
	 s_ens_h,s_ens_hv,create_hybens_localization_parameters,grd_loc,sp_loc,&
         readin_localization,nval_lenz_en,pwgtflg,betaflg
  use gridmod,only: regional
  use hybrid_ensemble_isotropic, only: init_rf_x,init_rf_y,init_rf_z,&
         normal_new_factorization_rf_z,normal_new_factorization_rf_x,&
         normal_new_factorization_rf_y,init_sf_xy
  use constants, only: one,zero
  use mpimod, only: mype
  implicit none

  character(len=40)  :: fname = 'hybens_locinfo'
  integer(i_kind) k,msig,istat
  integer(i_kind)    :: lunin = 47
  logical            :: lexist
  real(r_kind) s_ens_h_gu_x,s_ens_h_gu_y

! Allocate
  call create_hybens_localization_parameters

! Set up localization parameters as function of level

! if horizontal parameter is set <= 0, read in k-levels of localization parameters
  if (readin_localization) then

!   Check the status of input file
    inquire(file=trim(fname),exist=lexist)
    if ( lexist ) then
       open(lunin,file=trim(fname),form='formatted')
       rewind(lunin)
       read(lunin,100,iostat=istat) msig
       if ( msig /= grd_ens%nsig ) then
          write(6,*) 'HYBENS_LOCALIZATION_SETUP:  ***ERROR*** error in ',trim(fname)
          write(6,*) 'HYBENS_LOCALIZATION_SETUP:  levels do not match,msig[read in],nsig[defined] = ',msig,grd_ens%nsig
          close(lunin)
          call stop2(123)
       endif
       do k=1,grd_ens%nsig
         read(lunin,101) s_ens_hv(k),s_ens_vv(k)
       end do
       close(lunin)
    else 
      write(6,*) 'HYBENS_LOCALIZATION_SETUP:  ***ERROR*** INPUT FILE MISSING -- ',trim(fname)
    end if 
 100 format(I4)
 101 format(F8.1,3x,F5.1)

  else
!          assign all levels to same value, s_ens_h  (ran with this on 20100702 and reproduced results from
!                                                      rungsi62_hyb_dualres.sh)
     s_ens_hv=s_ens_h
     s_ens_vv=s_ens_v
  end if

! Set up localization filters

  call init_rf_z(s_ens_vv)
  call normal_new_factorization_rf_z

  if(regional) then
!     convert s_ens_h from km to grid units.
     call convert_km_to_grid_units(s_ens_h,s_ens_h_gu_x,s_ens_h_gu_y)
     call init_rf_x(s_ens_h_gu_x)
     call normal_new_factorization_rf_x
     call init_rf_y(s_ens_h_gu_y)
     call normal_new_factorization_rf_y
  else
     call init_sf_xy(jcap_ens)
  end if

  if (pwgtflg .or. betaflg) call setup_ens_pwgt

!  set value of nval_lenz_en here for now, but will need to rearrange so this can be set in control_vectors
!     and triggered by lsqrtb.
  if(regional) then
     nval_lenz_en=grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1)
  else
     nval_lenz_en=sp_loc%nc*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1)
  end if

end subroutine hybens_localization_setup

subroutine convert_km_to_grid_units(s_ens_h,s_ens_h_gu_x,s_ens_h_gu_y)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    convert_km_to_grid_units
!   prgmmr: parrish          org: np22                date: 2009-09-17
!
! abstract: change horizontal localization length scales from km to grid units.
!
! program history log:
!   2009-09-17  parrish
!
!   input argument list:
!     s_ens_h  - input localization length scale in km
!
!   output argument list:
!     s_ens_h_gu_x - output x localization length in grid units
!     s_ens_h_gu_y - output y localization length in grid units
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind,i_kind
  use hybrid_ensemble_isotropic,only: region_dx_ens,region_dy_ens
  implicit none

  real(r_kind),intent(in   ) ::s_ens_h
  real(r_kind),intent(  out) ::s_ens_h_gu_x,s_ens_h_gu_y
  logical,parameter:: debug=.false.

  if(debug) write(6,*)' in convert_km_to_grid_units, min, max region_dx_ens*.001=',&
                    .001_r_kind*minval(region_dx_ens),.001_r_kind*maxval(region_dx_ens)
  s_ens_h_gu_x=s_ens_h/(.001_r_kind*maxval(region_dx_ens))
  if(debug) write(6,*)' in convert_km_to_grid_units, min, max region_dy_ens*.001=',&
                    .001_r_kind*minval(region_dy_ens),.001_r_kind*maxval(region_dy_ens)
  s_ens_h_gu_y=s_ens_h/(.001_r_kind*maxval(region_dy_ens))

  if(debug) write(6,*)' in convert_km_to_grid_units, s_ens_h,s_ens_h_gu_x,y=', &
                    s_ens_h,s_ens_h_gu_x,s_ens_h_gu_y

end subroutine convert_km_to_grid_units

subroutine grads1(f,nvert,mype,fname)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grads1      generate grads output file for f on subdomains
!   prgmmr: parrish          org: np22                date: 2009-09-17
!
! abstract: create grads output file for variable f on subdomains.  used
!            to visualize the field f using grads software.
!
! program history log:
!   2009-09-17  parrish
!
!   input argument list:
!     f     - field to generate grads output file
!     nvert - number of vertical levels in f
!     mype  - local processor
!     fname - character string used to identify field f in grads output file
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_single,r_kind,i_kind
  use constants, only: one
  use gridmod, only: nlat,nlon,lon2,lat2
  implicit none

  integer(i_kind),intent(in   ) :: nvert,mype
  character(*)   ,intent(in   ) :: fname
  real(r_kind)   ,intent(in   ) :: f(lat2,lon2,nvert)

  real(r_kind),dimension(nlat,nlon)::work
  real(r_single) outfield(nlon,nlat)

  character(50) dsname,title,filename
! data dsname/'test.dat'/
  data title/'inmi'/
  character(112) datdes(50000)
  character(1) blank
  data blank/' '/
  data undef/-9.99e33_r_single/

  integer(i_kind) i,k,kend,kstart,next,np,ioutdes,ioutdat
  integer(i_kind) last,j,koutmax
  real(r_single) undef
  real(r_single) startp,pinc

  if(mype == 0) then
     np=nvert
     startp=1._r_single
     pinc=1._r_single
     ioutdes=98550
     ioutdat=98551
     write(filename,'(a,".des")')trim(fname)
     write(dsname,'(a,".dat")')trim(fname)
     open(unit=ioutdes,file=trim(filename),form='formatted')
     open(unit=ioutdat,file=trim(dsname),form='unformatted')
     rewind ioutdes
     rewind ioutdat
     do i=1,50000
        write(datdes(i),'(112a1)')(blank,k=1,112)
     end do
     write(datdes(1),'("DSET ",a50)')dsname
     write(datdes(2),'("options big_endian sequential")')
     write(datdes(3),'("TITLE ",a50)')title
     write(datdes(4),'("UNDEF ",e11.2)')undef
     next=5
     write(datdes(next),'("XDEF ",i5," LINEAR ",f7.2,f7.2)')nlon,startp,pinc
     next=next+1
     write(datdes(next),'("YDEF ",i5," LINEAR ",f7.2,f7.2)')nlat,startp,pinc
     next=next+1
     write(datdes(next),'("ZDEF ",i5," LINEAR ",f7.2,f7.2)')np,startp,pinc
     next=next+1
     koutmax=1
     write(datdes(next),'("TDEF ",i5," LINEAR 0Z23may1992 24hr")')koutmax
     next=next+1
     write(datdes(next),'("VARS 1")')
     next=next+1
     write(datdes(next),'("f   ",i5," 99 f   ")')nvert
     next=next+1
     write(datdes(next),'("ENDVARS")')
     last=next
     write(ioutdes,'(a112)')(datdes(i),i=1,last)
 
  end if

  do k=1,nvert
     call sub2grid_1(f(1,1,k),work,0,mype)
     if(mype == 0) then
        do j=1,nlon ; do i=1,nlat
           outfield(j,i)=work(i,j)
        end do ; end do
        write(ioutdat)outfield
     end if
  end do

  if(mype == 0) then
     close(ioutdes)
     close(ioutdat)
  end if

end subroutine grads1

subroutine sub2grid_1(sub,grid,gridpe,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grads1      generate grads output file for f on subdomains
!   prgmmr: parrish          org: np22                date: 2009-09-17
!
! abstract: straightforward, but inefficient code to convert a single variable
!            on subdomains to complete slab on one processor.
!
! program history log:
!   2009-09-17  parrish
!   2010-04-01  treadon - move strip to gridmod
!
!   input argument list:
!     sub   - field on subdomains
!     gridpe- processor that contains full slab representation of sub
!     mype  - local processor
!
!   output argument list:
!     grid  - field on complete slab
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$


  use kinds, only: r_kind,i_kind
  use constants, only: zero
  use gridmod, only: nlat,nlon,lat2,lon2,lat1,lon1,&
         ltosi,ltosj,iglobal,ijn,displs_g,itotsub,strip
  use mpimod, only: mpi_comm_world,ierror,mpi_rtype
  implicit none

  integer(i_kind)                  ,intent(in   ) :: gridpe,mype
  real(r_kind),dimension(lat2,lon2),intent(in   ) :: sub
  real(r_kind),dimension(nlat,nlon),intent(  out) :: grid

  real(r_kind),dimension(lat1*lon1):: zsm
  real(r_kind),dimension(itotsub):: work1
  integer(i_kind) mm1,i,j,k

  mm1=mype+1

  do j=1,lon1*lat1
     zsm(j)=zero
  end do
  call strip(sub,zsm,1)
  call mpi_gatherv(zsm,ijn(mm1),mpi_rtype, &
                 work1,ijn,displs_g,mpi_rtype, &
                 gridpe,mpi_comm_world,ierror)
  if(mype == gridpe) then
     do k=1,iglobal
        i=ltosi(k) ; j=ltosj(k)
        grid(i,j)=work1(k)
     end do
  end if

end subroutine sub2grid_1

subroutine grads1_ens(f,nvert,mype,fname)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grads1      generate grads output file for f on subdomains
!   prgmmr: parrish          org: np22                date: 2009-09-17
!
! abstract: create grads output file for variable f on subdomains.  used
!            to visualize the field f using grads software.
!
! program history log:
!   2009-09-17  parrish
!
!   input argument list:
!     f     - field to generate grads output file
!     nvert - number of vertical levels in f
!     mype  - local processor
!     fname - character string used to identify field f in grads output file
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_single,r_kind,i_kind
  use constants, only: one
  use gridmod, only: nlat,nlon,lon2,lat2
  use hybrid_ensemble_parameters, only: grd_ens
  implicit none

  integer(i_kind),intent(in   ) :: nvert,mype
  character(*)   ,intent(in   ) :: fname
  real(r_kind)   ,intent(in   ) :: f(lat2,lon2,nvert)

  real(r_kind),dimension(nlat,nlon)::work
  real(r_single) outfield(nlon,nlat)

  character(50) dsname,title,filename
! data dsname/'test.dat'/
  data title/'inmi'/
  character(112) datdes(50000)
  character(1) blank
  data blank/' '/
  data undef/-9.99e33_r_single/

  integer(i_kind) i,k,kend,kstart,next,np,ioutdes,ioutdat
  integer(i_kind) last,j,koutmax
  real(r_single) undef
  real(r_single) startp,pinc

  if(mype == 0) then
     np=nvert
     startp=1._r_single
     pinc=1._r_single
     ioutdes=98550
     ioutdat=98551
     write(filename,'(a,".des")')trim(fname)
     write(dsname,'(a,".dat")')trim(fname)
     open(unit=ioutdes,file=trim(filename),form='formatted')
     open(unit=ioutdat,file=trim(dsname),form='unformatted')
     rewind ioutdes
     rewind ioutdat
     do i=1,50000
        write(datdes(i),'(112a1)')(blank,k=1,112)
     end do
     write(datdes(1),'("DSET ",a50)')dsname
     write(datdes(2),'("options big_endian sequential")')
     write(datdes(3),'("TITLE ",a50)')title
     write(datdes(4),'("UNDEF ",e11.2)')undef
     next=5
     write(datdes(next),'("XDEF ",i5," LINEAR ",f7.2,f7.2)')nlon,startp,pinc
     next=next+1
     write(datdes(next),'("YDEF ",i5," LINEAR ",f7.2,f7.2)')nlat,startp,pinc
     next=next+1
     write(datdes(next),'("ZDEF ",i5," LINEAR ",f7.2,f7.2)')np,startp,pinc
     next=next+1
     koutmax=1
     write(datdes(next),'("TDEF ",i5," LINEAR 0Z23may1992 24hr")')koutmax
     next=next+1
     write(datdes(next),'("VARS 1")')
     next=next+1
     write(datdes(next),'("f   ",i5," 99 f   ")')nvert
     next=next+1
     write(datdes(next),'("ENDVARS")')
     last=next
     write(ioutdes,'(a112)')(datdes(i),i=1,last)
 
  end if

  do k=1,nvert
     call sub2grid_1(f(1,1,k),work,0,mype)
     if(mype == 0) then
        do j=1,nlon ; do i=1,nlat
           outfield(j,i)=work(i,j)
        end do ; end do
        write(ioutdat)outfield
     end if
  end do

  if(mype == 0) then
     close(ioutdes)
     close(ioutdat)
  end if

end subroutine grads1_ens

subroutine general_grads1(f,nvert,mype,fname,grd)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grads1      generate grads output file for f on subdomains
!   prgmmr: parrish          org: np22                date: 2009-09-17
!
! abstract: create grads output file for variable f on subdomains.  used
!            to visualize the field f using grads software.
!
! program history log:
!   2009-09-17  parrish
!
!   input argument list:
!     f     - field to generate grads output file
!     nvert - number of vertical levels in f
!     mype  - local processor
!     fname - character string used to identify field f in grads output file
!     grd   - contains communication info for f
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_single,r_kind,i_kind
  use constants, only: one
  use general_sub2grid_mod, only: sub2grid_info
  implicit none

  type(sub2grid_info),intent(in):: grd
  integer(i_kind),intent(in   ) :: nvert,mype
  character(*)   ,intent(in   ) :: fname
  real(r_kind)   ,intent(in   ) :: f(grd%lat2,grd%lon2,nvert)

  real(r_kind),dimension(grd%nlat,grd%nlon)::work
  real(r_single) outfield(grd%nlon,grd%nlat)

  character(50) dsname,title,filename
! data dsname/'test.dat'/
  data title/'inmi'/
  character(112) datdes(50000)
  character(1) blank
  data blank/' '/
  data undef/-9.99e33_r_single/

  integer(i_kind) i,k,kend,kstart,next,np,ioutdes,ioutdat
  integer(i_kind) last,j,koutmax
  real(r_single) undef
  real(r_single) startp,pinc

  if(mype == 0) then
     np=nvert
     startp=1._r_single
     pinc=1._r_single
     ioutdes=98550
     ioutdat=98551
     write(filename,'(a,".des")')trim(fname)
     write(dsname,'(a,".dat")')trim(fname)
     open(unit=ioutdes,file=trim(filename),form='formatted')
     open(unit=ioutdat,file=trim(dsname),form='unformatted')
     rewind ioutdes
     rewind ioutdat
     do i=1,50000
        write(datdes(i),'(112a1)')(blank,k=1,112)
     end do
     write(datdes(1),'("DSET ",a50)')dsname
     write(datdes(2),'("options big_endian sequential")')
     write(datdes(3),'("TITLE ",a50)')title
     write(datdes(4),'("UNDEF ",e11.2)')undef
     next=5
     write(datdes(next),'("XDEF ",i5," LINEAR ",f7.2,f7.2)')grd%nlon,startp,pinc
     next=next+1
     write(datdes(next),'("YDEF ",i5," LINEAR ",f7.2,f7.2)')grd%nlat,startp,pinc
     next=next+1
     write(datdes(next),'("ZDEF ",i5," LINEAR ",f7.2,f7.2)')np,startp,pinc
     next=next+1
     koutmax=1
     write(datdes(next),'("TDEF ",i5," LINEAR 0Z23may1992 24hr")')koutmax
     next=next+1
     write(datdes(next),'("VARS 1")')
     next=next+1
     write(datdes(next),'("f   ",i5," 99 f   ")')nvert
     next=next+1
     write(datdes(next),'("ENDVARS")')
     last=next
     write(ioutdes,'(a112)')(datdes(i),i=1,last)
 
  end if

  do k=1,nvert
     call general_sub2grid_1_ens(f(1,1,k),work,0,mype,grd)
     if(mype == 0) then
        do j=1,grd%nlon ; do i=1,grd%nlat
           outfield(j,i)=work(i,j)
        end do ; end do
        write(ioutdat)outfield
     end if
  end do

  if(mype == 0) then
     close(ioutdes)
     close(ioutdat)
  end if

end subroutine general_grads1

subroutine general_sub2grid_1_ens(sub,grid,gridpe,mype,grd)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grads1      generate grads output file for f on subdomains
!   prgmmr: parrish          org: np22                date: 2009-09-17
!
! abstract: straightforward, but inefficient code to convert a single variable
!            on subdomains to complete slab on one processor.
!
! program history log:
!   2009-09-17  parrish
!
!   input argument list:
!     sub   - field on subdomains
!     gridpe- processor that contains full slab representation of sub
!     mype  - local processor
!
!   output argument list:
!     grid  - field on complete slab
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$


  use kinds, only: r_kind,i_kind
  use constants, only: zero
  use mpimod, only: mpi_comm_world,ierror,mpi_rtype
  use general_sub2grid_mod, only: sub2grid_info
  implicit none

  type(sub2grid_info),intent(in):: grd
  integer(i_kind)                  ,intent(in   ) :: gridpe,mype
  real(r_kind),dimension(grd%lat2,grd%lon2),intent(in   ) :: sub
  real(r_kind),dimension(grd%nlat,grd%nlon),intent(  out) :: grid

  real(r_kind),dimension(grd%lat1,grd%lon1):: zsm
  real(r_kind),dimension(grd%itotsub):: work1
  integer(i_kind) mm1,i,i0,j,j0,k

  mm1=mype+1

  do j=2,grd%lon2-1
     j0=j-1
     do i=2,grd%lat2-1
        i0=i-1
        zsm(i0,j0)=sub(i,j)
     end do
  end do
  call mpi_gatherv(zsm,grd%ijn(mm1),mpi_rtype, &
                 work1,grd%ijn,grd%displs_g,mpi_rtype, &
                 gridpe,mpi_comm_world,ierror)
  if(mype == gridpe) then
     do k=1,grd%iglobal
        i=grd%ltosi(k) ; j=grd%ltosj(k)
        grid(i,j)=work1(k)
     end do
  end if

end subroutine general_sub2grid_1_ens

subroutine sub2grid_1_ens(sub,grid,gridpe,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    grads1      generate grads output file for f on subdomains
!   prgmmr: parrish          org: np22                date: 2009-09-17
!
! abstract: straightforward, but inefficient code to convert a single variable
!            on subdomains to complete slab on one processor.
!
! program history log:
!   2009-09-17  parrish
!
!   input argument list:
!     sub   - field on subdomains
!     gridpe- processor that contains full slab representation of sub
!     mype  - local processor
!
!   output argument list:
!     grid  - field on complete slab
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$


  use kinds, only: r_kind,i_kind
  use constants, only: zero
  use hybrid_ensemble_parameters, only: grd_ens
  use mpimod, only: mpi_comm_world,ierror,mpi_rtype
  implicit none

  integer(i_kind)                  ,intent(in   ) :: gridpe,mype
  real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2),intent(in   ) :: sub
  real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon),intent(  out) :: grid

  real(r_kind),dimension(grd_ens%lat1,grd_ens%lon1):: zsm
  real(r_kind),dimension(grd_ens%itotsub):: work1
  integer(i_kind) mm1,i,i0,j,j0,k

  mm1=mype+1

  do j=2,grd_ens%lon2-1
     j0=j-1
     do i=2,grd_ens%lat2-1
        i0=i-1
        zsm(i0,j0)=sub(i,j)
     end do
  end do
  call mpi_gatherv(zsm,grd_ens%ijn(mm1),mpi_rtype, &
                 work1,grd_ens%ijn,grd_ens%displs_g,mpi_rtype, &
                 gridpe,mpi_comm_world,ierror)
  if(mype == gridpe) then
     do k=1,grd_ens%iglobal
        i=grd_ens%ltosi(k) ; j=grd_ens%ltosj(k)
        grid(i,j)=work1(k)
     end do
  end if

end subroutine sub2grid_1_ens

subroutine get_region_lat_lon_ens(region_lat_ens,region_lon_ens,rlat_e,rlon_e,nlat_e,nlon_e)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    get_region_lat_lon_ens
!   prgmmr: parrish          org: np22                date: 2009-09-25
!
! abstract: compute earth lats and lons on regional ensemble grid, which has the same coordinates
!            as the regional analysis grid, but on a coarser grid with a halo zone to allow full
!            accuracy interpolation of ensemble perturbations to the entire domain of the analysis grid.
!
! program history log:
!   2010-09-25  parrish, initial documentation
!
!   input argument list:
!     rlat_e - regional ensemble grid latitude coordinates in analysis grid units
!     rlon_e - regional ensemble grid longitude coordinates in analysis grid units
!     nlat_e - number of ensemble grid latitudes
!     nlon_e - number of ensemble grid longitudes
!
!   output argument list:
!     region_lat_ens - earth lats at regional ensemble grid points
!     region_lon_ens - earth lons at regional ensemble grid points
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind,i_kind,r_single
  use constants, only: half,one,two,pi
  use gridmod, only: nlon,nlat,txy2ll
                               use constants, only: rad2deg     ! debug only
  use gridmod, only: region_lat,region_lon
  use mpimod, only: mype
  implicit none

  real(r_kind),intent(in):: rlat_e(nlat_e),rlon_e(nlon_e)
  integer(i_kind),intent(in):: nlon_e,nlat_e
  real(r_kind),intent(out):: region_lat_ens(nlat_e,nlon_e),region_lon_ens(nlat_e,nlon_e)

  integer(i_kind) i,j,k
  logical make_test_maps
  real(r_single),allocatable::out1e(:,:)
  real(r_single),allocatable::out1(:,:)
  real(r_kind),allocatable::region_lat_e(:,:),region_lon_e(:,:)
  real(r_kind) twopi

  twopi=two*pi

  do j=1,nlon_e
     do i=1,nlat_e
        call txy2ll(rlon_e(j),rlat_e(i),region_lon_ens(i,j),region_lat_ens(i,j))
        do k=-2,2
           if(region_lon_ens(i,j)<-pi) region_lon_ens(i,j)=region_lon_ens(i,j)+twopi
           if(region_lon_ens(i,j)> pi) region_lon_ens(i,j)=region_lon_ens(i,j)-twopi
        end do
     end do
  end do
                             if(mype==0) write(6,*)' min,max(region_lon_ens)=', &
                             rad2deg*minval(region_lon_ens),rad2deg*maxval(region_lon_ens)

  make_test_maps=.false.
  allocate(out1e(nlon_e,nlat_e))
  allocate(out1(nlon,nlat))
  if(make_test_maps.and.mype==0) then
         do j=1,nlon_e
            do i=1,nlat_e
               out1e(j,i)=region_lon_ens(i,j)
            end do
         end do
         call outgrads1(out1e,nlon_e,nlat_e,'region_lon_e')
                do j=1,nlon
                   do i=1,nlat
                      out1(j,i)=region_lon(i,j)
                   end do
                end do
                call outgrads1(out1,nlon,nlat,'region_lon')
         do j=1,nlon_e
            do i=1,nlat_e
               out1e(j,i)=region_lat_ens(i,j)
            end do
         end do
         call outgrads1(out1e,nlon_e,nlat_e,'region_lat_e')
                do j=1,nlon
                   do i=1,nlat
                      out1(j,i)=region_lat(i,j)
                   end do
                end do
                call outgrads1(out1,nlon,nlat,'region_lat')
  end if

end subroutine get_region_lat_lon_ens

subroutine get_region_dx_dy_ens(region_dx_ens,region_dy_ens)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    get_region_dx_dy_ens
!   prgmmr: parrish          org: np22                date: 2009-09-25
!
! abstract: compute earth grid dx, dy  on regional ensemble grid, which has the same coordinates
!            as the regional analysis grid, but on a coarser grid.
!
! program history log:
!   2010-09-25  parrish, initial documentation
!
!   input argument list:
!
!   output argument list:
!     region_dx_ens  - ensemble grid increment in x-direction (m)
!     region_dy_ens  - ensemble grid increment in y-direction (m)
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind,i_kind,r_single
  use hybrid_ensemble_parameters, only: nlon_ens,nlat_ens,p_e2a
  use egrid2agrid_mod, only: agrid2egrid
  use constants, only: half,one,two
  use gridmod, only: nlon,nlat
  use gridmod, only: region_dx,region_dy
  use mpimod, only: mype
  implicit none

  real(r_kind),intent(out):: region_dx_ens(nlat_ens,nlon_ens),region_dy_ens(nlat_ens,nlon_ens)

  integer(i_kind) i,j
  logical make_test_maps
  real(r_single),allocatable::out1ens(:,:),out1(:,:)
  real(r_kind) worka(1,nlat,nlon,1),worke(1,nlat_ens,nlon_ens,1)
  real(r_kind) ratio

  ratio=(nlon-one)/(nlon_ens-one)
  worka(1,:,:,1)=region_dx(:,:)
  call agrid2egrid(p_e2a,worka,worke,1,1)
  region_dx_ens(:,:)=worke(1,:,:,1)*ratio
  ratio=(nlat-one)/(nlat_ens-one)
  worka(1,:,:,1)=region_dy(:,:)
  call agrid2egrid(p_e2a,worka,worke,1,1)
  region_dy_ens(:,:)=worke(1,:,:,1)*ratio

  make_test_maps=.false.
  allocate(out1ens(nlon_ens,nlat_ens))
  allocate(out1(nlon,nlat))
  if(make_test_maps.and.mype==0) then
         do j=1,nlon_ens
            do i=1,nlat_ens
               out1ens(j,i)=region_dx_ens(i,j)
            end do
         end do
         call outgrads1(out1ens,nlon_ens,nlat_ens,'region_dx_ens')
                do j=1,nlon
                   do i=1,nlat
                      out1(j,i)=region_dx(i,j)
                   end do
                end do
                call outgrads1(out1,nlon,nlat,'region_dx')
         do j=1,nlon_ens
            do i=1,nlat_ens
               out1ens(j,i)=region_dy_ens(i,j)
            end do
         end do
         call outgrads1(out1ens,nlon_ens,nlat_ens,'region_dy_ens')
                do j=1,nlon
                   do i=1,nlat
                      out1(j,i)=region_dy(i,j)
                   end do
                end do
                call outgrads1(out1,nlon,nlat,'region_dy')
  end if
  deallocate(out1ens,out1)

end subroutine get_region_dx_dy_ens

subroutine get_regional_dual_res_grid(eps,r_e,n_a,n_e,x_a,x_e)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    get_regional_dual_res_grid
!   prgmmr: parrish          org: np22                date: 2009-09-25
!
! abstract: create ensemble grid centered on regional grid domain.
!           The ensemble grid extends beyond the edge of the analysis domain by
!           amount eps (in analysis grid units), and the resolution of the
!           ensemble grid is in ratio r_e to the analysis grid.
!
! program history log:
!   2010-09-21  parrish, initial documentation
!   2011-01-30  parrish -- Rename subroutine get_regional_gefs_grid to get_regional_dual_res_grid.
!
!   input argument list:
!     eps   - width of buffer zone--how far ensemble grid extends beyond analysis grid--
!               in analysis grid units.
!     r_e   - ratio of ensemble grid resolution to analysis grid resolution (approximate--
!               actual resolution can be slightly higher)
!     n_a   - number of analysis grid points in coordinate being processed
!
!   output argument list:
!     n_e   - number of ensemble grid points in coordinate being processed
!     x_a   - analysis coordinate in grid units (1 to n_a)
!     x_e   - ensemble grid coordinate in analysis grid units
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind,i_kind,r_single
  use constants, only: half,one,two
                               use mpimod, only: mype
  implicit none

  real(r_kind),intent(in):: eps      !  width of halo zone that ensemble grid extends beyond analysis grid
                                     !       in analysis grid units
  real(r_kind),intent(in):: r_e      !  ratio of ensemble grid increment to analysis grid increment.
                                     !     because analysis grid increment=1, bigl_e is ensemble grid
                                     !     increment in analysis grid units.
  integer(i_kind),intent(in):: n_a   !  number of analysis coordinate grid points in direction processed
  integer(i_kind),intent(out):: n_e  !  number of ensemble grid points in direction processed
  real(r_kind),intent(out)::x_a(n_a) !  analysis coordinate in analysis grid units
  real(r_kind),intent(out)::x_e(n_a) !  ensemble coordinate in analysis grid units

  real(r_kind) bigl_e,bigl_e_x,r_e_x
  integer(i_kind) i,j,n

!      schematic for one dimension of domain:
!
!       |       |                                |       |
!        < eps >                                  < eps >
!
!               1                               n_a

  bigl_e = n_a - one + two*eps
  n_e = ceiling(one + bigl_e/r_e)
  r_e_x = bigl_e/(n_e-one)
  do n=1,n_a
     x_a(n)=n
  end do
  do n=1,n_e
     x_e(n) =one - eps + r_e_x*(n-one)
  end do

end subroutine get_regional_dual_res_grid
