module ncepnems_io
!$$$ module documentation block
!           .      .    .                                       .
! module:   ncepnems_io
!   prgmmr: Huang       org: np23                date: 2010-02-22
!
! abstract: This module contains routines which handle input/output
!           operations for NCEP NEMS global atmospheric and surface files.
!
! program history log:
!   2010-02-22 Huang    Initial version.  Based on ncepgfs_io
!   2010-10-18 Huang    Remove subroutine reorder_gfsgrib for no longer been called in GSI
!                       For Now, subroutine sfc_interpolate is kept in ncepgfs_io.f90.
!                       When sigio and gfsio are both retired, i.e., remove ncepgfs_io.f90.
!                       move this routines back to this module
!   2011-03-03 Huang    Changes has been made to adopt to high resolution GSI run (T382 & T574)
!                       both for CPU and memory issues.
!                       Future development of nemsio need to consider a mapping routine be
!                       inserted between read-in forecast field and GSI first guess field,
!                       as well as GSI analysis field and write-out data field for forecast
!                       model. Due to computation resource, GSI may not be able to run at
!                       the same resolution as that of forecast model, e.g., run GSI at T382
!                       w/ T574 forecast model output.
!   2011-10-25 Huang    (1) Add unified error message routine to make the code cleaner
!                       (2) To reduce the memory allocation as low as possible, remove all
!                           reference to sfc_head and re-used the same local arrays.
!                           Remove unneeded nemsio_data & gfsdata.
!                       (3) Add parallel IO code to read_atm_
!   2011-11-01 Huang    (1) add integer nst_gsi to control the mode of NSST
!                       (2) add read_nemsnst to read ncep nst file
!                       (3) add subroutine write_nemssfc_nst to save sfc and nst files
!   2016-01-01 Li       (1) Move tran_gfssfc from ncepgfs_io.f90 to here
!                       (2) Modify write_sfc_nst_ to follows the update done in sfcio
!                       (3) Modify read_sfc_ to follows the update done in sfcio for more effective I/O
!   2016-04-20 Li       Modify to handle the updated nemsio sig file (P, DP & DPDT removed)
!   2016-08-18 li     - tic591: add read_sfc_anl & read_nemssfc_anl to read nemsio sfc file (isli only) with analysis resolution
!                               change/modify sfc_interpolate to be intrp22 to handle more general interpolation (2d to 2d)
!   2016-11-18 li     - tic615: change nst mask name from slmsk to land 
!   2017-08-30 li     - tic659: modify read_nems_sfc_ and read_sfc_  to read sfc file in 
!                               nemsio Gaussin grids generated by FV3 WriteComponent
!
! Subroutines Included:
!   sub read_nems       - driver to read ncep nems atmospheric and surface
!   sub read_nems_chem
!   sub read_nemsatm    - read ncep nems atmospheric file, scatter
!                         on grid to analysis subdomains
!   sub read_nemssfc    - read ncep nems surface file, scatter on grid to
!                         analysis subdomains
!   sub read_nemssfc_anl- read ncep EnKF nems surface file, scatter on grid to
!                         analysis subdomains
!   sub write_nems      - driver to write ncep nems atmospheric and surface
!                         analysis files
!   sub write_nemsatm   - gather on grid, write ncep nems atmospheric analysis file
!   sub write_nemssfc   - gather/write on grid ncep surface analysis file
!   sub read_nemsnst    - read ncep nst file, scatter on grid to analysis subdomains
!   sub write_nems_sfc_nst - gather/write on grid ncep surface & nst analysis file
!   sub intrp22         - interpolate from one grid to another grid (2D)
!   sub read_nems_sfcnst  - read sfc hist file, including sfc and nst vars, scatter on grid to analysis subdomains
!
! Variable Definitions:
!   The difference of time Info between operational GFS IO (gfshead%, sfc_head%),
!      analysis time (iadate), and NEMSIO (idate=)
!
!       gfshead & sfc_head            NEMSIO Header           Analysis time (obsmod)
!       ===================   ============================  ==========================
!         %idate(1)  Hour     idate(1)  Year                iadate(1)  Year
!         %idate(2)  Month    idate(2)  Month               iadate(2)  Month
!         %idate(3)  Day      idate(3)  Day                 iadate(3)  Day
!         %idate(4)  Year     idate(4)  Hour                iadate(4)  Hour
!                             idate(5)  Minute              iadate(5)  Minute
!                             idate(6)  Scaled seconds
!                             idate(7)  Seconds multiplier
!
!   The difference of header forecasting hour Info bewteen operational GFS IO
!      (gfshead%, sfc_head%) and NEMSIO
!
!           gfshead & sfc_head                NEMSIO Header
!       ==========================     ============================
!       %fhour  FCST Hour (r_kind)     nfhour     FCST Hour (i_kind)
!                                      nfminute   FCST Mins (i_kind)
!                                      nfsecondn  FCST Secs (i_kind) numerator
!                                      nfsecondd  FCST Secs (i_kind) denominator
!
!       %fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
!
!   nframe     - nframe is the number of grids extend outward from the
!                edge of modeling domain.
!
!                NEMSIO provides a more flexible read.  User can get the
!                size of record (1D) to be read from file header. The
!                normal record size should be delx*dely, i.e., total model
!                grid points.  However, some regional models also ouput
!                additional data of grids around the modeling domain
!                (buffer zone). For this type of output, nframe needs to
!                be know to calculate the size of record, i.e.,
!                   array size = (delx+2*nframe) * (dely+2*nframe)
!
!                However, nframe should always be zero for global model.
!                To simplify the code for reading and writing global model
!                files, we will not factor in the nframe for computing
!                array size or array index shift (by nframe) between
!                input/output array and internal GSI array.  The normal
!                size of I/O record remains as delx*dely.  Add a checking
!                routine to assure nframe=zero.
!
! attributes:
!   language: f90
!   machine:
!
! NOTE: When global meteorology switched to NEMS/GFS, all routines and
!       modules of old GFS (sigio) can be deactivated.  To keep the code
!       clean, all "nems" can be replaced by "gfs" for minimal changes
!       of GSI code structure.  For dual purpose, two distincit routine
!       names are used to accomodiate old and new systems.  It is now
!       controled by a namelist argument "use_gfs_nemsio"
!
!
!$$$ end documentation block

  use constants, only: zero,one,fv,r60,r3600
  implicit none

  private
  public read_nems
  public read_nems_chem
  public read_nemsatm
  public read_nemssfc
  public read_nemssfc_anl
  public write_nemsatm
  public write_nemssfc
  public read_nemsnst
  public write_nems_sfc_nst
  public intrp22
  public tran_gfssfc
  public error_msg

  interface read_nems
     module procedure read_
  end interface

  interface read_nems_chem
     module procedure read_chem_
  end interface

  interface read_nemsatm
     module procedure read_atm_
  end interface

  interface read_nemssfc
     module procedure read_nemssfc_
  end interface

  interface read_nemssfc_anl
     module procedure read_nemssfc_anl_
  end interface

  interface read_nemsnst
     module procedure read_nemsnst_
  end interface

  interface write_nemsatm
     module procedure write_atm_
  end interface

  interface write_nemssfc
     module procedure write_sfc_
  end interface

  interface write_nems_sfc_nst
     module procedure write_sfc_nst_
  end interface

  interface error_msg
     module procedure error_msg_
  end interface

  character(len=*),parameter::myname='ncepnems_io'

contains

  subroutine read_
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    read_nems
!
!   prgrmmr: Ho-Chun Huang
!
! abstract:
!
! program history log:
!   2010-03-31  Huang   - create routine based on read_gfs
!   2010-10-19  Huang   - remove spectral part for gridded NEMS/GFS
!   2011-05-01  todling - cwmr no longer in guess-grids; use metguess bundle now
!   2013-10-19  todling - metguess now holds background
!   2016-03-30  todling - update interface to general read (pass bundle)
!   2016-06-23  Li      - Add cloud partitioning, which was missed (based on GFS
!                         ticket #239, comment 18)
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

    use kinds, only: i_kind,r_kind
    use gridmod, only: sp_a,grd_a,lat2,lon2,nsig
    use guess_grids, only: ifilesig,nfldsig
    use gsi_metguess_mod, only: gsi_metguess_bundle
    use gsi_bundlemod, only: gsi_bundlegetpointer
    use gsi_bundlemod, only: gsi_bundlecreate
    use gsi_bundlemod, only: gsi_grid
    use gsi_bundlemod, only: gsi_gridcreate
    use gsi_bundlemod, only: gsi_bundle
    use gsi_bundlemod, only: gsi_bundledestroy
    use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info,general_sub2grid_destroy_info
    use mpimod, only: npe,mype
    use cloud_efr_mod, only: cloud_calc_gfs,set_cloud_lower_bound
    implicit none

    character(len=*),parameter::myname_=myname//'*read_'
    character(24) filename
    integer(i_kind):: it, istatus, inner_vars, num_fields
    integer(i_kind):: iret_ql,iret_qi

    real(r_kind),pointer,dimension(:,:  ):: ges_ps_it  =>NULL()
    real(r_kind),pointer,dimension(:,:  ):: ges_z_it   =>NULL()
    real(r_kind),pointer,dimension(:,:,:):: ges_u_it   =>NULL()
    real(r_kind),pointer,dimension(:,:,:):: ges_v_it   =>NULL()
    real(r_kind),pointer,dimension(:,:,:):: ges_div_it =>NULL()
    real(r_kind),pointer,dimension(:,:,:):: ges_vor_it =>NULL()
    real(r_kind),pointer,dimension(:,:,:):: ges_tv_it  =>NULL()
    real(r_kind),pointer,dimension(:,:,:):: ges_q_it   =>NULL()
    real(r_kind),pointer,dimension(:,:,:):: ges_oz_it  =>NULL()
    real(r_kind),pointer,dimension(:,:,:):: ges_cwmr_it=>NULL()
    real(r_kind),pointer,dimension(:,:,:):: ges_ql_it  => NULL()
    real(r_kind),pointer,dimension(:,:,:):: ges_qi_it  => NULL()

    type(sub2grid_info) :: grd_t
    logical regional
    logical:: l_cld_derived,zflag,inithead

    type(gsi_bundle) :: atm_bundle
    type(gsi_grid)   :: atm_grid
    integer(i_kind),parameter :: n2d=2
    integer(i_kind),parameter :: n3d=8
    character(len=4), parameter :: vars2d(n2d) = (/ 'z   ', 'ps  ' /)
    character(len=4), parameter :: vars3d(n3d) = (/ 'u   ', 'v   ', &
                                                    'vor ', 'div ', &
                                                    'tv  ', 'q   ', &
                                                    'cw  ', 'oz  ' /)
    real(r_kind),pointer,dimension(:,:):: ptr2d   =>NULL()
    real(r_kind),pointer,dimension(:,:,:):: ptr3d =>NULL()

    regional=.false.
    inner_vars=1
    num_fields=min(8*grd_a%nsig+2,npe)
!  Create temporary communication information fore read routines
    call general_sub2grid_create_info(grd_t,inner_vars,grd_a%nlat,grd_a%nlon, &
          grd_a%nsig,num_fields,regional)

!   Allocate bundle used for reading members
    call gsi_gridcreate(atm_grid,lat2,lon2,nsig)
    call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d)
    if(istatus/=0) then
      write(6,*) myname_,': trouble creating atm_bundle'
      call stop2(999)
    endif

    do it=1,nfldsig

       write(filename,'(''sigf'',i2.2)') ifilesig(it)

!      Read background fields into bundle
       call general_read_gfsatm_nems(grd_t,sp_a,filename,.true.,.true.,.true.,&
            atm_bundle,.true.,istatus)

       inithead=.false.
       zflag=.false.

!      Set values to actual MetGuess fields
       call set_guess_

       l_cld_derived = associated(ges_cwmr_it).and.&
                       associated(ges_q_it)   .and.&
                       associated(ges_ql_it)  .and.&
                       associated(ges_qi_it)  .and.&
                       associated(ges_tv_it)
!      call set_cloud_lower_bound(ges_cwmr_it)
       if (mype==0) write(6,*)'READ_GFS_NEMS: l_cld_derived = ', l_cld_derived

       if (l_cld_derived) then
          call cloud_calc_gfs(ges_ql_it,ges_qi_it,ges_cwmr_it,ges_q_it,ges_tv_it,.true.)
       end if

    end do
    call general_sub2grid_destroy_info(grd_t)
    call gsi_bundledestroy(atm_bundle,istatus)

    contains

    subroutine set_guess_

    call gsi_bundlegetpointer (atm_bundle,'ps',ptr2d,istatus)
    if (istatus==0) then
       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ps',ges_ps_it ,istatus)
       if(istatus==0) ges_ps_it = ptr2d
    endif
    call gsi_bundlegetpointer (atm_bundle,'z',ptr2d,istatus)
    if (istatus==0) then
       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'z' ,ges_z_it ,istatus)
       if(istatus==0) ges_z_it = ptr2d
    endif
    call gsi_bundlegetpointer (atm_bundle,'u',ptr3d,istatus)
    if (istatus==0) then
       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'u' ,ges_u_it ,istatus)
       if(istatus==0) ges_u_it = ptr3d
    endif
    call gsi_bundlegetpointer (atm_bundle,'v',ptr3d,istatus)
    if (istatus==0) then
       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'v' ,ges_v_it ,istatus)
       if(istatus==0) ges_v_it = ptr3d
    endif
    call gsi_bundlegetpointer (atm_bundle,'vor',ptr3d,istatus)
    if (istatus==0) then
       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'vor',ges_vor_it,istatus)
       if(istatus==0) ges_vor_it = ptr3d
    endif
    call gsi_bundlegetpointer (atm_bundle,'div',ptr3d,istatus)
    if (istatus==0) then
       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'div',ges_div_it,istatus)
       if(istatus==0) ges_div_it = ptr3d
    endif
    call gsi_bundlegetpointer (atm_bundle,'tv',ptr3d,istatus)
    if (istatus==0) then
       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'tv',ges_tv_it ,istatus)
       if(istatus==0) ges_tv_it = ptr3d
    endif
    call gsi_bundlegetpointer (atm_bundle,'q',ptr3d,istatus)
    if (istatus==0) then
       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'q' ,ges_q_it ,istatus)
       if(istatus==0) ges_q_it = ptr3d
    endif
    call gsi_bundlegetpointer (atm_bundle,'oz',ptr3d,istatus)
    if (istatus==0) then
       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'oz',ges_oz_it ,istatus)
       if(istatus==0) ges_oz_it = ptr3d
    endif
    call gsi_bundlegetpointer (atm_bundle,'cw',ptr3d,istatus)
    if (istatus==0) then
       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'cw',ges_cwmr_it,istatus)
       if(istatus==0) ges_cwmr_it = ptr3d
    endif
    call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ql',ges_ql_it,  iret_ql)
    call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qi',ges_qi_it,  iret_qi)
    if (iret_ql/=0) then
       if (mype==0) write(6,*)'READ_ NEMSIO: cannot get pointer to ql,iret_ql=',iret_ql
    endif
    if (iret_qi/=0) then
       if (mype==0) write(6,*)'READ_ NEMSIO: cannot get pointer to qi,iret_qi=',iret_qi
    endif

  end subroutine set_guess_

  end subroutine read_

  subroutine read_chem_ ( iyear, month,idd )
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    read_nems_chem
!
!   prgrmmr: todling
!
! abstract: fills chemguess_bundle with GFS chemistry.
!
! remarks:
!    1. Right now, only CO2 is done and even this is treated
!        as constant througout the assimialation window.
!    2. iyear and month could come from obsmod, but logically
!       this program should never depend on obsmod
!
!
! program history log:
!   2010-12-23  Huang   - initial code, based on read_gfs_chem
!   2011-06-29  todling - no explict reference to internal bundle arrays
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

    use kinds,   only: i_kind, r_kind
    use mpimod,  only: mype
    use gridmod, only: lat2,lon2,nsig,nlat,rlats,istart
    use ncepgfs_ghg, only: read_gfsco2
    use guess_grids, only: nfldsig
    use gsi_bundlemod, only: gsi_bundlegetpointer
    use gsi_chemguess_mod, only: gsi_chemguess_bundle
    use gsi_chemguess_mod, only: gsi_chemguess_get

    implicit none

!   Declared argument list
    integer(i_kind), intent(in):: iyear
    integer(i_kind), intent(in):: month
    integer(i_kind), intent(in):: idd

!   Declare local variables
    integer(i_kind) :: igfsco2, i, j, n, iret
    real(r_kind),dimension(lat2):: xlats
    real(r_kind),pointer,dimension(:,:,:)::p_co2=>NULL()
    real(r_kind),pointer,dimension(:,:,:)::ptr3d=>NULL()

    if(.not.associated(gsi_chemguess_bundle)) return
    call gsi_bundlegetpointer(gsi_chemguess_bundle(1),'co2',p_co2,iret)
    if(iret /= 0) return

!   Get subdomain latitude array
    j = mype + 1
    do i = 1, lat2
       n = min(max(1, istart(j)+i-2), nlat)
       xlats(i) = rlats(n)
    enddo

!   Read in CO2
    call gsi_chemguess_get ( 'i4crtm::co2', igfsco2, iret )
    call read_gfsco2 ( iyear,month,idd,igfsco2,xlats,&
                       lat2,lon2,nsig,mype, p_co2 )

! Approximation: setting all times co2 values equal to the daily co2 values

    do n = 2, nfldsig
       call gsi_bundlegetpointer(gsi_chemguess_bundle(n),'co2',ptr3d,iret)
       ptr3d = p_co2
    enddo

  end subroutine read_chem_

  subroutine read_atm_ (grd,filename,sp_a,uvflag,vordivflag,zflag, &
       g_z,g_ps,g_vor,g_div,g_u,g_v,&
       g_tv,g_q,g_cwmr,g_oz)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_nemsatm    read nems atm and send to all mpi tasks
!   prgmmr: Huang            org: np23                date: 2010-02-22
!
! abstract: read ncep nems/gfs atmospheric guess field and
!           scatter to subdomains
!
! program history log:
!   2010-02-22 Huang    Initial version.  Based on sub read_gfsatm
!   2011-02-28 Huang    Re-arrange the read sequence to be same as model
!                       write sequence.  Alsom allocate and deallocate
!                       temporary working array immediatelt before and after
!                       the processing and scattering first guess field to reduce
!                       maximum resident memory size.  Page fault can happen
!                       when running at high resolution GSI, e.g., T574.
!   2011-09-23 Huang    Add NEMS parallel IO capability
!   2013-10-25 todling  reposition fill_ns,filluv_ns to commvars
!
!   input argument list:
!     grd      - structure variable containing information about grid
!                    (initialized by general_sub2grid_create_info, located in
!                    general_sub2grid_mod.f90)
!     sp_a     - structure variable containing spectral information for analysis
!                    (initialized by general_init_spec_vars, located in
!                    general_specmod.f90)
!     uvflag   - logical to use u,v (.true.) or st,vp (.false.) perturbations
!     vordivflag - logical to determine if routine should output vorticity and
!                  divergence
!     zflag    - logical to determine if surface height field should be output
!
!   output argument list:
!     g_*      - guess fields
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use kinds, only: r_kind,i_kind
    use gridmod, only: ntracer,ncloud,reload,itotsub
    use general_commvars_mod, only: fill_ns,filluv_ns,fill2_ns,filluv2_ns,ltosj_s,ltosi_s
    use general_specmod, only: spec_vars
    use general_sub2grid_mod, only: sub2grid_info
    use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype,mype
    use nemsio_module, only: nemsio_init,nemsio_open,nemsio_close
    use nemsio_module, only: nemsio_gfile,nemsio_getfilehead,nemsio_readrecv
    use egrid2agrid_mod,only: g_egrid2agrid,g_create_egrid2agrid,egrid2agrid_parm,destroy_egrid2agrid
    use constants, only: two,pi,half,deg2rad
    use control_vectors, only: imp_physics
    implicit none

!   Declare local parameters
    real(r_kind),parameter:: r0_001 = 0.001_r_kind

!   Declare passed variables
    type(sub2grid_info)                   ,intent(in   ) :: grd
    character(len=24)                     ,intent(in   ) :: filename
    logical                               ,intent(in   ) :: uvflag,vordivflag,zflag
    real(r_kind),dimension(grd%lat2,grd%lon2)     ,intent(  out) :: g_z,g_ps
    real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent(  out) :: g_u,g_v,&
         g_vor,g_div,g_cwmr,g_q,g_oz,g_tv
    type(spec_vars)                       ,intent(in   ) :: sp_a

!   Declare local variables
    character(len=120) :: my_name = 'READ_NEMSATM'
    character(len=1)   :: null = ' '
    integer(i_kind),dimension(7):: idate
    integer(i_kind),dimension(4):: odate
    integer(i_kind) :: iret,nlatm2,nflds
    integer(i_kind) :: k,icount,icount_prev,mm1,i,j,kk
    integer(i_kind) :: mype_hs, mype_ps,nord_int
    integer(i_kind) :: latb, lonb, levs, nframe
    integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd
    integer(i_kind) :: istop = 101
    real(r_kind),allocatable,dimension(:,:) :: grid, grid_v, &
         grid_vor, grid_div, grid_b, grid_b2
    real(r_kind),allocatable,dimension(:,:,:) :: grid_c, grid2, grid_c2
    real(r_kind),allocatable,dimension(:)   :: work, work_vor, work_div, &
         work_v
    real(r_kind),allocatable,dimension(:,:) :: sub, sub_vor, sub_div, &
         sub_v
    real(r_kind),dimension(sp_a%nc):: spec_vor,spec_div
    real(r_kind),allocatable,dimension(:) :: rwork1d0, rwork1d1, rwork1d2
    real(r_kind),allocatable,dimension(:) :: rlats,rlons,clons,slons
    real(4),allocatable,dimension(:) :: r4lats,r4lons
    real(r_kind) :: fhour
    type(nemsio_gfile) :: gfile
    logical diff_res,eqspace
    logical,dimension(1) :: vector
    type(egrid2agrid_parm) :: p_high

!******************************************************************************
!   Initialize variables used below
    mm1=mype+1
    mype_hs=min(1,npe-1)
    mype_ps=0
    nlatm2=grd%nlat-2
    nflds=5*grd%nsig+1
    if(zflag) nflds=nflds+1
    if(vordivflag .or. .not. uvflag)nflds=nflds+2*grd%nsig
!   nflds=npe
    nflds=grd%nsig
    levs=grd%nsig

    allocate( work(grd%itotsub),work_v(grd%itotsub) )
    work=zero
    work_v=zero
    allocate( sub(grd%lat2*grd%lon2,max(grd%nsig,npe)),sub_v(grd%lat2*grd%lon2,max(grd%nsig,npe)) )
    allocate( sub_div(grd%lat2*grd%lon2,max(grd%nsig,npe)),sub_vor(grd%lat2*grd%lon2,max(grd%nsig,npe)) )
    if(mype < nflds)then

      call nemsio_init(iret=iret)
      if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'init',istop,iret)

      call nemsio_open(gfile,filename,'READ',iret=iret)
      if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'open',istop+1,iret)

      call nemsio_getfilehead(gfile,iret=iret, nframe=nframe, &
           nfhour=nfhour, nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd, &
           idate=idate, dimx=lonb, dimy=latb,dimz=levs)

      if( nframe /= 0 ) then
         if ( mype == 0 ) &
         write(6,*)trim(my_name),': ***ERROR***  nframe /= 0 for global model read, nframe = ', nframe
         call stop2(101)
      end if

      fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
      odate(1) = idate(4)  !hour
      odate(2) = idate(2)  !month
      odate(3) = idate(3)  !day
      odate(4) = idate(1)  !year
!
!  g_* array already pre-allocate as (lat2,lon2,<nsig>) => 2D and <3D> array
!
      diff_res=.false.
      if(latb /= nlatm2) then
         diff_res=.true.
         if ( mype == 0 ) write(6, &
            '(a,'': different spatial dimension nlatm2 = '',i4,tr1,''latb = '',i4)') &
            trim(my_name),nlatm2,latb
  !      call stop2(101)
      end if
      if(lonb /= grd%nlon) then
         diff_res=.true.
         if ( mype == 0 ) write(6, &
            '(a,'': different spatial dimension nlon   = '',i4,tr1,''lonb = '',i4)') &
            trim(my_name),grd%nlon,lonb
  !      call stop2(101)
      end if
      if(levs /= grd%nsig)then
         if ( mype == 0 ) write(6, &
            '(a,'': inconsistent spatial dimension nsig   = '',i4,tr1,''levs = '',i4)') &
            trim(my_name),grd%nsig,levs
         call stop2(101)
      end if
!
      allocate( grid(grd%nlon,nlatm2), grid_v(grd%nlon,nlatm2) )
      if(diff_res)then
         allocate( grid_b(lonb,latb),grid_c(latb+2,lonb,1),grid2(grd%nlat,grd%nlon,1))
         allocate( grid_b2(lonb,latb),grid_c2(latb+2,lonb,1))
      end if
      allocate( rwork1d0(latb*lonb) )
      allocate( rlats(latb+2),rlons(lonb),clons(lonb),slons(lonb),r4lats(lonb*latb),r4lons(lonb*latb))
      allocate(rwork1d1(latb*lonb))
      call nemsio_getfilehead(gfile,lat=r4lats,iret=iret)
      call nemsio_getfilehead(gfile,lon=r4lons,iret=iret)
      do j=1,latb
        rlats(latb+2-j)=deg2rad*r4lats(lonb/2+(j-1)*lonb)
      end do
      do j=1,lonb
        rlons(j)=deg2rad*r4lons(j)
      end do
      deallocate(r4lats,r4lons)
      rlats(1)=-half*pi
      rlats(latb+2)=half*pi
      do j=1,lonb
         clons(j)=cos(rlons(j))
         slons(j)=sin(rlons(j))
      end do

      nord_int=4
      eqspace=.false.
      call g_create_egrid2agrid(grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons, &
                              latb+2,rlats,lonb,rlons,&
                              nord_int,p_high,.true.,eqspace)
      deallocate(rlats,rlons)
    end if
!
!   Load values into rows for south and north pole before scattering
!
!   Terrain:  scatter to all mpi tasks
!
    if(zflag)then
       if (mype==mype_hs) then
          call nemsio_readrecv(gfile,'hgt', 'sfc',1,rwork1d0,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'hgt','read',istop+2,iret)
          if(diff_res)then
             grid_b=reshape(rwork1d0,(/size(grid_b,1),size(grid_b,2)/))
             vector(1)=.false.
             call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
             call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
             do kk=1,itotsub
               i=ltosi_s(kk)
               j=ltosj_s(kk)
               work(kk)=grid2(i,j,1)
             end do
          else
             grid=reshape(rwork1d0,(/size(grid,1),size(grid,2)/))
             call fill_ns(grid,work)
          end if
       endif
       call mpi_scatterv(work,grd%ijn_s,grd%displs_s,mpi_rtype,&
          g_z,grd%ijn_s(mm1),mpi_rtype,mype_hs,mpi_comm_world,ierror)
    end if

!   Surface pressure:  same procedure as terrain, but handled by task mype_ps
!
    if (mype==mype_ps) then
       call nemsio_readrecv(gfile,'pres','sfc',1,rwork1d0,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'pres','read',istop+3,iret)
       rwork1d1 = r0_001*rwork1d0
       if(diff_res)then
          vector(1)=.false.
          grid_b=reshape(rwork1d1,(/size(grid_b,1),size(grid_b,2)/))
          call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
          call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
          do kk=1,itotsub
            i=ltosi_s(kk)
            j=ltosj_s(kk)
            work(kk)=grid2(i,j,1)
          end do
       else
          grid=reshape(rwork1d1,(/size(grid,1),size(grid,2)/)) ! convert Pa to cb
          call fill_ns(grid,work)
       endif
    endif
    call mpi_scatterv(work,grd%ijn_s,grd%displs_s,mpi_rtype,&
       g_ps,grd%ijn_s(mm1),mpi_rtype,mype_ps,mpi_comm_world,ierror)

!   Divergence and voriticity.  Compute u and v from div and vor
    sub_vor=zero
    sub_div=zero
    sub    =zero
    sub_v  =zero
    icount     =0
    icount_prev=1
    allocate( work_vor(grd%itotsub),work_div(grd%itotsub) )
    do k=1,levs
       icount=icount+1
       if (mype==mod(icount-1,npe)) then
          ! Convert grid u,v to div and vor
          call nemsio_readrecv(gfile,'ugrd','mid layer',k,rwork1d0,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ugrd','read',istop+4,iret)
          call nemsio_readrecv(gfile,'vgrd','mid layer',k,rwork1d1,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'vgrd','read',istop+5,iret)
          if(diff_res)then
             grid_b=reshape(rwork1d0,(/size(grid_b,1),size(grid_b,2)/))
             grid_b2=reshape(rwork1d1,(/size(grid_b,1),size(grid_b,2)/))
             vector(1)=.true.
             call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons)
             call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
             do kk=1,itotsub
               i=ltosi_s(kk)
               j=ltosj_s(kk)
               work(kk)=grid2(i,j,1)
             end do
             do j=1,grd%nlon
               do i=2,grd%nlat-1
                 grid(j,grd%nlat-i)=grid2(i,j,1)
               end do
             end do
             call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector)
             do kk=1,itotsub
               i=ltosi_s(kk)
               j=ltosj_s(kk)
               work_v(kk)=grid2(i,j,1)
             end do
             do j=1,grd%nlon
               do i=2,grd%nlat-1
                 grid_v(j,grd%nlat-i)=grid2(i,j,1)
               end do
             end do
          else
             grid=reshape(rwork1d0,(/size(grid,1),size(grid,2)/))
             grid_v=reshape(rwork1d1,(/size(grid_v,1),size(grid_v,2)/))
             call filluv_ns(grid,grid_v,work,work_v)
          end if

          if(vordivflag .or. .not. uvflag)then

             allocate( grid_vor(grd%nlon,nlatm2), grid_div(grd%nlon,nlatm2) )
             call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1)
             call general_sptez_s_b(sp_a,sp_a,spec_div,grid_div,1)
             call general_sptez_s_b(sp_a,sp_a,spec_vor,grid_vor,1)

          ! Load values into rows for south and north pole
             call fill_ns(grid_div,work_div)
             call fill_ns(grid_vor,work_vor)
             deallocate(grid_vor,grid_div)
          end if
       endif
       ! Scatter to sub
       if (mod(icount,npe)==0 .or. icount==levs) then
          if(vordivflag .or. .not. uvflag)then
             call mpi_alltoallv(work_vor,grd%ijn_s,grd%displs_s,mpi_rtype,&
                sub_vor(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
                mpi_comm_world,ierror)
             call mpi_alltoallv(work_div,grd%ijn_s,grd%displs_s,mpi_rtype,&
                sub_div(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
                mpi_comm_world,ierror)
          end if
          if(uvflag)then
             call mpi_alltoallv(work,grd%ijn_s,grd%displs_s,mpi_rtype,&
                sub(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
                mpi_comm_world,ierror)
             call mpi_alltoallv(work_v,grd%ijn_s,grd%displs_s,mpi_rtype,&
                sub_v(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
                mpi_comm_world,ierror)
          end if
          icount_prev=icount+1
       endif
    end do
    deallocate(work_vor,work_div)

    ! Transfer vor,div,u,v into real(r_kind) guess arrays
    call reload(sub_vor,g_vor)
    call reload(sub_div,g_div)
    call reload(sub,g_u)
    call reload(sub_v,g_v)
    deallocate(sub_vor,sub_div)

!   Thermodynamic variable and Specific humidity:  communicate to all tasks
!
    sub=zero
    icount=0
    icount_prev=1
    do k=1,levs
       icount=icount+1
       if (mype==mod(icount-1,npe)) then

          call nemsio_readrecv(gfile,'spfh','mid layer',k,rwork1d0,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'spfh','read',istop+6,iret)
          if(diff_res)then
             grid_b=reshape(rwork1d0,(/size(grid_b,1),size(grid_b,2)/))
             vector(1)=.false.
             call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
             call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
             do kk=1,itotsub
               i=ltosi_s(kk)
               j=ltosj_s(kk)
               work(kk)=grid2(i,j,1)
             end do
          else
             grid=reshape(rwork1d0,(/size(grid,1),size(grid,2)/))
             call fill_ns(grid,work)
          end if

          call nemsio_readrecv(gfile,'tmp','mid layer',k,rwork1d1,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','read',istop+7,iret)
          allocate(rwork1d2(latb*lonb))
          rwork1d2 = rwork1d1*(one+fv*rwork1d0)
          if(diff_res)then
             grid_b=reshape(rwork1d2,(/size(grid_b,1),size(grid_b,2)/))
             vector(1)=.false.
             call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
             call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
             do kk=1,itotsub
               i=ltosi_s(kk)
               j=ltosj_s(kk)
               work_v(kk)=grid2(i,j,1)
             end do
          else
             grid_v=reshape(rwork1d2,(/size(grid_v,1),size(grid_v,2)/))
             call fill_ns(grid_v,work_v)
          end if

          deallocate(rwork1d2)
       endif

       if (mod(icount,npe)==0 .or. icount==levs) then
          call mpi_alltoallv(work_v,grd%ijn_s,grd%displs_s,mpi_rtype,&
             sub_v(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
             mpi_comm_world,ierror)
          call mpi_alltoallv(work,grd%ijn_s,grd%displs_s,mpi_rtype,&
             sub(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
             mpi_comm_world,ierror)
          icount_prev=icount+1
       endif
    end do
    call reload(sub_v,g_tv)
    call reload(sub,g_q)
    deallocate(sub_v,work_v)

    sub=zero
    icount=0
    icount_prev=1
    do k=1,levs
       icount=icount+1
       if (mype==mod(icount-1,npe)) then
          call nemsio_readrecv(gfile,'o3mr','mid layer',k,rwork1d0,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'o3mr','read',istop+8,iret)
          if(diff_res)then
             grid_b=reshape(rwork1d0,(/size(grid_b,1),size(grid_b,2)/))
             vector(1)=.false.
             call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
             call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
             do kk=1,itotsub
               i=ltosi_s(kk)
               j=ltosj_s(kk)
               work(kk)=grid2(i,j,1)
             end do
          else
             grid=reshape(rwork1d0,(/size(grid,1),size(grid,2)/))
             call fill_ns(grid,work)
          end if
       endif
       if (mod(icount,npe)==0 .or. icount==levs) then
          call mpi_alltoallv(work,grd%ijn_s,grd%displs_s,mpi_rtype,&
             sub(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
             mpi_comm_world,ierror)
          icount_prev=icount+1
       endif
    end do
    call reload(sub,g_oz)

!   Cloud condensate mixing ratio.

    if (ntracer>2 .or. ncloud>=1) then
       sub=zero
       icount=0
       icount_prev=1
       do k=1,levs
          icount=icount+1
          if (mype==mod(icount-1,npe)) then
             call nemsio_readrecv(gfile,'clwmr','mid layer',k,rwork1d0,iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','read',istop+9,iret)
             if (imp_physics == 11) then
                call nemsio_readrecv(gfile,'icmr','mid layer',k,rwork1d1,iret=iret)
                if (iret == 0) then
                   rwork1d0 = rwork1d0 + rwork1d1
                else
                   call error_msg(trim(my_name),trim(filename),'icmr','read',istop+10,iret)
                endif
             endif
             if(diff_res)then
                grid_b=reshape(rwork1d0,(/size(grid_b,1),size(grid_b,2)/))
                vector(1)=.false.
                call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
                call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
                do kk=1,itotsub
                  i=ltosi_s(kk)
                  j=ltosj_s(kk)
                  work(kk)=grid2(i,j,1)
                end do
             else
                grid=reshape(rwork1d0,(/size(grid,1),size(grid,2)/))
                call fill_ns(grid,work)
             end if

          endif
          if (mod(icount,npe)==0 .or. icount==levs) then
             call mpi_alltoallv(work,grd%ijn_s,grd%displs_s,mpi_rtype,&
                sub(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
                mpi_comm_world,ierror)
             icount_prev=icount+1
          endif
       end do
       call reload(sub,g_cwmr)
    else
       g_cwmr = zero
    endif

    if(mype < nflds)then
       if(diff_res) deallocate(grid_b,grid_b2,grid_c,grid_c2,grid2)
       call destroy_egrid2agrid(p_high)
       deallocate(rwork1d1,clons,slons)
       deallocate(rwork1d0)
       deallocate(grid,grid_v)
       call nemsio_close(gfile,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop+9,iret)
    end if
    deallocate(work,sub)

!   Print date/time stamp
    if ( mype == 0 ) write(6, &
       '(a,'': ges read/scatter,lonb,latb,levs= '',3i6,'',hour= '',f4.1,'',idate= '',4i5)') &
       trim(my_name),lonb,latb,levs,fhour,odate

  end subroutine read_atm_

  subroutine read_sfc_(sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, &
                           veg_type,soil_type,terrain,isli,use_sfc_any, &
                           tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d) 
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_sfc_     read nems hist file
!   prgmmr: Li            org: np23                date: 2017-08-30
!
! abstract: read nems sfc & nst combined file
!
! program history log:
!  
!   input argument list:
!     use_sfc_any - true if any processor uses extra surface fields
!
!   output argument list:
!     sfct      - surface temperature (skin temp)
!     soil_moi  - soil moisture of first layer
!     sno       - snow depth
!     soil_temp - soil temperature of first layer
!     veg_frac  - vegetation fraction
!     fact10    - 10 meter wind factor
!     sfc_rough - surface roughness
!     veg_type  - vegetation type
!     soil_type - soil type
!     terrain   - terrain height
!     isli      - sea/land/ice mask
!     tref      - optional, oceanic foundation temperature
!     dt_cool   - optional, sub-layer cooling amount at sub-skin layer
!     z_c       - optional, depth of sub-layer cooling layer
!     dt_warm   - optional, diurnal warming amount at sea surface 
!     z_w       - optional, depth of diurnal warming layer
!     c_0       - optional, coefficient to calculate d(Tz)/d(tr) in dimensionless
!     c_d       - optional, coefficient to calculate d(Tz)/d(tr) in m^-1
!     w_0       - optional, coefficient to calculate d(Tz)/d(tr) in dimensionless
!     w_d       - optional, coefficient to calculate d(Tz)/d(tr) in m^-1

!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use mpimod, only: mype
    use kinds, only: r_kind,i_kind,r_single
    use gridmod, only: nlat_sfc,nlon_sfc
    use guess_grids, only: nfldsfc,ifilesfc
    use constants, only: zero,two
    use nemsio_module, only:  nemsio_init,nemsio_open,nemsio_close
    use nemsio_module, only:  nemsio_gfile,nemsio_getfilehead,nemsio_readrecv
    implicit none

!   Declare passed variables
    logical,                                               intent(in)  :: use_sfc_any
    real(r_single),  dimension(nlat_sfc,nlon_sfc,nfldsfc), intent(out) :: sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough
    real(r_single),  dimension(nlat_sfc,nlon_sfc),         intent(out) :: veg_type,soil_type,terrain
    integer(i_kind), dimension(nlat_sfc,nlon_sfc),         intent(out) :: isli
    real(r_single),  optional, dimension(nlat_sfc,nlon_sfc,nfldsfc), intent(out) :: tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d
                                
!   Declare local parameters
    integer(i_kind), parameter :: nsfc_all=11
    integer(i_kind),dimension(7):: idate
    integer(i_kind),dimension(4):: odate
!   Declare local variables
    real(r_single),  dimension(nlat_sfc,nlon_sfc,nfldsfc) :: xt
    character(len=24)  :: filename
    character(len=120) :: my_name = 'READ_SFCNST'
    character(len=1)   :: null = ' '
    integer(i_kind) :: i,j,it,n,nsfc
    integer(i_kind) :: iret, nframe, lonb, latb
    integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd
    real(r_single)  :: fhour
    integer(i_kind) :: istop = 102
    real(r_single), allocatable, dimension(:)   :: rwork2d
    real(r_single), allocatable, dimension(:,:) :: work,outtmp
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!  Define read variable property   !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    type(nemsio_gfile) :: gfile
!-----------------------------------------------------------------------------
    call nemsio_init(iret=iret)
    if (iret /= 0) call error_msg(trim(my_name),null,null,'init',istop,iret)

    do it = 1, nfldsfc
! read a surface history file on the task
       write(filename,200)ifilesfc(it)
200    format('sfcf',i2.2)

       call nemsio_open(gfile,filename,'READ',iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret)

       call nemsio_getfilehead(gfile, idate=idate, iret=iret, nframe=nframe,   &
          nfhour=nfhour, nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd, &
          dimx=lonb, dimy=latb )

       if( nframe /= 0 ) then
          if ( mype == 0 ) &
          write(6,*)trim(my_name),': ***ERROR***  nframe /= 0 for global sfc hist read, nframe = ', nframe
          call stop2(102)
       end if

       fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
       odate(1) = idate(4)  !hour
       odate(2) = idate(2)  !month
       odate(3) = idate(3)  !day
       odate(4) = idate(1)  !year

       if ( (latb /= nlat_sfc-2) .or. (lonb /= nlon_sfc) ) then
          if ( mype == 0 ) write(6, &
             '(a,'': inconsistent spatial dimension '',''nlon,nlatm2 = '',2(i4,tr1),''-vs- sfc file lonb,latb = '',i4)') &
             trim(my_name),nlon_sfc,nlat_sfc-2,lonb,latb
          call stop2(102)
       endif
!
!      Read the surface records (lonb, latb)  and convert to GSI array pattern (nlat_sfc,nlon_sfc)
!      Follow the read order sfcio in ncepgfs_io
!
       allocate(work(lonb,latb))
       allocate(rwork2d(size(work,1)*size(work,2)))
       work    = zero
       rwork2d = zero

       if(it == 1)then
         nsfc=nsfc_all
       else
         nsfc=nsfc_all-4
       end if

       do n = 1, nsfc

          if (n == 1) then                               ! skin temperature

!            Tsea
             call nemsio_readrecv(gfile, 'tmp', 'sfc', 1, rwork2d, iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','read',istop,iret)
             work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
             call tran_gfssfc(work,sfct(1,1,it),lonb,latb)

          elseif(n == 2 .and. use_sfc_any) then          ! soil moisture

!            smc/soilw
             call nemsio_readrecv(gfile, 'smc', 'soil layer', 1, rwork2d, iret=iret)
             ! FV3 nemsio files use 'soilw 0-10cm down' insted of 'smc soil layer 1'
             if (iret /= 0) then
               if ( mype == 0 ) print *,'could not read smc, try to read soilw 0-10 cm down instead...'
               call nemsio_readrecv(gfile,'soilw','0-10 cm down',1,rwork2d,iret=iret)
               if (iret /= 0) &
               call error_msg(trim(my_name),trim(filename),'smc/soilw','read',istop,iret)
             endif
             work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
             call tran_gfssfc(work,soil_moi(1,1,it),lonb,latb)

          elseif(n == 3) then                            ! snow depth

!            sheleg
             call nemsio_readrecv(gfile, 'weasd','sfc', 1, rwork2d, iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'weasd','read',istop,iret)
             work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
             call tran_gfssfc(work,sno(1,1,it),lonb,latb)

          elseif(n == 4 .and. use_sfc_any) then          ! soil temperature

!            stc/tmp
             call nemsio_readrecv(gfile, 'stc', 'soil layer', 1, rwork2d, iret=iret)
             if (iret /= 0) then
               ! FV3 nemsio files use 'tmp 0-10cm down' insted of 'stc soil layer 1'
               if ( mype == 0 ) print *,'could not read stc, try to read tmp 0-10 cm down instead...'
               call nemsio_readrecv(gfile,'tmp','0-10 cm down',1,rwork2d,iret=iret)
               if (iret /= 0) &
               call error_msg(trim(my_name),trim(filename),'stc/tmp','read',istop,iret)
             endif
             work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
             call tran_gfssfc(work,soil_temp(1,1,it),lonb,latb)

          elseif(n == 5 .and. use_sfc_any) then          ! vegetation cover

!            vfrac
             call nemsio_readrecv(gfile, 'veg',  'sfc', 1, rwork2d, iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'veg','read',istop,iret)
             work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
             call tran_gfssfc(work,veg_frac(1,1,it),lonb,latb)

          elseif(n == 6) then                            ! 10m wind factor

!            f10m
             call nemsio_readrecv(gfile, 'f10m', '10 m above gnd', 1, rwork2d, iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'f10m','read',istop,iret)
             work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
             call tran_gfssfc(work,fact10(1,1,it),lonb,latb)

          elseif(n == 7) then                            ! suface roughness

!            zorl
             call nemsio_readrecv(gfile, 'sfcr', 'sfc', 1, rwork2d,iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'sfcr','read',istop,iret)
             work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
             call tran_gfssfc(work,sfc_rough(1,1,it),lonb,latb)

          elseif(n == 8 .and. use_sfc_any) then          ! vegetation type

!            vtype
             call nemsio_readrecv(gfile, 'vtype','sfc', 1, rwork2d, iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'vtype','read',istop,iret)
             work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
             call tran_gfssfc(work,veg_type,lonb,latb)

          elseif(n == 9 .and. use_sfc_any) then          ! soil type

!            stype
             call nemsio_readrecv(gfile, 'sotyp','sfc', 1, rwork2d, iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'sotyp','read',istop,iret)
             work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
             call tran_gfssfc(work,soil_type,lonb,latb)

          elseif(n == 10) then                           ! terrain

!            orog
             call nemsio_readrecv(gfile, 'orog', 'sfc', 1, rwork2d,iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'orog','read',istop,iret)
             work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
             call tran_gfssfc(work,terrain,lonb,latb)

          elseif(n == 11) then                           ! sea/land/ice flag

!            slmsk
             call nemsio_readrecv(gfile, 'land', 'sfc', 1, rwork2d, iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'land','read',istop,iret)
             work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
             allocate(outtmp(latb+2,lonb))
             call tran_gfssfc(work,outtmp,lonb,latb)
             do j=1,lonb
                do i=1,latb+2
                   isli(i,j) = nint(outtmp(i,j))
                end do
             end do
             deallocate(outtmp)

          endif
!      End of loop over data records
       enddo

       if( present(tref) ) then                         
          if ( mype == 0 ) write(6,*) ' read 9 optional NSST variables '

          call nemsio_readrecv(gfile, 'tref', 'sfc', 1, rwork2d, iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tref','read',istop,iret)
          work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
          call tran_gfssfc(work,tref(1,1,it),lonb,latb)

          call nemsio_readrecv(gfile, 'dtcool','sfc', 1, rwork2d, iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dtcool','read',istop,iret)
          work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
          call tran_gfssfc(work,dt_cool(1,1,it),lonb,latb)

          call nemsio_readrecv(gfile, 'zc',  'sfc', 1, rwork2d,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'z_c','read',istop,iret)
          work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
          call tran_gfssfc(work,z_c(1,1,it),lonb,latb)

          call nemsio_readrecv(gfile, 'xt',   'sfc', 1, rwork2d, iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'xt','read',istop,iret)
          work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
          call tran_gfssfc(work,xt(1,1,it),lonb,latb)

          call nemsio_readrecv(gfile, 'xz',   'sfc', 1, rwork2d, iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'xz','read',istop,iret)
          work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
          call tran_gfssfc(work,z_w(1,1,it),lonb,latb)
 
          call nemsio_readrecv(gfile, 'c0',  'sfc', 1, rwork2d, iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'c0','read',istop,iret)
          work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
          call tran_gfssfc(work,c_0(1,1,it),lonb,latb)

          call nemsio_readrecv(gfile, 'cd',  'sfc', 1, rwork2d, iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'cd','read',istop,iret)
          work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
          call tran_gfssfc(work,c_d(1,1,it),lonb,latb)

          call nemsio_readrecv(gfile, 'w0',  'sfc', 1, rwork2d, iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'w0','read',istop,iret)
          work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
          call tran_gfssfc(work,w_0(1,1,it),lonb,latb)

          call nemsio_readrecv(gfile, 'wd',  'sfc', 1, rwork2d,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'wd','read',istop,iret)
          work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
          call tran_gfssfc(work,w_d(1,1,it),lonb,latb)
!
!         Get diurnal warming amout at z=0
!
          do j = 1,nlon_sfc
             do i = 1,nlat_sfc
                if (z_w(i,j,it) > zero) then
                   dt_warm(i,j,it) = two*xt(i,j,it)/z_w(i,j,it)
                end if
             end do
          end do
       endif
!      Deallocate local work arrays
       deallocate(work,rwork2d)

       call nemsio_close(gfile,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret)
!
!      Print date/time stamp
       if ( mype == 0 ) write(6, &
          '(a,'': sfc read,nlon,nlat= '',2i6,'',hour= '',f4.1,'',idate= '',4i5)') &
          trim(my_name),lonb,latb,fhour,odate
!   End of loop over time levels
    end do
  end subroutine read_sfc_

  subroutine read_nemssfc_(iope,sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, &
                               veg_type,soil_type,terrain,isli,use_sfc_any, &
                               tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_nemssfc_     read nems hist file
!   prgmmr: xuli          org: np23                date: 2017-08-30
!
! abstract: read nems surface file
!
! program history log:
!   2017-08-30  li
!
!   input argument list:
!     iope        - mpi task handling i/o
!     use_sfc_any - true if any processor uses extra surface fields
!
!   output argument list:
!     sfct      - surface temperature (skin temp)
!     soil_moi  - soil moisture of first layer
!     sno       - snow depth
!     soil_temp - soil temperature of first layer
!     veg_frac  - vegetation fraction
!     fact10    - 10 meter wind factor
!     sfc_rough - surface roughness
!     veg_type  - vegetation type
!     soil_type - soil type
!     terrain   - terrain height
!     isli      - sea/land/ice mask
!     tref      - oceanic foundation temperature
!     dt_cool   - optional, sub-layer cooling amount at sub-skin layer
!     z_c       - optional, depth of sub-layer cooling layer
!     dt_warm   - optional, diurnal warming amount at sea surface 
!     z_w       - optional, depth of diurnal warming layer
!     c_0       - optional, coefficient to calculate d(Tz)/d(tf) 
!     c_d       - optional, coefficient to calculate d(Tz)/d(tf)
!     w_0       - optional, coefficient to calculate d(Tz)/d(tf) 
!     w_d       - optional, coefficient to calculate d(Tz)/d(tf) 
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use kinds, only: r_kind,i_kind,r_single
    use gridmod, only: nlat_sfc,nlon_sfc
    use guess_grids, only: nfldsfc,sfcmod_mm5,sfcmod_gfs
    use mpimod, only: mpi_itype,mpi_rtype4,mpi_comm_world,mype
    use constants, only: zero
    implicit none

!   Declare passed variables
    integer(i_kind),                                       intent(in)  :: iope
    logical,                                               intent(in)  :: use_sfc_any
    real(r_single),  dimension(nlat_sfc,nlon_sfc,nfldsfc), intent(out) :: sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough
    real(r_single),  dimension(nlat_sfc,nlon_sfc),         intent(out) :: veg_type,soil_type,terrain
    integer(i_kind), dimension(nlat_sfc,nlon_sfc),         intent(out) :: isli
    real(r_single), optional, dimension(nlat_sfc,nlon_sfc,nfldsfc), intent(out) :: tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d

!   Declare local variables
    integer(i_kind):: iret,npts,nptsall

!-----------------------------------------------------------------------------
!   Read surface history file on processor iope
    if(mype == iope)then
       if ( present(tref) ) then
          call read_sfc_(sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, &
                         veg_type,soil_type,terrain,isli,use_sfc_any, &
                         tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)
          write(*,*) 'read_sfc nemsio, with NSST variables'
       else
          call read_sfc_(sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, &
                         veg_type,soil_type,terrain,isli,use_sfc_any)
          write(*,*) 'read_sfc nemsio, without NSST variables'
       endif
    endif

!   Load onto all processors

    npts=nlat_sfc*nlon_sfc
    nptsall=npts*nfldsfc

    call mpi_bcast(sfct,      nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    call mpi_bcast(fact10,    nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    call mpi_bcast(sno,       nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    if(sfcmod_mm5 .or. sfcmod_gfs)then
       call mpi_bcast(sfc_rough, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    else
       sfc_rough = zero
    endif
    call mpi_bcast(terrain, npts,   mpi_rtype4,iope,mpi_comm_world,iret)
    call mpi_bcast(isli,    npts,   mpi_itype, iope,mpi_comm_world,iret)
    if(use_sfc_any)then
       call mpi_bcast(veg_frac, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
       call mpi_bcast(soil_temp,nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
       call mpi_bcast(soil_moi, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
       call mpi_bcast(veg_type, npts,   mpi_rtype4,iope,mpi_comm_world,iret)
       call mpi_bcast(soil_type,npts,   mpi_rtype4,iope,mpi_comm_world,iret)
    endif
    if ( present(tref) ) then
       call mpi_bcast(tref,    nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
       call mpi_bcast(dt_cool, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
       call mpi_bcast(z_c,     nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
       call mpi_bcast(dt_warm, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
       call mpi_bcast(z_w,     nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
       call mpi_bcast(c_0,     nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
       call mpi_bcast(c_d,     nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
       call mpi_bcast(w_0,     nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
       call mpi_bcast(w_d,     nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    endif

  end subroutine read_nemssfc_

  subroutine read_sfc_anl_(isli_anl)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_sfc_anl_     read nems surface file with analysis resolution
!
!   prgmmr: li            org: np23                date: 2016-08-18
!
! abstract: read nems surface file at analysis grids when nlon /= nlon_sfc or nlat /= nlat_sfc
!
! program history log:
!
!   input argument list:
!
!   output argument list:
!     isli      - sea/land/ice mask
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use mpimod, only: mype
    use kinds, only: r_kind,i_kind,r_single
    use gridmod, only: nlat,nlon
    use constants, only: zero
    use nemsio_module, only:  nemsio_init,nemsio_open,nemsio_close
    use nemsio_module, only:  nemsio_gfile,nemsio_getfilehead,nemsio_readrecv
    implicit none

!   Declare passed variables
    integer(i_kind), dimension(nlat,nlon),   intent(  out) :: isli_anl

!   Declare local parameters
    integer(i_kind),dimension(7):: idate
    integer(i_kind),dimension(4):: odate


!   Declare local variables
    character(len=24)  :: filename
    character(len=120) :: my_name = 'READ_NEMSSFC_ANL'
    character(len=1)   :: null = ' '
    integer(i_kind) :: i,j
    integer(i_kind) :: iret, nframe, lonb, latb
    integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd
    real(r_single) :: fhour
    integer(i_kind) :: istop = 102
    real(r_single), allocatable, dimension(:)   :: rwork2d
    real(r_single), allocatable, dimension(:,:) :: work,outtmp

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!  Define read variable property   !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
    type(nemsio_gfile) :: gfile
!-----------------------------------------------------------------------------

    call nemsio_init(iret=iret)
    if (iret /= 0) call error_msg(trim(my_name),null,null,'init',istop,iret)


    filename='sfcf06_anlgrid'
    call nemsio_open(gfile,trim(filename),'READ',iret=iret)
    if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret)

    call nemsio_getfilehead(gfile, idate=idate, iret=iret, nframe=nframe,   &
       nfhour=nfhour, nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd, &
       dimx=lonb, dimy=latb )

    if( nframe /= 0 ) then
       if ( mype == 0 ) &
       write(6,*)trim(my_name),': ***ERROR***  nframe /= 0 for global model read, nframe = ', nframe
       call stop2(102)
    end if

    fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
    odate(1) = idate(4)  !hour
    odate(2) = idate(2)  !month
    odate(3) = idate(3)  !day
    odate(4) = idate(1)  !year

    if ( (latb /= nlat-2) .or. (lonb /= nlon) ) then
       if ( mype == 0 ) write(6, &
          '(a,'': inconsistent spatial dimension '',''nlon,nlatm2 = '',2(i4,tr1),''-vs- sfc file lonb,latb = '',i4)') &
          trim(my_name),nlon,nlat-2,lonb,latb
       call stop2(102)
    endif
!
!   Read the surface records (lonb, latb)  and convert to GSI array pattern (nlat,nlon)
!   Follow the read order sfcio in ncepgfs_io
!
    allocate(work(lonb,latb))
    allocate(rwork2d(size(work,1)*size(work,2)))
    work    = zero
    rwork2d = zero

!   slmsk
    call nemsio_readrecv(gfile, 'land', 'sfc', 1, rwork2d, iret=iret)
    if (iret /= 0) call error_msg(trim(my_name),trim(filename),'land','read',istop,iret)
    work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
    allocate(outtmp(latb+2,lonb))
    call tran_gfssfc(work,outtmp,lonb,latb)
    do j=1,lonb
       do i=1,latb+2
          isli_anl(i,j) = nint(outtmp(i,j))
       end do
    end do
    deallocate(outtmp)

!   Deallocate local work arrays
    deallocate(work,rwork2d)

    call nemsio_close(gfile,iret=iret)
    if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret)
!
!   Print date/time stamp
    if ( mype == 0 ) write(6, &
       '(a,'': read_sfc_anl_ ,nlon,nlat= '',2i6,'',hour= '',f4.1,'',idate= '',4i5)') &
       trim(my_name),lonb,latb,fhour,odate
  end subroutine read_sfc_anl_

  subroutine read_nemssfc_anl_(iope,isli_anl)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_nemssfc_anl     read nems surface guess file with analysis resolution
!
!   prgmmr: xuli          org: np23                date: 2016-08-18
!
! abstract: read nems surface file at analysis grids
!
! program history log:
!
!   input argument list:
!     iope        - mpi task handling i/o
!
!   output argument list:
!     isli      - sea/land/ice mask
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use kinds, only: r_kind,i_kind,r_single
    use gridmod, only: nlat,nlon
    use mpimod, only: mpi_itype,mpi_comm_world,mype
    implicit none

!   Declare passed variables
    integer(i_kind),                               intent(in   ) :: iope
    integer(i_kind), dimension(nlat,nlon),         intent(  out) :: isli_anl


!   Declare local variables
    integer(i_kind):: iret,npts

!-----------------------------------------------------------------------------
!   Read surface file on processor iope
    if(mype == iope)then
       call read_sfc_anl_(isli_anl)
       write(*,*) 'read_sfc nemsio'
    end if

!   Load onto all processors
    npts=nlat*nlon
    call mpi_bcast(isli_anl,npts,mpi_itype,iope,mpi_comm_world,iret)

  end subroutine read_nemssfc_anl_

  subroutine read_nst_ (tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_nst_     read nems nst surface guess file (quadratic
!                                 Gaussin grids) without scattering to tasks
!   prgmmr: Huang            org: np23                date: 2011-11-01
!
! abstract: read nems surface NST file
!
! program history log:
!   2011-11-01  Huang    Initial version based on sub read_gfsnst
!   2016-03-13  Li       Modify for more effective I/O
!
!   input argument list:
!
!   output argument list:
!   tref     (:,:)                ! oceanic foundation temperature
!   dt_cool  (:,:)                ! sub-layer cooling amount at sub-skin layer
!   z_c      (:,:)                ! depth of sub-layer cooling layer
!   dt_warm  (:,:)                ! diurnal warming amount at sea surface (skin layer)
!   z_w      (:,:)                ! depth of diurnal warming layer
!   c_0      (:,:)                ! coefficient to calculate d(Tz)/d(tr) 
!   c_d      (:,:)                ! coefficient to calculate d(Tz)/d(tr) 
!   w_0      (:,:)                ! coefficient to calculate d(Tz)/d(tr) 
!   w_d      (:,:)                ! coefficient to calculate d(Tz)/d(tr) 
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use kinds, only: r_kind,i_kind,r_single
    use mpimod, only: mype
    use gridmod, only: nlat_sfc,nlon_sfc
    use constants, only: zero,two
    use guess_grids, only: nfldnst,ifilenst
    use nemsio_module, only:  nemsio_init,nemsio_open,nemsio_close
    use nemsio_module, only:  nemsio_gfile,nemsio_getfilehead,nemsio_readrecv
    implicit none

!   Declare passed variables
    real(r_single) , dimension(nlat_sfc,nlon_sfc,nfldnst), intent(  out) ::  &
         tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d
!   Declare local parameters
    integer(i_kind),parameter    :: n_nst=9
    integer(i_kind),dimension(7) :: idate
    integer(i_kind),dimension(4) :: odate

!   Declare local variables
    character(len=6)   :: filename
    character(len=120) :: my_name = 'READ_NEMSNST'
    character(len=1)   :: null = ' '
    integer(i_kind) :: i,j,it,latb,lonb
    integer(i_kind) :: iret, nframe
    integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd
    integer(i_kind) :: istop = 103
    real(r_single) :: fhour
    real(r_single), dimension(nlat_sfc,nlon_sfc,nfldnst) :: xt
    real(r_single), allocatable, dimension(:)   :: rwork2d
    real(r_single), allocatable, dimension(:,:) :: work

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!  Define read variable property   !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
    type(nemsio_gfile) :: gfile
!-----------------------------------------------------------------------------

    call nemsio_init(iret=iret)
    if (iret /= 0) call error_msg(trim(my_name),null,null,'init',istop,iret)


    do it=1,nfldnst
! read a nst file on the task
       write(filename,200)ifilenst(it)
200    format('nstf',i2.2)
       call nemsio_open(gfile,trim(filename),'READ',iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret)

       call nemsio_getfilehead(gfile, idate=idate, iret=iret, nframe=nframe,   &
          nfhour=nfhour, nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd, &
          dimx=lonb, dimy=latb )

       if( nframe /= 0 ) then
          if ( mype == 0 ) &
          write(6,*)trim(my_name),': ***ERROR***  nframe /= 0 for global model read, nframe = ', nframe
          call stop2(istop)
       end if

       fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
       odate(1) = idate(4)  !hour
       odate(2) = idate(2)  !month
       odate(3) = idate(3)  !day
       odate(4) = idate(1)  !year

       if ( (latb /= nlat_sfc-2) .or. (lonb /= nlon_sfc) ) then
          if ( mype == 0 ) &
             write(6,'(a,'': inconsistent spatial dimension nlon,nlatm2 = '',2(i4,tr1),''-vs- sfc file lonb,latb = '',i4)') &
             trim(my_name),nlon_sfc,nlat_sfc-2,lonb,latb
          call stop2(80)
       endif
!
!      Load surface fields into local work array
!      Follow NEMS/GFS sfcf read order
!
       allocate(work(lonb,latb))
       allocate(rwork2d(size(work,1)*size(work,2)))
       work    = zero
       rwork2d = zero

!      Tref
       call nemsio_readrecv(gfile, 'tref', 'sfc', 1, rwork2d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tref','read',istop,iret)
       work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
       call tran_gfssfc(work,tref(1,1,it),lonb,latb)

!      dt_cool
       call nemsio_readrecv(gfile, 'dtcool','sfc', 1, rwork2d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dt_cool','read',istop,iret)
       work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
       call tran_gfssfc(work,dt_cool(1,1,it),lonb,latb)

!      z_c
       call nemsio_readrecv(gfile, 'zc',  'sfc', 1, rwork2d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'z_c','read',istop,iret)
       work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
       call tran_gfssfc(work,z_c(1,1,it),lonb,latb)

!      xt
       call nemsio_readrecv(gfile, 'xt',   'sfc', 1, rwork2d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'xt','read',istop,iret)
       work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
       call tran_gfssfc(work,xt(1,1,it),lonb,latb)

!      xz
       call nemsio_readrecv(gfile, 'xz',   'sfc', 1, rwork2d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'xz','read',istop,iret)
       work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
       call tran_gfssfc(work,z_w(1,1,it),lonb,latb)
!
!      c_0
       call nemsio_readrecv(gfile, 'c0',  'sfc', 1, rwork2d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'c_0','read',istop,iret)
       work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
       call tran_gfssfc(work,c_0(1,1,it),lonb,latb)

!      c_d
       call nemsio_readrecv(gfile, 'cd',  'sfc', 1, rwork2d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'c_d','read',istop,iret)
       work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
       call tran_gfssfc(work,c_d(1,1,it),lonb,latb)

!      w_0
       call nemsio_readrecv(gfile, 'w0',  'sfc', 1, rwork2d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'w_0','read',istop,iret)
       work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
       call tran_gfssfc(work,w_0(1,1,it),lonb,latb)

!      w_d
       call nemsio_readrecv(gfile, 'wd',  'sfc', 1, rwork2d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'w_d','read',istop,iret)
       work(:,:)=reshape(rwork2d(:),(/size(work,1),size(work,2)/))
       call tran_gfssfc(work,w_d(1,1,it),lonb,latb)

!
!      Get diurnal warming amout at z=0
!
       do j = 1,nlon_sfc
          do i = 1,nlat_sfc
             if (z_w(i,j,it) > zero) then
                dt_warm(i,j,it) = two*xt(i,j,it)/z_w(i,j,it)
             end if
          end do
       end do

!      Deallocate local work arrays
       deallocate(work,rwork2d)

       call nemsio_close(gfile,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret)
!   End of loop over time levels
    end do
  end subroutine read_nst_


  subroutine read_nemsnst_ (iope,tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_nems_nst
!   prgmmr: li          org: np23                date: 2016-03-13
!
! abstract: read nems nst fields from a specific task and then broadcast to others
!
!   input argument list:
!     iope     - mpi task handling i/o
!
!   output argument list:
!   tref     (:,:)                        ! oceanic foundation temperature
!   dt_cool  (:,:)                        ! sub-layer cooling amount at sub-skin layer
!   z_c      (:,:)                        ! depth of sub-layer cooling layer
!   dt_warm  (:,:)                        ! diurnal warming amount at sea surface
!   z_w      (:,:)                        ! depth of diurnal warming layer
!   c_0      (:,:)                        ! coefficient to calculate d(Tz)/d(tr) in dimensionless
!   c_d      (:,:)                        ! coefficient to calculate d(Tz)/d(tr) in m^-1
!   w_0      (:,:)                        ! coefficient to calculate d(Tz)/d(tr) in dimensionless
!   w_d      (:,:)                        ! coefficient to calculate d(Tz)/d(tr) in m^-1
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use kinds, only: r_kind,i_kind,r_single
    use gridmod, only: nlat_sfc,nlon_sfc
    use guess_grids, only: nfldnst
    use mpimod, only: mpi_itype,mpi_rtype4,mpi_comm_world
    use mpimod, only: mype
    use constants, only: zero
    implicit none

!   Declare passed variables
    integer(i_kind),                                      intent(in   ) :: iope
    real(r_single), dimension(nlat_sfc,nlon_sfc,nfldnst), intent(  out) :: &
                    tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d

!   Declare local variables
    integer(i_kind):: iret,npts,nptsall

!-----------------------------------------------------------------------------
!   Read nst file on processor iope
    if(mype == iope)then
       write(*,*) 'read_nst nemsio'
       call read_nst_(tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)
    end if

!   Load onto all processors

    npts=nlat_sfc*nlon_sfc
    nptsall=npts*nfldnst

    call mpi_bcast(tref,    nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    call mpi_bcast(dt_cool, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    call mpi_bcast(z_c,     nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    call mpi_bcast(dt_warm, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    call mpi_bcast(z_w,     nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    call mpi_bcast(c_0,     nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    call mpi_bcast(c_d,     nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    call mpi_bcast(w_0,     nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
    call mpi_bcast(w_d,     nptsall,mpi_rtype4,iope,mpi_comm_world,iret)

  end subroutine read_nemsnst_


  subroutine write_atm_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin)

!$$$  subprogram documentation block
!                .      .    .
! subprogram:    write_nemsatm --- Gather, transform, and write out
!
!   prgmmr: Huang            org: np23                date: 2010-02-22
!
! abstract: This routine gathers fields needed for the GSI analysis
!           file from subdomains and then transforms the fields from
!           analysis grid to model guess grid, then written to an
!           atmospheric analysis file.
!
! program history log:
!   2010-02-22  Huang    Initial version.  Based on write_gfsatm
!   2011-02-14  Huang    Re-arrange the write sequence to be same as model
!                        read/rite sequence.
!   2013-10-25  todling  reposition load_grid to commvars
!   2016-07-28  mahajan  update with bundling ability
!
!   input argument list:
!     filename  - file to open and write to
!     mype_out  - mpi task to write output file
!    gfs_bundle - bundle containing fields on subdomains
!     ibin      - time bin
!
!   output argument list:
!
! attributes:
!   language: f90
!   machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP
!
!$$$ end documentation block

! !USES:
    use kinds, only: r_kind,i_kind

    use constants, only: r1000,fv,one,zero,qcmin,r0_05,t0c

    use mpimod, only: mpi_rtype
    use mpimod, only: mpi_comm_world
    use mpimod, only: ierror
    use mpimod, only: mype

    use guess_grids, only: ifilesig
    use guess_grids, only: ges_prsl,ges_prsi
    use guess_grids, only: load_geop_hgt,geop_hgti,ges_geopi

    use gridmod, only: ntracer
    use gridmod, only: ncloud
    use gridmod, only: strip,jcap_b,bk5

    use general_commvars_mod, only: load_grid,fill2_ns,filluv2_ns
    use general_specmod, only: spec_vars

    use obsmod, only: iadate

    use nemsio_module, only: nemsio_gfile,nemsio_open,nemsio_init,&
         nemsio_getfilehead,nemsio_close,nemsio_writerecv,nemsio_readrecv
    use gsi_4dvar, only: ibdate,nhr_obsbin,lwrite4danl
    use general_sub2grid_mod, only: sub2grid_info
    use egrid2agrid_mod,only: g_egrid2agrid,g_create_egrid2agrid,egrid2agrid_parm,destroy_egrid2agrid
    use constants, only: two,pi,half,deg2rad
    use gsi_bundlemod, only: gsi_bundle
    use gsi_bundlemod, only: gsi_bundlegetpointer
    use control_vectors, only: imp_physics,lupp
    use cloud_efr_mod, only: cloud_calc_gfs

    implicit none

! !INPUT PARAMETERS:

    type(sub2grid_info), intent(in) :: grd
    type(spec_vars),     intent(in) :: sp_a
    character(len=24),   intent(in) :: filename  ! file to open and write to
    integer(i_kind),     intent(in) :: mype_out  ! mpi task to write output file
    type(gsi_bundle),    intent(in) :: gfs_bundle
    integer(i_kind),     intent(in) :: ibin      ! time bin

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

    real(r_kind),parameter:: r0_001 = 0.001_r_kind
    character(6):: fname_ges
    character(len=120) :: my_name = 'WRITE_NEMSATM'
    character(len=1)   :: null = ' '
    integer(i_kind),dimension(7):: idate, jdate
    integer(i_kind),dimension(4):: odate
    integer(i_kind) :: k, mm1, nlatm2, nord_int, i, j, kk
    integer(i_kind) :: iret, lonb, latb, levs, istatus
    integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd
    integer(i_kind) :: istop = 104
    integer(i_kind),dimension(5):: mydate
    integer(i_kind),dimension(8) :: ida,jda
    real(r_kind),dimension(5)    :: fha
    real(r_kind)    :: fhour

    real(r_kind),pointer,dimension(:,:) :: sub_ps
    real(r_kind),pointer,dimension(:,:,:) :: sub_u,sub_v,sub_tv
    real(r_kind),pointer,dimension(:,:,:) :: sub_q,sub_oz,sub_cwmr

    real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig) :: sub_dzb,sub_dza
    real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig) :: sub_prsl
    real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig+1) :: sub_prsi

    real(r_kind),dimension(grd%lat1*grd%lon1)     :: psm
    real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig):: sub_dp
    real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: tvsm,prslm, usm, vsm
    real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: dpsm, qsm, ozsm
    real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: cwsm, dzsm
    real(r_kind),dimension(max(grd%iglobal,grd%itotsub))     :: work1,work2
    real(r_kind),dimension(grd%nlon,grd%nlat-2):: grid
    real(r_kind),allocatable,dimension(:) :: rwork1d,rwork1d1,rlats,rlons,clons,slons
    real(4),allocatable,dimension(:) :: r4lats,r4lons
    real(r_kind),allocatable,dimension(:,:) :: grid_b,grid_b2
    real(r_kind),allocatable,dimension(:,:,:) :: grid_c, grid3, grid_c2, grid3b

    type(nemsio_gfile) :: gfile,gfileo
    logical diff_res,eqspace
    logical,dimension(1) :: vector
    type(egrid2agrid_parm) :: p_low,p_high

!*************************************************************************
!   Initialize local variables
    mm1=mype+1
    nlatm2=grd%nlat-2
    diff_res=.false.

    istatus=0
    call gsi_bundlegetpointer(gfs_bundle,'ps', sub_ps,  iret); istatus=istatus+iret
    call gsi_bundlegetpointer(gfs_bundle,'u',  sub_u,   iret); istatus=istatus+iret
    call gsi_bundlegetpointer(gfs_bundle,'v',  sub_v,   iret); istatus=istatus+iret
    call gsi_bundlegetpointer(gfs_bundle,'tv', sub_tv,  iret); istatus=istatus+iret
    call gsi_bundlegetpointer(gfs_bundle,'q',  sub_q,   iret); istatus=istatus+iret
    call gsi_bundlegetpointer(gfs_bundle,'oz', sub_oz,  iret); istatus=istatus+iret
    call gsi_bundlegetpointer(gfs_bundle,'cw', sub_cwmr,iret); istatus=istatus+iret
    if ( istatus /= 0 ) then
       if ( mype == 0 ) then
         write(6,*) 'write_atm_: ERROR'
         write(6,*) 'Missing some of the required fields'
         write(6,*) 'Aborting ... '
      endif
      call stop2(999)
    endif

    if ( sp_a%jcap /= jcap_b ) then
        if ( mype == 0 ) write(6, &
            '('' dual resolution for nems sp_a%jcap,jcap_b = '',2i6)') &
            sp_a%jcap,jcap_b
        diff_res = .true.
    endif


    ! Single task writes analysis data to analysis file
    if ( mype == mype_out ) then
       write(fname_ges,'(''sigf'',i2.2)') ifilesig(ibin)

       ! Read header information from first guess file.
       call nemsio_init(iret)
       if ( iret /= 0 ) call error_msg(trim(my_name),null,null,'init',istop,iret)

       call nemsio_open(gfile,trim(fname_ges),'read',iret)
       if ( iret /= 0 ) call error_msg(trim(my_name),trim(fname_ges),null,'open',istop,iret)

       call nemsio_getfilehead(gfile, iret=iret, nfhour=nfhour, &
            nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd, &
            idate=idate, dimx=lonb, dimy=latb, dimz=levs)
       if ( iret /= 0 ) then
          write(6,*) trim(my_name),': problem with nemsio_getfilehead, Status = ',iret
          call stop2(103)
       endif
       if ( levs /= grd%nsig ) then
          write(6,*) trim(my_name),': problem in data dimension background levs = ',levs,' nsig = ',grd%nsig
          call stop2(103)
       endif

       ! copy input header info to output header info
       gfileo=gfile

       ! Update header information (with ibdate) and write it to analysis file (w/ _open statement).
       mydate=ibdate
       fha(:)=zero ; ida=0; jda=0
       fha(2)=real(nhr_obsbin*(ibin-1))  ! relative time interval in hours
       ida(1)=mydate(1) ! year
       ida(2)=mydate(2) ! month
       ida(3)=mydate(3) ! day
       ida(4)=0         ! time zone
       ida(5)=mydate(4) ! hour

       ! Move date-time forward by nhr_assimilation hours
       call w3movdat(fha,ida,jda)

       jdate(1) = jda(1)     ! analysis year
       jdate(2) = jda(2)     ! analysis month
       jdate(3) = jda(3)     ! analysis day
       jdate(4) = jda(5)     ! analysis hour
       jdate(5) = iadate(5)  ! analysis minute
       jdate(6) = 0          ! analysis scaled seconds
       jdate(7) = idate(7)   ! analysis seconds multiplier

       nfhour   =0       !  new forecast hour, zero at analysis time
       nfminute =0
       nfsecondn=0
       nfsecondd=100      ! default for denominator

       fhour = zero
       odate(1) = jdate(4)  !hour
       odate(2) = jdate(2)  !month
       odate(3) = jdate(3)  !day
       odate(4) = jdate(1)  !year

       ! open new output file with new header gfileo with "write" access.
       ! Use this call to update header as well

       call nemsio_open(gfileo,trim(filename),'write',iret=iret, &
          idate=jdate, nfhour=nfhour, nfminute=nfminute, &
          nfsecondn=nfsecondn, nfsecondd=nfsecondd)
       if ( iret /= 0 ) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret)

       ! Allocate structure arrays to hold data
       allocate(rwork1d(latb*lonb),rwork1d1(latb*lonb))
       if (imp_physics == 11) allocate(grid3b(grd%nlat,grd%nlon,1))
       if ( diff_res .or. imp_physics == 11 .or. lupp) then
          allocate( grid_b(lonb,latb),grid_c(latb+2,lonb,1),grid3(grd%nlat,grd%nlon,1))
          allocate( grid_b2(lonb,latb),grid_c2(latb+2,lonb,1))
          allocate( rlats(latb+2),rlons(lonb),clons(lonb),slons(lonb),r4lats(lonb*latb),r4lons(lonb*latb))
          call nemsio_getfilehead(gfile,lat=r4lats,iret=iret)
          call nemsio_getfilehead(gfile,lon=r4lons,iret=iret)
          do j=1,latb
            rlats(latb+2-j)=deg2rad*r4lats(lonb/2+(j-1)*lonb)
          enddo
          rlats(1)=-half*pi
          rlats(latb+2)=half*pi
          do j=1,lonb
            rlons(j)=deg2rad*r4lons(j)
          enddo
          do j=1,lonb
             clons(j)=cos(rlons(j))
             slons(j)=sin(rlons(j))
          enddo

          nord_int=4
          eqspace=.false.
          call g_create_egrid2agrid(grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons, &
                                latb+2,rlats,lonb,rlons,&
                                nord_int,p_low,.false.,eqspace=eqspace)
          call g_create_egrid2agrid(latb+2,rlats,lonb,rlons, &
                                grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons,&
                                nord_int,p_high,.false.,eqspace=eqspace)

          deallocate(rlats,rlons,r4lats,r4lons)
       endif ! if ( diff_res )

       !   Terrain
       !   Write out input file surface height

       call nemsio_readrecv(gfile,'hgt', 'sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'hgt','writeread',istop,iret)
       call nemsio_writerecv(gfileo,'hgt','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'hgt','write',istop,iret)
    endif ! if ( mype == mype_out )

    sub_prsl = ges_prsl(:,:,:,ibin)
    sub_prsi = ges_prsi(:,:,:,ibin)

    do k=1,grd%nsig
       sub_dp(:,:,k) = sub_prsi(:,:,k) - sub_prsi(:,:,k+1)
    end do

    ! Calculate delz increment for UPP
    if (lupp) then
       if ((.not. lwrite4danl) .or. ibin == 1) ges_geopi = geop_hgti
       do k=1,grd%nsig
          sub_dzb(:,:,k) = ges_geopi(:,:,k+1,ibin) - ges_geopi(:,:,k,ibin)
       enddo

       if ((.not. lwrite4danl) .or. ibin == 1) call load_geop_hgt
       do k=1,grd%nsig
          sub_dza(:,:,k) = geop_hgti(:,:,k+1,ibin) - geop_hgti(:,:,k,ibin)
       enddo

       sub_dza = sub_dza - sub_dzb !sub_dza is increment
    endif
    
    ! Strip off boundary points from subdomains
    call strip(sub_ps  ,psm)
    call strip(sub_tv  ,tvsm  ,grd%nsig)
    call strip(sub_q   ,qsm   ,grd%nsig)
    call strip(sub_oz  ,ozsm  ,grd%nsig)
    call strip(sub_cwmr,cwsm  ,grd%nsig)
    call strip(sub_dp  ,dpsm  ,grd%nsig)
    call strip(sub_prsl,prslm ,grd%nsig)
    call strip(sub_u   ,usm   ,grd%nsig)
    call strip(sub_v   ,vsm   ,grd%nsig)
    if (lupp) call strip(sub_dza ,dzsm  ,grd%nsig)

    ! Thermodynamic variable
    ! The GSI analysis variable is virtual temperature (Tv).   For NEMSIO
    ! output we need the sensible temperature.

    ! Convert Tv to T
    tvsm = tvsm/(one+fv*qsm)

    ! Generate and write analysis fields

    ! Surface pressure.
    call mpi_gatherv(psm,grd%ijn(mm1),mpi_rtype,&
         work1,grd%ijn,grd%displs_g,mpi_rtype,&
         mype_out,mpi_comm_world,ierror)
    if (mype==mype_out) then
       if(diff_res .or. lupp)then
          call nemsio_readrecv(gfile,'pres','sfc',1,rwork1d,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'pres','read',istop,iret)
          rwork1d1 = r0_001*rwork1d
          grid_b=reshape(rwork1d1,(/size(grid_b,1),size(grid_b,2)/))
          vector(1)=.false.
          call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
          call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector)
          do kk=1,grd%iglobal
             i=grd%ltosi(kk)
             j=grd%ltosj(kk)
             grid3(i,j,1)=work1(kk)-grid3(i,j,1)
             if (lupp) work1(kk)=grid3(i,j,1)
          end do
          if (lupp) then
             do k=1,grd%nsig
                do kk=1,grd%iglobal
                   i=grd%ltosi(kk)
                   j=grd%ltosj(kk)
                   grid3(i,j,1)=work1(kk)*(bk5(k)-bk5(k+1))
                enddo
                call g_egrid2agrid(p_high,grid3,grid_c2,1,1,vector)
                call nemsio_readrecv(gfile,'dpres','mid layer',k,rwork1d,iret=iret)
                if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dpres','read',istop,iret)
                grid_b2=reshape(rwork1d,(/size(grid_b2,1),size(grid_b2,2)/))
                do j=1,latb
                   do i=1,lonb
                      grid_b2(i,j)=grid_b2(i,j)+r1000*(grid_c2(latb-j+2,i,1))
                   enddo
                enddo
                rwork1d = reshape(grid_b2,(/size(rwork1d)/))
                call nemsio_writerecv(gfileo,'dpres','mid layer',k,rwork1d,iret=iret)
                if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dpres','write',istop,iret)
             enddo
             do kk=1,grd%iglobal
                i=grd%ltosi(kk)
                j=grd%ltosj(kk)
                grid3(i,j,1)=work1(kk)
             enddo
          endif
          call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
          do j=1,latb
             do i=1,lonb
                grid_b(i,j)=r1000*(grid_b(i,j)+grid_c(latb-j+2,i,1))
             end do
          end do
          rwork1d = reshape(grid_b,(/size(rwork1d)/))
       else
          call load_grid(work1,grid)
          grid = grid*r1000
          rwork1d = reshape(grid,(/size(rwork1d)/))
       end if
       call nemsio_writerecv(gfileo,'pres','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'psfc','write',istop,iret)
    endif

!   u, v
    do k=1,grd%nsig
       call mpi_gatherv(usm(1,k),grd%ijn(mm1),mpi_rtype,&
            work1,grd%ijn,grd%displs_g,mpi_rtype,&
            mype_out,mpi_comm_world,ierror)
       call mpi_gatherv(vsm(1,k),grd%ijn(mm1),mpi_rtype,&
            work2,grd%ijn,grd%displs_g,mpi_rtype,&
            mype_out,mpi_comm_world,ierror)
       if (mype==mype_out) then
          if(diff_res)then
             call nemsio_readrecv(gfile,'ugrd','mid layer',k,rwork1d,iret=iret)
             call nemsio_readrecv(gfile,'vgrd','mid layer',k,rwork1d1,iret=iret)
             grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/))
             grid_b2=reshape(rwork1d1,(/size(grid_b,1),size(grid_b,2)/))
             vector(1)=.true.
             call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons)
             call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector)
             do kk=1,grd%iglobal
                i=grd%ltosi(kk)
                j=grd%ltosj(kk)
                grid3(i,j,1)=work1(kk)-grid3(i,j,1)
             end do
             call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
             do j=1,latb
                do i=1,lonb
                   grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1)
                end do
             end do
             call g_egrid2agrid(p_low,grid_c2,grid3,1,1,vector)
             do kk=1,grd%iglobal
                i=grd%ltosi(kk)
                j=grd%ltosj(kk)
                grid3(i,j,1)=work2(kk)-grid3(i,j,1)
             end do
             call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
             do j=1,latb
                do i=1,lonb
                   grid_b2(i,j)=grid_b2(i,j)+grid_c(latb-j+2,i,1)
                end do
             end do
             rwork1d = reshape(grid_b,(/size(rwork1d)/))
             rwork1d1 = reshape(grid_b2,(/size(rwork1d1)/))

          else
             call load_grid(work1,grid)
             rwork1d = reshape(grid,(/size(rwork1d)/))
             call load_grid(work2,grid)
             rwork1d1 = reshape(grid,(/size(rwork1d1)/))
          end if

          ! Zonal wind
          call nemsio_writerecv(gfileo,'ugrd','mid layer',k,rwork1d,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ugrd','write',istop,iret)
          ! Meridional wind
          call nemsio_writerecv(gfileo,'vgrd','mid layer',k,rwork1d1,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'vgrd','write',istop,iret)
       endif
    end do

!   Thermodynamic variable
    do k=1,grd%nsig
       call mpi_gatherv(tvsm(1,k),grd%ijn(mm1),mpi_rtype,&
            work1,grd%ijn,grd%displs_g,mpi_rtype,&
            mype_out,mpi_comm_world,ierror)
       if (mype == mype_out) then
          if(diff_res)then
             call nemsio_readrecv(gfile,'tmp','mid layer',k,rwork1d,iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','read',istop,iret)
             grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/))
             vector(1)=.false.
             call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
             call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector)
             do kk=1,grd%iglobal
                i=grd%ltosi(kk)
                j=grd%ltosj(kk)
                grid3(i,j,1)=work1(kk)-grid3(i,j,1)
             end do
             call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
             do j=1,latb
                do i=1,lonb
                   grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1)
                end do
             end do
             rwork1d = reshape(grid_b,(/size(rwork1d)/))
          else
             call load_grid(work1,grid)
             rwork1d = reshape(grid,(/size(rwork1d)/))
          end if
          call nemsio_writerecv(gfileo,'tmp','mid layer',k,rwork1d,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','write',istop,iret)
       endif
    end do

!   Specific humidity
    do k=1,grd%nsig
       call mpi_gatherv(qsm(1,k),grd%ijn(mm1),mpi_rtype,&
            work1,grd%ijn,grd%displs_g,mpi_rtype,&
            mype_out,mpi_comm_world,ierror)
       if (mype == mype_out) then
          if(diff_res)then
             call nemsio_readrecv(gfile,'spfh','mid layer',k,rwork1d,iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'spfh','read',istop,iret)
             grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/))
             vector(1)=.false.
             call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
             call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector)
             do kk=1,grd%iglobal
                i=grd%ltosi(kk)
                j=grd%ltosj(kk)
                grid3(i,j,1)=work1(kk)-grid3(i,j,1)
             end do
             call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
             do j=1,latb
                do i=1,lonb
                   grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1)
                end do
             end do
             rwork1d = reshape(grid_b,(/size(rwork1d)/))
          else
             call load_grid(work1,grid)
             rwork1d = reshape(grid,(/size(rwork1d)/))
          end if
          call nemsio_writerecv(gfileo,'spfh','mid layer',k,rwork1d,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'spfh','write',istop,iret)
       endif
    end do

!   Ozone
    do k=1,grd%nsig
       call mpi_gatherv(ozsm(1,k),grd%ijn(mm1),mpi_rtype,&
            work1,grd%ijn,grd%displs_g,mpi_rtype,&
            mype_out,mpi_comm_world,ierror)
       if (mype == mype_out) then
          if(diff_res)then
             call nemsio_readrecv(gfile,'o3mr','mid layer',k,rwork1d,iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'o3mr','read',istop,iret)
             grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/))
             vector(1)=.false.
             call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
             call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector)
             do kk=1,grd%iglobal
                i=grd%ltosi(kk)
                j=grd%ltosj(kk)
                grid3(i,j,1)=work1(kk)-grid3(i,j,1)
             end do
             call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
             do j=1,latb
                do i=1,lonb
                   grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1)
                end do
             end do
             rwork1d = reshape(grid_b,(/size(rwork1d)/))
          else
             call load_grid(work1,grid)
             rwork1d = reshape(grid,(/size(rwork1d)/))
          end if
          call nemsio_writerecv(gfileo,'o3mr','mid layer',k,rwork1d,iret=iret)
          if (iret /= 0) call error_msg(trim(my_name),trim(filename),'o3mr','write',istop,iret)
       endif
    end do

!   Cloud condensate mixing ratio
    if (ntracer>2 .or. ncloud>=1) then

       do k=1,grd%nsig
          call mpi_gatherv(cwsm(1,k),grd%ijn(mm1),mpi_rtype,&
               work1,grd%ijn,grd%displs_g,mpi_rtype,&
               mype_out,mpi_comm_world,ierror)
          if (imp_physics == 11) then
             call mpi_gatherv(tvsm(1,k),grd%ijn(mm1),mpi_rtype,&
                  work2,grd%ijn,grd%displs_g,mpi_rtype,&
                  mype_out,mpi_comm_world,ierror)
          endif
          if (mype == mype_out) then
             if(diff_res .or. imp_physics == 11)then
                call nemsio_readrecv(gfile,'clwmr','mid layer',k,rwork1d,iret=iret)
                if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','read',istop,iret)
                grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/))
                if (imp_physics == 11) then
                   call nemsio_readrecv(gfile,'icmr','mid layer',k,rwork1d1,iret=iret)
                   if (iret /= 0) call error_msg(trim(my_name),trim(filename),'icmr','read',istop,iret)
                   grid_b2=reshape(rwork1d1,(/size(grid_b2,1),size(grid_b2,2)/))
                   grid_b = grid_b + grid_b2
                endif
                vector(1)=.false.
                call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
                call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector)
                do kk=1,grd%iglobal
                   i=grd%ltosi(kk)
                   j=grd%ltosj(kk)
                   grid3(i,j,1)=work1(kk)-max(grid3(i,j,1),qcmin)
                   if (imp_physics == 11) then
                      work2(kk) = -r0_05*(work2(kk) - t0c)
                      work2(kk) = max(zero,work2(kk))
                      work2(kk) = min(one,work2(kk))
                      grid3b(i,j,1)=grid3(i,j,1)
                      grid3(i,j,1)=grid3b(i,j,1)*(one - work2(kk))
                   endif
                end do
                call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
                if (imp_physics == 11) grid_b = grid_b - grid_b2
                do j=1,latb
                   do i=1,lonb
                      grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1)
                   end do
                end do
                rwork1d = reshape(grid_b,(/size(rwork1d)/))
                if (imp_physics == 11) then
                   do kk=1,grd%iglobal
                      i=grd%ltosi(kk)
                      j=grd%ltosj(kk)
                      grid3(i,j,1)=grid3b(i,j,1)*work2(kk)
                   end do
                   call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
                   do j=1,latb
                      do i=1,lonb
                         grid_b2(i,j)=grid_b2(i,j)+grid_c(latb-j+2,i,1)
                      end do
                   end do
                   rwork1d1 = reshape(grid_b2,(/size(rwork1d1)/))
                endif
             else
                call load_grid(work1,grid)
                rwork1d = reshape(grid,(/size(rwork1d)/))
             endif
             call nemsio_writerecv(gfileo,'clwmr','mid layer',k,rwork1d,iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','write',istop,iret)
             if (imp_physics == 11) then
                call nemsio_writerecv(gfileo,'icmr','mid layer',k,rwork1d1,iret=iret)
                if (iret /= 0) call error_msg(trim(my_name),trim(filename),'icmr','write',istop,iret)

                if (lupp) then
                   call nemsio_readrecv(gfile,'rwmr','mid layer',k,rwork1d,iret=iret)
                   if (iret /= 0) call error_msg(trim(my_name),trim(filename),'rwmr','read',istop,iret)
                   call nemsio_writerecv(gfileo,'rwmr','mid layer',k,rwork1d,iret=iret)
                   if (iret /= 0) call error_msg(trim(my_name),trim(filename),'rwmr','write',istop,iret)

                   call nemsio_readrecv(gfile,'snmr','mid layer',k,rwork1d,iret=iret)
                   if (iret /= 0) call error_msg(trim(my_name),trim(filename),'snmr','read',istop,iret)
                   call nemsio_writerecv(gfileo,'snmr','mid layer',k,rwork1d,iret=iret)
                   if (iret /= 0) call error_msg(trim(my_name),trim(filename),'snmr','write',istop,iret)

                   call nemsio_readrecv(gfile,'grle','mid layer',k,rwork1d,iret=iret)
                   if (iret /= 0) call error_msg(trim(my_name),trim(filename),'grle','read',istop,iret)
                   call nemsio_writerecv(gfileo,'grle','mid layer',k,rwork1d,iret=iret)
                   if (iret /= 0) call error_msg(trim(my_name),trim(filename),'grle','write',istop,iret)
                endif

                call nemsio_readrecv(gfile,'cld_amt','mid layer',k,rwork1d,iret=iret)
                if (iret == 0) then
                   call nemsio_writerecv(gfileo,'cld_amt','mid layer',k,rwork1d,iret=iret)
                   if (iret /= 0) call error_msg(trim(my_name),trim(filename),'cld_amt','write',istop,iret)
                endif

             endif
          endif !mype == mype_out
       end do
    endif !ntracer

! Variables needed by the Unified Post Processor (dzdt, delz, delp)
    if (lupp) then
       if (mype == mype_out) then
          do k=1,grd%nsig
             call nemsio_readrecv(gfile,'dzdt','mid layer',k,rwork1d,iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dzdt','read',istop,iret)
             call nemsio_writerecv(gfileo,'dzdt','mid layer',k,rwork1d,iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dzdt','write',istop,iret)
          enddo
       endif
       do k=1,grd%nsig
          call mpi_gatherv(dzsm(1,k),grd%ijn(mm1),mpi_rtype,&
               work1,grd%ijn,grd%displs_g,mpi_rtype,&
               mype_out,mpi_comm_world,ierror)
          if (mype == mype_out) then
             call nemsio_readrecv(gfile,'delz','mid layer',k,rwork1d,iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'delz','read',istop,iret)
             if(diff_res)then
                grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/))
                do kk=1,grd%iglobal
                   i=grd%ltosi(kk)
                   j=grd%ltosj(kk)
                   grid3(i,j,1)=work1(kk)
                end do
                call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
                do j=1,latb
                   do i=1,lonb
                      grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1)
                   end do
                end do
                rwork1d = reshape(grid_b,(/size(rwork1d)/))
             else
                call load_grid(work1,grid)
                rwork1d = rwork1d + reshape(grid,(/size(rwork1d)/))
             end if
             call nemsio_writerecv(gfileo,'delz','mid layer',k,rwork1d,iret=iret)
             if (iret /= 0) call error_msg(trim(my_name),trim(filename),'delz','write',istop,iret)
          endif
       end do
    endif
    
!
! Deallocate local array
!
    if (mype==mype_out) then
       if (diff_res .or. lupp .or. imp_physics == 11) deallocate(grid_b,grid_b2,grid_c,grid_c2,grid3,clons,slons)
       if (imp_physics == 11) deallocate(grid3b)

       call nemsio_close(gfile,iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_ges),null,'close',istop,iret)

       call nemsio_close(gfileo,iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret)
!
! Deallocate local array
!
       deallocate(rwork1d,rwork1d1)
!
       write(6,'(a,'': atm anal written for lonb,latb,levs= '',3i6,'',valid hour= '',f4.1,'',idate= '',4i5)') &
          trim(my_name),lonb,latb,levs,fhour,odate
    endif

  end subroutine write_atm_

  subroutine write_sfc_ (filename,mype_sfc,dsfct)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    write_nemssfc --- Write surface analysis to file
!
!   prgmmr: Huang            org: np23                date: 2010-02-22
!
! abstract:     This routine writes the updated surface analysis.  At
!               this point (20101020) the only surface field update by
!               the gsi is the skin temperature.  The current (20101020)
!               GDAS setup does use the updated surface file.  Rather,
!               the output from surface cycle is used as the surface
!               analysis for subsequent NEMS/GFS runs.
!
!               The routine gathers surface fields from subdomains,
!               reformats the data records, and then writes each record
!               to the output file.
!
!               Since the gsi only update the skin temperature, all
!               other surface fields are simply read from the guess
!               surface file and written to the analysis file.
!
! program history log:
!   2010-02-22  Huang    Initial version.  Based on write_gfssfc
!   2011-04-01  Huang    change type of buffer2, grid2 from single to r_kind
!   2013-10-25  todling - reposition ltosi and others to commvars
!
!   input argument list:
!     filename  - file to open and write to
!     dsfct     - delta skin temperature
!     mype_sfc  - mpi task to write output file
!
!   output argument list:
!
! attributes:
!   language: f90
!   machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP
!
!$$$ end documentation block

! !USES:
    use kinds, only: r_kind,i_kind,r_single

    use mpimod, only: mpi_rtype
    use mpimod, only: mpi_comm_world
    use mpimod, only: ierror
    use mpimod, only: mype

    use gridmod, only: nlat,nlon
    use gridmod, only: lat1,lon1
    use gridmod, only: lat2,lon2
    use gridmod, only: iglobal
    use gridmod, only: ijn
    use gridmod, only: displs_g
    use gridmod, only: itotsub
    use gridmod, only: rlats,rlons,rlats_sfc,rlons_sfc

    use general_commvars_mod, only: ltosi,ltosj

    use obsmod, only: iadate

    use constants, only: zero

    use nemsio_module, only:  nemsio_init,nemsio_open,nemsio_close,nemsio_readrecv
    use nemsio_module, only:  nemsio_gfile,nemsio_getfilehead
    use nemsio_module, only:  nemsio_readrec, nemsio_writerec, nemsio_writerecv

    implicit none

! !INPUT PARAMETERS:
    character(24)                    ,intent(in   ) :: filename  ! file to open and write to

    real(r_kind),dimension(lat2,lon2),intent(in   ) :: dsfct   ! delta skin temperature

    integer(i_kind)                  ,intent(in   ) :: mype_sfc ! mpi task to write output file

! !OUTPUT PARAMETERS:

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

!   Declare local parameters
    character( 6),parameter:: fname_ges='sfcf06'
!   Declare local variables
    character(len=120) :: my_name = 'WRITE_NEMSSFC'
    character(len=1)   :: null = ' '
    integer(i_kind),dimension(7):: idate, jdate
    integer(i_kind),dimension(4):: odate
    integer(i_kind) :: i, j, ip1, jp1, ilat, ilon, jj, mm1
    integer(i_kind) :: nlatm2, n, nrec, lonb, latb, iret
    integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd
    integer(i_kind) :: istop = 105
    real(r_kind)    :: fhour

    real(r_kind),dimension(lat1,lon1):: sfcsub
    real(r_kind),dimension(nlon,nlat):: grid
    real(r_kind),dimension(max(iglobal,itotsub)):: sfcall
    real(r_kind),allocatable,dimension(:,:) :: tsea
    real(r_kind),allocatable,dimension(:)   :: rwork1d
    real(r_single),dimension(nlon,nlat):: buffer
    real(r_single),allocatable,dimension(:,:) :: buffer2,grid2

    type(nemsio_gfile) :: gfile, gfileo
!*****************************************************************************

!   Initialize local variables
    mm1=mype+1
    nlatm2=nlat-2

!   Gather skin temperature information from all tasks.
    do j=1,lon1
       jp1 = j+1
       do i=1,lat1
          ip1 = i+1
          sfcsub(i,j)=dsfct(ip1,jp1)
       end do
    end do
    call mpi_gatherv(sfcsub,ijn(mm1),mpi_rtype,&
         sfcall,ijn,displs_g,mpi_rtype,mype_sfc,&
         mpi_comm_world,ierror)

!   Only MPI task mype_sfc writes the surface file.
    if (mype==mype_sfc) then

!      Reorder updated skin temperature to output format
       do i=1,iglobal
          ilon=ltosj(i)
          ilat=ltosi(i)
          grid(ilon,ilat)=sfcall(i)
       end do
       do j=1,nlat
          jj=nlat-j+1
          do i=1,nlon
             buffer(i,j)=grid(i,jj)
          end do
       end do

!      Read surface guess file
       call nemsio_init(iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),null,null,'init',istop,iret)

       call nemsio_open(gfile,fname_ges,'read',iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_ges),null,'open',istop,iret)
!
       call nemsio_getfilehead(gfile, nrec=nrec, idate=idate, dimx=lonb, &
          dimy=latb, nfhour=nfhour, nfminute=nfminute, nfsecondn=nfsecondn, &
          nfsecondd=nfsecondd, iret=iret)
!
!      Replace header record date with analysis time from iadate
!
       jdate(1) = iadate(1)  ! analysis year
       jdate(2) = iadate(2)  ! analysis month
       jdate(3) = iadate(3)  ! analysis day
       jdate(4) = iadate(4)  ! analysis hour
       jdate(5) = iadate(5)  ! analysis minute
       jdate(5) = 0          ! analysis minute
       jdate(6) = 0          ! analysis scaled seconds
       jdate(7) = idate(7)   ! analysis seconds multiplier

       nfhour=0       !  new forecast hour, zero at analysis time
       nfminute=0
       nfsecondn=0
       nfsecondd=100      ! default for denominator

       fhour    = zero
       odate(1) = jdate(4)  !hour
       odate(2) = jdate(2)  !month
       odate(3) = jdate(3)  !day
       odate(4) = jdate(1)  !year
!
! Start to write output sfc file : filename
!      open new output file with new header gfileo with "write" access.
!      Use this call to update header as well
!
!
       gfileo=gfile      ! copy input header info to output header info
                         ! need to do this before nemsio_close(gfile)
       call nemsio_open(gfileo,filename,'write',iret=iret, idate=jdate, nfhour=nfhour,&
          nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd )
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret)
!
!      First copy entire data from fname_ges to filename, then do selective update
!
       allocate(rwork1d(lonb*latb))
       allocate(buffer2(lonb,latb))
       allocate(grid2(lonb,latb))
       allocate(tsea(lonb,latb))

       do n = 1, nrec
         call nemsio_readrec (gfile, n,rwork1d,iret=iret)
         if ( iret /= 0 ) write(6,*) 'readrec  nrec = ', n, '  Status = ', iret
         call nemsio_writerec(gfileo,n,rwork1d,iret=iret)
         if ( iret /= 0 ) write(6,*) 'writerec nrec = ', n, '  Status = ', iret
       end do
!
! Only sea surface temperature will be updated in the SFC files
!

       call nemsio_readrecv(gfile,'tmp','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_ges),'tmp','read',istop,iret)
       tsea=reshape(rwork1d,(/size(tsea,1),size(tsea,2)/))

       if ( (latb /= nlatm2) .or. (lonb /= nlon) ) then
          write(6,*)trim(my_name),':  different grid dimensions analysis', &
             ' vs sfc. interpolating sfc temperature nlon,nlat-2=',nlon,  &
             nlatm2,' -vs- sfc file lonb,latb=',lonb,latb
          call intrp22(buffer, rlons,rlats,nlon,nlat, &
                       buffer2,rlons_sfc,rlats_sfc,lonb,latb)
       else
          do j=1,latb
             do i=1,lonb
                buffer2(i,j)=buffer(i,j+1)
             end do
          end do
       endif

       grid2 = tsea + buffer2
       rwork1d = reshape( grid2,(/size(rwork1d)/) )

       deallocate(buffer2)

!      update tsea record
       call nemsio_writerecv(gfileo,'tmp','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','write',istop,iret)
       deallocate(rwork1d)

       call nemsio_close(gfile, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_ges),null,'close',istop,iret)

       call nemsio_close(gfileo,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret)

       write(6,'(a,'': sfc anal written for lonb,latb= '',2i6,'',valid hour= '',f4.1,'',idate= '',4i5)') &
          trim(my_name),lonb,latb,fhour,odate
    endif
  end subroutine write_sfc_

  subroutine write_sfc_nst_ (mype_so,dsfct)

!$$$  subprogram documentation block
!                .      .    .
! subprogram:    write_sfc_nst --- Write both sfc and nst surface analysis to file
!
!   prgmmr: Huang            org: np23                date: 2011-11-01
!
! abstract:     This routine writes the sfc & nst analysis files and is nst_gsi dependent.
!               Tr (foundation temperature), instead of skin temperature, is the analysis variable.
!               nst_gsi >  2: Tr analysis is on
!               nst_gsi <= 2: Tr analysis is off
!
!               The routine gathers Tr field from subdomains,
!               reformats the data records, and then writes each record
!               to the output files.
!
!               Since the gsi only update the Tr temperature, all
!               other fields in surface are simply read from the guess
!               files and written to the analysis file.
!
! program history log:
!   2011-11-01  Huang    initial version based on routine write_gfs_sfc_nst
!   2013-10-25  todling - reposition ltosi and others to commvars
!   2016-01-01  li      - update write_sfc_nst_ (nemsio) as for write_gfs_sfc_nst (sfcio)
!
!   input argument list:
!     dsfct     - delta skin temperature
!     mype_so   - mpi task to write output file
!
!   output argument list:
!
! attributes:
!   language: f90
!   machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP
!
!$$$ end documentation block

! !USES:
    use kinds, only: r_kind,i_kind,r_single

    use mpimod, only: mpi_rtype,mpi_itype
    use mpimod, only: mpi_comm_world
    use mpimod, only: ierror
    use mpimod, only: mype

    use gridmod, only: nlat,nlon
    use gridmod, only: lat1,lon1
    use gridmod, only: lat2,lon2
    use gridmod, only: nlat_sfc,nlon_sfc
    use gridmod, only: iglobal
    use gridmod, only: ijn
    use gridmod, only: displs_g
    use gridmod, only: itotsub

    use general_commvars_mod, only: ltosi,ltosj

    use obsmod, only: iadate

    use constants, only: zero,two,tfrozen,z_w_max
    use constants, only: zero_single

    use guess_grids, only: isli2
    use gsi_nstcouplermod, only: nst_gsi,zsea1,zsea2
    use gridmod, only: rlats,rlons,rlats_sfc,rlons_sfc

    use nemsio_module, only:  nemsio_init,nemsio_open,nemsio_close,nemsio_readrecv
    use nemsio_module, only:  nemsio_gfile,nemsio_getfilehead
    use nemsio_module, only:  nemsio_readrec, nemsio_writerec, nemsio_writerecv

    implicit none

! !INPUT PARAMETERS:

    real(r_kind),dimension(lat2,lon2),intent(in   ) :: dsfct     ! delta skin temperature
    integer(i_kind)                  ,intent(in   ) :: mype_so   ! mpi task to write output file

! !OUTPUT PARAMETERS:

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

!   Declare local parameters
    character(6), parameter:: fname_sfcges = 'sfcf06'
    character(6), parameter:: fname_sfcgcy = 'sfcgcy'
    character(6), parameter:: fname_sfctsk = 'sfctsk'
    character(6), parameter:: fname_sfcanl = 'sfcanl'
    character(6), parameter:: fname_nstges = 'nstf06'
    character(6), parameter:: fname_nstanl = 'nstanl'
    character(6), parameter:: fname_dtfanl = 'dtfanl'

!   Declare local variables
    integer(i_kind), parameter:: io_dtfanl = 54
    integer(i_kind), parameter:: nprep=15
    real(r_kind),parameter :: houra = zero_single
    character(len=120) :: my_name = 'WRITE_SFC_NST'
    character(len=1)   :: null = ' '
    integer(i_kind),dimension(7):: idate, jdate
    integer(i_kind),dimension(4):: odate
    integer(i_kind) :: i, j, ip1, jp1, ilat, ilon, mm1
    integer(i_kind) :: lonb, latb, nlatm2, n, nrec_sfc, nrec_nst, iret
    integer(i_kind) :: lonb_nst, latb_nst
    integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd
    integer(i_kind) :: istop = 106
    real(r_kind)    :: fhour
    real(r_single)  :: r_zsea1,r_zsea2

    real(r_kind),    dimension(lat1,lon1):: dsfct_sub
    integer(i_kind), dimension(lat1,lon1):: isli_sub

    real(r_kind),    dimension(max(iglobal,itotsub)):: dsfct_all
    integer(i_kind), dimension(max(iglobal,itotsub)):: isli_all

    real(r_kind),    dimension(nlat,nlon):: dsfct_glb,dsfct_tmp
    integer(i_kind), dimension(nlat,nlon):: isli_glb,isli_tmp

    real(r_kind),    dimension(nlat_sfc,nlon_sfc)  :: dsfct_gsi
    integer(i_kind), dimension(nlat_sfc,nlon_sfc)  :: isli_gsi

    real(r_kind),    dimension(nlon_sfc,nlat_sfc-2):: dsfct_anl
    real(r_single),  dimension(nlon_sfc,nlat_sfc-2):: dtzm
    real(r_single),  dimension(nlat_sfc,nlon_sfc)  :: work

    real(r_single),   allocatable, dimension(:,:) :: tsea,xt,xs,xu,xv,xz,zm,xtts,xzts,dt_cool,z_c, &
                                                     c_0,c_d,w_0,w_d,d_conv,ifd,tref,qrain
    real(r_single),   allocatable, dimension(:,:) :: slmsk_ges,slmsk_anl
    real(r_single),   allocatable, dimension(:)   :: rwork1d

    type(nemsio_gfile) :: gfile_sfcges,gfile_sfcgcy,gfile_nstges,gfile_sfctsk,gfile_sfcanl,gfile_nstanl

!*****************************************************************************

!   Initialize local variables
    mm1=mype+1
    nlatm2=nlat-2
!
!   Extract the analysis increment and surface mask in subdomain without the buffer
!
    do j=1,lon1
       jp1 = j+1
       do i=1,lat1
          ip1 = i+1
          dsfct_sub(i,j) = dsfct(ip1,jp1)
          isli_sub (i,j) = isli2(ip1,jp1)
       end do
    end do
!
!   Gather global analysis increment and surface mask info from subdomains
!
    call mpi_gatherv(dsfct_sub,ijn(mm1),mpi_rtype,&
         dsfct_all,ijn,displs_g,mpi_rtype,mype_so ,&
         mpi_comm_world,ierror)

    call mpi_gatherv(isli_sub,ijn(mm1),mpi_itype,&
         isli_all,ijn,displs_g,mpi_itype,mype_so ,&
         mpi_comm_world,ierror)

!   Only MPI task mype_so  writes the surface file.
    if (mype==mype_so ) then

      write(*,'(a,5(1x,a6))') 'write_nems_sfc_nst:',fname_sfcges,fname_nstges,fname_sfctsk,fname_sfcanl,fname_nstanl
!
!     get Tf analysis increment and surface mask at analysis (lower resolution) grids
!
      do i=1,iglobal
         ilon=ltosj(i)
         ilat=ltosi(i)
         dsfct_glb(ilat,ilon) = dsfct_all(i)
         isli_glb (ilat,ilon) = isli_all (i)
      end do
!
!      write dsfct_anl to a data file for later use (at eupd step at present)
!
       open(io_dtfanl,file=fname_dtfanl,form='unformatted')
       write(io_dtfanl) nlon,nlat
       write(io_dtfanl) dsfct_glb
       write(io_dtfanl) isli_glb

!      Initiate nemsio
       call nemsio_init(iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),null,null,'init',istop,iret)

!      open nsst guess file
       call nemsio_open(gfile_nstges,trim(fname_nstges),'read',iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),null,'open',istop,iret)
!      open surface guess file
       call nemsio_open(gfile_sfcges,trim(fname_sfcges),'read',iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcges),null,'open',istop,iret)
!      open surface gcycle file
       call nemsio_open(gfile_sfcgcy,trim(fname_sfcgcy),'read',iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcgcy),null,'open',istop,iret)

!      read a few surface guess file header records
       call nemsio_getfilehead(gfile_sfcges, nrec=nrec_sfc, idate=idate, &
          dimx=lonb, dimy=latb, nfhour=nfhour, nfminute=nfminute, &
          nfsecondn=nfsecondn, nfsecondd=nfsecondd, iret=iret)

!      read some nsst guess file header records (dimensions)
       call nemsio_getfilehead(gfile_nstges, nrec=nrec_nst, dimx=lonb_nst,dimy=latb_nst,iret=iret)

       write(6,*) 'nrec_sfc, nrec_nst = ',nrec_sfc, nrec_nst

!      check the dimensions consistency in sfc, nst files and the used.
       if ( latb /= latb_nst .or. lonb /= lonb_nst ) then
          write(6,*) 'Inconsistent dimension for sfc & nst files. latb,lonb : ',latb,lonb, &
                     'latb_nst,lonb_nst : ',latb_nst,lonb_nst
          call stop2(80)
       endif

       if ( nlat_sfc /= latb+2 .or. nlon_sfc /= lonb ) then
          write(6,*) 'Inconsistent dimension for used and read. nlat_sfc,nlon_sfc : ',nlat_sfc,nlon_sfc, &
                     'latb+2,lonb :',latb+2,lonb
          call stop2(81)
       endif
!
       allocate(slmsk_ges(lonb,latb),slmsk_anl(lonb,latb))
       allocate(rwork1d(lonb*latb))

!      read slmsk in fname_sfcges to get slmsk_ges
       call nemsio_readrecv(gfile_sfcges, 'land', 'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcges),'land','read',istop,iret)
       slmsk_ges=reshape(rwork1d,(/size(slmsk_ges,1),size(slmsk_ges,2)/))

!      read slmsk in fname_sfcgcy to get slmsk_anl
       call nemsio_readrecv(gfile_sfcgcy, 'land', 'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcgcy),'land','read',istop,iret)
       slmsk_anl=reshape(rwork1d,(/size(slmsk_anl,1),size(slmsk_anl,2)/))
!
!      Replace header record date with analysis time from iadate
!
       jdate(1) = iadate(1)  ! analysis year
       jdate(2) = iadate(2)  ! analysis month
       jdate(3) = iadate(3)  ! analysis day
       jdate(4) = iadate(4)  ! analysis hour
       jdate(5) = iadate(5)  ! analysis minute
       jdate(5) = 0          ! analysis minute
       jdate(6) = 0          ! analysis scaled seconds
       jdate(7) = idate(7)   ! analysis seconds multiplier

       nfhour=0              !  new forecast hour, zero at analysis time
       nfminute=0
       nfsecondn=0
       nfsecondd=100         ! default for denominator

       fhour    = zero
       odate(1) = jdate(4)   !hour
       odate(2) = jdate(2)   !month
       odate(3) = jdate(3)   !day
       odate(4) = jdate(1)   !year

       if ( (latb /= nlatm2) .or. (lonb /= nlon) ) then
          write(6,*)'WRITE_NEMSIO_SFC_NST:  different grid dimensions analysis vs sfc. interpolating sfc temperature  ',&
               ', nlon,nlat-2=',nlon,nlatm2,' -vs- sfc file lonb,latb=',lonb,latb
          write(6,*) ' WRITE_NEMSIO_SFC_NST, nlon_sfc,nlat_sfc : ',  nlon_sfc,nlat_sfc
!
!         Get the expanded values for a surface type (0 = water now) and the new mask
!
          call int2_msk_glb_prep(dsfct_glb,isli_glb,dsfct_tmp,isli_tmp,nlat,nlon,0,nprep)
!
!         Get updated/analysis surface mask info from sfcgcy file
!
          call tran_gfssfc(slmsk_anl,work,lonb,latb)
          do j=1,lonb
             do i=1,latb+2
                isli_gsi(i,j) = nint(work(i,j))
             end do
          end do
!
!         Interpolate dsfct_tmp(nlat,nlon) to dsfct_gsi(nlat_sfc,nlon_sfc) with surface mask accounted
!
          call int22_msk_glb(dsfct_tmp,isli_tmp,rlats,rlons,nlat,nlon, &
                             dsfct_gsi,isli_gsi,rlats_sfc,rlons_sfc,nlat_sfc,nlon_sfc,0)
!
!         transform the dsfct_gsi(latb+2,lonb) to dsfct_anl(lonb,latb) for sfc file format
!
          do j = 1, latb
             do i = 1, lonb
                dsfct_anl(i,j) = dsfct_gsi(latb+2-j,i)
             end do
          end do

       else
!
!         transform the dsfct_glb(nlat,nlon) to dsfct_anl(lonb,latb) for sfc file
!         format when nlat == latb-2 & nlon = lonb
!
          do j=1,latb
             do i=1,lonb
                dsfct_anl(i,j)=dsfct_glb(latb+1-j,i)
             end do
          end do
       endif                 ! if ( (latb /= nlatm2) .or. (lonb /= nlon) ) then

!
! Start to write output sfc file : fname_sfcanl & fname_nstanl
!      open new output file with new header gfile_sfcanl and gfile_nstanl with "write" access.
!      Use this call to update header as well
!
!      copy input header info to output header info for sfcanl, need to do this before nemsio_close(gfile)
!
       gfile_sfcanl=gfile_sfcgcy
!      open nemsio sfcanl
       call nemsio_open(gfile_sfcanl,trim(fname_sfcanl),'write',iret=iret, idate=jdate, nfhour=nfhour,&
          nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd )
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcanl),null,'open',istop,iret)

       gfile_sfctsk=gfile_sfcgcy
!      open nemsio sfctsk
       call nemsio_open(gfile_sfctsk,trim(fname_sfctsk),'write',iret=iret, idate=jdate, nfhour=nfhour,&
          nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd )
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfctsk),null,'open',istop,iret)
!
!      copy input header info to output header info for nstanl, need to do this before nemsio_close(gfile)
!
       gfile_nstanl=gfile_nstges
!      open nemsio nstanl
       call nemsio_open(gfile_nstanl,trim(fname_nstanl),'write',iret=iret, idate=jdate, nfhour=nfhour,&
          nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd )
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),null,'open',istop,iret)
! Allocate work array (rwork1d) and tsea in sfc file
       allocate(tsea(lonb,latb))

! Allocate nsst variables
       allocate(xt(lonb,latb))
       allocate(xs(lonb,latb))
       allocate(xu(lonb,latb))
       allocate(xv(lonb,latb))
       allocate(xz(lonb,latb))
       allocate(zm(lonb,latb))
       allocate(xtts(lonb,latb))
       allocate(xzts(lonb,latb))
       allocate(dt_cool(lonb,latb))
       allocate(z_c(lonb,latb))
       allocate(c_0(lonb,latb))
       allocate(c_d(lonb,latb))
       allocate(w_0(lonb,latb))
       allocate(w_d(lonb,latb))
       allocate(d_conv(lonb,latb))
       allocate(ifd(lonb,latb))
       allocate(tref(lonb,latb))
       allocate(qrain(lonb,latb))
!
!      First copy entire data from sfcgcy to fname_anl, then do selective update
!
!      read the nrec_sfc variables from sfcgcy and then write then to sfcanl
!
       do n = 1, nrec_sfc
          call nemsio_readrec(gfile_sfcgcy,n,rwork1d,iret=iret)
          if ( iret /= 0 ) write(6,*) 'readrec for gfile_sfcgcy,  nrec_sfc = ', n, '  Status = ', iret
          call nemsio_writerec(gfile_sfcanl,n,rwork1d,iret=iret)
          if ( iret /= 0 ) write(6,*) 'writerec for gfile_sfcanl, nrec_sfc = ', n, '  Status = ', iret
          call nemsio_writerec(gfile_sfctsk,n,rwork1d,iret=iret)
          if ( iret /= 0 ) write(6,*) 'writerec for gfile_sfctsk, nrec_sfc = ', n, '  Status = ', iret
       end do

       write(*,*) 'read gfile_sfcgcy, and the write to gfile_sfcanl, gfile_sfctsk'
!
!      For sfcanl, Only tsea (sea surface temperature) will be updated in the SFC
!                  Need values from nstges for tref update
!      read tsea from sfcges
       call nemsio_readrecv(gfile_sfcges,'tmp','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcges),'tmp','read',istop,iret)
       tsea=reshape(rwork1d,(/size(tsea,1),size(tsea,2)/))

!      For nstanl, Only tref (foundation temperature) is updated by analysis
!                  others are updated for snow melting case
!      read 18 nsst variables from nstges
! xt
       call nemsio_readrecv(gfile_nstges, 'xt',    'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xt','read',istop,iret)
       xt=reshape(rwork1d,(/size(xt,1),size(xt,2)/))
! xs
       call nemsio_readrecv(gfile_nstges, 'xs',    'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xs','read',istop,iret)
       xs=reshape(rwork1d,(/size(xs,1),size(xs,2)/))
! xu
       call nemsio_readrecv(gfile_nstges, 'xu',    'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xu','read',istop,iret)
       xu=reshape(rwork1d,(/size(xu,1),size(xu,2)/))
! xv
       call nemsio_readrecv(gfile_nstges, 'xv',    'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xv','read',istop,iret)
       xv=reshape(rwork1d,(/size(xv,1),size(xv,2)/))
! xz
       call nemsio_readrecv(gfile_nstges, 'xz',    'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xz','read',istop,iret)
       xz=reshape(rwork1d,(/size(xz,1),size(xz,2)/))
! zm
       call nemsio_readrecv(gfile_nstges, 'zm',    'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'zm','read',istop,iret)
       zm=reshape(rwork1d,(/size(zm,1),size(zm,2)/))
! xtts
       call nemsio_readrecv(gfile_nstges, 'xtts',    'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xtts','read',istop,iret)
       xtts=reshape(rwork1d,(/size(xtts,1),size(xtts,2)/))
! xzts
       call nemsio_readrecv(gfile_nstges, 'xzts',    'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'xzts','read',istop,iret)
       xzts=reshape(rwork1d,(/size(xzts,1),size(xzts,2)/))
! dt_cool
       call nemsio_readrecv(gfile_nstges, 'dtcool','sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'dt_cool','read',istop,iret)
       dt_cool=reshape(rwork1d,(/size(dt_cool,1),size(dt_cool,2)/))
! z_c
       call nemsio_readrecv(gfile_nstges, 'zc','sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'zc','read',istop,iret)
       z_c=reshape(rwork1d,(/size(z_c,1),size(z_c,2)/))
! c_0
       call nemsio_readrecv(gfile_nstges, 'c0','sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'c0','read',istop,iret)
       c_0=reshape(rwork1d,(/size(c_0,1),size(c_0,2)/))
! c_d
       call nemsio_readrecv(gfile_nstges, 'cd','sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'cd','read',istop,iret)
       c_d=reshape(rwork1d,(/size(c_d,1),size(c_d,2)/))
! w_0
       call nemsio_readrecv(gfile_nstges, 'w0','sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'w0','read',istop,iret)
       w_0=reshape(rwork1d,(/size(w_0,1),size(w_0,2)/))
! w_d
       call nemsio_readrecv(gfile_nstges, 'wd','sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'wd','read',istop,iret)
       w_d=reshape(rwork1d,(/size(w_d,1),size(w_d,2)/))
! tref
       call nemsio_readrecv(gfile_nstges, 'tref',  'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'tref','read',istop,iret)
       tref=reshape(rwork1d,(/size(tref,1),size(tref,2)/))
! d_conv
       call nemsio_readrecv(gfile_nstges, 'dconv',  'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'dconv','read',istop,iret)
       d_conv=reshape(rwork1d,(/size(d_conv,1),size(d_conv,2)/))
! ifd
       call nemsio_readrecv(gfile_nstges, 'ifd',  'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'ifd','read',istop,iret)
       ifd=reshape(rwork1d,(/size(ifd,1),size(ifd,2)/))
! qrain
       call nemsio_readrecv(gfile_nstges, 'qrain',  'sfc', 1, rwork1d, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),'qrain','read',istop,iret)
       qrain=reshape(rwork1d,(/size(qrain,1),size(qrain,2)/))
!
!      update tref (in nst file) & tsea (in the surface file) when Tr analysis is on
!      reset NSSTM variables for new open water grids
!
       if ( nst_gsi > 2 ) then
!
!         For the new open water (sea ice just melted) grids, (1) set dsfct_anl = zero; (2) reset the NSSTM variables
!
!         Notes: slmsk_ges is the mask of the background
!                slmsk_anl is the mask of the analysis
!
          where ( slmsk_anl(:,:) == zero .and. slmsk_ges(:,:) == two )

            dsfct_anl(:,:)        = zero

            xt(:,:)      = zero
            xs(:,:)      = zero
            xu(:,:)      = zero
            xv(:,:)      = zero
            xz(:,:)      = z_w_max
            zm(:,:)      = zero
            xtts(:,:)    = zero
            xzts(:,:)    = zero
            dt_cool(:,:) = zero
            z_c(:,:)     = zero
            c_0(:,:)     = zero
            c_d(:,:)     = zero
            w_0(:,:)     = zero
            w_d(:,:)     = zero
            d_conv(:,:)  = zero
            ifd(:,:)     = zero
            tref(:,:)    = tfrozen
            qrain(:,:)   = zero
          end where
!
!         update analysis variable: Tref (foundation temperature) for nst file
!
          where ( slmsk_anl(:,:) == zero )
             tref(:,:) = max(tref(:,:) + dsfct_anl(:,:),tfrozen)
          elsewhere
             tref(:,:) = tsea(:,:)
          end where
!
!         update SST: tsea for sfc file with NSST profile
!
          r_zsea1 = 0.001_r_single*real(zsea1)
          r_zsea2 = 0.001_r_single*real(zsea2)
          call dtzm_2d(xt,xz,dt_cool,z_c,slmsk_anl,r_zsea1,r_zsea2,lonb,latb,dtzm)

          where ( slmsk_anl(:,:) == zero )
             tsea(:,:) = max(tref(:,:) + dtzm(:,:), tfrozen)
          end where

       else          ! when (nst_gsi <= 2)

          do j=1,latb
             do i=1,lonb
                tref(i,j) = tsea(i,j)  ! keep tref as tsea before analysis
             end do
          end do
!
!         For the new open water (sea ice just melted) grids, reset the NSSTM variables
!
          where ( slmsk_anl(:,:) == zero .and. slmsk_ges(:,:) == two )

            xt(:,:)      = zero
            xs(:,:)      = zero
            xu(:,:)      = zero
            xv(:,:)      = zero
            xz(:,:)      = z_w_max
            zm(:,:)      = zero
            xtts(:,:)    = zero
            xzts(:,:)    = zero
            dt_cool(:,:) = zero
            z_c(:,:)     = zero
            c_0(:,:)     = zero
            c_d(:,:)     = zero
            w_0(:,:)     = zero
            w_d(:,:)     = zero
            d_conv(:,:)  = zero
            ifd(:,:)     = zero
            tref(:,:)    = tfrozen
            qrain(:,:)   = zero
          end where
!
!         update tsea when NO Tf analysis
!
          do j=1,latb
             do i=1,lonb
                tsea(i,j) = max(tsea(i,j) + dsfct_anl(i,j),tfrozen)
             end do
          end do

       endif                   ! if ( nst_gsi > 2 ) then
!
!      update tsea record in sfcanl
!
       rwork1d = reshape(tsea, (/size(rwork1d)/) )
       call nemsio_writerecv(gfile_sfcanl,'tmp','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcanl),'tmp','write',istop,iret)
       write(6,100) fname_sfcanl,lonb,latb,houra,iadate(1:4),iret
100    format(' WRITE_NEMSIO_SFC_NST:  update tsea in ',a6,2i6,1x,f4.1,4(i4,1x),' with iret=',i2)
!
!      update tsea record in sfctsk
!
       rwork1d = reshape(tsea, (/size(rwork1d)/) )
       call nemsio_writerecv(gfile_sfctsk,'tmp','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfctsk),'tmp','write',istop,iret)
       write(6,101) fname_sfctsk,lonb,latb,houra,iadate(1:4),iret
101    format(' WRITE_NEMSIO_SFC_NST:  update tsea in ',a6,2i6,1x,f4.1,4(i4,1x),' with iret=',i2)
!
!      update nsst records in nstanl
!
! slmsk
       rwork1d = reshape( slmsk_anl,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'land','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'land','write',istop,iret)
! xt
       rwork1d = reshape( xt,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'xt','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xt','write',istop,iret)
! xs
       rwork1d = reshape( xs,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'xs','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xs','write',istop,iret)
! xu
       rwork1d = reshape( xu,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'xu','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xu','write',istop,iret)
! xv
       rwork1d = reshape( xv,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'xv','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xv','write',istop,iret)
! xz
       rwork1d = reshape( xz,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'xz','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xz','write',istop,iret)
! zm
       rwork1d = reshape( zm,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'zm','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'zm','write',istop,iret)
! xtts
       rwork1d = reshape( xtts,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'xtts','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xtts','write',istop,iret)
! xzts
       rwork1d = reshape( xzts,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'xzts','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'xzts','write',istop,iret)
! z_0
       rwork1d = reshape( dt_cool,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'dtcool','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'dtcool','write',istop,iret)
! z_c
       rwork1d = reshape( z_c,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'zc','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'zc','write',istop,iret)
! c_0
       rwork1d = reshape( c_0,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'c0','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'c0','write',istop,iret)
! c_d
       rwork1d = reshape( c_d,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'cd','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'cd','write',istop,iret)
! w_0
       rwork1d = reshape( w_0,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'w0','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'w0','write',istop,iret)
! w_d
       rwork1d = reshape( w_d,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'wd','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'wd','write',istop,iret)
! d_conv
       rwork1d = reshape( d_conv,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'dconv','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'dconv','write',istop,iret)
! ifd
       rwork1d = reshape( ifd,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'ifd','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'ifd','write',istop,iret)
! tref
       rwork1d = reshape( tref,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'tref','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'tref','write',istop,iret)
! qrain
       rwork1d = reshape( qrain,(/size(rwork1d)/) )
       call nemsio_writerecv(gfile_nstanl,'qrain','sfc',1,rwork1d,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),'qrain','write',istop,iret)

       write(6,200) fname_nstanl,lonb,latb,houra,iadate(1:4),iret
200    format(' WRITE_NEMSIO_SFC_NST:  update variables in ',a6,2i6,1x,f4.1,4(i4,1x),' with iret=',i2)

       deallocate(xt,xs,xu,xv,xz,zm,xtts,xzts,dt_cool,z_c,c_0,c_d,w_0,w_d,d_conv,ifd,tref,qrain)
       deallocate(rwork1d)

       call nemsio_close(gfile_sfcges, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcges),null,'close',istop,iret)

       call nemsio_close(gfile_sfcgcy, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcgcy),null,'close',istop,iret)

       call nemsio_close(gfile_nstges, iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstges),null,'close',istop,iret)

       call nemsio_close(gfile_sfcanl,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfcanl),null,'close',istop,iret)

       call nemsio_close(gfile_nstanl,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_nstanl),null,'close',istop,iret)

       call nemsio_close(gfile_sfctsk,iret=iret)
       if (iret /= 0) call error_msg(trim(my_name),trim(fname_sfctsk),null,'close',istop,iret)

       write(6,'(a,'': nemsio sfc_nst anal written for lonb,latb= '',2i6,'',valid hour= '',f4.1,'',idate= '',4i5)') &
          trim(my_name),lonb,latb,fhour,odate
    endif
  end subroutine write_sfc_nst_

  subroutine error_msg_(sub_name,file_name,var_name,action,stop_code,error_code,lprint)
    use mpimod, only: mype
    use kinds, only: i_kind
    implicit none

    character(len=*), intent(in) :: sub_name,file_name,var_name,action
    integer(i_kind),  intent(in) :: stop_code, error_code
    logical, optional,intent(in) :: lprint
    
    if ( mype == 0 .or. present(lprint) ) then
       select case (trim(action))
       case('init')
          write(6,'(a,'':  PROBLEM with nemsio_init, Status = '', i3)') &
             trim(sub_name), error_code
       case('open')
          write(6,'(a,'':  ***ERROR*** problem opening file '',a,'', Status = '', i3)') &
             trim(sub_name), trim(file_name), error_code
       case('close')
          write(6,'(a,'':  ***ERROR*** problem closing file '',a,'', Status = '', i3)') &
             trim(sub_name), trim(file_name), error_code
       case default
          write(6,'(a,'':  ***ERROR*** '',a,tr1,a,'',variable = '',a,'',Status = '',i3)') &
             trim(sub_name),trim(action),trim(file_name),trim(var_name),error_code
       end select
     end if
     if ( stop_code /= 0 ) call stop2(stop_code)
  end subroutine error_msg_

  subroutine intrp22(a,rlons_a,rlats_a,nlon_a,nlat_a, &
                     b,rlons_b,rlats_b,nlon_b,nlat_b)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    intrp22 --- interpolates from one 2-d grid to another 2-d grid
!                            like analysis to surface grid or vice versa
!   prgrmmr:     li -  initial version; org: np2
!
! abstract:      This routine interpolates a grid to b grid
!
! program history log:
!
!   input argument list:
!     rlons_a - longitudes of input array
!     rlats_a - latitudes of input array
!     nlon_a  - number of longitude of input array
!     nlat_a  - number of latitude of input array
!     rlons_b - longitudes of output array
!     rlats_b - latitudes of output array
!     nlon_b  - number of longitude of output array
!     nlat_b  - number of latitude of output array
!     a       - input values
!
!   output argument list:
!     b       - output values
!
! attributes:
!   language: f90
!   machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP
!
!$$$ end documentation block

! !USES:
    use kinds, only: r_kind,i_kind,r_single
    use constants, only: zero,one

    implicit none

! !INPUT PARAMETERS:
    integer(i_kind)                 ,intent(in   ) :: nlon_a,nlat_a,nlon_b,nlat_b
    real(r_kind), dimension(nlon_a) ,intent(in   ) :: rlons_a
    real(r_kind), dimension(nlat_a) ,intent(in   ) :: rlats_a
    real(r_kind), dimension(nlon_b) ,intent(in   ) :: rlons_b
    real(r_kind), dimension(nlat_b) ,intent(in   ) :: rlats_b

    real(r_single), dimension(nlon_a,nlat_a),intent(in   ) :: a

! !OUTPUT PARAMETERS:
    real(r_single), dimension(nlon_b,nlat_b),intent(  out) :: b

!   Declare local variables
    integer(i_kind) i,j,ix,iy,ixp,iyp
    real(r_kind) dx1,dy1,dx,dy,w00,w01,w10,w11,bout,dlat,dlon

!*****************************************************************************

    b=zero
!   Loop over all points to get interpolated value
    do j=1,nlat_b
       dlat=rlats_b(j)
       call grdcrd1(dlat,rlats_a,nlat_a,1)
       iy=int(dlat)
       iy=min(max(1,iy),nlat_a)
       dy  =dlat-iy
       dy1 =one-dy
       iyp=min(nlat_a,iy+1)

       do i=1,nlon_b
          dlon=rlons_b(i)
          call grdcrd1(dlon,rlons_a,nlon_a,1)
          ix=int(dlon)
          dx  =dlon-ix
          dx=max(zero,min(dx,one))
          dx1 =one-dx
          w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy

          ix=min(max(0,ix),nlon_a)
          ixp=ix+1
          if(ix==0) ix=nlon_a
          if(ixp==nlon_a+1) ixp=1
          bout=w00*a(ix,iy)+w01*a(ix,iyp)+w10*a(ixp,iy)+w11*a(ixp,iyp)
          b(i,j)=bout

       end do
    end do


!   End of routine
    return
  end subroutine intrp22

  subroutine tran_gfssfc(ain,aout,lonb,latb)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    tran_gfssfc     transform gfs surface file to analysis grid
!   prgmmr: derber          org: np2                date: 2003-04-10
!
! abstract: transform gfs surface file to analysis grid
!
! program history log:
!   2012-31-38  derber  - initial routine
!
!   input argument list:
!     ain      - input surface record on processor iope
!     lonb     - input number of longitudes
!     latb     - input number of latitudes
!
!   output argument list:
!     aout     - output transposed surface record
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind,r_single
  use constants, only: zero
  use sfcio_module, only: sfcio_realkind
  implicit none

! Declare passed variables
  integer(i_kind)                  ,intent(in ) :: lonb,latb
  real(sfcio_realkind),dimension(lonb,latb),intent(in ) :: ain
  real(r_single),dimension(latb+2,lonb),intent(out) :: aout

! Declare local variables
  integer(i_kind) i,j
  real(r_kind) sumn,sums
! of surface guess array
  sumn = zero
  sums = zero
  do i=1,lonb
     sumn = ain(i,1)    + sumn
     sums = ain(i,latb) + sums
  end do
  sumn = sumn/float(lonb)
  sums = sums/float(lonb)
!  Transfer from local work array to surface guess array
  do j = 1,lonb
     aout(1,j)=sums
     do i=2,latb+1
        aout(i,j) = ain(j,latb+2-i)
     end do
     aout(latb+2,j)=sumn
  end do

  return
  end subroutine tran_gfssfc

end module ncepnems_io