subroutine get_gefs_ensperts_dualres
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    get_gefs_ensperts_dualres copy of get_gefs_ensperts for dual resolution
!   prgmmr: kleist           org: np22                date: 2010-01-05
!
! abstract: read ensemble members, and construct ensemble perturbations, for use
!             with hybrid ensemble option.  ensemble spread is also written out as
!             a byproduct for diagnostic purposes.
!
!
! program history log:
!   2010-01-05  kleist, initial documentation
!   2010-02-17  parrish - make changes to allow dual resolution capability
!   2010-03-24  derber - use generalized genqsat rather than specialized for this resolution
!   2010-04-06  parrish - comment out call to ensemble_spread_dualres for now--there is a bug
!                           and daryl has a newer version which is faster and uses less memory.
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

  use gridmod, only: idsl5
  use hybrid_ensemble_parameters, only: n_ens
  use hybrid_ensemble_isotropic_regional, only: st_en,vp_en,t_en,rh_en,oz_en,cw_en,p_en,sst_en
  use constants,only: zero,ione,izero,half,fv,rd_over_cp,one
  use mpimod, only: mpi_comm_world,ierror,mype
  use kinds, only: r_kind,i_kind
  use hybrid_ensemble_parameters, only: grd_ens,nlat_ens,nlon_ens,sp_ens
  implicit none

  real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig+1) :: pri
  real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig):: vor,div,u,v,tv,q,cwmr,oz,qs,rh,tsen,prsl
  real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2):: z,ps,sst2
  real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon):: sst_full,dum
  real(r_kind),dimension(grd_ens%latlon1n):: stbar,vpbar,tbar,rhbar,ozbar,cwbar
  real(r_kind),dimension(grd_ens%latlon11):: pbar,sstbar
  real(r_kind) bar_norm,sig_norm,kapr,kap1

  integer(i_kind),dimension(grd_ens%nlat,grd_ens%nlon):: idum
  integer(i_kind) iret,i,j,k,m,n,il,jl,mm1,iderivative
  character(24) filename
  logical ice

  stbar=zero ; vpbar=zero ; tbar=zero ; rhbar=zero ; ozbar=zero ; cwbar=zero 
  pbar=zero ; sstbar =zero

  mm1=mype+ione
  kap1=rd_over_cp+one
  kapr=one/rd_over_cp

  do n=1,n_ens
    write(filename,100) n
100        format('sigf06_ens_mem',i3.3)

! Use read_gfs routine to get ensemble members
! *** NOTE ***
! For now, everything is assumed to be at the same resolution
    if (mype==0)write(6,*) 'CALL READ_GFSATM FOR ENS FILE : ',filename
    call general_read_gfsatm(grd_ens,sp_ens,filename,mype,z,ps,vor,div,u,v,tv,q,cwmr,oz,iret)

! Temporarily, try to use the u,v option (i.e. uv_hyb_ens)...eventually two
! options exist:
!   1. Create modified version of read_gfsatm that only passes back out the things
!      we need, i.e. SF, VP, TV, RH, OZ, CW, and PS
!   2. Create subroutine to go from gridded values of vor/div on subdomain to get
!      SF/VP on subdomain (I can't believe this doesn't exist yet!)
  
! Compute RH
! First step is go get sensible temperature and 3d pressure
    do k=1,grd_ens%nsig
      do j=1,grd_ens%lon2
        do i=1,grd_ens%lat2
          tsen(i,j,k)= tv(i,j,k)/(one+fv*max(zero,q(i,j,k)))
        end do
      end do
    end do
! Get 3d pressure field now on interfaces
    call general_getprs_glb(ps,tv,pri)
    if (idsl5/=2_i_kind) then
      do j=1,grd_ens%lon2
        do i=1,grd_ens%lat2
          do k=1,grd_ens%nsig
            prsl(i,j,k)=((pri(i,j,k)**kap1-pri(i,j,k+ione)**kap1)/&
                           (kap1*(pri(i,j,k)-pri(i,j,k+ione))))**kapr
          end do
        end do
      end do
    else 
      do j=1,grd_ens%lon2
        do i=1,grd_ens%lat2
          do k=1,grd_ens%nsig
            prsl(i,j,k)=(pri(i,j,k)+pri(i,j,k+ione))*half
          end do
        end do
      end do
    end if

    ice=.true.
    iderivative=0
    call genqsat(qs,tsen,prsl,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig,ice,iderivative)

    do k=1,grd_ens%nsig
      do j=1,grd_ens%lon2
        do i=1,grd_ens%lat2
          rh(i,j,k) = q(i,j,k)/qs(i,j,k)
        end do
      end do
    end do

    m=izero
    do k=1,grd_ens%nsig
      do j=1,grd_ens%lon2
        do i=1,grd_ens%lat2
          m=m+ione
          st_en(m,n) = u(i,j,k)
          vp_en(m,n) = v(i,j,k)
          t_en(m,n)  = tv(i,j,k)
          rh_en(m,n) = rh(i,j,k)
          oz_en(m,n) = oz(i,j,k)
          cw_en(m,n) = cwmr(i,j,k)

          stbar(m)=stbar(m)+u(i,j,k)
          vpbar(m)=vpbar(m)+v(i,j,k)
           tbar(m)= tbar(m)+ tv(i,j,k)
          rhbar(m)=rhbar(m)+rh(i,j,k)
          ozbar(m)=ozbar(m)+oz(i,j,k)
          cwbar(m)=cwbar(m)+cwmr(i,j,k)
        end do
      end do
    end do
    m=izero
    do j=1,grd_ens%lon2
      do i=1,grd_ens%lat2
        m=m+ione
        p_en(m,n) = ps(i,j)
        pbar(m)=pbar(m)+ ps(i,j)
      end do
    end do
  end do ! end do over ensemble

! Convert to mean
  bar_norm = one/float(n_ens)
  do i=1,grd_ens%latlon1n
    stbar(i)=stbar(i)*bar_norm
    vpbar(i)=vpbar(i)*bar_norm
    tbar(i) = tbar(i)*bar_norm
    rhbar(i)=rhbar(i)*bar_norm
    ozbar(i)=ozbar(i)*bar_norm
    cwbar(i)=cwbar(i)*bar_norm
  end do
  do i=1,grd_ens%latlon11
    pbar(i)=pbar(i)*bar_norm
  end do

  call mpi_barrier(mpi_comm_world,ierror)
! Before converting to perturbations, get ensemble spread
! call ens_spread_dualres(stbar,vpbar,tbar,rhbar,ozbar,cwbar,pbar,mype)
! call mpi_barrier(mpi_comm_world,ierror)

! Convert ensemble members to perturbations
   sig_norm=sqrt(one/max(one,n_ens-one))

  do n=1,n_ens
    do i=1,grd_ens%latlon1n
      st_en(i,n)=(st_en(i,n)-stbar(i))*sig_norm
      vp_en(i,n)=(vp_en(i,n)-vpbar(i))*sig_norm
      t_en(i,n) =( t_en(i,n)- tbar(i))*sig_norm
      rh_en(i,n)=(rh_en(i,n)-rhbar(i))*sig_norm
      oz_en(i,n)=(oz_en(i,n)-ozbar(i))*sig_norm
      cw_en(i,n)=(cw_en(i,n)-cwbar(i))*sig_norm
    end do
    do i=1,grd_ens%latlon11
      p_en(i,n)=(p_en(i,n)- pbar(i))*sig_norm
    end do
  end do

!  since initial version is ignoring sst perturbations, skip following code for now.  revisit
!   later--creating general_read_gfssfc, analogous to general_read_gfsatm above.
!! GET SST PERTURBATIONS HERE
!  do n=1,n_ens
!    write(filename,105) n
!105        format('sfcf06_ens_mem',i3.3)
!
!! This will get full 2d nlat x nlon sst field
!    if(mype==0)write(6,*) 'CALL READ_GFSSFC FOR ENS FILE : ',filename
!    call read_gfssfc(filename,mype,&
!         dum,sst_full,dum, &
!         dum,dum,dum,dum,dum,idum,dum,dum)
!
!    call mpi_barrier(mpi_comm_world,ierror)
!
!! mpi barrier here?
!    do j=1,grd_ens%lon2
!      jl=j+grd_ens%jstart(mm1)-2
!      jl=min0(max0(1,jl),grd_ens%nlon)
!      do i=1,grd_ens%lat2
!        il=i+grd_ens%istart(mm1)-2
!        il=min0(max0(1,il),grd_ens%nlat)
!        sst2(i,j)=sst_full(il,jl)
!      end do
!    end do
!  
!    m=izero
!    do j=1,grd_ens%lon2
!      do i=1,grd_ens%lat2
!        m=m+ione
!        sst_en(m,n) = sst2(i,j)
!        sstbar(m)=sstbar(m)+ sst2(i,j)
!      end do
!    end do
!  end do ! end do over ensemble
!
!  do n=1,n_ens
!    do i=1,grd_ens%latlon11
!      sst_en(i,n)=(sst_en(i,n)- sstbar(i)*bar_norm)
!    end do
!  end do

! dtk: temporarily ignore sst perturbations in hybrid
  do n=1,n_ens
    do i=1,grd_ens%latlon11
      sst_en(i,n)=zero
    end do
  end do


  return
end subroutine get_gefs_ensperts_dualres

subroutine ens_spread_dualres(stbar,vpbar,tbar,rhbar,ozbar,cwbar,pbar,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ens_spread_dualres  output ensemble spread for diagnostics
!   prgmmr: kleist           org: np22                date: 2010-01-05
!
! abstract: compute ensemble spread on ensemble grid, interpolate to analysis grid
!             and write out for diagnostic purposes.
!
!
! program history log:
!   2010-01-05  kleist, initial documentation
!   2010-02-28  parrish - make changes to allow dual resolution capability
!
!   input argument list:
!     stbar  - ensemble mean stream function (or u wind component)
!     vpbar  - ensemble mean stream function (or v wind component)
!      tbar  - ensemble mean Tv
!     rhbar  - ensemble mean rh
!     ozbar  - ensemble mean ozone
!     cwbar  - ensemble mean cloud water
!      pbar  - ensemble mean surface pressure
!      mype  - current processor number
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
  use kinds, only: r_single,r_kind,i_kind
  use hybrid_ensemble_parameters, only: n_ens,grd_ens,grd_anl,p_e2a,uv_hyb_ens
  use hybrid_ensemble_isotropic_regional, only: st_en,vp_en,t_en,rh_en,oz_en,cw_en,p_en
  use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info,general_sube2suba_r_double
  use constants, only:  zero,two,half,one
  implicit none

  real(r_kind),dimension(grd_ens%latlon1n),intent(in):: stbar,vpbar,tbar,rhbar,ozbar,cwbar
  real(r_kind),dimension(grd_ens%latlon11),intent(in):: pbar
  integer(i_kind),intent(in):: mype

  real(r_kind),dimension(grd_ens%latlon11*(grd_ens%nsig*6+1)):: sube
  real(r_kind),dimension(grd_anl%latlon11*(grd_anl%nsig*6+1)):: suba
  real(r_kind) sp_norm
  type(sub2grid_info)::se,sa

  integer(i_kind) i,ii,n
  integer(i_kind) ist,ivp,it,irh,ioz,icw,ip
  logical regional
  integer(i_kind) num_fields,inner_vars
  logical,allocatable::vector(:)

  sp_norm=(one/float(n_ens))

  sube=zero
  do n=1,n_ens
    ii=0
    do i=1,grd_ens%latlon1n
       ii=ii+1
       sube(ii) =sube(ii)  + (st_en(i,n)-stbar(i))*(st_en(i,n)-stbar(i)) 
    end do
    do i=1,grd_ens%latlon1n
       ii=ii+1
       sube(ii) =sube(ii)  + (vp_en(i,n)-vpbar(i))*(vp_en(i,n)-vpbar(i)) 
    end do
    do i=1,grd_ens%latlon1n
       ii=ii+1
       sube(ii) =sube(ii)  + ( t_en(i,n)- tbar(i))*( t_en(i,n)- tbar(i)) 
    end do
    do i=1,grd_ens%latlon1n
       ii=ii+1
       sube(ii) =sube(ii)  + (rh_en(i,n)-rhbar(i))*(rh_en(i,n)-rhbar(i)) 
    end do
    do i=1,grd_ens%latlon1n
       ii=ii+1
       sube(ii) =sube(ii)  + (oz_en(i,n)-ozbar(i))*(oz_en(i,n)-ozbar(i)) 
    end do
    do i=1,grd_ens%latlon1n
       ii=ii+1
       sube(ii) =sube(ii)  + (cw_en(i,n)-cwbar(i))*(cw_en(i,n)-cwbar(i)) 
    end do
    do i=1,grd_ens%latlon11
       ii=ii+1
       sube(ii) =sube(ii)  + ( p_en(i,n)- pbar(i))*( p_en(i,n)- pbar(i)) 
    end do
  end do
 
  do i=1,grd_ens%latlon11*(grd_ens%nsig*6+1)
    sube(i) = sqrt(sp_norm*sube(i))
  end do

  if(grd_ens%latlon1n == grd_anl%latlon1n) then
     suba=sube
  else
     regional=.false.
     inner_vars=1
     num_fields=6*grd_ens%nsig+1
     allocate(vector(num_fields))
     vector=.false.
     vector(1:2*grd_ens%nsig)=uv_hyb_ens   !  assume here that 1st two 3d variables are either u,v or psi,chi
     call general_sub2grid_create_info(se,inner_vars,grd_ens%nlat,grd_ens%nlon,grd_ens%nsig,num_fields, &
                                       regional,vector)
     call general_sub2grid_create_info(sa,inner_vars,grd_anl%nlat,grd_anl%nlon,grd_anl%nsig,num_fields, &
                                       regional,vector)
     deallocate(vector)
     call general_sube2suba_r_double(se,sa,p_e2a,sube,suba,regional)
  end if

  ist=1
  ivp=grd_anl%latlon1n+ist
  it =grd_anl%latlon1n+ivp
  irh=grd_anl%latlon1n+it
  ioz=grd_anl%latlon1n+irh
  icw=grd_anl%latlon1n+ioz
  ip =grd_anl%latlon1n+icw
  call write_spread_dualres(suba(ist),suba(ivp),suba(it),suba(irh),suba(ioz),suba(icw),suba(ip),mype)

  return
end subroutine ens_spread_dualres


subroutine write_spread_dualres(a,b,c,d,e,f,g,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    write_spread_dualres   write ensemble spread for diagnostics
!   prgmmr: kleist           org: np22                date: 2010-01-05
!
! abstract: write ensemble spread (previously interpolated to analysis grid)
!             for diagnostic purposes.
!
!
! program history log:
!   2010-01-05  kleist, initial documentation
!   2010-02-28  parrish - make changes to allow dual resolution capability
!
!   input argument list:
!     a    -  spread variable 1
!     b    -  spread variable 2
!     c    -  spread variable 3
!     d    -  spread variable 4
!     e    -  spread variable 5
!     f    -  spread variable 6
!     g    -  spread variable 7
!     mype -  current processor number
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
  use kinds, only: r_kind,i_kind,r_single
  use hybrid_ensemble_parameters, only: grd_anl
  implicit none

  integer(i_kind),intent(in):: mype
  character(255):: grdfile

  real(r_kind),dimension(grd_anl%lat2,grd_anl%lon2,grd_anl%nsig),intent(in):: a,b,c,d,e,f
  real(r_kind),dimension(grd_anl%lat2,grd_anl%lon2),intent(in):: g

  real(r_kind),dimension(grd_anl%nlat,grd_anl%nlon,grd_anl%nsig):: ag,bg,cg,dg,eg,fg
  real(r_kind),dimension(grd_anl%nlat,grd_anl%nlon):: gg

  real(r_single),dimension(grd_anl%nlon,grd_anl%nlat,grd_anl%nsig):: a4,b4,c4,d4,e4,f4
  real(r_single),dimension(grd_anl%nlon,grd_anl%nlat):: g4

  integer(i_kind) ncfggg,iret,i,j,k

! gather stuff to processor 0
  do k=1,grd_anl%nsig
    call gather_stuff2(a(1,1,k),ag(1,1,k),mype,0)
    call gather_stuff2(b(1,1,k),bg(1,1,k),mype,0)
    call gather_stuff2(c(1,1,k),cg(1,1,k),mype,0)
    call gather_stuff2(d(1,1,k),dg(1,1,k),mype,0)
    call gather_stuff2(e(1,1,k),eg(1,1,k),mype,0)
    call gather_stuff2(f(1,1,k),fg(1,1,k),mype,0)
  end do
  call gather_stuff2(g,gg,mype,0)

  if (mype==0) then
    write(6,*) 'WRITE OUT NEW VARIANCES'
! load single precision arrays
    do k=1,grd_anl%nsig
      do j=1,grd_anl%nlon
        do i=1,grd_anl%nlat
          a4(j,i,k)=ag(i,j,k)
          b4(j,i,k)=bg(i,j,k)
          c4(j,i,k)=cg(i,j,k)
          d4(j,i,k)=dg(i,j,k)
          e4(j,i,k)=eg(i,j,k)
          f4(j,i,k)=fg(i,j,k)
        end do
      end do
    end do
    do j=1,grd_anl%nlon
      do i=1,grd_anl%nlat
        g4(j,i)=gg(i,j)
      end do
    end do

! Create byte-addressable binary file for grads
    grdfile='ens_spread.grd'
    ncfggg=len_trim(grdfile)
    call baopenwt(22,grdfile(1:ncfggg),iret)
    call wryte(22,4*grd_anl%nlat*grd_anl%nlon*grd_anl%nsig,a4)
    call wryte(22,4*grd_anl%nlat*grd_anl%nlon*grd_anl%nsig,b4)
    call wryte(22,4*grd_anl%nlat*grd_anl%nlon*grd_anl%nsig,c4)
    call wryte(22,4*grd_anl%nlat*grd_anl%nlon*grd_anl%nsig,d4)
    call wryte(22,4*grd_anl%nlat*grd_anl%nlon*grd_anl%nsig,e4)
    call wryte(22,4*grd_anl%nlat*grd_anl%nlon*grd_anl%nsig,f4)
    call wryte(22,4*grd_anl%nlat*grd_anl%nlon,g4)
    call baclose(22,iret)
  end if

  return
end subroutine write_spread_dualres

subroutine general_getprs_glb(ps,tv,prs)
! subprogram:    getprs       get 3d pressure or 3d pressure deriv
!   prgmmr: kleist           org: np20                date: 2005-09-29
!
! abstract: calculate 3d pressure and its horizontal derivatives
!
! program history log:
!   2005-09-29  kleist
!   2006-04-12  treadon - replace sigi with bk5
!   2006-07-31  kleist  - analysis variable change from ln(ps) to ps
!   2007-05-08  kleist  - add generalized vert coord and derivative call
!   2007-07-26  cucurull- compute 3d pressure and derivatives in different subroutines
!                       - remove gues_tv from argument list; clean up code
!   2008-06-04  safford - rm unused uses
!   2008-09-05  lueken  - merged ed's changes into q1fy09 code
!   2010-02-23  parrish - copy getprs and generalize for dual resolution.
!
!   input argument list:
!     ps       - surface pressure
!
!   output argument list:
!     prs        - 3d pressure
!
! attributes:
!   language:  f90
!   machine:   ibm/RS6000 SP
!
!$$$ end documentation block

  use kinds,only: r_kind,i_kind
  use constants,only: ione,zero,half,one_tenth,rd_over_cp,one
  use gridmod,only: nsig,lat2,lon2,ak5,bk5,ck5,tref5,idvc5
  use gridmod,only: wrf_nmm_regional,nems_nmmb_regional,eta1_ll,eta2_ll,pdtop_ll,pt_ll,&
       regional,wrf_mass_regional,twodvar_regional
  use hybrid_ensemble_parameters, only: grd_ens
! use guess_grids, only: ges_tv,ntguessig
                                                                   use mpimod, only: mype
  implicit none

! Declare passed variables
  real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2)          ,intent(in   ) :: ps
  real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,nsig)     ,intent(in   ) :: tv
  real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,nsig+ione),intent(  out) :: prs

! Declare local variables
  real(r_kind) kapr,trk
  integer(i_kind) i,j,k,k2    ! ,it

! Declare local parameter
  real(r_kind),parameter:: ten = 10.0_r_kind

                                     
  kapr=one/rd_over_cp
  prs=zero 
! it=ntguessig

  if (regional) then
     if(wrf_nmm_regional.or.nems_nmmb_regional) then
        do k=1,nsig+ione
           do j=1,grd_ens%lon2
              do i=1,grd_ens%lat2
                 prs(i,j,k)=one_tenth* &
                      (eta1_ll(k)*pdtop_ll + &
                      eta2_ll(k)*(ten*ps(i,j)-pdtop_ll-pt_ll) + &
                      pt_ll)
              end do
           end do
        end do
     elseif (wrf_mass_regional .or. twodvar_regional) then
        do k=1,nsig+ione
           do j=1,grd_ens%lon2
              do i=1,grd_ens%lat2
                 prs(i,j,k)=one_tenth*(eta1_ll(k)*(ten*ps(i,j)-pt_ll) + pt_ll)
              end do
           end do
        end do
     endif
  else
     k=ione
     k2=nsig+ione
     do j=1,grd_ens%lon2
        do i=1,grd_ens%lat2
           prs(i,j,k)=ps(i,j)
           prs(i,j,k2)=zero
        end do
     end do
     if (idvc5 /= 3_i_kind) then
        do k=2,nsig
           do j=1,grd_ens%lon2
              do i=1,grd_ens%lat2
                 prs(i,j,k)=ak5(k)+bk5(k)*ps(i,j)
              end do
           end do
        end do
     else
        do k=2,nsig
           do j=1,grd_ens%lon2
              do i=1,grd_ens%lat2
                 trk=(half*(tv(i,j,k-ione)+tv(i,j,k))/tref5(k))**kapr
                 prs(i,j,k)=ak5(k)+(bk5(k)*ps(i,j))+(ck5(k)*trk)
              end do
           end do
        end do
     end if
  end if

  return
end subroutine general_getprs_glb