!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !MODULE:  guess_grids --- Guess-related grid definitions
!
! !INTERFACE:
!

module guess_grids

! !USES:
 
  use kinds, only: r_single,r_kind,i_kind
  use gridmod, only: regional
  use gridmod, only: wrf_nmm_regional,nems_nmmb_regional
  use gridmod, only: eta1_ll
  use gridmod, only: eta2_ll
  use gridmod, only: aeta1_ll
  use gridmod, only: aeta2_ll
  use gridmod, only: pdtop_ll
  use gridmod, only: pt_ll
  use regional_io, only: update_pint

  ! meteorological guess (beyond standard ones)
  use gsi_metguess_mod, only: gsi_metguess_create_grids
  use gsi_metguess_mod, only: gsi_metguess_destroy_grids
  use gsi_metguess_mod, only: gsi_metguess_get
  use gsi_metguess_mod, only: gsi_metguess_bundle

  ! chem trace gases
  use gsi_chemguess_mod, only: gsi_chemguess_create_grids
  use gsi_chemguess_mod, only: gsi_chemguess_destroy_grids
  use gsi_chemguess_mod, only: gsi_chemguess_get

  implicit none

! !DESCRIPTION: module containing variables related to the guess fields
!
! !REVISION HISTORY:
!
!   2003-12-01  kleist
!   2004-05-14  kleist, documentation
!   2004-07-15  todling, protex-compliant prologue; added onlys
!   2005-06-01  treadon - add routine add_rtm_layers
!   2005-06-03  parrish - add horizontal derivatives of guess fields
!   2005-06-10  devenyi/treadon - initialize nfldsig and nfldsfc
!   2005-07-06  parrish - add update_pint, arrays ges_pint, ges_pd
!   2005-08-03  parrish - add array to hold roughness length
!   2005-09-29  kleist - add derivatives of terrain, move prsi allocation
!   2005-11-21  kleist - add tendency arrays
!   2005-11-29  derber - add ozmz remove psfcg
!   2005-11-30  derber - combine create_atm_grids and create_pcp_grids (and destroys)
!   2006-02-02  treadon - prefix prsi,prsl,lnprsl,prslk with "ges_"
!   2006-03-07  treadon - remove ges_prslk (no longer needed)
!   2006-04-14  treadon - add bias_tskin
!   2006-04-17  treadon - add ges_psfcavg and ges_prslavg
!   2006-04-21  kleist - add ges divt and agvt arrays
!   2006-07-28  derber  - clean up add_rtm_layers routine 
!   2006-07-28  derber  - add ges_tsen (sensible temperature) array
!   2006-07-31  kleist - use ges_ps instead of ln(ps)
!   2006-09-20  cucurull - add ln(ges_prsi) array
!   2006-09-29  treadon - add flags to control 10m wind factor recompute
!   2007-05-30  h.liu - remove ozmz
!   2007-06-21  rancic - add pbl (ges_teta)
!   2006-12-01  todling - remove bias stuff; merging GMAO bias correction scheme
!   2006-12-15  todling - add _initialized parameters to control allocations
!   2007-03-15  todling - merged in da Silva/Cruz ESMF changes 
!   2008-02-07  eliu    - fixed the unit difference between prsitmp
!                         (kPa) and toa_pressure (hPa).
!   2009-08-19  guo     - added sfc_grids_allocated_, ges_grids_allocated_,
!			  and gesfinfo_created_ to track the state of the data.
!			  for multi-pass observer.
!			- merged destroy_sfc_grids() and destroy_sfct().
!   2008-08-25  hu    - add array definitions for hydrometeor fields
!                     - add subroutine create_cld_grids and destroy_cld_grids
!   2010-04-16  hou   - add array definitions ges_co2 (co2 mixing ratio) and
!                       control variable igfsco2
!   2010-04-22  todling - remove tracers,vtid,pdryini,xncld
!   2010-05-19  todling - add chem init and destroy (revamp Hou's implementation)
!   2010-08-31  cucurull - add logical use_compress
!   2010-09-15  pagowski - add cmaq
!   2010-12-20  cucurull - add integer nsig_ext 
!   2011-01-05  cucurull - add real gpstop
!   2011-02-11  zhu      - add ges_gust,ges_vis,ges_pblh
!   2011-03-13  li      - add for nst FCST file
!   2011-04-29  todling  - some of cloud fields move to wrf_guess_mod; some to met_guess
!   2011-05-01  todling - cwmr no longer in guess-grids; use metguess bundle now
!   2011-11-01  eliu    - modified condition to allocate/deallocate arrays related to 
!                         cloud water tendencies and derivatives 
!   2011-12-27  kleist  - add 4d guess array for saturation specific humidity
!   2012-01-11  Hu      - add GSD PBL height
!   2013-02-22  Carley  - Add NMMB to GSD PBL height calc
!   2014-03-12  Hu      - Add ges_q2 
!
! !AUTHOR: 
!   kleist           org: np20                date: 2003-12-01
!
!EOP
!-------------------------------------------------------------------------

! set default to private
  private
! set subroutines to public
  public :: create_sfc_grids
  public :: create_ges_grids
  public :: destroy_ges_grids
#ifdef TO_BE_REMOVED
  public :: destroy_sfct
#endif
  public :: destroy_sfc_grids
  public :: create_gesfinfo
  public :: destroy_gesfinfo
  public :: load_prsges
  public :: load_geop_hgt
  public :: load_gsdpbl_hgt
  public :: add_rtm_layers
  public :: load_fact10
  public :: comp_fact10
  public :: guess_grids_print
  public :: guess_grids_stats
  public :: create_metguess_grids
  public :: destroy_metguess_grids
  public :: create_chemges_grids
  public :: destroy_chemges_grids
! set passed variables to public
  public :: ntguessig,ges_ps,ges_tv,ges_prsi,ges_oz,ges_psfcavg,ges_prslavg
  public :: isli2,ges_prsl,ges_z,ges_q,ges_v,ges_u,nfldsig,ges_vor,ges_div
  public :: ges_ozlat,ges_ozlon,ges_qlat,ges_teta,ges_cwmr_lat
  public :: ges_cwmr_lon,ges_v_lon,ges_u_lat,ges_u_lon,ges_v_lat,ges_qlon
  public :: ges_tvlon,ges_tvlat,ges_prs_ten,ges_tv_ten,ges_v_ten,ges_cwmr_ten
  public :: ges_oz_ten,ges_q_ten,fact_tv,tropprs,sfct,ges_u_ten,ges_ps_lat
  public :: ges_ps_lon,ntguessfc,ntguesnst,dsfct,ifilesig,veg_frac,soil_type,veg_type
  public :: sno2,ifilesfc,ifilenst,sfc_rough,fact10,sno,isli,soil_temp,soil_moi
  public :: nfldsfc,nfldnst,hrdifsig,ges_tsen,sfcmod_mm5,sfcmod_gfs,ifact10,hrdifsfc,hrdifnst
  public :: ges_pd,ges_pint,geop_hgti,ges_lnprsi,ges_lnprsl,geop_hgtl,pt_ll,pbl_height
  public :: ges_gust,ges_vis,ges_pblh,ges_qsat
  public :: use_compress,nsig_ext,gpstop
  public :: ges_th2,ges_soilt1,ges_tslb,ges_smois,ges_tsk,ges_q2
  public :: efr_ql,efr_qi,efr_qr,efr_qs,efr_qg,efr_qh
! ajl
  public :: ges_no2,ges_no2lat,ges_no2lon,ges_no2_ten

  public :: ges_initialized
  public :: tnd_initialized
  public :: drv_initialized

  public :: nfldsig_all,nfldsig_now,hrdifsig_all
  public :: nfldsfc_all,nfldsfc_now,hrdifsfc_all
  public :: nfldnst_all,nfldnst_now,hrdifnst_all
  public :: extrap_intime
  public :: ntguessig_ref
  public :: ntguessfc_ref
  public :: ntguesnst_ref

  logical:: sfcmod_gfs = .false.    ! .true. = recompute 10m wind factor using gfs physics
  logical:: sfcmod_mm5 = .false.    ! .true. = recompute 10m wind factor using mm5 physics

  logical:: use_compress = .false. ! true to turn on compressibility factor in geopotential heights

  logical, save :: ges_initialized = .false.
  logical, save :: tnd_initialized = .false.
  logical, save :: drv_initialized = .false.

  integer(i_kind) ntguessig         ! location of actual guess time for sigma fields
  integer(i_kind) ntguessfc         ! location of actual guess time for sfc fields
  integer(i_kind) ntguesnst         ! location of actual guess time for nst FCST fields

  integer(i_kind), save:: ntguessig_ref	! replace ntguessig as the storage for its original value
  integer(i_kind), save:: ntguessfc_ref	! replace ntguessfc as the storage for its original value
  integer(i_kind), save:: ntguesnst_ref ! replace ntguesnst as the storage for its original value

  integer(i_kind):: ifact10 = 0     ! 0 = use 10m wind factor from guess
  integer(i_kind):: nsig_ext = 13   ! use 13 layers above model top to compute the bending angle for gpsro

  ! number of guess sigma/surface times are set in GSI_gridComp.rc

  real(r_kind), allocatable, dimension(:), save:: hrdifsig_all  ! a list of all times
  real(r_kind), allocatable, dimension(:), save:: hrdifsfc_all  ! a list of all times
  real(r_kind), allocatable, dimension(:), save:: hrdifnst_all  ! a list of all times

  integer(i_kind), save:: nfldsig_all	! expected total count of time slots
  integer(i_kind), save:: nfldsfc_all
  integer(i_kind), save:: nfldnst_all

  integer(i_kind), save:: nfldsig	! actual count of in-cache time slots
  integer(i_kind), save:: nfldsfc
  integer(i_kind), save:: nfldnst       ! actual count of in-cache time slots for NST file

  integer(i_kind), save:: nfldsig_now	! current count of filled time slots
  integer(i_kind), save:: nfldsfc_now
  integer(i_kind), save:: nfldnst_now

  logical, save:: extrap_intime		! compute o-f interpolate within the time ranges of guess_grids,
  					! or also extrapolate outside the time ranges.

  real(r_kind), allocatable, dimension(:):: hrdifsig  ! times for cached sigma guess_grid
  real(r_kind), allocatable, dimension(:):: hrdifsfc  ! times for cached surface guess_grid
  real(r_kind), allocatable, dimension(:):: hrdifnst  ! times for cached nst guess_grid

  integer(i_kind),allocatable, dimension(:)::ifilesfc  ! array used to open the correct surface guess files
  integer(i_kind),allocatable, dimension(:)::ifilesig  ! array used to open the correct sigma guess files
  integer(i_kind),allocatable, dimension(:)::ifilenst  ! array used to open the correct nst guess files

  integer(i_kind),allocatable,dimension(:,:,:):: isli    ! snow/land/ice mask
  integer(i_kind),allocatable,dimension(:,:,:):: isli_g  ! isli on horiz/global grid
  integer(i_kind),allocatable,dimension(:,:):: isli2     ! snow/land/ice mask at analysis time

  real(r_kind),allocatable,dimension(:,:,:):: sno2  ! sno depth on subdomain


  real(r_kind):: gpstop=30.0_r_kind   ! maximum gpsro height used in km 
                                      ! geometric height for ref, impact height for bnd

  real(r_kind):: ges_psfcavg                            ! average guess surface pressure 
  real(r_kind),allocatable,dimension(:):: ges_prslavg   ! average guess pressure profile

  real(r_kind),allocatable,dimension(:,:):: tropprs     ! guess tropopause pressure
  real(r_kind),allocatable,dimension(:,:,:):: dsfct     ! delta skin temperature

  real(r_kind),allocatable,dimension(:,:,:):: fact10    ! 10 meter wind factor
  real(r_kind),allocatable,dimension(:,:,:):: sfct      ! guess skin temperature
  real(r_kind),allocatable,dimension(:,:,:):: sno       ! snow-ice mask
  real(r_kind),allocatable,dimension(:,:,:):: veg_type  ! vegetation type
  real(r_kind),allocatable,dimension(:,:,:):: veg_frac  ! vegetation fraction(0-1.0)
  real(r_kind),allocatable,dimension(:,:,:):: sfc_rough ! sfc roughness length
  real(r_kind),allocatable,dimension(:,:,:):: soil_type ! soil type
  real(r_kind),allocatable,dimension(:,:,:):: soil_temp ! soil temperature of first layer
  real(r_kind),allocatable,dimension(:,:,:):: soil_moi  ! soil moisture of first layer
  

  real(r_kind),allocatable,dimension(:,:,:,:):: geop_hgtl ! guess geopotential height at mid-layers
  real(r_kind),allocatable,dimension(:,:,:,:):: geop_hgti ! guess geopotential height at level interfaces

  real(r_kind),allocatable,dimension(:,:,:):: ges_z      ! topography
  real(r_kind),allocatable,dimension(:,:,:):: ges_ps     ! log(surface pressure)
  real(r_kind),allocatable,dimension(:,:,:):: ges_ps_lat ! log(ps)/lat for pcp routine
  real(r_kind),allocatable,dimension(:,:,:):: ges_ps_lon ! log(ps)/lon for pcp routine
  real(r_kind),allocatable,dimension(:,:,:):: ges_gust   ! wind gust speed
  real(r_kind),allocatable,dimension(:,:,:):: ges_vis    ! visibility
  real(r_kind),allocatable,dimension(:,:,:):: ges_pblh   ! pbl height
  real(r_kind),allocatable,dimension(:,:,:):: pbl_height  !  GSD PBL height in hPa
                                                         ! Guess Fields ...
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_prsi  ! interface pressure
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_prsl  ! layer midpoint pressure
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_lnprsl! log(layer midpoint pressure)
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_lnprsi! log(interface pressure)
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_u     ! zonal wind
  real(r_kind),allocatable,dimension(:,:,:):: ges_u_lat ! zonal wind/lat
  real(r_kind),allocatable,dimension(:,:,:):: ges_u_lon ! zonal wind/lon
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_v     ! meridional wind
  real(r_kind),allocatable,dimension(:,:,:):: ges_v_lat ! meridional wind/lat
  real(r_kind),allocatable,dimension(:,:,:):: ges_v_lon ! meridional wind/lon
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_vor   ! vorticity
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_div   ! divergence
  real(r_kind),allocatable,dimension(:,:,:):: ges_cwmr_lat  ! cloud condensate mixing ratio/lat
  real(r_kind),allocatable,dimension(:,:,:):: ges_cwmr_lon  ! cloud condensate mixing ratio/lon
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_q     ! specific humidity
  real(r_kind),allocatable,dimension(:,:,:):: ges_qlon  ! q/lat for pcp routine advection calc
  real(r_kind),allocatable,dimension(:,:,:):: ges_qlat  ! q/lon for pcp routine advection calc
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_oz    ! ozone mixing ratio
  real(r_kind),allocatable,dimension(:,:,:):: ges_ozlat ! ozone mixing ratio/lat
  real(r_kind),allocatable,dimension(:,:,:):: ges_ozlon ! ozone mixing ratio/lon
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_pint  ! pint variable (nmm only)
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_tv    ! virtual temperature
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_tsen  ! sensible temperature
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_teta  ! potential temperature
  real(r_kind),allocatable,dimension(:,:,:):: ges_tvlat ! tv/lat for pcp routine advection calc
  real(r_kind),allocatable,dimension(:,:,:):: ges_tvlon ! tv/lon for pcp routine advection calc 

! ajl
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_no2    ! no2 mixing ratio
  real(r_kind),allocatable,dimension(:,:,:):: ges_no2lat ! no2 mixing ratio/lat
  real(r_kind),allocatable,dimension(:,:,:):: ges_no2lon ! no2 mixing ratio/lon
  real(r_kind),allocatable,dimension(:,:,:):: ges_no2_ten   ! no2 tendency

  real(r_kind),allocatable,dimension(:,:,:)::ges_pd        ! pdges (for nmm only)
  real(r_kind),allocatable,dimension(:,:,:):: ges_prs_ten  ! 3d pressure tendency
  real(r_kind),allocatable,dimension(:,:,:):: ges_u_ten    ! u tendency
  real(r_kind),allocatable,dimension(:,:,:):: ges_v_ten    ! v tendency
  real(r_kind),allocatable,dimension(:,:,:):: ges_tv_ten   ! Tv tendency
  real(r_kind),allocatable,dimension(:,:,:):: ges_q_ten    ! q tendency
  real(r_kind),allocatable,dimension(:,:,:):: ges_oz_ten   ! ozone tendency
  real(r_kind),allocatable,dimension(:,:,:):: ges_cwmr_ten ! cloud water tendency
  real(r_kind),allocatable,dimension(:,:,:):: fact_tv      ! 1./(one+fv*ges_q) for virt to sen calc.
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_qsat      ! 4d qsat array
! for GSD soil nudging
  real(r_kind),allocatable,dimension(:,:,:):: ges_th2       ! 2-m potential temperature
  real(r_kind),allocatable,dimension(:,:,:):: ges_q2        ! 2-m moisture
  real(r_kind),allocatable,dimension(:,:,:):: ges_tsk       ! skin temperature
  real(r_kind),allocatable,dimension(:,:,:):: ges_soilt1    ! TEMPERATURE INSIDE SNOW
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_tslb    ! SOIL TEMPERATURE
  real(r_kind),allocatable,dimension(:,:,:,:):: ges_smois   ! SOIL MOISTURE   

  real(r_kind),allocatable,dimension(:,:,:,:):: efr_ql     ! effective radius for cloud liquid water
  real(r_kind),allocatable,dimension(:,:,:,:):: efr_qi     ! effective radius for cloud ice
  real(r_kind),allocatable,dimension(:,:,:,:):: efr_qr     ! effective radius for rain
  real(r_kind),allocatable,dimension(:,:,:,:):: efr_qs     ! effective radius for snow
  real(r_kind),allocatable,dimension(:,:,:,:):: efr_qg     ! effective radius for graupel
  real(r_kind),allocatable,dimension(:,:,:,:):: efr_qh     ! effective radius for hail

  interface guess_grids_print
     module procedure print1r8_
     module procedure print2r8_
     module procedure print3r8_
     module procedure print4r8_
  end interface
  interface guess_grids_stats
     module procedure guess_grids_stats3d_
     module procedure guess_grids_stats2d_
  end interface


  logical,save:: sfc_grids_allocated_=.false.
  logical,save:: ges_grids_allocated_=.false.
  logical,save:: gesfinfo_created_=.false.

contains

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: create_sfc_grids --- Allocate memory for surface related grids
!
! !INTERFACE:
!
  subroutine create_sfc_grids

! !USES:

   use gridmod, only: lat2,lon2,nlat,nlon
   use constants, only: zero

   use mpeu_util, only: die,tell
   implicit none

! !DESCRIPTION: allocate memory for surface related grids
!
! !REVISION HISTORY:
!   2003-12-01  kleist
!   2004-05-14  kleist, documentation
!   2004-07-15  todling, protex-compliant prologue
!   2004-07-28  treadon - remove subroutine call list, pass variables via modules
!   2005-06-03  parrish - allocate and initialize sfct_lat and sfct_lon
!   2007-03-15  todling - merged in da Silva/Cruz ESMF changes
!   2008-12-5   todling - add time dimension to dsfct
!   2009-01-23  todling - zero out arrays
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
!   kleist          org: w/nmc20     date: 2003-12-01
!
!EOP
!-------------------------------------------------------------------------

    integer(i_kind) :: i,j,it,istatus
    if(sfc_grids_allocated_) call die('create_sfc_grids','alread allocated')
    sfc_grids_allocated_=.true.

    allocate( isli_g(nlat,nlon,nfldsfc),&
         isli2(lat2,lon2),sno2(lat2,lon2,nfldsfc),&
         stat=istatus)
    if (istatus/=0) write(6,*)'CREATE_SFC_GRIDS(1):  allocate error, istatus=',&
         istatus,lat2,lon2,nlat,nlon,nfldsfc

#ifndef HAVE_ESMF
    allocate( isli(lat2,lon2,nfldsfc),fact10(lat2,lon2,nfldsfc),&
         dsfct(lat2,lon2,nfldsfc),sfct(lat2,lon2,nfldsfc),sno(lat2,lon2,nfldsfc),&
         veg_type(lat2,lon2,nfldsfc),veg_frac(lat2,lon2,nfldsfc),&
         sfc_rough(lat2,lon2,nfldsfc),&
         soil_type(lat2,lon2,nfldsfc),soil_temp(lat2,lon2,nfldsfc),&
         soil_moi(lat2,lon2,nfldsfc), &
         stat=istatus)
    if (istatus/=0) write(6,*)'CREATE_SFC_GRIDS(2):  allocate error, istatus=',&
         istatus,lat2,lon2,nlat,nlon,nfldsfc
#endif /* HAVE_ESMF */

    do it=1,nfldsfc
       do j=1,nlon
          do i=1,nlat
             isli_g(i,j,it)=0
          end do
       end do
    end do

#ifndef HAVE_ESMF
    do it=1,nfldsfc
       do j=1,lon2
          do i=1,lat2
             isli(i,j,it)=0
             fact10(i,j,it)=zero
             sfct(i,j,it)=zero
             dsfct(i,j,it)=zero
             sno(i,j,it)=zero
             veg_type(i,j,it)=zero
             veg_frac(i,j,it)=zero
             soil_type(i,j,it)=zero
             soil_temp(i,j,it)=zero
             soil_moi(i,j,it)=zero
          end do
       end do
    end do
#endif

    do it=1,nfldsfc
       do j=1,lon2
          do i=1,lat2
             sno2(i,j,it)=zero
          end do
       end do
    end do

    do j=1,lon2
       do i=1,lat2
          isli2(i,j)=0
       end do
    end do

    return
  end subroutine create_sfc_grids

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE:  create_ges_grids --- Alloc grid for guess and bias
!
! !INTERFACE:
!
  subroutine create_ges_grids(switch_on_derivatives,tendsflag)

! !USES:

    use constants,only: zero,one
    use gridmod, only: lat2,lon2,nsig,regional,nsig_soil
    use control_vectors, only: cvars3d
    use mpeu_util, only: die, tell, getindex
    implicit none

! !INPUT PARAMETERS:

    logical,intent(in   ) :: switch_on_derivatives    ! for for horizontal derivatives
    logical,intent(in   ) :: tendsflag                ! for time tendencies


! !OUTPUT PARAMETERS:

! !DESCRIPTION: allocate grids to hold guess and bias correction fields
!
! !REVISION HISTORY:
!   2004-06-03  treadon, original code
!   2004-07-15  todling, protex-compliant prologue; added onlys
!   2004-07-28  treadon - remove subroutine call list, pass variables via modules
!   2005-06-03  parrish - allocate/initialize _lat,_lon derivatives for u,v,cwmr,oz
!   2005-06-08  treadon - pass switch_on_derivatives via argument list
!   2005-07-06  parrish - add update_pint, arrays ges_pint, ges_pd
!   2005-07-27  kleist  - modified to include some shared arrays
!   2006-01-10  treadon - remove mype from calling list (not used)
!   2006-07-31  kleist  - use ges_ps arrays instead of ln(ps)
!   2006-06-08  zhang,b - change "biascor>0" to "biascor>=0" for debug purposes
!   2006-12-04  todling - remove bias initialization; rename routine
!   2006-12-15  todling - protection to allow initializing ges/tnd/drv at will
!   2007-03-15  todling - merged in da Silva/Cruz ESMF changes
!   2011-02-09  zhu     - add ges_gust,ges_vis,ges_pblh
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
!   treadon          org: w/nmc20     date: 2004-06-03
!
!EOP
!-------------------------------------------------------------------------

    integer(i_kind) i,j,k,n,ivar,istatus
    if(ges_grids_allocated_) call die('create_ges_grids','already allocated')
    ges_grids_allocated_=.true.

    if ( .not. ges_initialized ) then

#ifndef HAVE_ESMF
       nfldsig_all=nfldsig
       nfldsfc_all=nfldsfc
       nfldnst_all=nfldnst
       nfldsig_now=0 ! _now variables are not used if not for ESMF
       nfldsfc_now=0
       nfldnst_now=0
       extrap_intime=.true.
#endif /* HAVE_ESMF */

!      Allocate and zero guess grids
!       write(6,*)'nfldsig',nfldsig ! ajl
!       call flush(6) 
       allocate ( ges_prsi(lat2,lon2,nsig+1,nfldsig),ges_prsl(lat2,lon2,nsig,nfldsig),&
            ges_lnprsl(lat2,lon2,nsig,nfldsig),ges_lnprsi(lat2,lon2,nsig+1,nfldsig),&
            ges_tsen(lat2,lon2,nsig,nfldsig),&
            ges_teta(lat2,lon2,nsig,nfldsig),&
            geop_hgtl(lat2,lon2,nsig,nfldsig), &
            geop_hgti(lat2,lon2,nsig+1,nfldsig),ges_prslavg(nsig),&
            tropprs(lat2,lon2),fact_tv(lat2,lon2,nsig),&
            ges_qsat(lat2,lon2,nsig,nfldsig),stat=istatus)
       if (istatus/=0) write(6,*)'CREATE_GES_GRIDS(ges_prsi,..):  allocate error, istatus=',&
            istatus,lat2,lon2,nsig,nfldsig
#ifndef HAVE_ESMF
       allocate (ges_z(lat2,lon2,nfldsig),ges_ps(lat2,lon2,nfldsig),&
            ges_u(lat2,lon2,nsig,nfldsig),ges_v(lat2,lon2,nsig,nfldsig),&
            ges_vor(lat2,lon2,nsig,nfldsig),ges_div(lat2,lon2,nsig,nfldsig),&
            ges_q(lat2,lon2,nsig,nfldsig),&
            ges_oz(lat2,lon2,nsig,nfldsig),ges_tv(lat2,lon2,nsig,nfldsig),&
            ges_gust(lat2,lon2,nfldsig),ges_vis(lat2,lon2,nfldsig),&
            ges_pblh(lat2,lon2,nfldsig), &
            pbl_height(lat2,lon2,nfldsig),stat=istatus)
       if (istatus/=0) write(6,*)'CREATE_GES_GRIDS(ges_z,..):  allocate error, istatus=',&
            istatus,lat2,lon2,nsig,nfldsig
       allocate (ges_no2(lat2,lon2,nsig,nfldsig)) ! ajl
#endif /* HAVE_ESMF */
       if(update_pint) then
          allocate(ges_pint(lat2,lon2,nsig+1,nfldsig),ges_pd(lat2,lon2,nfldsig),&
               stat=istatus)
          if (istatus/=0) write(6,*)'CREATE_GES_GRIDS(ges_pint,..):  allocate error, istatus=',&
            istatus,lat2,lon2,nsig,nfldsig
       endif

       allocate ( ges_th2(lat2,lon2,nfldsig), ges_q2(lat2,lon2,nfldsig),&
         ges_soilt1(lat2,lon2,nfldsig),ges_tslb(lat2,lon2,nsig_soil,nfldsig),&
         ges_smois(lat2,lon2,nsig_soil,nfldsig), ges_tsk(lat2,lon2,nfldsig),&
         stat=istatus)
       if (istatus/=0) write(6,*)'CREATE_GES_GRIDS(ges_th2,..):  allocate error, istatus=',&
            istatus,lat2,lon2,nsig,nfldsig

       ges_initialized = .true.

!  Default for ges_psfcavg
       ges_psfcavg=zero
       do i=1,nsig
          ges_prslavg(i)=zero
       end do

       do j=1,lon2
          do i=1,lat2
             tropprs(i,j)=zero
          end do
       end do

       do k=1,nsig
          do j=1,lon2
             do i=1,lat2
                fact_tv(i,j,k)=one
             end do
          end do
       end do

#ifndef HAVE_ESMF
       do n=1,nfldsig
          do j=1,lon2
             do i=1,lat2
                ges_z(i,j,n)=zero
                ges_ps(i,j,n)=zero
                ges_gust(i,j,n)=zero
                ges_vis(i,j,n)=zero
                ges_pblh(i,j,n)=zero
                pbl_height(i,j,n)=zero
             end do
          end do
       end do
       do n=1,nfldsig
          do k=1,nsig
             do j=1,lon2
                do i=1,lat2
                   ges_u(i,j,k,n)=zero
                   ges_v(i,j,k,n)=zero
                   ges_vor(i,j,k,n)=zero
                   ges_div(i,j,k,n)=zero
                   ges_q(i,j,k,n)=zero
                   ges_oz(i,j,k,n)=zero
                   ges_tv(i,j,k,n)=zero
                   ges_qsat(i,j,k,n)=zero
!                  ges_pint(i,j,k,n)=zero
                   ges_no2(i,j,k,n)=zero ! ajl
                end do
             end do
          end do
       end do
#endif /* HAVE_ESMF */
       do n=1,nfldsig
          do k=1,nsig
             do j=1,lon2
                do i=1,lat2
                   ges_prsl(i,j,k,n)=zero
                   ges_lnprsl(i,j,k,n)=zero
                   ges_tsen(i,j,k,n)=zero
                   ges_teta(i,j,k,n)=zero
                   geop_hgtl(i,j,k,n)=zero
                end do
             end do
          end do
          do k=1,nsig+1
             do j=1,lon2
                do i=1,lat2
                   ges_prsi(i,j,k,n)=zero
                   ges_lnprsi(i,j,k,n)=zero
                   geop_hgti(i,j,k,n)=zero
                end do
             end do
          end do
       end do
       if(update_pint) then
          do n=1,nfldsig
             do k=1,nsig+1
                do j=1,lon2
                   do i=1,lat2
                      ges_pint(i,j,k,n)=zero
                   end do
                end do
             end do
             do j=1,lon2
                do i=1,lat2
                   ges_pd(i,j,n)=zero
                end do
             end do
          end do
       end if

       allocate (efr_ql(lat2,lon2,nsig,nfldsig),efr_qi(lat2,lon2,nsig,nfldsig), &
                 efr_qr(lat2,lon2,nsig,nfldsig),efr_qs(lat2,lon2,nsig,nfldsig), &
                 efr_qg(lat2,lon2,nsig,nfldsig),efr_qh(lat2,lon2,nsig,nfldsig))
       do n=1,nfldsig
          do k=1,nsig
             do j=1,lon2
                do i=1,lat2
                   efr_ql(i,j,k,n)=zero
                   efr_qi(i,j,k,n)=zero
                   efr_qr(i,j,k,n)=zero
                   efr_qs(i,j,k,n)=zero
                   efr_qg(i,j,k,n)=zero
                   efr_qh(i,j,k,n)=zero
                end do
             end do
          end do
       end do

! for GSD  soil nudging
       do n=1,nfldsig
          do k=1,nsig_soil
             do j=1,lon2
                do i=1,lat2
                   ges_tslb(i,j,k,n)=zero
                   ges_smois(i,j,k,n)=zero
                end do
             end do
          end do
          do j=1,lon2
             do i=1,lat2
                ges_th2(i,j,n)=zero
                ges_q2(i,j,n)=zero
                ges_tsk(i,j,n)=zero
                ges_soilt1(i,j,n)=zero
             end do
          end do
       end do

    end if ! ges_initialized
    
!   If tendencies option on, allocate/initialize _ten arrays to zero
    if (.not.tnd_initialized .and. tendsflag) then
       allocate(ges_prs_ten(lat2,lon2,nsig+1),ges_u_ten(lat2,lon2,nsig),&
                ges_v_ten(lat2,lon2,nsig),ges_tv_ten(lat2,lon2,nsig),&
                ges_q_ten(lat2,lon2,nsig),ges_oz_ten(lat2,lon2,nsig),&
                stat=istatus)
       if (istatus/=0) write(6,*)'CREATE_GES_GRIDS(ges_prs_ten,..):  allocate error, istatus=',&
            istatus,lat2,lon2,nsig
       allocate(ges_no2_ten(lat2,lon2,nsig)) ! ajl
       tnd_initialized = .true.
       do k=1,nsig
          do j=1,lon2
             do i=1,lat2
                ges_u_ten(i,j,k)=zero
                ges_v_ten(i,j,k)=zero
                ges_tv_ten(i,j,k)=zero
                ges_q_ten(i,j,k)=zero
                ges_oz_ten(i,j,k)=zero
                ges_prs_ten(i,j,k)=zero
                ges_no2_ten(i,j,k)=zero ! ajl
             end do
          end do
       end do
!      Get pointer to could water mixing ratio, and alloc tendency if cwmr present in guess
!      call gsi_metguess_get ( 'var::cw', ivar, istatus )
!      if (regional .and. nems_nmmb_regional) ivar=getindex(cvars3d,'cw')
       ivar=getindex(cvars3d,'cw')
       if (ivar>0) then
           allocate(ges_cwmr_ten(lat2,lon2,nsig),stat=istatus)
           if (istatus/=0) write(6,*)'CREATE_GES_GRIDS(ges_cwmr_ten):  allocate error, istatus=',&
                istatus,lat2,lon2,nsig
            do k=1,nsig
               do j=1,lon2
                  do i=1,lat2
                     ges_cwmr_ten(i,j,k)=zero
                  end do
               end do
            end do
       endif

       do j=1,lon2
          do i=1,lat2
             ges_prs_ten(i,j,nsig+1)=zero
          end do
       end do
    end if

!   If derivatives option on, allocate and initialize derivatives arrays to 0.0
    if (.not.drv_initialized .and. switch_on_derivatives) then
       allocate(ges_u_lat(lat2,lon2,nsig),ges_u_lon(lat2,lon2,nsig),&
            ges_v_lat(lat2,lon2,nsig),ges_v_lon(lat2,lon2,nsig),&
            ges_ozlat(lat2,lon2,nsig),ges_ozlon(lat2,lon2,nsig),&
            ges_ps_lat(lat2,lon2,nfldsig),ges_ps_lon(lat2,lon2,nfldsig),&
            ges_tvlat(lat2,lon2,nsig),ges_tvlon(lat2,lon2,nsig),&
            ges_qlat(lat2,lon2,nsig),ges_qlon(lat2,lon2,nsig),&
            stat=istatus)
       if (istatus/=0) write(6,*)'CREATE_GES_GRIDS(ges_u_lat,..):  allocate error, istatus=',&
            istatus,lat2,lon2,nsig,nfldsig
       allocate(ges_no2lat(lat2,lon2,nsig),ges_no2lon(lat2,lon2,nsig))
       drv_initialized = .true.
       do k=1,nsig
          do j=1,lon2
             do i=1,lat2
                ges_u_lat(i,j,k)=zero
                ges_u_lon(i,j,k)=zero
                ges_v_lat(i,j,k)=zero
                ges_v_lon(i,j,k)=zero
                ges_ozlat(i,j,k)=zero
                ges_ozlon(i,j,k)=zero
                ges_tvlat(i,j,k)=zero
                ges_tvlon(i,j,k)=zero
                ges_qlat(i,j,k)=zero
                ges_qlon(i,j,k)=zero
                ges_no2lat(i,j,k)=zero ! ajl
                ges_no2lon(i,j,k)=zero
             end do
          end do
       end do
!      Get pointer to could water mixing ratio, and alloc directional derivatives if cwmr present in guess
!      call gsi_metguess_get ( 'var::cw', ivar, istatus )
!      if (regional .and. nems_nmmb_regional) ivar=getindex(cvars3d,'cw')
       ivar=getindex(cvars3d,'cw')
       if (ivar>0) then
           allocate(ges_cwmr_lat(lat2,lon2,nsig),ges_cwmr_lon(lat2,lon2,nsig),&
            stat=istatus)
            if (istatus/=0) write(6,*)'CREATE_GES_GRIDS(ges_cwmr_lat,..):  allocate error, istatus=',&
                istatus,lat2,lon2,nsig,nfldsig
            do k=1,nsig
               do j=1,lon2
                  do i=1,lat2
                     ges_cwmr_lat(i,j,k)=zero
                     ges_cwmr_lon(i,j,k)=zero
                  end do
               end do
            end do
       endif
       do n=1,nfldsig
          do j=1,lon2
             do i=1,lat2
                ges_ps_lat(i,j,n)=zero
                ges_ps_lon(i,j,n)=zero
             end do
          end do
       end do
    endif  ! end if switch_derivatives block

    return
  end subroutine create_ges_grids

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: create_metguess_grids --- initialize extra meterological guess
!
! !INTERFACE:
!
  subroutine create_metguess_grids(istatus)

! !USES:
  use gridmod, only: lat2,lon2,nsig
  implicit none

! !OUTPUT PARAMETERS:

  integer(i_kind), intent(out) :: istatus

! !DESCRIPTION: initialize extra meteorological background fields beyond 
!               the standard ones - wired-in this module.
!
! !REVISION HISTORY:
!   2011-04-29  todling
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; Linux Cluster
!
! !AUTHOR: 
!   todling         org: w/nmc20     date: 2011-04-29
!
!EOP
!-------------------------------------------------------------------------
   character(len=*),parameter::myname_='create_metguess_grids'
   integer(i_kind) :: nmguess                   ! number of meteorol. fields (namelist)
   character(len=256),allocatable:: mguess(:)   ! names of meterol. fields

   istatus=0
  
!  When proper connection to ESMF is complete,
!  the following will not be needed here
!  ------------------------------------------
   call gsi_metguess_get('dim',nmguess,istatus)
   if(istatus/=0) then
      write(6,*) myname_, ': trouble getting number of met-guess fields'
      return
   endif
   if(nmguess==0) return
   if (nmguess>0) then
       allocate (mguess(nmguess))
       call gsi_metguess_get('gsinames',mguess,istatus)
       if(istatus/=0) then
          write(6,*) myname_, ': trouble getting name of met-guess fields'
          return
       endif

!      Allocate memory for guess files for trace gases
!      ------------------------------------------------
       call gsi_metguess_create_grids(lat2,lon2,nsig,nfldsig,istatus)
       if(istatus/=0) then
          write(6,*) myname_, ': trouble allocating mem for extra met-guess'
          return
       endif
   endif

  end subroutine create_metguess_grids

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: destroy_metguess_grids --- destroy extra meterological background
!
! !INTERFACE:
!
  subroutine destroy_metguess_grids(istatus)
! !USES:
  implicit none
! !OUTPUT PARAMETERS:
  integer(i_kind),intent(out)::istatus
! !DESCRIPTION: destroy extra meterological background
!
! !REVISION HISTORY:
!   2011-04-29  todling
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; Linux Cluster
!
! !AUTHOR: 
!   todling         org: w/nmc20     date: 2011-04-29
!
!EOP
  istatus=0
  call gsi_metguess_destroy_grids(istatus)
  end subroutine destroy_metguess_grids

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: create_chemges_grids --- initialize chem component
!
! !INTERFACE:
!
  subroutine create_chemges_grids(istatus)

! !USES:
  use gridmod, only: lat2,lon2,nsig
  implicit none

! !OUTPUT PARAMETERS:

  integer(i_kind), intent(out) :: istatus

! !DESCRIPTION: initialize chem background
!
! !REVISION HISTORY:
!   2010-05-19  todling
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; Linux Cluster
!
! !AUTHOR: 
!   todling         org: w/nmc20     date: 2010-05-19
!
!EOP
!-------------------------------------------------------------------------
  character(len=*),parameter::myname_='create_chemges_grids'
   integer(i_kind) :: ntgases                   ! number of tracer gases (namelist)
   character(len=256),allocatable:: tgases(:)   ! names of tracer gases

  istatus=0
  
!  When proper connection to ESMF is complete,
!  the following will not be needed here
!  ------------------------------------------
   call gsi_chemguess_get('dim',ntgases,istatus)
   if(istatus/=0) then
      write(6,*) myname_, ': trouble getting number of chem/gases'
      return
   endif
   if(ntgases==0) return
   if (ntgases>0) then
      allocate (tgases(ntgases))
      call gsi_chemguess_get('gsinames',tgases,istatus)
      if(istatus/=0) then
         write(6,*) myname_, ': trouble getting name of chem/gases'
         return
      endif

!     Allocate memory for guess files for trace gases
!     ------------------------------------------------
      call gsi_chemguess_create_grids(lat2,lon2,nsig,nfldsig,istatus)
      if(istatus/=0) then
         write(6,*) myname_, ': trouble allocating mem for chem/gases'
         return
      endif
   endif

  end subroutine create_chemges_grids

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: destroy_chemges_grids --- destroy chem background
!
! !INTERFACE:
!
  subroutine destroy_chemges_grids(istatus)
! !USES:
  implicit none
! !OUTPUT PARAMETERS:
  integer(i_kind),intent(out)::istatus
! !DESCRIPTION: destroy chem background
!
! !REVISION HISTORY:
!   2010-05-19  todling
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; Linux Cluster
!
! !AUTHOR: 
!   todling         org: w/nmc20     date: 2010-05-19
!
!EOP
  istatus=0
  call gsi_chemguess_destroy_grids(istatus)
  end subroutine destroy_chemges_grids

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: destroy_ges_grids --- Dealloc guess and bias fields
!
! !INTERFACE:
!
  subroutine destroy_ges_grids(switch_on_derivatives,tendsflag)

! !USES:

    use mpeu_util, only: die, tell,getindex
    use control_vectors, only: cvars3d
    implicit none

! !INPUT PARAMETERS:
    logical,intent(in   ) :: switch_on_derivatives    ! flag for horizontal derivatives
    logical,intent(in   ) :: tendsflag                ! flag for tendency
    
! !DESCRIPTION: deallocate guess and bias grids
!
! !REVISION HISTORY:
!   2003-12-01  kleist
!   2004-05-14  kleist, documentation
!   2004-07-15  todling, protex-compliant prologue; added onlys
!   2005-06-03  parrish - deallocate _lat,_lon arrays for u,v,cwmr,oz
!   2005-06-08  treadon - check flag to see if need to deallocate derivatives
!   2005-07-06  parrish - add update_pint, arrays ges_pint, ges_pd
!   2005-07-27  kleist  - modified to include some shared arrays
!   2006-07-31  kleist  - use ges_ps arrays instead of ln(ps)
!   2006-12-04  todling - remove bias destroy; rename routine
!   2006-12-15  todling - using internal switches to deallc(tnds/drvs)
!   2007-03-15  todling - merged in da Silva/Cruz ESMF changes
!   2011-11-01  eliu    - add access to cvars3d & getindex 
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR: 
!   kleist          org: w/nmc20     date: 2003-12-01
!
!EOP
!-------------------------------------------------------------------------
    integer(i_kind):: ivar,istatus

    deallocate(ges_prsi,ges_prsl,ges_lnprsl,ges_lnprsi,&
         ges_tsen,ges_teta,geop_hgtl,geop_hgti,ges_prslavg,&
         tropprs,fact_tv,ges_qsat,stat=istatus)
    if (istatus/=0) &
         write(6,*)'DESTROY_GES_GRIDS(ges_prsi,..):  deallocate error, istatus=',&
         istatus
#ifndef HAVE_ESMF
    deallocate(ges_z,ges_ps,&
         ges_u,ges_v,ges_vor,ges_div,ges_q,&
         ges_oz,ges_tv,&
         stat=istatus)
    if (istatus/=0) &
         write(6,*)'DESTROY_GES_GRIDS(ges_z,..):  deallocate error, istatus=',&
         istatus
       deallocate(pbl_height,stat=istatus)
    if (istatus/=0) &
         write(6,*)'DESTROY_GES_GRIDS(pbl_height,..):  deallocate error, istatus=',&
         istatus
    deallocate(ges_no2)
#endif /* HAVE_ESMF */
    if(update_pint) then
       deallocate(ges_pint,ges_pd,stat=istatus)
       if (istatus/=0) &
            write(6,*)'DESTROY_GES_GRIDS(ges_pint,..):  deallocate error, istatus=',&
            istatus
    endif
    deallocate(efr_ql,efr_qi,efr_qr,efr_qs,efr_qg,efr_qh)
! GSD soil nudging
    deallocate(ges_th2,ges_q2,ges_soilt1,ges_tslb,ges_smois,ges_tsk,stat=istatus)
!
    if (drv_initialized .and.switch_on_derivatives) then
!      Get pointer to could water mixing ratio, and alloc tendency if cwmr present in guess
!      call gsi_metguess_get ( 'var::cw', ivar, istatus )
       ivar=getindex(cvars3d,'cw')
       if (ivar>0) then
           deallocate(ges_cwmr_lat,ges_cwmr_lon,&
            stat=istatus)
            if (istatus/=0) &
                 write(6,*)'DESTROY_GES_GRIDS(ges_cwmr_lat,..):  deallocate error, istatus=',&
                 istatus
       endif
       deallocate(ges_u_lat,ges_u_lon,ges_v_lat,ges_v_lon,&
            ges_ozlat,ges_ozlon,&
            ges_ps_lat,ges_ps_lon,ges_tvlat,ges_tvlon,&
            ges_qlat,ges_qlon,stat=istatus)
       if (istatus/=0) &
            write(6,*)'DESTROY_GES_GRIDS(ges_u_lat,..):  deallocate error, istatus=',&
            istatus
       deallocate(ges_no2lat,ges_no2lon)
    endif
    if (tnd_initialized .and. tendsflag) then
!      Get pointer to could water mixing ratio, and alloc tendency if cwmr present in guess
!      call gsi_metguess_get ( 'var::cw', ivar, istatus )
       ivar=getindex(cvars3d,'cw')
       if (ivar>0) then
           deallocate(ges_cwmr_ten,stat=istatus)
           if (istatus/=0) &
                write(6,*)'DESTROY_GES_GRIDS(ges_cwmr_ten,..):  deallocate error, istatus=',&
                istatus
       endif
       deallocate(ges_u_ten,ges_v_ten,ges_tv_ten,ges_prs_ten,ges_q_ten,&
            ges_oz_ten,stat=istatus)
       if (istatus/=0) &
            write(6,*)'DESTROY_GES_GRIDS(ges_u_ten,..):  deallocate error, istatus=',&
            istatus
       deallocate(ges_no2_ten)
    endif
    return
  end subroutine destroy_ges_grids

#ifdef TO_BE_REMOVED
!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: destroy_sfct --- Deallocate sfct only
!
! !INTERFACE:
!
  subroutine destroy_sfct

! !USES:

   use mpeu_util, only: die, tell
   implicit none

! !DESCRIPTION: deallocate surface temperature field
!
! !REVISION HISTORY:
!   2008-06-30  derber
!   2008-09-05  lueken - add subprogram doc block
!   2009-01-02  todling - replaced doc block with protex-prologue
!   2009-01-17  todling - dealloc isli2,sno2 was misplaced
!
! !REMARKS:
!   language: f90
!   machine:ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR: 
!   derber          org: w/nmc2     date: 2008-06-30
!
!EOP
!-------------------------------------------------------------------------

    integer(i_kind):: istatus

    deallocate(isli2,sno2,stat=istatus)
    if (istatus/=0) &
         write(6,*)'DESTROY_SFCT:  deallocate error, istatus=',&
         istatus
#ifndef HAVE_ESMF
    deallocate(sfct,dsfct,stat=istatus)
    if (istatus/=0) &
         write(6,*)'DESTROY_SFCT:  deallocate error, istatus=',&
         istatus
#endif /* HAVE_ESMF */

    return
  end subroutine destroy_sfct
#endif

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: destroy_sfc_grids --- Deallocate surface fields
!
! !INTERFACE:
!
  subroutine destroy_sfc_grids

! !USES:

   use mpeu_util, only: die, tell
   implicit none
   
! !DESCRIPTION: deallocate surface related grids
!
! !REVISION HISTORY:
!   2003-12-01  kleist
!   2004-05-14  kleist, documentation
!   2004-07-15  todling, protex-compliant prologue
!   2005-06-03  parrish - deallocate sfct_lat and sfct_lon
!   2007-03-15  todling - merged in da Silva/Cruz ESMF changes
!   2008-06-30  derber - remove sfct deallocate to allow earlier call
!   2009-01-17  todling - move isli2,sno2 into destroy_sfct
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR: 
!   kleist          org: w/nmc20     date: 2003-12-01
!
!EOP
!-------------------------------------------------------------------------

    integer(i_kind):: istatus
    if(.not.sfc_grids_allocated_) call die('destroy_sfc_grids_','not allocated')
    sfc_grids_allocated_=.false.

    deallocate(isli_g,isli2,sno2,stat=istatus)
    if (istatus/=0) &
         write(6,*)'DESTROY_SFC_GRIDS:  deallocate error, istatus=',&
         istatus
#ifndef HAVE_ESMF
    deallocate(isli,fact10,dsfct,sfct,sno,veg_type,veg_frac,soil_type,&
         sfc_rough,soil_temp,soil_moi,stat=istatus)
    if (istatus/=0) &
         write(6,*)'DESTROY_SFC_GRIDS:  deallocate error, istatus=',&
         istatus
#endif /* HAVE_ESMF */

    return
  end subroutine destroy_sfc_grids

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: create_gesfinfo --- Allocate guess-files information arrays
!
! !INTERFACE:
!
  subroutine create_gesfinfo

! !USES:

   use mpeu_util, only: die, tell
   implicit none

! !DESCRIPTION: allocate guess-files information arrays
!
! !REVISION HISTORY:
!   2009-01-08  todling
!
! !REMARKS:
!   language: f90
!   machine:ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR: 
!   todling          org: w/nmc2     date: 2009-01-08
!
!EOP
!-------------------------------------------------------------------------

    integer(i_kind):: istatus
    if(gesfinfo_created_) call die('create_gesfinfo','already created')
    gesfinfo_created_=.true.

#ifndef HAVE_ESMF
    nfldsig_all=nfldsig
    nfldsfc_all=nfldsfc
    nfldnst_all=nfldnst
    nfldsig_now=0	! _now variables are not used if not for ESMF
    nfldsfc_now=0
    nfldnst_now=0
    extrap_intime=.true.
    allocate(hrdifsfc(nfldsfc),ifilesfc(nfldsfc), &
             hrdifnst(nfldnst),ifilenst(nfldnst), &
             hrdifsig(nfldsig),ifilesig(nfldsig), &
	     hrdifsfc_all(nfldsfc_all), &
             hrdifnst_all(nfldnst_all), &
	     hrdifsig_all(nfldsig_all), &
	     stat=istatus)
    if (istatus/=0) &
         write(6,*)'CREATE_GESFINFO(hrdifsfc,..):  allocate error, istatus=',&
         istatus
#endif /* HAVE_ESMF */

    return
  end subroutine create_gesfinfo

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: destroy_gesfinfo --- Deallocate guess-files information
!
! !INTERFACE:
!
  subroutine destroy_gesfinfo

! !USES:

   use mpeu_util, only: die
   implicit none

! !DESCRIPTION: deallocate guess-files information
!
! !REVISION HISTORY:
!   2009-01-08  todling
!
! !REMARKS:
!   language: f90
!   machine:ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR: 
!   todling          org: w/nmc2     date: 2009-01-08
!
!EOP
!-------------------------------------------------------------------------

    integer(i_kind):: istatus
    if(.not.gesfinfo_created_) call die('destroy_gesfinfo','not created')
    gesfinfo_created_=.false.

#ifndef HAVE_ESMF
    deallocate(hrdifsfc,ifilesfc,hrdifnst,ifilenst,hrdifsig,ifilesig, &
    	hrdifsfc_all,hrdifnst_all,hrdifsig_all,stat=istatus)
    if (istatus/=0) &
         write(6,*)'DESTROY_GESFINFO:  deallocate error, istatus=',&
         istatus

    nfldsfc_all=0
    nfldnst_all=0
    nfldsig_all=0
    nfldsfc    =0
    nfldnst    =0
    nfldsig    =0
#endif /* HAVE_ESMF */

    return
  end subroutine destroy_gesfinfo


!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: load_prsges --- Populate guess pressure arrays
!
! !INTERFACE:
!
  subroutine load_prsges

! !USES:

    use constants,only: zero,one,rd_over_cp,one_tenth,half,ten
    use gridmod, only: lat2,lon2,nsig,ak5,bk5,ck5,tref5,idvc5,&
         regional,wrf_nmm_regional,nems_nmmb_regional,wrf_mass_regional,&
         cmaq_regional,pt_ll,aeta2_ll,&
         aeta1_ll,eta2_ll,pdtop_ll,eta1_ll,twodvar_regional,idsl5
    use gridmod, only : raqms ! ajl
    use raqmsmod, only : raqpsol,raqpdash ! ajl
    implicit none

! !DESCRIPTION: populate guess pressure arrays
!
! !REVISION HISTORY:
!   2003-10-15  kleist
!   2004-03-22  parrish, regional capability added
!   2004-05-14  kleist, documentation
!   2004-07-15  todling, protex-compliant prologue; added onlys
!   2004-07-28  treadon - remove subroutine call list, pass variables via modules
!   2005-05-24  pondeca - add regional surface analysis option
!   2006-04-14  treadon - unify global calculations to use ak5,bk5
!   2006-04-17  treadon - add ges_psfcavg and ges_prslavg for regional
!   2006-07-31  kleist  - use ges_ps instead of ln(ps)
!   2007-05-08  kleist  - add fully generalized coordinate for pressure calculation
!   2011-07-07  todling - add cap for log(pressure) calculation
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
!   kleist          org: w/nmc20     date: 2003-10-15
!
!EOP
!-------------------------------------------------------------------------

!   Declare local parameter
    real(r_kind),parameter:: r1013=1013.0_r_kind

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

    kap1=rd_over_cp+one
    kapr=one/rd_over_cp

    do jj=1,nfldsig
       do k=1,nsig+1
          do j=1,lon2
             do i=1,lat2
                if(regional) then
                   if (wrf_nmm_regional.or.nems_nmmb_regional.or.&
                        cmaq_regional ) &
                      ges_prsi(i,j,k,jj)=one_tenth* &
                             (eta1_ll(k)*pdtop_ll + &
                              eta2_ll(k)*(ten*ges_ps(i,j,jj)-pdtop_ll-pt_ll) + &
                              pt_ll)

                   if (wrf_mass_regional .or. twodvar_regional) &
                      ges_prsi(i,j,k,jj)=one_tenth*(eta1_ll(k)*(ten*ges_ps(i,j,jj)-pt_ll) + pt_ll)
                elseif  (raqms)then ! ajl
                  ges_prsi(i,j,k,jj)=raqpsol(i,j,k,jj) ! ajl
!                  end ajl
                else
                   if (idvc5==1 .or. idvc5==2) then
                      ges_prsi(i,j,k,jj)=ak5(k)+(bk5(k)*ges_ps(i,j,jj))
                   else if (idvc5==3) then
                      if (k==1) then
                         ges_prsi(i,j,k,jj)=ges_ps(i,j,jj)
                      else if (k==nsig+1) then
                         ges_prsi(i,j,k,jj)=zero
                      else
                         trk=(half*(ges_tv(i,j,k-1,jj)+ges_tv(i,j,k,jj))/tref5(k))**kapr
                         ges_prsi(i,j,k,jj)=ak5(k)+(bk5(k)*ges_ps(i,j,jj))+(ck5(k)*trk)
                      end if
                   end if
                endif
                ges_prsi(i,j,k,jj)=max(ges_prsi(i,j,k,jj),zero)
                ges_lnprsi(i,j,k,jj)=log(max(ges_prsi(i,j,k,jj),0.0001_r_kind))
             end do
          end do
       end do
    end do

    if(regional) then
       if (wrf_nmm_regional.or.nems_nmmb_regional.or.cmaq_regional) then
! load using aeta coefficients
          do jj=1,nfldsig
             do k=1,nsig
                do j=1,lon2
                   do i=1,lat2
                      ges_prsl(i,j,k,jj)=one_tenth* &
                                  (aeta1_ll(k)*pdtop_ll + &
                                   aeta2_ll(k)*(ten*ges_ps(i,j,jj)-pdtop_ll-pt_ll) + &
                                   pt_ll)
                      ges_lnprsl(i,j,k,jj)=log(ges_prsl(i,j,k,jj))

                   end do
                end do
             end do
          end do
       end if   ! end if wrf_nmm regional block
       if (wrf_mass_regional .or. twodvar_regional) then
! load using aeta coefficients
          do jj=1,nfldsig
             do k=1,nsig
                do j=1,lon2
                   do i=1,lat2
                      ges_prsl(i,j,k,jj)=one_tenth*(aeta1_ll(k)*(ten*ges_ps(i,j,jj)-pt_ll)+pt_ll)
                      ges_lnprsl(i,j,k,jj)=log(ges_prsl(i,j,k,jj))
                   end do
                end do
             end do
          end do
       end if   ! end if wrf_mass regional block


    elseif  (raqms)then ! ajl
      do jj=1,nfldsig
        do k=1,nsig
          do j=1,lon2
            do i=1,lat2
              ges_prsl(i,j,k,jj)=raqpdash(i,j,k,jj)
              ges_lnprsl(i,j,k,jj)=log(ges_prsl(i,j,k,jj))
            end do
          end do
        end do
      end do
     
    else
!      load mid-layer pressure by using phillips vertical interpolation
       if (idsl5/=2) then
          do jj=1,nfldsig
             do j=1,lon2
                do i=1,lat2
                   do k=1,nsig
                      ges_prsl(i,j,k,jj)=((ges_prsi(i,j,k,jj)**kap1-ges_prsi(i,j,k+1,jj)**kap1)/&
                           (kap1*(ges_prsi(i,j,k,jj)-ges_prsi(i,j,k+1,jj))))**kapr
                      ges_lnprsl(i,j,k,jj)=log(ges_prsl(i,j,k,jj))
                   end do
                end do
             end do
          end do

!      load mid-layer pressure by simple averaging
       else
          do jj=1,nfldsig
             do j=1,lon2
                do i=1,lat2
                   do k=1,nsig
                      ges_prsl(i,j,k,jj)=(ges_prsi(i,j,k,jj)+ges_prsi(i,j,k+1,jj))*half
                      ges_lnprsl(i,j,k,jj)=log(ges_prsl(i,j,k,jj))
                   end do
                end do
             end do
          end do
       endif

    end if  !  end regional/global block

! For regional applications only, load variables containing mean
! surface pressure and pressure profile at the layer midpoints
    if (regional) then
       ges_psfcavg = r1013
       if (wrf_nmm_regional.or.nems_nmmb_regional.or.cmaq_regional) then
          do k=1,nsig
             ges_prslavg(k)=aeta1_ll(k)*pdtop_ll+aeta2_ll(k)*(r1013-pdtop_ll-pt_ll)+pt_ll
          end do
       else
          do k=1,nsig
             ges_prslavg(k)=aeta1_ll(k)*(r1013-pt_ll)+pt_ll
          end do
       endif
    endif


    return
  end subroutine load_prsges

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: load_geop_hgt --- Populate guess geopotential height
!
! !INTERFACE:
!
  subroutine load_geop_hgt

! !USES:

    use constants, only: one,eps, rd, grav, half, t0c, fv
    use constants, only: cpf_a0, cpf_a1, cpf_a2, cpf_b0, cpf_b1, cpf_c0, cpf_c1, cpf_d, cpf_e
    use constants, only: psv_a, psv_b, psv_c, psv_d
    use constants, only: ef_alpha, ef_beta, ef_gamma
    use gridmod, only: lat2, lon2, nsig, twodvar_regional

    implicit none

! !INPUT PARAMETERS:


! !DESCRIPTION: populate guess geopotential height
!
! !REVISION HISTORY:
!   2003-10-15  treadon
!   2004-05-14  kleist, documentation
!   2004-07-15  todling, protex-compliant prologue
!   2004-10-28  treadon - replace "tiny" with "tiny_r_kind"
!   2004-12-15  treadon - replace use of Paul van Delst's Geopotential
!                         function with simple integration of hydrostatic
!                         equation (done to be consistent with Lidia
!                         Cucurull's GPS work)
!   2005-05-24  pondeca - add regional surface analysis option
!   2010-08-27  cucurull - add option to compute and use compressibility factors in geopot heights
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
!   treadon          org: w/nmc20      date: 2003-10-15
!
!EOP
!-------------------------------------------------------------------------

    integer(i_kind) i,j,k,jj
    real(r_kind) h,dz,rdog
    real(r_kind),dimension(nsig+1):: height
    real(r_kind) cmpr, x_v, rl_hm, fact, pw, tmp_K, tmp_C, prs_sv, prs_a, ehn_fct, prs_v
    real(r_kind),parameter:: thousand = 1000.0_r_kind

    if (twodvar_regional) return

    rdog = rd/grav

    if (use_compress) then

!     Compute compressibility factor (Picard et al 2008) and geopotential heights at midpoint 
!     of each layer

       do jj=1,nfldsig
          do j=1,lon2
             do i=1,lat2
                k  = 1
                fact    = one + fv * ges_q(i,j,k,jj)
                pw      = eps + ges_q(i,j,k,jj)*( one - eps )
                tmp_K   = ges_tv(i,j,k,jj) / fact
                tmp_C   = tmp_K - t0c
                prs_sv  = exp(psv_a*tmp_K**2 + psv_b*tmp_K + psv_c + psv_d/tmp_K)  ! Pvap sat, eq A1.1 (Pa)
                prs_a   = thousand * exp(half*(log(ges_prsi(i,j,k,jj)) + log(ges_prsl(i,j,k,jj))))     ! (Pa) 
                ehn_fct = ef_alpha + ef_beta*prs_a + ef_gamma*tmp_C**2 ! enhancement factor (eq. A1.2)
                prs_v   = ges_q(i,j,k,jj) * prs_a / pw   ! vapor pressure (Pa)
                rl_hm   = prs_v / prs_sv    ! relative humidity
                x_v     = rl_hm * ehn_fct * prs_sv / prs_a     ! molar fraction of water vapor (eq. A1.3)
 
                ! Compressibility factor (eq A1.4 from Picard et al 2008)
                cmpr = one - (prs_a/tmp_K) * (cpf_a0 + cpf_a1*tmp_C + cpf_a2*tmp_C**2 &
                           + (cpf_b0 + cpf_b1*tmp_C)*x_v + (cpf_c0 + cpf_c1*tmp_C)*x_v**2 ) &
                           + (prs_a**2/tmp_K**2) * (cpf_d + cpf_e*x_v**2)

                h  = rdog * ges_tv(i,j,k,jj)
                dz = h * cmpr * log(ges_prsi(i,j,k,jj)/ges_prsl(i,j,k,jj))
                height(k) = ges_z(i,j,jj) + dz   

                do k=2,nsig
                   fact    = one + fv * half * (ges_q(i,j,k-1,jj)+ges_q(i,j,k,jj))
                   pw      = eps + half * (ges_q(i,j,k-1,jj)+ges_q(i,j,k,jj)) * (one - eps)
                   tmp_K   = half * (ges_tv(i,j,k-1,jj)+ges_tv(i,j,k,jj)) / fact
                   tmp_C   = tmp_K - t0c
                   prs_sv  = exp(psv_a*tmp_K**2 + psv_b*tmp_K + psv_c + psv_d/tmp_K)  ! eq A1.1 (Pa)
                   prs_a   = thousand * exp(half*(log(ges_prsl(i,j,k-1,jj))+log(ges_prsl(i,j,k,jj))))   ! (Pa)
                   ehn_fct = ef_alpha + ef_beta*prs_a + ef_gamma*tmp_C**2 ! enhancement factor (eq. A1.2)
                   prs_v   = half*(ges_q(i,j,k-1,jj)+ges_q(i,j,k,jj) ) * prs_a / pw   ! (Pa)
                   rl_hm   = prs_v / prs_sv    ! relative humidity
                   x_v     = rl_hm * ehn_fct * prs_sv / prs_a     ! molar fraction of water vapor (eq. A1.3)
                   cmpr    = one - (prs_a/tmp_K) * ( cpf_a0 + cpf_a1*tmp_C + cpf_a2*tmp_C**2 &
                             + (cpf_b0 + cpf_b1*tmp_C)*x_v + (cpf_c0 + cpf_c1*tmp_C)*x_v**2 ) &
                             + (prs_a**2/tmp_K**2) * (cpf_d + cpf_e*x_v**2)
                   h       = rdog * half * (ges_tv(i,j,k-1,jj)+ges_tv(i,j,k,jj))
                   dz      = h * cmpr * log(ges_prsl(i,j,k-1,jj)/ges_prsl(i,j,k,jj))
                   height(k) = height(k-1) + dz
                end do

                do k=1,nsig
                   geop_hgtl(i,j,k,jj)=height(k) - ges_z(i,j,jj)
                end do
             enddo
          enddo
       enddo

!      Compute compressibility factor (Picard et al 2008) and geopotential heights at interface
!      between layers

       do jj=1,nfldsig
          do j=1,lon2
             do i=1,lat2
                k=1
                height(k) = ges_z(i,j,jj)

                do k=2,nsig
                   fact    = one + fv * ges_q(i,j,k-1,jj)
                   pw      = eps + ges_q(i,j,k-1,jj)*(one - eps)
                   tmp_K   = ges_tv(i,j,k-1,jj) / fact
                   tmp_C   = tmp_K - t0c
                   prs_sv  = exp(psv_a*tmp_K**2 + psv_b*tmp_K + psv_c + psv_d/tmp_K)  ! eq A1.1 (Pa)
                   prs_a   = thousand * exp(half*(log(ges_prsi(i,j,k-1,jj))+log(ges_prsi(i,j,k,jj)))) 
                   ehn_fct = ef_alpha + ef_beta*prs_a + ef_gamma*tmp_C**2 ! enhancement factor (eq. A1.2)
                   prs_v   = ges_q(i,j,k-1,jj) * prs_a / pw   ! vapor pressure (Pa)
                   rl_hm   = prs_v / prs_sv    ! relative humidity
                   x_v     = rl_hm * ehn_fct * prs_sv / prs_a     ! molar fraction of water vapor (eq. A1.3)
                   cmpr    = one - (prs_a/tmp_K) * ( cpf_a0 + cpf_a1*tmp_C + cpf_a2*tmp_C**2 &
                            + (cpf_b0 + cpf_b1*tmp_C)*x_v + (cpf_c0 + cpf_c1*tmp_C)*x_v**2 ) &
                            + (prs_a**2/tmp_K**2) * (cpf_d + cpf_e*x_v**2)
                   h       = rdog * ges_tv(i,j,k-1,jj)
                   dz      = h * cmpr * log(ges_prsi(i,j,k-1,jj)/ges_prsi(i,j,k,jj))
                   height(k) = height(k-1) + dz
                enddo

                k=nsig+1
                fact    = one + fv* ges_q(i,j,k-1,jj)
                pw      = eps + ges_q(i,j,k-1,jj)*(one - eps)
                tmp_K   = ges_tv(i,j,k-1,jj) / fact
                tmp_C   = tmp_K - t0c
                prs_sv  = exp(psv_a*tmp_K**2 + psv_b*tmp_K + psv_c + psv_d/tmp_K)  ! eq A1.1 (Pa)
                prs_a   = thousand * exp(half*(log(ges_prsi(i,j,k-1,jj))+log(ges_prsl(i,j,k-1,jj))))     ! (Pa)
                ehn_fct = ef_alpha + ef_beta*prs_a + ef_gamma*tmp_C**2 ! enhancement factor (eq. A1.2)
                prs_v   = ges_q(i,j,k-1,jj) * prs_a / pw  
                rl_hm   = prs_v / prs_sv    ! relative humidity
                x_v     = rl_hm * ehn_fct * prs_sv / prs_a     ! molar fraction of water vapor (eq. A1.3)
                cmpr    = one - (prs_a/tmp_K) * ( cpf_a0 + cpf_a1*tmp_C + cpf_a2*tmp_C**2 &
                          + (cpf_b0 + cpf_b1*tmp_C)*x_v + (cpf_c0 + cpf_c1*tmp_C)*x_v**2 ) &
                          + (prs_a**2/tmp_K**2) * (cpf_d + cpf_e*x_v**2)
                h       = rdog * ges_tv(i,j,k-1,jj)
                dz      = h * cmpr * log(ges_prsi(i,j,k-1,jj)/ges_prsl(i,j,k-1,jj))
                height(k) = height(k-1) + dz
 
                do k=1,nsig+1
                   geop_hgti(i,j,k,jj)=height(k) - ges_z(i,j,jj)
                end do
             enddo
          enddo
       enddo

    else

!      Compute geopotential height at midpoint of each layer
       do jj=1,nfldsig
          do j=1,lon2
             do i=1,lat2
                k  = 1
                h  = rdog * ges_tv(i,j,k,jj)
                dz = h * log(ges_prsi(i,j,k,jj)/ges_prsl(i,j,k,jj))
                height(k) = ges_z(i,j,jj) + dz
 
                do k=2,nsig
                   h  = rdog * half * (ges_tv(i,j,k-1,jj)+ges_tv(i,j,k,jj))
                   dz = h * log(ges_prsl(i,j,k-1,jj)/ges_prsl(i,j,k,jj))
                   height(k) = height(k-1) + dz
                end do

                do k=1,nsig
                   geop_hgtl(i,j,k,jj)=height(k) - ges_z(i,j,jj)
                end do
             end do
          end do
       end do

!      Compute geopotential height at interface between layers
       do jj=1,nfldsig
          do j=1,lon2
             do i=1,lat2
                k=1
                height(k) = ges_z(i,j,jj)

                do k=2,nsig
                   h  = rdog * ges_tv(i,j,k-1,jj)
                   dz = h * log(ges_prsi(i,j,k-1,jj)/ges_prsi(i,j,k,jj))
                   height(k) = height(k-1) + dz
                end do

                k=nsig+1
                h = rdog * ges_tv(i,j,k-1,jj)
                dz = h * log(ges_prsi(i,j,k-1,jj)/ges_prsl(i,j,k-1,jj))
                height(k) = height(k-1) + dz

                do k=1,nsig+1
                   geop_hgti(i,j,k,jj)=height(k) - ges_z(i,j,jj)
                end do
             end do
          end do
       end do

    endif

    return
  end subroutine load_geop_hgt

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: load_gsdpbl_hgt --- Populate PBL height
!
! !INTERFACE:
!
  subroutine load_gsdpbl_hgt(mype)

! !USES:

    use constants, only: one,rd_over_cp_mass,r1000,ten,zero,two
    use gridmod, only: lat2, lon2, nsig,wrf_mass_regional, &
         twodvar_regional,nems_nmmb_regional

    implicit none

! !INPUT PARAMETERS:


! !DESCRIPTION: populate guess geopotential height in millibars
!
!
! !REVISION HISTORY:
!   2011-06-06  Ming Hu
!   2013-02-22  Jacob Carley - Added NMMB
!
! !REMARKS:
!   language: f90
!   machine:  JET
!
! !AUTHOR:
!
!EOP
!-------------------------------------------------------------------------

    integer(i_kind)              , intent(in   ) :: mype
    integer(i_kind) i,j,k,jj
    real(r_kind),dimension(nsig):: thetav
    real(r_kind),dimension(nsig):: pbk
    real(r_kind) :: thsfc, d

    if (twodvar_regional) return

!   Compute geopotential height at midpoint of each layer
    do jj=1,nfldsig
       do j=1,lon2
          do i=1,lat2

             do k=1,nsig

                if (wrf_mass_regional)  pbk(k) = aeta1_ll(k)*(ges_ps(i,j,1)*ten-pt_ll)+pt_ll
		if (nems_nmmb_regional) then
		   pbk(k) = aeta1_ll(k)*pdtop_ll + aeta2_ll(k)*(ten*ges_ps(i,j,jj) & 
		            -pdtop_ll-pt_ll) + pt_ll   			    			    
		end if
				
		thetav(k)  = ges_tv(i,j,k,jj)*(r1000/pbk(k))**rd_over_cp_mass
             end do

             pbl_height(i,j,jj) = zero
             thsfc = thetav(1)
             k=1
             DO while (abs(pbl_height(i,j,jj)) < 0.0001_r_kind)
               if( thetav(k) > thsfc + 1.0_r_kind ) then
                 pbl_height(i,j,jj) = float(k) - (thetav(k) - (thsfc + 1.0_r_kind))/   &
                             max((thetav(k)-thetav(k-1)),0.01_r_kind)
               endif
               k=k+1
             ENDDO
             if(abs(pbl_height(i,j,jj)) < 0.0001_r_kind) pbl_height(i,j,jj)=two
             k=int(pbl_height(i,j,jj))
             if( k < 1 .or. k > nsig-1) then
                write(6,*) ' Error in PBL height calculation ',mype,i,j,pbl_height(i,j,jj)
             endif
             d=pbl_height(i,j,jj) - k
             pbl_height(i,j,jj) = pbk(k) * (one-d) + pbk(k+1) * d

          end do
       end do
    end do

    return
  end subroutine load_gsdpbl_hgt

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: add_rtm_layers --- Add pressure layers for RTM use
!
! !INTERFACE:
!
  subroutine add_rtm_layers(prsitmp,prsltmp,prsitmp_ext,prsltmp_ext,klevel)

! !USES:

    use constants, only: half,one_tenth
    use gridmod, only: nsig,msig,nlayers
    use crtm_module, only: toa_pressure

    implicit none

! !INPUT PARAMETERS:
    integer(i_kind),dimension(msig)  ,intent(  out) :: klevel

    real(r_kind)   ,dimension(nsig+1),intent(in   ) :: prsitmp
    real(r_kind)   ,dimension(nsig)  ,intent(in   ) :: prsltmp

    real(r_kind)   ,dimension(msig+1),intent(  out) :: prsitmp_ext
    real(r_kind)   ,dimension(msig)  ,intent(  out) :: prsltmp_ext


! !DESCRIPTION:  Add pressure layers for use in RTM
!
! !REVISION HISTORY:
!   2005-06-01  treadon
!   2006-05-10  derber modify how levels are added above model top
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
!   treadon          org: w/nmc20      date: 2005-06-01
!
!EOP
!-------------------------------------------------------------------------

!   Declare local variables
    integer(i_kind) k,kk,l
    real(r_kind) dprs,toa_pressure01

    toa_pressure01=toa_pressure*one_tenth

!   Check if model top pressure above rtm top pressure, where prsitmp
!   is in kPa and toa_pressure is in hPa.
    if (prsitmp(nsig) < toa_pressure01)then
       write(6,*)'ADD_RTM_LAYERS:  model top pressure(hPa)=', &
            prsitmp(nsig),&
            ' above rtm top pressure(hPa)=',toa_pressure01
       call stop2(35)
    end if

!   Linear in pressure sub-divsions
    kk=0
    do k = 1,nsig
       if (nlayers(k)<=1) then
          kk = kk + 1
          prsltmp_ext(kk) = prsltmp(k)
          prsitmp_ext(kk) = prsitmp(k)
          klevel(kk) = k
       else
          if (k/=nsig) then
             dprs = (prsitmp(k+1)-prsitmp(k))/nlayers(k)
          else
             dprs = (toa_pressure01-prsitmp(k))/nlayers(k)
          end if
          prsitmp_ext(kk+1) = prsitmp(k)
          do l=1,nlayers(k)
             kk=kk + 1
             prsitmp_ext(kk+1) = prsitmp(k) + dprs*l
             prsltmp_ext(kk) = half*(prsitmp_ext(kk+1)+prsitmp_ext(kk))
             klevel(kk) = k
          end do
       endif
    end do

!   Set top of atmosphere pressure
    prsitmp_ext(msig+1) = toa_pressure01

  end subroutine add_rtm_layers

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: load_fact10 --- Compute 10m wind factor
!
! !INTERFACE:
!
  subroutine load_fact10

! !USES:

    use gridmod, only: lat2,lon2
    implicit none

! !INPUT PARAMETERS:

! !DESCRIPTION: compute 10m wind factor
!
! !REVISION HISTORY:
!   2006-09-26  treadon
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
!   treadon          org: w/nmc20      date: 2006-09-26
!
!EOP
!-------------------------------------------------------------------------

!   Declare local variables
    logical iqtflg
    integer(i_kind):: i,j,it,itt,nt,regime
    integer(i_kind),dimension(nfldsfc):: indx
    real(r_kind):: u10ges,v10ges,t2ges,q2ges

    nt=0
    indx=1
    do i=1,nfldsfc
       if(abs(hrdifsfc(i)-hrdifsig(i))<0.001_r_kind) then
          nt=nt+1
          indx(nt) = i
       endif
    end do

    if (sfcmod_gfs) then
       do it=1,nt
          itt=indx(it)
          do j=1,lon2
             do i=1,lat2
                call compute_fact10(ges_u(i,j,1,itt),ges_v(i,j,1,itt),&
                     ges_tsen(i,j,1,itt),ges_q(i,j,1,itt),&
                     ges_ps(i,j,itt),ges_prsi(i,j,1,itt), &
                     ges_prsi(i,j,2,itt),sfct(i,j,itt), &
                     sfc_rough(i,j,itt),isli(i,j,itt),fact10(i,j,itt))
             end do
          end do
       end do
    endif

    if (sfcmod_mm5) then
       iqtflg=.true.
       do it=1,nt
          itt=indx(it)
          do j=1,lon2
             do i=1,lat2
                call SFC_WTQ_FWD (&
                     ges_ps(i,j,itt),&
                     sfct(i,j,itt),&
                     ges_lnprsl(i,j,1,itt),&
                     ges_tv(i,j,1,itt),&
                     ges_q(i,j,1,itt),&
                     ges_u(i,j,1,itt),&
                     ges_v(i,j,1,itt),&
                     ges_lnprsl(i,j,2,itt),&
                     ges_tv(i,j,2,itt),&
                     ges_q(i,j,2,itt),&
                     geop_hgtl(i,j,1,itt),&
                     sfc_rough(i,j,itt),&
                     isli(i,j,itt),&
                     fact10(i,j,itt),&
                     u10ges,v10ges,t2ges,q2ges,regime,iqtflg)
             end do
          end do
       end do
    endif

    return
  end subroutine load_fact10

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE: comp_fact10  ---  compute 10m wind factor
!
! !INTERFACE:
!
  subroutine comp_fact10(dlat,dlon,dtime,skint,sfcrough,islimsk,mype,factw)

! !USES:

    use gridmod, only: nlat,nlon,&
         lon1,istart,jstart
    use constants, only: zero,one
    implicit none

! !INPUT PARAMETERS:

    real(r_kind)   ,intent(in   ) :: dlat,dlon,dtime,skint,sfcrough
    real(r_kind)   ,intent(inout) :: factw
    integer(i_kind),intent(in   ) :: mype,islimsk

! !DESCRIPTION: compute 10m wind factor
!
! !REVISION HISTORY:
!   2006-09-26  treadon
!   2008-12-05  todling - use dsfct(:,:,ntguessfc) for calculation
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
!   treadon          org: w/nmc20      date: 2006-09-26
!
!EOP
!-------------------------------------------------------------------------

!   Declare local parameters

!   Declare local variables
    logical iqtflg
    integer(i_kind) ix,ix1,ixp,iy,iy1,iyp,regime
    integer(i_kind) itsig,itsigp,j,m1,islimsk2
    real(r_kind) w00,w01,w10,w11
    real(r_kind) delx,dely,delx1,dely1,dtsig,dtsigp
    real(r_kind):: u10ges,v10ges,t2ges,q2ges
    real(r_kind):: pgesin,ugesin,vgesin,qgesin,tgesin,prsigesin1
    real(r_kind):: prsigesin2,lnpgesin1,lnpgesin2,tgesin2,qgesin2,geopgesin,ts

    islimsk2=islimsk
    if(islimsk2 > 2)islimsk2=islimsk2-3
    m1=mype+1
!   Set spatial interpolation indices and weights
    ix1=dlat
    ix1=max(1,min(ix1,nlat))
    delx=dlat-ix1
    delx=max(zero,min(delx,one))
    ix=ix1-istart(m1)+2
    ixp=ix+1
    if(ix1==nlat) then
       ixp=ix
    end if
    delx1=one-delx

    iy1=dlon
    dely=dlon-iy1
    iy=iy1-jstart(m1)+2
    if(iy<1) then
       iy1=iy1+nlon
       iy=iy1-jstart(m1)+2
    end if
    if(iy>lon1+1) then
       iy1=iy1-nlon
       iy=iy1-jstart(m1)+2
    end if
    iyp=iy+1
    dely1=one-dely


    w00=delx1*dely1; w10=delx*dely1; w01=delx1*dely; w11=delx*dely
!   Get time interpolation factors for sigma files
    if(dtime > hrdifsig(1) .and. dtime < hrdifsig(nfldsig))then
       do j=1,nfldsig-1
          if(dtime > hrdifsig(j) .and. dtime <= hrdifsig(j+1))then
             itsig=j
             itsigp=j+1
             dtsig=((hrdifsig(j+1)-dtime)/(hrdifsig(j+1)-hrdifsig(j)))
          end if
       end do
    else if(dtime <=hrdifsig(1))then
       itsig=1
       itsigp=1
       dtsig=one
    else
       itsig=nfldsig
       itsigp=nfldsig
       dtsig=one
    end if
    dtsigp=one-dtsig

    ts =(dsfct(ix,iy ,ntguessfc)*w00 + dsfct(ixp,iy ,ntguessfc)*w10 +          &
         dsfct(ix,iyp,ntguessfc)*w01 + dsfct(ixp,iyp,ntguessfc)*w11) + skint

    pgesin=(ges_ps(ix,iy ,itsig )*w00+ges_ps(ixp,iy ,itsig )*w10+ &
            ges_ps(ix,iyp,itsig )*w01+ges_ps(ixp,iyp,itsig )*w11)*dtsig + &
           (ges_ps(ix,iy ,itsigp)*w00+ges_ps(ixp,iy ,itsigp)*w10+ &
            ges_ps(ix,iyp,itsigp)*w01+ges_ps(ixp,iyp,itsigp)*w11)*dtsigp
    ugesin=(ges_u(ix,iy ,1,itsig )*w00+ges_u(ixp,iy ,1,itsig )*w10+ &
            ges_u(ix,iyp,1,itsig )*w01+ges_u(ixp,iyp,1,itsig )*w11)*dtsig + &
           (ges_u(ix,iy ,1,itsigp)*w00+ges_u(ixp,iy ,1,itsigp)*w10+ &
            ges_u(ix,iyp,1,itsigp)*w01+ges_u(ixp,iyp,1,itsigp)*w11)*dtsigp
    vgesin=(ges_v(ix,iy ,1,itsig )*w00+ges_v(ixp,iy ,1,itsig )*w10+ &
            ges_v(ix,iyp,1,itsig )*w01+ges_v(ixp,iyp,1,itsig )*w11)*dtsig + &
           (ges_v(ix,iy ,1,itsigp)*w00+ges_v(ixp,iy ,1,itsigp)*w10+ &
            ges_v(ix,iyp,1,itsigp)*w01+ges_v(ixp,iyp,1,itsigp)*w11)*dtsigp
    qgesin=(ges_q(ix,iy ,1,itsig )*w00+ges_q(ixp,iy ,1,itsig )*w10+ &
            ges_q(ix,iyp,1,itsig )*w01+ges_q(ixp,iyp,1,itsig )*w11)*dtsig + &
           (ges_q(ix,iy ,1,itsigp)*w00+ges_q(ixp,iy ,1,itsigp)*w10+ &
            ges_q(ix,iyp,1,itsigp)*w01+ges_q(ixp,iyp,1,itsigp)*w11)*dtsigp


    if (sfcmod_gfs) then
       tgesin    =(ges_tsen(ix ,iy ,1,itsig )*w00+ &
                   ges_tsen(ixp,iy ,1,itsig )*w10+ &
                   ges_tsen(ix ,iyp,1,itsig )*w01+ &
                   ges_tsen(ixp,iyp,1,itsig )*w11)*dtsig + &
                  (ges_tsen(ix ,iy ,1,itsigp)*w00+ &
                   ges_tsen(ixp,iy ,1,itsigp)*w10+ &
                   ges_tsen(ix ,iyp,1,itsigp)*w01+ &
                   ges_tsen(ixp,iyp,1,itsigp)*w11)*dtsigp
       prsigesin1=(ges_prsi(ix ,iy ,1,itsig )*w00+ &
                   ges_prsi(ixp,iy ,1,itsig )*w10+ &
                   ges_prsi(ix ,iyp,1,itsig )*w01+ &
                   ges_prsi(ixp,iyp,1,itsig )*w11)*dtsig + &
                  (ges_prsi(ix ,iy ,1,itsigp)*w00+ &
                   ges_prsi(ixp,iy ,1,itsigp)*w10+ &
                   ges_prsi(ix ,iyp,1,itsigp)*w01+ &
                   ges_prsi(ixp,iyp,1,itsigp)*w11)*dtsigp
       prsigesin2=(ges_prsi(ix ,iy ,2,itsig )*w00+ &
                   ges_prsi(ixp,iy ,2,itsig )*w10+ &
                   ges_prsi(ix ,iyp,2,itsig )*w01+ &
                   ges_prsi(ixp,iyp,2,itsig )*w11)*dtsig + &
                  (ges_prsi(ix ,iy ,2,itsigp)*w00+ &
                   ges_prsi(ixp,iy ,2,itsigp)*w10+ &
                   ges_prsi(ix ,iyp,2,itsigp)*w01+ &
                   ges_prsi(ixp,iyp,2,itsigp)*w11)*dtsigp
       call compute_fact10(ugesin,vgesin,tgesin,qgesin,pgesin, &
            prsigesin1,prsigesin2,ts,sfcrough,islimsk,factw)
    else if (sfcmod_mm5)then
       iqtflg=.true.
       lnpgesin1 =(ges_lnprsl(ix ,iy ,1,itsig )*w00+ &
                   ges_lnprsl(ixp,iy ,1,itsig )*w10+ &
                   ges_lnprsl(ix ,iyp,1,itsig )*w01+ &
                   ges_lnprsl(ixp,iyp,1,itsig )*w11)*dtsig + &
                  (ges_lnprsl(ix ,iy ,1,itsigp)*w00+ &
                   ges_lnprsl(ixp,iy ,1,itsigp)*w10+ &
                   ges_lnprsl(ix ,iyp,1,itsigp)*w01+ &
                   ges_lnprsl(ixp,iyp,1,itsigp)*w11)*dtsigp
       lnpgesin2 =(ges_lnprsl(ix ,iy ,2,itsig )*w00+ &
                   ges_lnprsl(ixp,iy ,2,itsig )*w10+ &
                   ges_lnprsl(ix ,iyp,2,itsig )*w01+ &
                   ges_lnprsl(ixp,iyp,2,itsig )*w11)*dtsig + &
                  (ges_lnprsl(ix ,iy ,2,itsigp)*w00+ &
                   ges_lnprsl(ixp,iy ,2,itsigp)*w10+ &
                   ges_lnprsl(ix ,iyp,2,itsigp)*w01+ &
                   ges_lnprsl(ixp,iyp,2,itsigp)*w11)*dtsigp
       tgesin    =(ges_tv(ix ,iy ,1,itsig )*w00+ &
                   ges_tv(ixp,iy ,1,itsig )*w10+ &
                   ges_tv(ix ,iyp,1,itsig )*w01+ &
                   ges_tv(ixp,iyp,1,itsig )*w11)*dtsig + &
                  (ges_tv(ix ,iy ,1,itsigp)*w00+ &
                   ges_tv(ixp,iy ,1,itsigp)*w10+ &
                   ges_tv(ix ,iyp,1,itsigp)*w01+ &
                   ges_tv(ixp,iyp,1,itsigp)*w11)*dtsigp
       tgesin2   =(ges_tv(ix ,iy ,2,itsig )*w00+ &
                   ges_tv(ixp,iy ,2,itsig )*w10+ &
                   ges_tv(ix ,iyp,2,itsig )*w01+ &
                   ges_tv(ixp,iyp,2,itsig )*w11)*dtsig + &
                  (ges_tv(ix ,iy ,2,itsigp)*w00+ &
                   ges_tv(ixp,iy ,2,itsigp)*w10+ &
                   ges_tv(ix ,iyp,2,itsigp)*w01+ &
                   ges_tv(ixp,iyp,2,itsigp)*w11)*dtsigp
       qgesin2   =(ges_q(ix ,iy ,2,itsig )*w00+ &
                   ges_q(ixp,iy ,2,itsig )*w10+ &
                   ges_q(ix ,iyp,2,itsig )*w01+ &
                   ges_q(ixp,iyp,2,itsig )*w11)*dtsig + &
                  (ges_q(ix ,iy ,2,itsigp)*w00+ &
                   ges_q(ixp,iy ,2,itsigp)*w10+ &
                   ges_q(ix ,iyp,2,itsigp)*w01+ &
                   ges_q(ixp,iyp,2,itsigp)*w11)*dtsigp
       geopgesin =(geop_hgtl(ix ,iy ,1,itsig )*w00+ &
                   geop_hgtl(ixp,iy ,1,itsig )*w10+ &
                   geop_hgtl(ix ,iyp,1,itsig )*w01+ &
                   geop_hgtl(ixp,iyp,1,itsig )*w11)*dtsig + &
                  (geop_hgtl(ix ,iy ,1,itsigp)*w00+ &
                   geop_hgtl(ixp,iy ,1,itsigp)*w10+ &
                   geop_hgtl(ix ,iyp,1,itsigp)*w01+ &
                   geop_hgtl(ixp,iyp,1,itsigp)*w11)*dtsigp
       call SFC_WTQ_FWD (pgesin,ts,lnpgesin1,tgesin,qgesin,ugesin,vgesin, &
                lnpgesin2,tgesin2,qgesin2,geopgesin,sfcrough,islimsk, &
                factw,u10ges,v10ges,t2ges,q2ges,regime,iqtflg)
    endif

    return
  end subroutine comp_fact10


!-------------------------------------------------------------------------
   subroutine guess_grids_stats3d_(name,a,mype)
!-------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    guess_grids_stats3d_
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-04  lueken - added subprogram doc block
!
!   input argument list:
!    name
!    a
!    mype     - mpi task id
!
!   output argument log:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

   use constants, only: zero
   use mpimod, only: ierror,mpi_rtype,mpi_sum,mpi_comm_world
   use gridmod, only: lon1,lat1,nsig

   implicit none

   character(len=*)             , intent(in   ) :: name
   real(r_kind),dimension(:,:,:), intent(in   ) :: a
   integer(i_kind)              , intent(in   ) :: mype


! local variables
   integer(i_kind) :: i,j,k
   real(r_kind),dimension(nsig+1):: work_a,work_a1
   real(r_kind),dimension(nsig):: amz ! global mean profile of a
   real(r_kind) :: rms

! start

!  Calculate global means for a

!  Calculate sums for a to estimate variance.
   work_a = zero
   do k = 1,nsig
      do j = 2,lon1+1
         do i = 2,lat1+1
            work_a(k) = work_a(k) + a(i,j,k)
         end do
      end do
   end do
   work_a(nsig+1)=float(lon1*lat1)

   call mpi_allreduce(work_a,work_a1,nsig+1,mpi_rtype,mpi_sum,&
       mpi_comm_world,ierror)

   amz=zero
   do k=1,nsig
      if (work_a1(nsig+1)>zero) amz(k)=work_a1(k)/work_a1(nsig+1)
      rms=sqrt(amz(k)**2/work_a1(nsig+1))
      if (mype==0) write(*,100) trim(name),k,amz(k),rms
   enddo
100 format(a,': Level, Global mean, RMS = ',i3,1P2E16.8)

   end subroutine guess_grids_stats3d_

!-------------------------------------------------------------------------
   subroutine guess_grids_stats2d_(name,a,mype)
!-------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    guess_grids_stats2d_
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-04  lueken - added subprogram doc block
!
!   input argument list:
!    name
!    a
!    mype     - mpi task id
!
!   output argument log:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

   use constants, only: zero
   use mpimod, only: ierror,mpi_rtype,mpi_sum,mpi_comm_world
   use gridmod, only: lon1,lat1

   implicit none

   character(len=*)           , intent(in   ) :: name
   real(r_kind),dimension(:,:), intent(in   ) :: a
   integer(i_kind)            , intent(in   ) :: mype


! local variables
   integer(i_kind) :: i,j
   real(r_kind),dimension(2):: work_a,work_a1
   real(r_kind) :: amz, rms

! start

!  Calculate global means for a

!  Calculate sums for a to estimate variance.
   work_a = zero
   do j = 2,lon1+1
      do i = 2,lat1+1
         work_a(1) = work_a(1) + a(i,j)
      end do
   end do
   work_a(2)=float(lon1*lat1)

   call mpi_allreduce(work_a,work_a1,2,mpi_rtype,mpi_sum,&
       mpi_comm_world,ierror)

   amz=zero
   if (work_a1(2)>zero) amz=work_a1(1)/work_a1(2)
   rms=sqrt(amz**2/work_a1(2))      
   if (mype==0) write(*,100) trim(name),amz,rms
100 format(a,': Global mean, RMS = ',1P2E16.8)

   end subroutine guess_grids_stats2d_

!-------------------------------------------------------------------------
   subroutine pstats_(a,amiss,avg,rms)
!-------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    pstats_
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-04  lueken - added subprogram doc block
!
!   input argument list:
!    amiss
!    a
!
!   output argument log:
!    avg,rms
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
   use constants, only: zero
   implicit none

   real(r_kind),dimension(:,:), intent(in   ) :: a      ! array var
   real(r_kind)               , intent(in   ) :: amiss  ! undef
   real(r_kind)               , intent(  out) :: avg,rms


! local variables
   integer(i_kind) :: i,j
   integer(i_kind) :: allcnt,cnt

! start

   allcnt=0
   cnt=0
   avg=zero
   rms=zero
   do i=1,size(a,1)
      do j=1,size(a,2)
         if(a(i,j)/=amiss) then
            cnt=cnt+1
            avg=avg+a(i,j)
         endif
         allcnt = allcnt+1
      end do
   end do
   avg=avg/max(1,cnt)
   rms=sqrt(avg*avg/max(1,cnt))      

   end subroutine pstats_

!-------------------------------------------------------------------------
   subroutine print1r8_(name,fld,undef)
!-------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    print1r8_
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-04  lueken - added subprogram doc block
!
!   input argument list:
!    name
!    fld 
!    undef
!
!   output argument log:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
   implicit none
   character(len=*)         , intent(in   ) :: name
   real(r_kind),dimension(:), intent(in   ) :: fld
   real(r_kind)             , intent(in   ) :: undef
! 
   write(6,100) trim(name),minval(fld),maxval(fld),sum(fld),undef
100 format(a,': range,sum = ',1P3E16.4)
   end subroutine print1r8_

!-------------------------------------------------------------------------
   subroutine print2r8_(name,fld,undef)
!-------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    print2r8_
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-04  lueken - added subprogram doc block
!
!   input argument list:
!    name
!    fld
!    undef
!
!   output argument log:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

   implicit none
   character(len=*)           , intent(in   ) :: name
   real(r_kind),dimension(:,:), intent(in   ) :: fld
   real(r_kind)               , intent(in   ) :: undef
! 
   real(r_kind) avg,rms
   write(6,100) trim(name),minval(fld),maxval(fld),sum(fld)
   call pstats_(fld,UNDEF,avg,rms)
   write(6,99) trim(name),avg,rms
100 format(a,': range,sum = ',1P3E16.4)
99  format(a,': avg, rms = ',1P2E16.4)
   end subroutine print2r8_

!-------------------------------------------------------------------------
   subroutine print3r8_(name,fld,undef,allk)
!-------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    print3r8_
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-04  lueken - added subprogram doc block
!
!   input argument list:
!    name
!    fld
!    undef
!    allk
!
!   output argument log:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
   implicit none
   character(len=*)             , intent(in   ) :: name
   real(r_kind),dimension(:,:,:), intent(in   ) :: fld
   real(r_kind)                 , intent(in   ) :: undef
   logical,optional             , intent(in   ) :: allk
! 
   logical prntlevs
   integer(i_kind) k
   real(r_kind) avg,rms

   if(present(allk)) prntlevs=allk
   if(prntlevs) then
      do k=1,size(fld,3)
         write(6,101) trim(name),k,minval(fld(:,:,k)),maxval(fld(:,:,k)),sum(fld(:,:,k))
         call pstats_(fld(:,:,k),UNDEF,avg,rms)
         write(6,99) trim(name),avg,rms
      end do
   else
      write(6,100) trim(name),minval(fld),maxval(fld),sum(fld)
   end if
101 format(a,': time or lev,range,sum = ',i3,1P3E16.4)
100 format(a,': range,sum = ',1P3E16.4)
99  format(a,': avg, rms = ',1P2E16.4)
   end subroutine print3r8_

!-------------------------------------------------------------------------
   subroutine print4r8_(name,fld,undef,allk)
!-------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    print4r8_
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-04  lueken - added subprogram doc block
!
!   input argument list:
!    name
!    fld
!    undef
!    allk
!
!   output argument log:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

   implicit none
   character(len=*)               , intent(in   ) :: name
   real(r_kind),dimension(:,:,:,:), intent(in   ) :: fld
   real(r_kind)                   , intent(in   ) :: undef
   logical,optional               , intent(in   ) :: allk
! 
   logical prntlevs
   integer(i_kind) k,it
   real(r_kind) avg,rms

   if(present(allk)) prntlevs=allk
   if(prntlevs) then
      do it=1,size(fld,4)
         do k=1,size(fld,3)
            write(6,101) trim(name),it,k,minval(fld(:,:,k,it)),maxval(fld(:,:,k,it)),sum(fld(:,:,k,it))
            call pstats_(fld(:,:,k,it),UNDEF,avg,rms)
            write(6,99) trim(name),avg,rms
         end do
      end do
   else
      write(6,100) trim(name),minval(fld),maxval(fld),sum(fld)
   end if
101 format(a,': time,lev,range,sum = ',i3,i3,1P3E16.4)
100 format(a,': range,sum = ',1P3E16.4)
99  format(a,': avg, rms = ',1P2E16.4)
   end subroutine print4r8_
    
end module guess_grids
