INITPOST_NETCDF.f

Go to the documentation of this file.

!         and 2D diag. output (d2d_chem) for GEFS-Aerosols and CCPP-Chem model.
      SUBROUTINE initpost_netcdf(ncid2d,ncid3d)


      use netcdf
      use vrbls4d, only: dust, salt, suso, soot, waso, smoke, fv3dust, coarsepm,                &
              no3,nh4, pp25, pp10 
      use vrbls3d, only: t, q, uh, vh, pmid, pint, alpint, dpres, zint, zmid, o3,               &
              qqr, qqs, cwm, qqi, qqw, omga, rhomid, q2, cfr, rlwtt, rswtt, tcucn,              &
              tcucns, train, el_pbl, exch_h, vdifftt, vdiffmois, dconvmois, nradtt,             &
              o3vdiff, o3prod, o3tndy, mwpv, unknown, vdiffzacce, zgdrag,cnvctummixing,         &
              vdiffmacce, mgdrag, cnvctvmmixing, ncnvctcfrac, cnvctumflx, cnvctdmflx,           &
              cnvctzgdrag, sconvmois, cnvctmgdrag, cnvctdetmflx, duwt, duem, dusd, dudp,        &
              dusv,ssem,sssd,ssdp,sswt,sssv,bcem,bcsd,bcdp,bcwt,bcsv,ocem,ocsd,ocdp,ocwt,ocsv, &
              wh, qqg, ref_10cm, qqnifa, qqnwfa, avgpmtf, avgozcon, aextc55, taod5503d,         &
              effri, effrl, effrs

      use vrbls2d, only: f, pd, fis, pblh, ustar, z0, ths, qs, twbs, qwbs, avgcprate,           &
              cprate, avgprec, prec, lspa, sno, sndepac, si, cldefi, th10, q10, tshltr, pshltr, &
              tshltr, albase, avgalbedo, avgtcdc, czen, czmean, mxsnal, landfrac, radot, sigt4, &
              cfrach, cfracl, cfracm, avgcfrach, qshltr, avgcfracl, avgcfracm, cnvcfr,          &
              islope, cmc, grnflx, vegfrc, acfrcv, ncfrcv, acfrst, ncfrst, ssroff,              &
              bgroff, rlwin, rlwtoa, cldwork, alwin, alwout, alwtoa, rswin, rswinc,             &
              rswout, aswin, auvbin, auvbinc, aswout, aswtoa, sfcshx, sfclhx, subshx,           &
              snopcx, sfcux, sfcvx, sfcuxi, sfcvxi, sfcuvx, gtaux, gtauy, potevp, u10, v10, smstav,&
              smstot, ivgtyp, isltyp, sfcevp, sfcexc, acsnow, acsnom, sst, thz0, qz0,           &
              uz0, vz0, ptop, htop, pbot, hbot, ptopl, pbotl, ttopl, ptopm, pbotm, ttopm,       &
              ptoph, pboth, pblcfr, ttoph, runoff, tecan, tetran, tedir, twa, maxtshltr,        &
              mintshltr, maxrhshltr, fdnsst, acgraup, graup_bucket, acfrain, frzrn_bucket,      &
              snow_acm, snow_bkt, snownc, graupelnc, qrmax,                                     &
              minrhshltr, dzice, smcwlt, suntime, fieldcapa, htopd, hbotd, htops, hbots,        &
              cuppt, dusmass, ducmass, dusmass25, ducmass25, aswintoa,rel_vort_maxhy1,          &
              maxqshltr, minqshltr, acond, sr, u10h, v10h,refd_max, w_up_max, w_dn_max,         &
              up_heli_max,up_heli_min,up_heli_max03,up_heli_min03,rel_vort_max01,u10max, v10max,  &
              avgedir,avgecan,paha,pahi,avgetrans,avgesnow,avgprec_cont,avgcprate_cont,rel_vort_max, &
              avisbeamswin,avisdiffswin,airbeamswin,airdiffswin,refdm10c_max,wspd10max, &
              alwoutc,alwtoac,aswoutc,aswtoac,alwinc,aswinc,avgpotevp,snoavg, &
              ti,aod550,du_aod550,ss_aod550,su_aod550,oc_aod550,bc_aod550,prate_max,maod,dustpm10, &
              dustcb,bccb,occb,sulfcb,sscb,dustallcb,ssallcb,dustpm,sspm,pp25cb,pp10cb,no3cb,nh4cb,&
              pwat, ebb, hwp, aqm_aod550, ltg1_max,ltg2_max,ltg3_max
      use soil,  only: sldpth, sllevel, sh2o, smc, stc
      use masks, only: lmv, lmh, htm, vtm, gdlat, gdlon, dx, dy, hbm2, sm, sice
      use physcons_post, only: grav => con_g, fv => con_fvirt, rgas => con_rd,                     &
                            eps => con_eps, epsm1 => con_epsm1
      use params_mod, only: erad, dtr, tfrz, h1, d608, rd, p1000, capa,pi
      use lookup_mod, only: thl, plq, ptbl, ttbl, rdq, rdth, rdp, rdthe, pl, qs0, sqs, sthe,    &
                            ttblq, rdpq, rdtheq, stheq, the0q, the0
      use ctlblk_mod, only: me, mpi_comm_comp, icnt, idsp, jsta, jend, ihrst, idat, sdat, ifhr, &
              ifmin, filename, tprec, tclod, trdlw, trdsw, tsrfc, tmaxmin, td3d, restrt, sdat,  &
              jend_m, imin, imp_physics, dt, spval, pdtop, pt, qmin, nbin_du, nphs, dtq2, ardlw,&
              ardsw, asrfc, avrain, avcnvc, theat, gdsdegr, spl, lsm, alsl, im, jm, im_jm, lm,  &
              jsta_2l, jend_2u, nsoil, lp1, icu_physics, ivegsrc, novegtype, nbin_ss, nbin_bc,  &
              nbin_oc, nbin_su, nbin_no3, nbin_nh4, gocart_on,gccpp_on, nasa_on,pt_tbl,hyb_sigp,&
              filenameflux, filenameaer,                                               &
              isf_surface_physics,rdaod, d2d_chem, modelname, aqf_on,                         &
              ista, iend, ista_2l, iend_2u,iend_m
      use gridspec_mod, only: maptype, gridtype, latstart, latlast, lonstart, lonlast, cenlon,  &
              dxval, dyval, truelat2, truelat1, psmapf, cenlat,lonstartv, lonlastv, cenlonv,    &
              latstartv, latlastv,cenlatv,latstart_r,latlast_r,lonstart_r,lonlast_r, standlon,  &
              latse,lonse,latnw,lonnw
      use upp_physics, only: fpvsnew
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
      implicit none
!
!      type(nemsio_gfile) :: nfile,ffile,rfile
      integer,parameter          :: nvar2d=48
!      character(nemsio_charkind) :: name2d(nvar2d)
      integer                    :: nvar3d, numdims
!      character(nemsio_charkind), allocatable :: name3din(:), name3dout(:)
!      character(nemsio_charkind) :: varname,levtype
!
!     INCLUDE/SET PARAMETERS.
!     
      include "mpif.h"

!     integer,parameter:: MAXPTS=1000000 ! max im*jm points
!
!     real,parameter:: con_g       =9.80665e+0! gravity
!     real,parameter:: con_rv      =4.6150e+2 ! gas constant H2O
!     real,parameter:: con_rd      =2.8705e+2 ! gas constant air
!     real,parameter:: con_fvirt   =con_rv/con_rd-1.
!     real,parameter:: con_eps     =con_rd/con_rv
!     real,parameter:: con_epsm1   =con_rd/con_rv-1
!
! This version of INITPOST shows how to initialize, open, read from, and
! close a NetCDF dataset. In order to change it to read an internal (binary)
! dataset, do a global replacement of _ncd_ with _int_. 

      real, parameter    :: gravi = 1.0/grav
      character(len=20)  :: varname, vcoordname
      integer            :: status, fldsize, fldst, recn, recn_vvel
      character             startdate*19,sysdepinfo*80,cgar*1
      character             startdate2(19)*4, flatlon*40
      logical            :: read_lonlat=.true.
! 
!     NOTE: SOME INTEGER VARIABLES ARE READ INTO DUMMY ( A REAL ). THIS IS OK
!     AS LONG AS REALS AND INTEGERS ARE THE SAME SIZE.
!
!     ALSO, EXTRACT IS CALLED WITH DUMMY ( A REAL ) EVEN WHEN THE NUMBERS ARE
!     INTEGERS - THIS IS OK AS LONG AS INTEGERS AND REALS ARE THE SAME SIZE.
      LOGICAL runb,singlrst,subpost,nest,hydro,ioomg,ioall
      logical, parameter :: debugprint = .false., zerout = .false.
!     logical, parameter :: debugprint = .true.,  zerout = .false.
      logical :: convert_rad_to_deg=.false.
      CHARACTER*32 varcharval 
!      CHARACTER*40 CONTRL,FILALL,FILMST,FILTMP,FILTKE,FILUNV,FILCLD,FILRAD,FILSFC
      CHARACTER*4  resthr
      CHARACTER    fname*255,envar*50
      INTEGER      idate(8),jdate(8),jpds(200),jgds(200),kpds(200),kgds(200)
!     LOGICAL*1    LB(IM,JM)
!     
!     INCLUDE COMMON BLOCKS.
!
!     DECLARE VARIABLES.
!     
!      REAL fhour
      integer nfhour ! forecast hour from nems io file
      integer fhzero !bucket
      real dtp !physics time step
      real dz
      REAL rinc(5)

      REAL dummy(im,jm)
!jw
      integer ii,jj,js,je,iyear,imn,iday,itmp,ioutcount,istatus,       &
              i,j,l,ll,k,kf,irtn,igdout,n,index,nframe,                &
              nframed2,iunitd3d,ierr,idum,iret,nrec,idrt
      integer ncid3d,ncid2d,varid,nhcas,varid_bl,iret_bl
      real    tstart,tlmh,tsph,es,fact,soilayert,soilayerb,zhour,dum,  &
              tvll,pmll,tv, tx1, tx2

      character*20,allocatable :: recname(:)
      integer,     allocatable :: reclev(:), kmsk(:,:)
      real,        allocatable :: glat1d(:), glon1d(:), qstl(:)
      real,        allocatable :: wrk1(:,:), wrk2(:,:)
      real,        allocatable :: p2d(:,:),  t2d(:,:),  q2d(:,:),      &
                                  qs2d(:,:), cw2d(:,:), cfr2d(:,:)
      real, dimension(lm+1)    :: ak5, bk5
      real*8, allocatable      :: pm2d(:,:), pi2d(:,:)
      real,   allocatable      :: tmp(:)
      real                     :: buf(ista_2l:iend_2u,jsta_2l:jend_2u)
      real                     :: buf2(ista_2l:iend_2u,jsta_2l:jend_2u)
      real                     :: buf3d(ista_2l:iend_2u,jsta_2l:jend_2u,lm)
      real                     :: chem_2d(ista_2l:iend_2u,jsta_2l:jend_2u)
      real                     :: chemt(ista_2l:iend_2u,jsta_2l:jend_2u,lm)
      real                     :: dt1(ista_2l:iend_2u,jsta_2l:jend_2u,lm)
      real                     :: dt2(ista_2l:iend_2u,jsta_2l:jend_2u,lm)
      real                     :: dt3(ista_2l:iend_2u,jsta_2l:jend_2u,lm)
      real                     :: dt4(ista_2l:iend_2u,jsta_2l:jend_2u,lm)
      real                     :: dt5(ista_2l:iend_2u,jsta_2l:jend_2u,lm)

!     real buf(ista_2l:iend_2u,jsta_2l:jend_2u),bufsoil(im,nsoil,jsta_2l:jend_2u)   &
!         ,buf3d(ista_2l:iend_2u,jsta_2l:jend_2u,lm),buf3d2(im,lp1,jsta_2l:jend_2u)

      real lat
      integer isa, jsa, latghf, jtem, idvc, idsl, nvcoord, ip1, nn, npass

      integer, parameter    :: npass2=5, npass3=30
      real, parameter       :: third=1.0/3.0
      INTEGER, DIMENSION(2) :: ij4min, ij4max
      REAL                  :: omgmin, omgmax
      real, allocatable :: d2d(:,:), u2d(:,:), v2d(:,:), omga2d(:,:)
      REAL, ALLOCATABLE :: ps2d(:,:),psx2d(:,:),psy2d(:,:)
      real, allocatable :: div3d(:,:,:)
      real(kind=4),allocatable :: vcrd(:,:)
      real                     :: dum_const 
      real, allocatable :: ext550(:,:,:)

      if (modelname == 'FV3R') then
         allocate(ext550(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
      endif

!***********************************************************************
!     START INIT HERE.
!
      WRITE(6,*)'INITPOST:  ENTER INITPOST_NETCDF'
      WRITE(6,*)'me=',me,  &
           'jsta_2l=',jsta_2l,'jend_2u=', &
           jend_2u,'im=',im, &
           'ista_2l=',ista_2l,'iend_2u=',iend_2u, &
           'ista=',ista,'iend=',iend, &
           'iend_m=',iend_m
!     
      isa = (ista+iend) / 2
      jsa = (jsta+jend) / 2

!$omp parallel do private(i,j)
      do j = jsta_2l, jend_2u
        do i= ista_2l, iend_2u
          buf(i,j) = spval
        enddo
      enddo

      status=nf90_get_att(ncid3d,nf90_global,'ak',ak5)
      if(status/=0)then
       print*,'ak not found; assigning missing value'
       ak5=spval
      else
       if(me==0)print*,'ak5= ',ak5
      end if 
      status=nf90_get_att(ncid3d,nf90_global,'idrt',idrt)
      if(status/=0)then
       print*,'idrt not in netcdf file,reading grid'
       status=nf90_get_att(ncid3d,nf90_global,'grid',varcharval)
       if(status/=0)then
        print*,'idrt and grid not in netcdf file, set default to latlon'
        idrt=0
        maptype=0
       else
        if(trim(varcharval)=='rotated_latlon')then
         maptype=207
         idrt=207
         status=nf90_get_att(ncid3d,nf90_global,'cen_lon',dum_const)
         if(status/=0)then
          print*,'cen_lon not found; assigning missing value'
          cenlon=spval
         else
          if(dum_const<0.)then
           cenlon=nint((dum_const+360.)*gdsdegr)
          else
           cenlon=dum_const*gdsdegr
          end if 
         end if
         status=nf90_get_att(ncid3d,nf90_global,'cen_lat',dum_const)
         if(status/=0)then
          print*,'cen_lat not found; assigning missing value'
          cenlat=spval
         else
          cenlat=dum_const*gdsdegr
         end if

         status=nf90_get_att(ncid3d,nf90_global,'lon1',dum_const)
         if(status/=0)then
          print*,'lonstart_r not found; assigning missing value'
          lonstart_r=spval
         else
          if(dum_const<0.)then
           lonstart_r=nint((dum_const+360.)*gdsdegr)
          else
           lonstart_r=dum_const*gdsdegr
          end if
         end if
         status=nf90_get_att(ncid3d,nf90_global,'lat1',dum_const)
         if(status/=0)then
          print*,'latstart_r not found; assigning missing value'
          latstart_r=spval
         else
          latstart_r=dum_const*gdsdegr
         end if

         status=nf90_get_att(ncid3d,nf90_global,'lon2',dum_const)
         if(status/=0)then
          print*,'lonlast_r not found; assigning missing value'
          lonlast_r=spval
         else
          if(dum_const<0.)then
           lonlast_r=nint((dum_const+360.)*gdsdegr)
          else
           lonlast_r=dum_const*gdsdegr
          end if
         end if
         status=nf90_get_att(ncid3d,nf90_global,'lat2',dum_const)
         if(status/=0)then
          print*,'latlast_r not found; assigning missing value'
          latlast_r=spval
         else
          latlast_r=dum_const*gdsdegr
         end if

         status=nf90_get_att(ncid3d,nf90_global,'dlon',dum_const)
         if(status/=0)then
          print*,'dlmd not found; assigning missing value'
          dxval=spval
         else
          dxval=dum_const*gdsdegr
         end if
         status=nf90_get_att(ncid3d,nf90_global,'dlat',dum_const)
         if(status/=0)then
          print*,'dphd not found; assigning missing value'
          dyval=spval
         else
          dyval=dum_const*gdsdegr
         end if

!         print*,'lonstart,latstart,cenlon,cenlat,dyval,dxval', &
!         lonstart,latstart,cenlon,cenlat,dyval,dxval

! Jili Dong add support for regular lat lon (2019/03/22) start
        else if(trim(varcharval)=='latlon')then
         maptype=0
         idrt=0

         status=nf90_get_att(ncid3d,nf90_global,'lon1',dum_const)
         if(status/=0)then
          print*,'lonstart not found; assigning missing value'
          lonstart=spval
         else
          if(dum_const<0.)then
           lonstart=nint((dum_const+360.)*gdsdegr)
          else
           lonstart=dum_const*gdsdegr
          end if
         end if
         status=nf90_get_att(ncid3d,nf90_global,'lat1',dum_const)
         if(status/=0)then
          print*,'latstart not found; assigning missing value'
          latstart=spval
         else
          latstart=dum_const*gdsdegr
         end if

         status=nf90_get_att(ncid3d,nf90_global,'lon2',dum_const)
         if(status/=0)then
          print*,'lonlast not found; assigning missing value'
          lonlast=spval
         else
          if(dum_const<0.)then
           lonlast=nint((dum_const+360.)*gdsdegr)
          else
           lonlast=dum_const*gdsdegr
          end if
         end if
         status=nf90_get_att(ncid3d,nf90_global,'lat2',dum_const)
         if(status/=0)then
          print*,'latlast not found; assigning missing value'
          latlast=spval
         else
          latlast=dum_const*gdsdegr
         end if

         status=nf90_get_att(ncid3d,nf90_global,'dlon',dum_const)
         if(status/=0)then
          print*,'dlmd not found; assigning missing value'
          dxval=spval
         else
          dxval=dum_const*gdsdegr
         end if
         status=nf90_get_att(ncid3d,nf90_global,'dlat',dum_const)
         if(status/=0)then
          print*,'dphd not found; assigning missing value'
          dyval=spval
         else
          dyval=dum_const*gdsdegr
         end if

         print*,'lonstart,latstart,dyval,dxval', &
         lonstart,lonlast,latstart,latlast,dyval,dxval

! Jili Dong add support for regular lat lon (2019/03/22) end 
 
        ELSE IF (trim(varcharval)=='lambert_conformal')then

          maptype=1
          idrt=1
          status=nf90_get_att(ncid3d,nf90_global,'cen_lon',dum_const)
          if(status/=0)then
            print*,'cen_lon not found; assigning missing value'
            cenlon=spval
          else
            if(dum_const<0.)then
              cenlon=nint((dum_const+360.)*gdsdegr)
            else
              cenlon=dum_const*gdsdegr
            end if
          end if
          status=nf90_get_att(ncid3d,nf90_global,'cen_lat',dum_const)
          if(status/=0)then
            print*,'cen_lat not found; assigning missing value'
            cenlat=spval
          else
            cenlat=dum_const*gdsdegr
          end if

          status=nf90_get_att(ncid3d,nf90_global,'lon1',dum_const)
          if(status/=0)then
            print*,'lonstart not found; assigning missing value'
            lonstart=spval
          else
            if(dum_const<0.)then
              lonstart=nint((dum_const+360.)*gdsdegr)
            else
              lonstart=dum_const*gdsdegr
            end if
          end if
          status=nf90_get_att(ncid3d,nf90_global,'lat1',dum_const)
          if(status/=0)then
            print*,'latstart not found; assigning missing value'
            latstart=spval
          else
            latstart=dum_const*gdsdegr
          end if

          status=nf90_get_att(ncid3d,nf90_global,'stdlat1',dum_const)
          if(status/=0)then
            print*,'stdlat1 not found; assigning missing value'
            truelat1=spval
          else
            truelat1=dum_const*gdsdegr
          end if
          status=nf90_get_att(ncid3d,nf90_global,'stdlat2',dum_const)
          if(status/=0)then
            print*,'stdlat2 not found; assigning missing value'
            truelat2=spval
          else
            truelat2=dum_const*gdsdegr
          end if

          status=nf90_get_att(ncid3d,nf90_global,'dx',dum_const)
          if(status/=0)then
            print*,'dx not found; assigning missing value'
            dxval=spval
          else
            dxval=dum_const*1.e3
          end if
          status=nf90_get_att(ncid3d,nf90_global,'dy',dum_const)
          if(status/=0)then
            print*,'dphd not found; assigning missing value'
            dyval=spval
          else
            dyval=dum_const*1.e3
          end if

          standlon = cenlon
          print*,'lonstart,latstart,cenlon,cenlat,truelat1,truelat2, &
                  stadlon,dyval,dxval', &
          lonstart,latstart,cenlon,cenlat,truelat1,truelat2,standlon,dyval,dxval

        else if(trim(varcharval)=='gaussian')then
         maptype=4
         idrt=4
        else ! setting default maptype
         maptype=0
         idrt=0
        end if
       end if
      end if
      if(me==0)print*,'idrt MAPTYPE= ',idrt,maptype
!     STEP 1.  READ MODEL OUTPUT FILE
!
!
!***
!
! LMH and LMV  always = LM for sigma-type vert coord

!$omp parallel do private(i,j)
      do j = jsta_2l, jend_2u
        do i = ista_2l, iend_2u
          lmv(i,j) = lm
          lmh(i,j) = lm
        end do
      end do

! HTM VTM all 1 for sigma-type vert coord

!$omp parallel do private(i,j,l)
      do l = 1, lm
        do j = jsta_2l, jend_2u
          do i = ista_2l, iend_2u
            htm(i,j,l) = 1.0
            vtm(i,j,l) = 1.0
          end do
        end do
      end do

      status=nf90_get_att(ncid3d,nf90_global,'nhcas',nhcas)
      if(status/=0)then
       print*,'nhcas not in netcdf file, set default to nonhydro'
       nhcas=0
      end if
      if(me==0)print*,'nhcas= ',nhcas
      if (nhcas == 0 ) then  !non-hydrostatic case
       nrec=18
       allocate (recname(nrec))
       recname=[character(len=20) :: 'ugrd','vgrd','spfh','tmp','o3mr', &
                                     'presnh','dzdt', 'clwmr','dpres',  &
                                     'delz','icmr','rwmr',              &
                                     'snmr','grle','smoke','dust',      &
                                     'coarsepm','ext550']
      else
       nrec=8
       allocate (recname(nrec))
       recname=[character(len=20) :: 'ugrd','vgrd','tmp','spfh','o3mr', &
                                     'hypres', 'clwmr','dpres']
      endif

!     write(*,*)'nrec=',nrec
      !allocate(recname(nrec),reclevtyp(nrec),reclev(nrec))
      allocate(glat1d(jm),glon1d(im))

! hardwire idate for now
!      idate=(/2017,08,07,00,0,0,0,0/)
! get cycle start time
      status=nf90_inq_varid(ncid3d,'time',varid)
      if(status/=0)then
       print*,'time not in netcdf file, stopping'
       stop 1
      else
       status=nf90_get_att(ncid3d,varid,'units',varcharval)
       if(status/=0)then
         print*,'time unit not available'
       else
         print*,'time unit read from netcdf file= ',varcharval
! assume use hours as unit
!       idate_loc=index(varcharval,'since')+6
         read(varcharval,101)idate(1),idate(2),idate(3),idate(4),idate(5)
       end if
!       Status=nf90_inquire_dimension(ncid3d,varid,len=ntimes)
!       allocate(fhours(ntimes))
!       status = nf90_inq_varid(ncid3d,varid,fhours)
!       Status=nf90_get_var(ncid3d,varid,nfhour,start=(/1/), &
!              count=(/1/))
!       if(Status/=0)then
!        print*,'forecast hour not in netcdf file, stopping'
!        STOP 1
!       end if
      end if
 101  format(t13,i4,1x,i2,1x,i2,1x,i2,1x,i2)
      print*,'idate= ',idate(1:5)

! Jili Dong check output format for coordinate reading
      status=nf90_inq_varid(ncid3d,'grid_xt',varid)
      status=nf90_inquire_variable(ncid3d,varid,ndims = numdims)
      if(numdims==1.and.modelname=="FV3R") then
        read_lonlat=.true.
      else
        read_lonlat=.false.
      end if
      

! Jili Dong add support for new write component output
! get longitude 
      if (read_lonlat) then
        status=nf90_inq_varid(ncid3d,'lon',varid)
        status=nf90_inquire_variable(ncid3d,varid,ndims = numdims)
        if(debugprint)print*,'number of dim for gdlon ',numdims
      else
        status=nf90_inq_varid(ncid3d,'grid_xt',varid)
        status=nf90_inquire_variable(ncid3d,varid,ndims = numdims)
        if(debugprint)print*,'number of dim for gdlon ',numdims
      end if
      if(numdims==1)then
        status=nf90_get_var(ncid3d,varid,glon1d)  
        do j=jsta,jend
          do i=ista,iend
            gdlon(i,j) = real(glon1d(i),kind=4)
          end do
        end do
        lonstart = nint(glon1d(1)*gdsdegr)
        lonlast  = nint(glon1d(im)*gdsdegr)

! Jili Dong add support for regular lat lon (2019/03/22) start
       if (maptype == 0) then
        if(lonstart<0.)then
         lonstart=lonstart+360.*gdsdegr
        end if
        if(lonlast<0.)then
         lonlast=lonlast+360.*gdsdegr
        end if
       end if
! Jili Dong add support for regular lat lon (2019/03/22) end 

      else if(numdims==2)then
        status=nf90_get_var(ncid3d,varid,dummy)
        if(maxval(abs(dummy))<2.0*pi)convert_rad_to_deg=.true. 
        if(convert_rad_to_deg)then
         do j=jsta,jend
          do i=ista,iend
            gdlon(i,j) = real(dummy(i,j),kind=4)*180./pi
          end do
         end do
        else
         do j=jsta,jend
          do i=ista,iend
            gdlon(i,j) = real(dummy(i,j),kind=4)
          end do
         end do
        end if
        if(convert_rad_to_deg)then
         lonstart = nint(dummy(1,1)*gdsdegr)*180./pi
         lonlast  = nint(dummy(im,jm)*gdsdegr)*180./pi
         lonse    = nint(dummy(im,1)*gdsdegr)*180./pi
         lonnw    = nint(dummy(1,jm)*gdsdegr)*180./pi
        else
         lonstart = nint(dummy(1,1)*gdsdegr)
         lonlast  = nint(dummy(im,jm)*gdsdegr)
         lonse    = nint(dummy(im,1)*gdsdegr)
         lonnw    = nint(dummy(1,jm)*gdsdegr)
        end if

! Jili Dong add support for regular lat lon (2019/03/22) start
       if (maptype == 0) then
        if(lonstart<0.)then
         lonstart=lonstart+360.*gdsdegr
        end if
        if(lonlast<0.)then
         lonlast=lonlast+360.*gdsdegr
        end if
       end if
! Jili Dong add support for regular lat lon (2019/03/22) end 

      end if
      print*,'lonstart,lonlast ',lonstart,lonlast 
! Jili Dong add support for new write component output
! get latitude
      if (read_lonlat) then
        status=nf90_inq_varid(ncid3d,'lat',varid)
        status=nf90_inquire_variable(ncid3d,varid,ndims = numdims)
        if(debugprint)print*,'number of dim for gdlat ',numdims
      else
        status=nf90_inq_varid(ncid3d,'grid_yt',varid)
        status=nf90_inquire_variable(ncid3d,varid,ndims = numdims)
        if(debugprint)print*,'number of dim for gdlat ',numdims
      end if
      if(numdims==1)then
        status=nf90_get_var(ncid3d,varid,glat1d)
        do j=jsta,jend
          do i=ista,iend
            gdlat(i,j) = real(glat1d(j),kind=4)
          end do
        end do
        latstart = nint(glat1d(1)*gdsdegr)
        latlast  = nint(glat1d(jm)*gdsdegr)
      else if(numdims==2)then
        status=nf90_get_var(ncid3d,varid,dummy)
        if(maxval(abs(dummy))<pi)convert_rad_to_deg=.true.
        if(convert_rad_to_deg)then
         do j=jsta,jend
          do i=ista,iend
            gdlat(i,j) = real(dummy(i,j),kind=4)*180./pi
          end do
         end do
        else
         do j=jsta,jend
          do i=ista,iend
            gdlat(i,j) = real(dummy(i,j),kind=4)
          end do
         end do
        end if
        if(convert_rad_to_deg)then
         latstart = nint(dummy(1,1)*gdsdegr)*180./pi
         latlast  = nint(dummy(im,jm)*gdsdegr)*180./pi
         latse    = nint(dummy(im,1)*gdsdegr)*180./pi
         latnw    = nint(dummy(1,jm)*gdsdegr)*180./pi
        else
         latstart = nint(dummy(1,1)*gdsdegr)
         latlast  = nint(dummy(im,jm)*gdsdegr)
         latse    = nint(dummy(im,1)*gdsdegr)
         latnw    = nint(dummy(1,jm)*gdsdegr)
        end if
      end if
      print*,'laststart,latlast = ',latstart,latlast
      if(debugprint)print*,'me sample gdlon gdlat= ' &
     ,me,gdlon(isa,jsa),gdlat(isa,jsa)

! Specify grid staggering type
      gridtype = 'A'
      maptype=idrt
      if (me == 0) print *, 'maptype and gridtype is ', &
      maptype,gridtype
 
      if(gridtype == 'A')then
        lonstartv=lonstart
        lonlastv=lonlast
        latstartv=latstart
        latlastv=latlast
        cenlatv=cenlat
        cenlonv=cenlon
      end if

      if(debugprint)then
        if (me == 0)then
          do i=1,nrec
            print *,'recname=',trim(recname(i))
          end do
        end if
      end if


      deallocate(glat1d,glon1d)

      print*,'idate = ',(idate(i),i=1,7)
!      print*,'nfhour = ',nfhour
      
! sample print point
      ii = im/2
      jj = jm/2
      
      print *,me,'max(gdlat)=', maxval(gdlat),  &
                 'max(gdlon)=', maxval(gdlon)
      CALL exch(gdlat(ista_2l,jsta_2l))
      CALL exch(gdlon(ista_2l,jsta_2l))
      print *,'after call EXCH,me=',me

!$omp parallel do private(i,j,ip1)
      do j = jsta, jend_m
        do i = ista, iend_m
          ip1 = i + 1
!          if (ip1 > im) ip1 = ip1 - im
          if(maptype==207)then
            dx(i,j) = erad*dxval*dtr/gdsdegr
          else
            dx(i,j) = erad*cos(gdlat(i,j)*dtr) *(gdlon(ip1,j)-gdlon(i,j))*dtr
          endif
          if(maptype==207)then
            dy(i,j)= erad*dyval*dtr/gdsdegr
          else
            dy(i,j) = erad*(gdlat(i,j+1)-gdlat(i,j))*dtr  ! like A*DPH
          endif
!     F(I,J)=1.454441e-4*sin(gdlat(i,j)*DTR)         ! 2*omeg*sin(phi)
!     if (i == ii .and. j == jj) print*,'sample LATLON, DY, DY='    &
!           ,i,j,GDLAT(I,J),GDLON(I,J),DX(I,J),DY(I,J)
        end do
      end do
      if(debugprint)print*,'me sample dx dy= ' &
     ,me,dx(isa,jsa),dy(isa,jsa)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          f(i,j) = 1.454441e-4*sin(gdlat(i,j)*dtr)   ! 2*omeg*sin(phi)
        end do
      end do
      
      iyear = idate(1)
      imn   = idate(2)  
      iday  = idate(3) 
      ihrst = idate(4)
      imin  = idate(5)
      jdate = 0
      idate = 0 
!
      print*,'start yr mo day hr min =',iyear,imn,iday,ihrst,imin
      print*,'processing yr mo day hr min='                            &
             ,idat(3),idat(1),idat(2),idat(4),idat(5)
!
      idate(1) = iyear
      idate(2) = imn
      idate(3) = iday
      idate(5) = ihrst
      idate(6) = imin
      sdat(1)  = imn
      sdat(2)  = iday
      sdat(3)  = iyear
      jdate(1) = idat(3)
      jdate(2) = idat(1)
      jdate(3) = idat(2)
      jdate(5) = idat(4)
      jdate(6) = idat(5)
!
      print *,' idate=',idate
      print *,' jdate=',jdate
!
      CALL w3difdat(jdate,idate,0,rinc)
!
      print *,' rinc=',rinc
      ifhr = nint(rinc(2)+rinc(1)*24.)
      print *,' ifhr=',ifhr
      ifmin = nint(rinc(3))
!      if(ifhr /= nint(fhour))print*,'find wrong Grib file';stop
      print*,' in INITPOST ifhr ifmin fileName=',ifhr,ifmin,filename
      
! Getting tstart
      tstart = 0.
      print*,'tstart= ',tstart
      
! Getiing restart
      
      restrt = .true.  ! set RESTRT as default
            
      IF(tstart > 1.0e-2)THEN
        ifhr    = ifhr+nint(tstart)
        rinc    = 0
        idate   = 0
        rinc(2) = -1.0*ifhr
        call w3movdat(rinc,jdate,idate)
        sdat(1) = idate(2)
        sdat(2) = idate(3)
        sdat(3) = idate(1)
        ihrst   = idate(5)       
        print*,'new forecast hours for restrt run= ',ifhr
        print*,'new start yr mo day hr min =',sdat(3),sdat(1)           &
               ,sdat(2),ihrst,imin
      END IF 
      
! GFS does not need DT to compute accumulated fields, set it to one
!      VarName='DT'
      dt   = 1

      hbm2 = 1.0

! start reading 3d netcdf output
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(1),uh(ista_2l,jsta_2l,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(2),vh(ista_2l,jsta_2l,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(3),q(ista_2l,jsta_2l,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(4),t(ista_2l,jsta_2l,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(5),o3(ista_2l,jsta_2l,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(7),wh(ista_2l,jsta_2l,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(8),qqw(ista_2l,jsta_2l,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(9),dpres(ista_2l,jsta_2l,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(10),buf3d(ista_2l,jsta_2l,1),lm)
       do l=1,lm
       do j=jsta,jend
         do i=ista,iend
            cwm(i,j,l)=spval
! dong add missing value
           if (wh(i,j,l) < spval) then
            omga(i,j,l)=(-1.)*wh(i,j,l)*dpres(i,j,l)/abs(buf3d(i,j,l))
           else
            omga(i,j,l) = spval
           end if
!           if(t(i,j,l)>1000.)print*,'bad T ',t(i,j,l)
         enddo
       enddo
       enddo
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(11),qqi(ista_2l,jsta_2l,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(12),qqr(ista_2l,jsta_2l,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(13),qqs(ista_2l,jsta_2l,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(14),qqg(ista_2l,jsta_2l,1),lm)
! read for regional FV3
       if (modelname == 'FV3R') then
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(15),smoke(ista_2l,jsta_2l,1,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(16),fv3dust(ista_2l,jsta_2l,1,1),lm)
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(17),coarsepm(ista_2l,jsta_2l,1,1),lm)
       call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,recname(18),ext550(ista_2l,jsta_2l,1),lm)
       endif

! Compute max QRAIN in the column to be used later in precip type computation
       do j = jsta_2l, jend_2u
        do i = 1, im
           qrmax(i,j)=0.
        end do
       end do

! calculate CWM from FV3 output
       do l=1,lm
       do j=jsta,jend
         do i=ista,iend
            qrmax(i,j)=max(qrmax(i,j),qqr(i,j,l))
            cwm(i,j,l)=qqg(i,j,l)+qqs(i,j,l)+qqr(i,j,l)+qqi(i,j,l)+qqw(i,j,l)
         enddo
       enddo
       if(debugprint)print*,'sample l,t,q,u,v,w= ',isa,jsa,l &
       ,t(isa,jsa,l),q(isa,jsa,l),uh(isa,jsa,l),vh(isa,jsa,l) &
       ,wh(isa,jsa,l)
       if(debugprint)print*,'sample l cwm for FV3',l, &
          cwm(isa,jsa,l)
      end do 

! instantaneous 3D cloud fraction
      if ( imp_physics==11) then !GFDL MP
        varname='cld_amt'
        call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
        spval,varname,cfr(ista_2l,jsta_2l,1),lm)
      else

        iret_bl = nf90_inq_varid(ncid2d,'cldfra_bl',varid_bl)
        iret = nf90_inq_varid(ncid2d,'cldfra',varid)

        if(iret_bl==nf90_noerr .and. iret==nf90_noerr) then
          write(*,*) 'WARNING: BOTH cldfra_bl AND cldfra ARE AVAILABLE. USING cldfra.'
          varname='cldfra'
        else if(iret_bl==nf90_noerr) then
          varname='cldfra_bl'
        else if(iret==nf90_noerr) then
          varname='cldfra'
        else
          varname='nope'
        endif
          
        if(varname /= 'nope') then
          call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
               spval,varname,cfr(ista_2l,jsta_2l,1),lm)
        endif
      endif
!      do l=1,lm
!       if(debugprint)print*,'sample ',VarName,'isa,jsa,l =' &
!          ,cfr(isa,jsa,l),isa,jsa,l
!      enddo

!     WL add cieffr for Thompson scheme cloud ice effective radius
      varname='cieffr'
      call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,effri(ista_2l,jsta_2l,1),lm)

!     WL add cleffr for Thompson scheme cloud water effective radius
      varname='cleffr'
      call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,effrl(ista_2l,jsta_2l,1),lm)

!     WL add cseffr for Thompson scheme snow effective radius
      varname='cseffr'
      call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,effrs(ista_2l,jsta_2l,1),lm)

!=====================================
! For AQF Hourly average field PM2.5
!=====================================

      if (aqf_on) then

       ! *********** VarName need to be in lower case ************
       ! === It will cause problem if not use the lower case =====
       ! *********************************************************

       !-- rename input o3_ave and pm25_ave to NCO grib2 name OZCON and PMTF

       varname='o3_ave'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,varname,avgozcon(ista_2l,jsta_2l,1),lm)

       varname='pm25_ave'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,varname,avgpmtf(ista_2l,jsta_2l,1),lm)

       varname='aod'
       call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
       spval,varname,aqm_aod550(ista_2l,jsta_2l))

      endif     ! -- aqf_on
!============================

! read for regional FV3
      if (modelname == 'FV3R') then
! max hourly updraft velocity
      varname='upvvelmax'
      call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,w_up_max(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' = ',w_up_max(isa,jsa)
! max hourly downdraft velocity
      varname='dnvvelmax'
      call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,w_dn_max(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' = ',w_dn_max(isa,jsa)
! max hourly updraft helicity
      varname='uhmax25'
      call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,up_heli_max(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' = ',up_heli_max(isa,jsa)
! min hourly updraft helicity
      varname='uhmin25'
      call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,up_heli_min(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' = ',up_heli_min(isa,jsa)
! max hourly 0-3km updraft helicity
      varname='uhmax03'
      call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,up_heli_max03(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' = ',up_heli_max03(isa,jsa)
! min hourly 0-3km updraft helicity
      varname='uhmin03'
      call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,up_heli_min03(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' = ',up_heli_min03(isa,jsa)

! max 0-1km relative vorticity max 
      varname='maxvort01'
      call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,rel_vort_max01(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' = ',rel_vort_max01(isa,jsa)
! max 0-2km relative vorticity max
      varname='maxvort02'
      call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,rel_vort_max(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' =',rel_vort_max(isa,jsa)
! max hybrid lev 1 relative vorticity max
      varname='maxvorthy1'
      call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,rel_vort_maxhy1(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' =',rel_vort_maxhy1(isa,jsa)
! biomass burning emissions
      varname='ebb_smoke_hr'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ebb(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' =',ebb(isa,jsa)
! hourly wildfire potential
      varname='hwp'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,hwp(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' =',hwp(isa,jsa)
      endif

! lightning threat index 1
      varname='ltg1_max'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ltg1_max(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' =',ltg1_max(isa,jsa)

! lightning threat index 2
      varname='ltg2_max'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ltg2_max(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' =',ltg2_max(isa,jsa)

! lightning threat index 3
      varname='ltg3_max'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ltg3_max(ista_2l,jsta_2l))
     if(debugprint)print*,'sample ',varname,' =',ltg3_max(isa,jsa)

! surface pressure
      varname='pressfc'
      call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,pint(ista_2l,jsta_2l,lp1))
      do j=jsta,jend
        do i=ista,iend
!          if(pint(i,j,lp1)>1.0E6 .or. pint(ista_2l,jsta_2l,lp1)<50000.) &
!           print*,'bad psfc ',i,j,pint(i,j,lp1)
        end do
      end do
      if(debugprint)print*,'sample ',varname,' =',pint(isa,jsa,lp1)

      pt = ak5(1)

      do j=jsta,jend
        do i=ista,iend
          pint(i,j,1)= pt
        end do
      end do

      do l=2,lp1
        do j=jsta,jend
          do i=ista,iend
            if (dpres(i,j,l-1)<spval .and. pint(i,j,l-1)<spval) then
              pint(i,j,l)= pint(i,j,l-1) + dpres(i,j,l-1)
            else
              pint(i,j,l)=spval
            endif
          enddo
        enddo
!        if (me == 0) print*,'sample model pint,pmid' ,ii,jj,l &
!          ,pint(ii,jj,l),pmid(ii,jj,l)
      end do

!compute pmid from averaged two layer pint
      do l=lm,1,-1
        do j=jsta,jend
          do i=ista,iend
            if (pint(i,j,l)<spval .and. pint(i,j,l+1)<spval) then
              pmid(i,j,l)=0.5*(pint(i,j,l)+pint(i,j,l+1))
            else
              pmid(i,j,l)=spval
            endif
          enddo
        enddo
      enddo

! surface height from FV3 
! dong set missing value for zint
!      zint=spval
      varname='hgtsfc'
      call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,zint(ista_2l,jsta_2l,lp1))
      if(debugprint)print*,'sample ',varname,' =',zint(isa,jsa,lp1)
      do j=jsta,jend
        do i=ista,iend
          if (zint(i,j,lp1) /= spval) then
            fis(i,j)      = zint(i,j,lp1) * grav
          else
            fis(i,j)      = spval
          endif
        enddo
      enddo

      do l=lm,1,-1
        do j=jsta,jend
          do i=ista,iend
            if(zint(i,j,l+1)/=spval .and. buf3d(i,j,l)/=spval)then
!make sure delz is positive
             zint(i,j,l)=zint(i,j,l+1)+abs(buf3d(i,j,l))
!             if(zint(i,j,l)>1.0E6)print*,'bad H ',i,j,l,zint(i,j,l)
            else
             zint(i,j,l)=spval
            end if
          end do
        end do
        if(debugprint)print*,'sample zint= ',isa,jsa,l,zint(isa,jsa,l)
      end do

      do l=lp1,1,-1
        do j=jsta,jend
          do i=ista,iend
            alpint(i,j,l)=log(pint(i,j,l))
          end do
        end do
      end do

      do l=lm,1,-1
        do j=jsta,jend
          do i=ista,iend
            if(zint(i,j,l+1)/=spval .and. zint(i,j,l)/=spval &
            .and. pmid(i,j,l)/=spval)then
             zmid(i,j,l)=zint(i,j,l+1)+(zint(i,j,l)-zint(i,j,l+1))* &
                    (log(pmid(i,j,l))-alpint(i,j,l+1))/ &
                    (alpint(i,j,l)-alpint(i,j,l+1))
             if(zmid(i,j,l)>1.0e6)print*,'bad Hmid ',i,j,l,zmid(i,j,l)
            else
             zmid(i,j,l)=spval
            endif
          end do
        end do
      end do

      
      print *, 'gocart_on=',gocart_on
      print *, 'gccpp_on=',gccpp_on
      print *, 'nasa_on=',nasa_on
      if (gocart_on .or.gccpp_on .or. nasa_on) then

! GFS output dust in nemsio (GOCART)
        dustcb=0.0
        dustallcb=0.0
!       DUST = SPVAL

        varname='dust1'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt1(ista_2l,jsta_2l,1),lm)
        varname='dust2'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt2(ista_2l,jsta_2l,1),lm)
        varname='dust3'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt3(ista_2l,jsta_2l,1),lm)
        varname='dust4'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt4(ista_2l,jsta_2l,1),lm)
        varname='dust5'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt5(ista_2l,jsta_2l,1),lm)

      
        do l=1,lm

!$omp parallel do private(i,j)
          do j=jsta,jend
          do i=ista,iend
          dust(i,j,l,1)=dt1(i,j,l)
          dust(i,j,l,2)=dt2(i,j,l)
          dust(i,j,l,3)=dt3(i,j,l)
          dust(i,j,l,4)=dt4(i,j,l)
          dust(i,j,l,5)=dt5(i,j,l)
          

           dustcb(i,j)=dustcb(i,j)+&
           (dust(i,j,l,1)+0.38*dust(i,j,l,2))* &
           dpres(i,j,l)/grav


           dustallcb(i,j)=dustallcb(i,j)+ &
           (dust(i,j,l,1)+dust(i,j,l,2)+ &
           dust(i,j,l,3)+0.74*dust(i,j,l,4))* &
           dpres(i,j,l)/grav
           enddo
           enddo
        end do ! do loop for l

 
! GFS output sea salt in nemsio (GOCART)
        sscb=0.0
        ssallcb=0.0
        
        varname='seas1'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt1(ista_2l,jsta_2l,1),lm)

        varname='seas2'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt2(ista_2l,jsta_2l,1),lm)

        varname='seas3'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt3(ista_2l,jsta_2l,1),lm)

        varname='seas4'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt4(ista_2l,jsta_2l,1),lm)

        varname='seas5'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt5(ista_2l,jsta_2l,1),lm)



        do l=1,lm
!$omp parallel do private(i,j)
          do j=jsta,jend
          do i=ista,iend
          salt(i,j,l,1)=dt1(i,j,l)
          salt(i,j,l,2)=dt2(i,j,l)
          salt(i,j,l,3)=dt3(i,j,l)
          salt(i,j,l,4)=dt4(i,j,l)
          salt(i,j,l,5)=dt5(i,j,l)

            sscb(i,j)=sscb(i,j)+ &
         (salt(i,j,l,1)+salt(i,j,l,2)+0.83*salt(i,j,l,3))*  &
           dpres(i,j,l)/grav


          ssallcb(i,j)=ssallcb(i,j)+ &
         (salt(i,j,l,1)+salt(i,j,l,2)+salt(i,j,l,3)+salt(i,j,l,4))* &
           dpres(i,j,l)/grav
           enddo
           enddo
        end do ! do loop for l
! GFS output black carbon in nemsio (GOCART)
        bccb=0.0


        varname='bc1'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt1(ista_2l,jsta_2l,1),lm)

        varname='bc2'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt2(ista_2l,jsta_2l,1),lm)

        do l=1,lm
!$omp parallel do private(i,j)
          do j=jsta,jend
          do i=ista,iend

          soot(i,j,l,1)=dt1(i,j,l)
          soot(i,j,l,2)=dt2(i,j,l)

            bccb(i,j)=bccb(i,j)+ &
        (soot(i,j,l,1)+soot(i,j,l,2))* &
           dpres(i,j,l)/grav
           enddo
           enddo
        end do ! do loop for l

        occb=0.0
! GFS output organic carbon in nemsio (GOCART)

        varname='oc1'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt1(ista_2l,jsta_2l,1),lm)

        varname='oc2'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt2(ista_2l,jsta_2l,1),lm)
        do l=1,lm
!$omp parallel do private(i,j)
          do j=jsta,jend
          do i=ista,iend
          waso(i,j,l,1)=dt1(i,j,l)
          waso(i,j,l,2)=dt2(i,j,l)

            occb(i,j)=occb(i,j)+ &
        (waso(i,j,l,1)+waso(i,j,l,2))* &
           dpres(i,j,l)/grav
           enddo
           enddo
        end do ! do loop for l

! GFS output sulfate in netcdf (GOCART)
        sulfcb=0.0

!       SUSO = SPVAL
        if (gocart_on .or. gccpp_on) then
        varname='sulf'
        endif

        if (nasa_on) then
        varname='so4'
        endif
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt1(ista_2l,jsta_2l,1),lm)
        
         do l=1,lm
!$omp parallel do private(i,j)
          do j=jsta,jend
          do i=ista,iend
          suso(i,j,l,1)=dt1(i,j,l)

            sulfcb(i,j)=sulfcb(i,j)+ &
        suso(i,j,l,1)* &
           dpres(i,j,l)/grav
           enddo
           enddo
        end do ! do loop for l

        if (nasa_on) then
! GFS output nitrate in netcdf (GOCART)
        no3cb=0.0
        varname='no3an1'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt1(ista_2l,jsta_2l,1),lm)

        varname='no3an2'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt2(ista_2l,jsta_2l,1),lm)

        varname='no3an3'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt3(ista_2l,jsta_2l,1),lm)

         do l=1,lm
!$omp parallel do private(i,j)
          do j=jsta,jend
          do i=ista,iend
          no3(i,j,l,1)=dt1(i,j,l)
          no3(i,j,l,2)=dt2(i,j,l)
          no3(i,j,l,3)=dt3(i,j,l)

            no3cb(i,j)=no3cb(i,j)+ &
        (no3(i,j,l,1)+no3(i,j,l,2)+no3(i,j,l,3))* &
           dpres(i,j,l)/grav
           enddo
           enddo
        end do ! do loop for l

! GFS output NH4 in netcdf (GOCART)
        nh4cb=0.0
        varname='nh4a'
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt1(ista_2l,jsta_2l,1),lm)

         do l=1,lm
!$omp parallel do private(i,j)
          do j=jsta,jend
          do i=ista,iend
          nh4(i,j,l,1)=dt1(i,j,l)

            nh4cb(i,j)=nh4cb(i,j)+ &
        nh4(i,j,l,1)* &
           dpres(i,j,l)/grav
           enddo
           enddo
        end do ! do loop for l


        endif !nasa_on

! GFS output pp25 in nemsio (GOCART)
        pp25cb=0.0

        if (gocart_on .or. gccpp_on) then
        varname='pp25'
        endif

        if (nasa_on) then
        varname='pm25'
        endif
       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt1(ista_2l,jsta_2l,1),lm)


! GFS output pp10 in nemsio (GOCART)
        pp10cb=0.0
        if (gocart_on .or. gccpp_on) then
        varname='pp10'
        endif

        if (nasa_on) then
        varname='pm10'
        endif

       call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
       spval,varname,dt2(ista_2l,jsta_2l,1),lm)

        do l=1,lm
!$omp parallel do private(i,j)
          do j=jsta,jend
          do i=ista,iend
          pp25(i,j,l,1)=dt1(i,j,l)
          pp10(i,j,l,1)=dt2(i,j,l)


            pp25cb(i,j)=pp25cb(i,j)+ &
        pp25(i,j,l,1)* dpres(i,j,l)/grav

            pp10cb(i,j)=pp10cb(i,j)+ &
        pp10(i,j,l,1)* dpres(i,j,l)/grav
           enddo
           enddo
        end do ! do loop for l
! -- compute air density RHOMID and remove negative tracer values
        do l=1,lm
!$omp parallel do private(i,j,tv)
          do j=jsta,jend
            do i=ista,iend

              tv = t(i,j,l) * (h1+d608*max(q(i,j,l),qmin))
              rhomid(i,j,l) = pmid(i,j,l) / (rd*tv)
              do n = 1,  nbin_du
                IF ( dust(i,j,l,n) < spval) THEN
                  dust(i,j,l,n) = max(dust(i,j,l,n), 0.0)
                ENDIF
              enddo
              do n = 1,  nbin_ss
                IF ( salt(i,j,l,n) < spval) THEN
                  salt(i,j,l,n) = max(salt(i,j,l,n), 0.0)
                ENDIF
              enddo
              do n = 1,  nbin_oc
                IF ( waso(i,j,l,n) < spval) THEN
                  waso(i,j,l,n) = max(waso(i,j,l,n), 0.0)
                ENDIF
              enddo
              do n = 1,  nbin_bc
                IF ( soot(i,j,l,n) < spval) THEN
                  soot(i,j,l,n) = max(soot(i,j,l,n), 0.0)
                ENDIF
              enddo
              do n = 1,  nbin_su
                IF ( suso(i,j,l,n) < spval) THEN
                  suso(i,j,l,n) = max(suso(i,j,l,n), 0.0)
                ENDIF
              enddo
          if (nasa_on) then
              do n = 1,  nbin_no3
                IF ( no3(i,j,l,n) < spval) THEN
                  no3(i,j,l,n) = max(no3(i,j,l,n), 0.0)
                ENDIF
              enddo
              do n = 1,  nbin_nh4
                IF ( nh4(i,j,l,n) < spval) THEN
                  nh4(i,j,l,n) = max(nh4(i,j,l,n), 0.0)
                ENDIF
              enddo
          endif !nasa_on
            end do
          end do
        end do
             l=lm
!$omp parallel do private(i,j)
          do j=jsta,jend
            do i=ista,iend
            dustcb(i,j) = max(dustcb(i,j), 0.0)
            dustallcb(i,j) = max(dustallcb(i,j), 0.0)
            sscb(i,j) = max(sscb(i,j), 0.0)
            ssallcb(i,j) = max(ssallcb(i,j), 0.0)
            bccb(i,j) = max(bccb(i,j), 0.0)
            occb(i,j) = max(occb(i,j), 0.0)
            sulfcb(i,j) = max(sulfcb(i,j), 0.0)
            if (nasa_on) then
            no3cb(i,j) = max(no3cb(i,j), 0.0)
            nh4cb(i,j) = max(nh4cb(i,j), 0.0)
            endif
            pp25cb(i,j) = max(pp25cb(i,j), 0.0)
            pp10cb(i,j) = max(pp10cb(i,j), 0.0)

!      PM25 dust and seasalt      
       dustpm(i,j)=(dust(i,j,l,1)+0.38*dust(i,j,l,2))*rhomid(i,j,l) !ug/m3
       dustpm10(i,j)=(dust(i,j,l,1)+dust(i,j,l,2)+dust(i,j,l,3)+ &
       0.74*dust(i,j,l,4))*rhomid(i,j,l) !ug/m3
       sspm(i,j)=(salt(i,j,l,1)+salt(i,j,l,2)+ &
       0.83*salt(i,j,l,3))*rhomid(i,j,l)  !ug/m3 

       if (gocart_on .or. gccpp_on) then
!      PM10 concentration
       dusmass(i,j)=(dust(i,j,l,1)+dust(i,j,l,2)+dust(i,j,l,3)+ &
       0.74*dust(i,j,l,4)+salt(i,j,l,1)+salt(i,j,l,2)+salt(i,j,l,3)+ &
       salt(i,j,l,4) + soot(i,j,l,1)+soot(i,j,l,2)+waso(i,j,l,1)+ &
       waso(i,j,l,2) +suso(i,j,l,1)+pp25(i,j,l,1)+pp10(i,j,l,1)) &
       *rhomid(i,j,l)  !ug/m3
!      PM25 concentration       
       dusmass25(i,j)=(dust(i,j,l,1)+0.38*dust(i,j,l,2)+ &
       salt(i,j,l,1)+salt(i,j,l,2)+0.83*salt(i,j,l,3) + &
       soot(i,j,l,1)+soot(i,j,l,2)+waso(i,j,l,1)+ &
       waso(i,j,l,2) +suso(i,j,l,1)+pp25(i,j,l,1))*rhomid(i,j,l)  !ug/m3

!      PM10 column
        ducmass(i,j)=dustallcb(i,j)+ssallcb(i,j)+bccb(i,j)+ &
         occb(i,j)+sulfcb(i,j)+pp25cb(i,j)+pp10cb(i,j)
!      PM25 column
        ducmass25(i,j)=dustcb(i,j)+sscb(i,j)+bccb(i,j)+occb(i,j) &
         +sulfcb(i,j)+pp25cb(i,j)
       endif !gocart_on

       if (nasa_on) then
!      PM10 concentration
       dusmass(i,j)=pp10(i,j,l,1)*rhomid(i,j,l)  !ug/m3
!      PM25 concentration       
       dusmass25(i,j)=pp25(i,j,l,1)*rhomid(i,j,l)  !ug/m3

!      PM10 column
        ducmass(i,j)=pp10cb(i,j)
!      PM25 column
        ducmass25(i,j)=pp25cb(i,j)
       endif !nasa_on


            end do
          end do
       

      endif                     ! endif for gocart_on & nasa_on

!         ',ll,waso(isa,jsa,ll,2)
!         ',ll,waso(isa,jsa,ll,2)


!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!

! done with 3d file, close it for now
      status=nf90_close(ncid3d)
      deallocate(recname)

! IVEGSRC=1 for IGBP, 0 for USGS, 2 for UMD
      varname='IVEGSRC'
      status=nf90_get_att(ncid2d,nf90_global,'IVEGSRC',ivegsrc)
      if (status /= 0) then
       print*,varname,' not found-Assigned 1 for IGBP as default'
       ivegsrc=1
      end if
      if (me == 0) print*,'IVEGSRC= ',ivegsrc

! set novegtype based on vegetation classification
      if(ivegsrc==2)then
       novegtype=13
      else if(ivegsrc==1)then
       novegtype=20
      else if(ivegsrc==0)then
       novegtype=24
      end if
      if (me == 0) print*,'novegtype= ',novegtype

      status=nf90_get_att(ncid2d,nf90_global,'fhzero',fhzero)
      if (status /= 0) then
       print*,'fhzero not found-Assigned 3 hours as default'
       fhzero=3
      end if
      if (me == 0) print*,'fhzero= ',fhzero
!
      status=nf90_get_att(ncid2d,nf90_global,'dtp',dtp)
      if (status /= 0) then
       print*,'dtp not found-Assigned 90s as default'
       dtp=90
      end if
      if (me == 0) print*,'dtp= ',dtp
! Initializes constants for Ferrier microphysics
      if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95) then
        CALL microinit(imp_physics)
      end if

        tprec   = float(fhzero)
        if(ifhr>240)tprec=12.
        tclod   = tprec
        trdlw   = tprec
        trdsw   = tprec
        tsrfc   = tprec
        tmaxmin = tprec
        td3d    = tprec
        print*,'tprec = ',tprec


      varname='refl_10cm'
      call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ref_10cm(ista_2l,jsta_2l,1),lm)
!      do l=1,lm
!       if(debugprint)print*,'sample ',VarName,'isa,jsa,l =' &
!          ,REF_10CM(isa,jsa,l),isa,jsa,l
!      enddo

! turbulence kinetic energy (QKE = 2*TKE)
      varname='qke'
      call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,q2(ista_2l,jsta_2l,1),lm)
      do l=1,lm
      do j=jsta,jend
      do i=ista,iend
        q2(i,j,l)=q2(i,j,l)/2.0
      enddo
      enddo
      enddo

! ice-friendly aerosol number concentration
      varname='nifa'
      call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,qqnifa(ista_2l,jsta_2l,1),lm)

! water-friendly aerosol number concentration
      varname='nwfa'
      call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,qqnwfa(ista_2l,jsta_2l,1),lm)

      varname='land' 
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sm)
      if(debugprint)print*,'sample ',varname,' =',sm((ista+iend)/2,(jsta+jend)/2)

!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= spval) sm(i,j) = 1.0 - sm(i,j)
        enddo
      enddo

! sea ice mask 

      varname    = 'icec'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sice)
     if(debugprint)print*,'sample ',varname,' = ',sice(isa,jsa)

!      where(sice /=spval .and. sice >=1.0)sm=0.0 !sea ice has sea
!      mask=0
! GFS flux files have land points with non-zero sea ice, per Iredell,
! these
! points have sea ice changed to zero, i.e., trust land mask more than
! sea ice
!     where(sm/=spval .and. sm==0.0)sice=0.0 !specify sea ice=0 at land

!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= spval .and. sm(i,j) == 0.0) sice(i,j) = 0.0
        enddo
      enddo


! PBL height using nemsio
      varname    = 'hpbl'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,pblh)
      if(debugprint)print*,'sample ',varname,' = ',pblh(isa,jsa)

! frictional velocity using nemsio
      varname='fricv'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ustar)
!     if(debugprint)print*,'sample ',VarName,' = ',ustar(isa,jsa)

! roughness length using getgb
      varname='sfcr'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,z0)
!     if(debugprint)print*,'sample ',VarName,' = ',z0(isa,jsa)

! sfc exchange coeff
      varname='sfexc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sfcexc)

! accumulated snowfall
      varname='tsnowp'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,snow_acm)
! snowfall bucket
      varname='tsnowpb'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,snow_bkt)

! accumulated graupel/sleet
      varname='frozr'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,acgraup)

! graupel/sleet bucket
      varname='frozrb'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,graup_bucket)

! accumulated freezing rain
      varname='frzr'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,acfrain)

! freezing rain bucket
      varname='frzrb'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,frzrn_bucket)

! time step snow (in m)
      varname='snow'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,snownc)

! time step graupel (in m)
      varname='graupel'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,graupelnc)

! aerodynamic conductance
      varname='acond'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,acond)
      if(debugprint)print*,'sample ',varname,' = ',acond(isa,jsa)
! mid day avg albedo
      varname='albdo_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgalbedo)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (avgalbedo(i,j) /= spval) avgalbedo(i,j) = avgalbedo(i,j) * 0.01
        enddo
      enddo
      if(debugprint)print*,'sample ',varname,' = ',avgalbedo(isa,jsa)

! surface potential T  using getgb
      varname='tmpsfc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ths)

!     where(ths/=spval)ths=ths*(p1000/pint(:,:,lp1))**CAPA ! convert to THS

!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (ths(i,j) /= spval) then
!    write(*,*)' i=',i,' j=',j,' ths=',ths(i,j),' pint=',pint(i,j,lp1)
            ths(i,j) = ths(i,j) * (p1000/pint(i,j,lp1))**capa
          endif
          qs(i,j)    = spval ! GFS does not have surface specific humidity
          twbs(i,j)  = spval ! GFS does not have inst sensible heat flux
          qwbs(i,j)  = spval ! GFS does not have inst latent heat flux
!assign sst
          if (sm(i,j) /= 0.0 .and. ths(i,j) < spval ) then
            if (sice(i,j) >= 0.15) then
              sst(i,j) = 271.4
            else
              sst(i,j) = ths(i,j) * (pint(i,j,lp1)/p1000)**capa
            endif
          else
              sst(i,j) = spval
          endif
        enddo
      enddo
     if(debugprint)print*,'sample ',varname,' = ',ths(isa,jsa)

! foundation temperature
      varname='tref'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,fdnsst)
      if(debugprint)print*,'sample ',varname,' = ',fdnsst(isa,jsa)

          
!  GFS does not have time step and physics time step, make up ones since they
! are not really used anyway
!      NPHS=1.
!      DT=90.
!      DTQ2 = DT * NPHS  !MEB need to get physics DT
      dtq2 = dtp !MEB need to get physics DT
      nphs=1
      dt = dtq2/nphs   !MEB need to get DT
      tsph = 3600./dt

! convective precip in m per physics time step using getgb
! read 3 hour bucket
      varname='cpratb_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgcprate)
!     where(avgcprate /= spval)avgcprate=avgcprate*dtq2/1000. ! convert to m
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (avgcprate(i,j) /= spval) avgcprate(i,j) = avgcprate(i,j) * (dtq2*0.001)
        enddo
      enddo
!     if(debugprint)print*,'sample ',VarName,' = ',avgcprate(isa,jsa)
      
!      print*,'maxval CPRATE: ', maxval(CPRATE)

! read continuous bucket
      varname='cprat_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgcprate_cont)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (avgcprate_cont(i,j) /= spval) avgcprate_cont(i,j) = &
            avgcprate_cont(i,j) * (dtq2*0.001)
        enddo
      enddo
!     if(debugprint)print*,'sample ',VarName,' = ',avgcprate_cont(isa,jsa)

!      print*,'maxval CPRATE: ', maxval(CPRATE)

! precip rate in m per physics time step using getgb
      varname='prateb_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgprec)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if(avgprec(i,j) /= spval)avgprec(i,j)=avgprec(i,j)*(dtq2*0.001)
        enddo
      enddo

     if(debugprint)print*,'sample ',varname,' = ',avgprec(isa,jsa)
      
!      prec = avgprec !set avg cprate to inst one to derive other fields

      varname='prate_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgprec_cont)
!     where(avgprec /= spval)avgprec=avgprec*dtq2/1000. ! convert to m
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (avgprec_cont(i,j) /=spval)avgprec_cont(i,j)=avgprec_cont(i,j) &
               * (dtq2*0.001)
        enddo
      enddo

     if(debugprint)print*,'sample ',varname,' = ',avgprec_cont(isa,jsa)
! precip rate in m per physics time step
      varname='tprcp'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,prec)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (prec(i,j) /= spval) prec(i,j)=prec(i,j)* (dtq2*0.001) &
              * 1000. / dtp
        enddo
      enddo

! convective precip rate in m per physics time step
      varname='cnvprcp'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,cprate)
!set cprate as 0.
      do j=jsta,jend
        do i=ista,iend
           if (cprate(i,j) /= spval) then
            cprate(i,j) = max(0.,cprate(i,j)) * (dtq2*0.001) &
                 * 1000. / dtp
          else
            cprate(i,j) = 0.
          endif
        enddo
      enddo

! GFS does not have accumulated total, gridscale, and convective precip, will use inst precip to derive in SURFCE.f

! max hourly surface precipitation rate
      varname='pratemax'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,prate_max)
     if(debugprint)print*,'sample ',varname,' = ',prate_max(isa,jsa)
! max hourly 1-km agl reflectivity
      varname='refdmax'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,refd_max)
     if(debugprint)print*,'sample ',varname,' = ',refd_max(isa,jsa)
! max hourly -10C reflectivity
      varname='refdmax263k'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,refdm10c_max)
     if(debugprint)print*,'sample ',varname,' = ',refdm10c_max(isa,jsa)

! max hourly u comp of 10m agl wind
      varname='u10max'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,u10max)
     if(debugprint)print*,'sample ',varname,' = ',u10max(isa,jsa)
! max hourly v comp of 10m agl wind
      varname='v10max'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,v10max)
     if(debugprint)print*,'sample ',varname,' = ',v10max(isa,jsa)
! max hourly 10m agl wind speed
      varname='spd10max'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,wspd10max)
     if(debugprint)print*,'sample ',varname,' = ',wspd10max(isa,jsa)

! inst snow water eqivalent using nemsio
      varname='weasd'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sno)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j)==0.) sno(i,j) = spval
        enddo
      enddo
     if(debugprint)print*,'sample ',varname,' = ',sno(isa,jsa)

! ave snow cover 
      varname='snowc_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,snoavg)
! snow cover is multipled by 100 in SURFCE before writing it out
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j)==1.0 .and. sice(i,j)==0.) snoavg(i,j)=spval
          if(snoavg(i,j)/=spval)snoavg(i,j)=snoavg(i,j)/100.
        end do
      end do

! snow depth in mm using nemsio
      varname='snod'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,si)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j)==1.0 .and. sice(i,j)==0.) si(i,j)=spval
          if (si(i,j) /= spval) si(i,j) = si(i,j) * 1000.0
          cldefi(i,j) = spval ! GFS does not have convective cloud efficiency
          lspa(i,j)   = spval ! GFS does not have similated precip
          th10(i,j)   = spval ! GFS does not have 10 m theta
          th10(i,j)   = spval ! GFS does not have 10 m theta
          q10(i,j)    = spval ! GFS does not have 10 m humidity
          albase(i,j) = spval ! GFS does not have snow free albedo
        enddo
      enddo
     if(debugprint)print*,'sample ',varname,' = ',si(isa,jsa)

! 2m T using nemsio
      varname='tmp2m'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,tshltr)
     if(debugprint)print*,'sample ',varname,' = ',tshltr(isa,jsa)

! GFS does not have 2m pres, estimate it, also convert t to theta 
      Do j=jsta,jend
        Do i=ista,iend
          pshltr(i,j)=pint(i,j,lm+1)*exp(-0.068283/tshltr(i,j))
          tshltr(i,j)= tshltr(i,j)*(p1000/pshltr(i,j))**capa ! convert to theta
!          if (j == jm/2 .and. mod(i,50) == 0)
!     +   print*,'sample 2m T and P after scatter= '
!     +   ,i,j,tshltr(i,j),pshltr(i,j)
        end do
      end do

! 2m specific humidity using nemsio
      varname='spfh2m'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,qshltr)
     if(debugprint)print*,'sample ',varname,' = ',qshltr(isa,jsa)
      
! time averaged column cloud fractionusing nemsio
      varname='tcdc_aveclm'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgtcdc)
!     where(avgtcdc /= spval)avgtcdc=avgtcdc/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (avgtcdc(i,j) /= spval) avgtcdc(i,j) = avgtcdc(i,j) * 0.01
        enddo
      enddo
     if(debugprint)print*,'sample ',varname,' = ',avgtcdc(isa,jsa)

! GFS probably does not use zenith angle
!$omp parallel do private(i,j)
      do j=jsta_2l,jend_2u
        do i=ista_2l,iend_2u
          czen(i,j)   = spval
          czmean(i,j) = spval      
        enddo
      enddo

! maximum snow albedo in fraction using nemsio
      varname='snoalb'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,mxsnal)

! land fraction
      varname='lfrac'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,landfrac)
     
! GFS probably does not use sigt4, set it to sig*t^4
!$omp parallel do private(i,j,tlmh)
      Do j=jsta,jend
        Do i=ista,iend
          tlmh = t(i,j,lm) * t(i,j,lm)
          sigt4(i,j) = 5.67e-8 * tlmh * tlmh
        End do
      End do

! TG is not used, skip it for now

! GFS does not have inst cloud fraction for high, middle, and low cloud
!$omp parallel do private(i,j)
      do j=jsta_2l,jend_2u
        do i=ista_2l,iend_2u
          cfrach(i,j) = spval
          cfracl(i,j) = spval
          cfracm(i,j) = spval
        enddo
      enddo

! ave high cloud fraction using nemsio
      varname='tcdc_avehcl'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgcfrach)
!     where(avgcfrach /= spval)avgcfrach=avgcfrach/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (avgcfrach(i,j) /= spval) avgcfrach(i,j) = avgcfrach(i,j) * 0.01
        enddo
      enddo
     if(debugprint)print*,'sample ',varname,' = ',avgcfrach(isa,jsa)

! ave low cloud fraction using nemsio
      varname='tcdc_avelcl'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgcfracl)
!     where(avgcfracl /= spval)avgcfracl=avgcfracl/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (avgcfracl(i,j) /= spval) avgcfracl(i,j) = avgcfracl(i,j) * 0.01
        enddo
      enddo
     if(debugprint)print*,'sample ',varname,' = ',avgcfracl(isa,jsa)
      
! ave middle cloud fraction using nemsio
      varname='tcdc_avemcl'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgcfracm)
!     where(avgcfracm /= spval)avgcfracm=avgcfracm/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (avgcfracm(i,j) /= spval) avgcfracm(i,j) = avgcfracm(i,j) * 0.01
        enddo
      enddo
     if(debugprint)print*,'sample ',varname,' = ',avgcfracm(isa,jsa)
      
! inst convective cloud fraction using nemsio
      varname='tcdccnvcl'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,cnvcfr)
!     where(cnvcfr /= spval)cnvcfr=cnvcfr/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (cnvcfr(i,j) /= spval) cnvcfr(i,j)= cnvcfr(i,j) * 0.01
        enddo
      enddo
!     if(debugprint)print*,'sample ',VarName,' = ',cnvcfr(isa,jsa)
      
! slope type using nemsio
      varname='sltyp'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,buf)
!$omp parallel do private(i,j)
      do j = jsta_2l, jend_2u
        do i=ista,iend
          if (buf(i,j) < spval) then
             islope(i,j) = nint(buf(i,j))
          else
             islope(i,j) = 0
          endif
        enddo
      enddo
!     if(debugprint)print*,'sample ',VarName,' = ',islope(isa,jsa)

! plant canopy sfc wtr in m 
      varname='cnwat'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,cmc)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (cmc(i,j) /= spval) cmc(i,j) = cmc(i,j) * 0.001
          if (sm(i,j) /= 0.0) cmc(i,j) = spval
        enddo
      enddo
!     if(debugprint)print*,'sample ',VarName,' = ',cmc(isa,jsa)
      
!$omp parallel do private(i,j)
      do j=jsta_2l,jend_2u
        do i=ista_2l,iend_2u
          grnflx(i,j) = spval ! GFS does not have inst ground heat flux
        enddo
      enddo

! frozen precip fraction using nemsio
      varname='cpofp'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sr)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if(sr(i,j) /= spval) then
!set range within (0,1)
            sr(i,j)=min(1.,max(0.,sr(i,j)))
          endif
        enddo
      enddo

! sea ice skin temperature
      varname='tisfc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ti)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sice(i,j) == spval .or. sice(i,j) == 0.) ti(i,j)=spval
        enddo
      enddo

! vegetation fraction in fraction. using nemsio
      varname='veg'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,vegfrc)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (vegfrc(i,j) /= spval) then
            vegfrc(i,j) = vegfrc(i,j) * 0.01
          else
            vegfrc(i,j) = 0.0
          endif
        enddo
      enddo
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) vegfrc(i,j) = spval
        enddo
      enddo
!     if(debugprint)print*,'sample ',VarName,' = ',vegfrc(isa,jsa)
      
! GFS doesn not yet output soil layer thickness, assign SLDPTH to be the same as nam

         sldpth(1) = 0.10
         sldpth(2) = 0.3
         sldpth(3) = 0.6
         sldpth(4) = 1.0

! Eric James, 1 Oct 2021: Because FV3 does not have 1d var "zs", used to
! assign soil depths for RUC LSM, hard wire 9 soil depths here
! so they aren't missing.

       IF (nsoil==9) THEN
         sllevel(1) = 0.0
         sllevel(2) = 0.01
         sllevel(3) = 0.04
         sllevel(4) = 0.1
         sllevel(5) = 0.3
         sllevel(6) = 0.6
         sllevel(7) = 1.0
         sllevel(8) = 1.6
         sllevel(9) = 3.0
       END IF
 
! liquid volumetric soil mpisture in fraction using nemsio
      varname='soill1'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sh2o(ista_2l,jsta_2l,1))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) sh2o(i,j,1) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,sh2o(isa,jsa,1)

      varname='soill2'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sh2o(ista_2l,jsta_2l,2))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) sh2o(i,j,2) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,sh2o(isa,jsa,2)

      varname='soill3'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sh2o(ista_2l,jsta_2l,3))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) sh2o(i,j,3) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,sh2o(isa,jsa,3)

      varname='soill4' 
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sh2o(ista_2l,jsta_2l,4))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) sh2o(i,j,4) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,sh2o(isa,jsa,4)

! volumetric soil moisture using nemsio
      varname='soilw1'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,smc(ista_2l,jsta_2l,1))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) smc(i,j,1) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,1)
      
      varname='soilw2'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,smc(ista_2l,jsta_2l,2))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) smc(i,j,2) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,2)
      
      varname='soilw3'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,smc(ista_2l,jsta_2l,3))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) smc(i,j,3) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,3)
      
      varname='soilw4'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,smc(ista_2l,jsta_2l,4))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) smc(i,j,4) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,4)

      IF (nsoil==9) THEN

      varname='soilw5'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,smc(ista_2l,jsta_2l,5))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) smc(i,j,5) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,5)

      varname='soilw6'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,smc(ista_2l,jsta_2l,6))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) smc(i,j,6) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,6)

      varname='soilw7'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,smc(ista_2l,jsta_2l,7))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) smc(i,j,7) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,7)

      varname='soilw8'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,smc(ista_2l,jsta_2l,8))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) smc(i,j,8) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,8)

      varname='soilw9'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,smc(ista_2l,jsta_2l,9))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) smc(i,j,9) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,9)

      END IF

! soil temperature using nemsio
      varname='soilt1'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,stc(ista_2l,jsta_2l,1))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,1) = spval
          !if (sm(i,j) /= 0.0) stc(i,j,1) = spval
        enddo
      enddo
     if(debugprint)print*,'sample l','stc',' = ',1,stc(isa,jsa,1)
      
      varname='soilt2'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,stc(ista_2l,jsta_2l,2))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,2) = spval
          !if (sm(i,j) /= 0.0) stc(i,j,2) = spval
        enddo
      enddo
     if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,2)
      
      varname='soilt3'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,stc(ista_2l,jsta_2l,3))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,3) = spval 
          !if (sm(i,j) /= 0.0) stc(i,j,3) = spval
        enddo
      enddo
     if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,3)
      
      varname='soilt4'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,stc(ista_2l,jsta_2l,4))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,4) = spval
          !if (sm(i,j) /= 0.0) stc(i,j,4) = spval
        enddo
      enddo
     if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,4)

      IF (nsoil==9) THEN

      varname='soilt5'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,stc(ista_2l,jsta_2l,5))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,5) = spval
          !if (sm(i,j) /= 0.0) stc(i,j,5) = spval
        enddo
      enddo
     if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,5)

      varname='soilt6'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,stc(ista_2l,jsta_2l,6))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,6) = spval
          !if (sm(i,j) /= 0.0) stc(i,j,6) = spval
        enddo
      enddo
     if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,6)

      varname='soilt7'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,stc(ista_2l,jsta_2l,7))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,7) = spval
          !if (sm(i,j) /= 0.0) stc(i,j,7) = spval         enddo
        enddo
      enddo
     if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,7)

      varname='soilt8'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,stc(ista_2l,jsta_2l,8))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,8) = spval
          !if (sm(i,j) /= 0.0) stc(i,j,8) = spval
        enddo
      enddo
     if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,8)

      varname='soilt9'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,stc(ista_2l,jsta_2l,9))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,9) = spval
          !if (sm(i,j) /= 0.0) stc(i,j,9) = spval
        enddo
      enddo
     if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,9)

      END IF
!
! E. James - 27 Sep 2022: this is for RRFS, adding smoke and dust
! extinction; it needs to be after ZINT is defined.
!
      if (modelname == 'FV3R') then
       do l = 1, lm
        do j = jsta_2l, jend_2u
         do i = ista_2l, iend_2u
          if(ext550(i,j,l)<spval)then
            taod5503d( i, j, l) = ext550( i, j, l )
            dz = zint( i, j, l ) - zint( i, j, l+1 )
            aextc55( i, j, l ) = taod5503d( i, j, l ) / dz
          endif
         if(debugprint.and.i==im/2.and.j==(jsta+jend)/2)print*,'sample taod5503d= ',   &
           i,j,l,taod5503d( i, j, l )
         if(debugprint.and.i==im/2.and.j==(jsta+jend)/2)print*,'sample dz= ',          &
           dz
         if(debugprint.and.i==im/2.and.j==(jsta+jend)/2)print*,'sample AEXTC55= ',     &
           i,j,l,aextc55( i, j, l )
         end do
        end do
       end do
       deallocate(ext550)
      end if

!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          acfrcv(i,j) = spval ! GFS does not output time averaged convective and strat cloud fraction, set acfrcv to spval, ncfrcv to 1
          ncfrcv(i,j) = 1.0
          acfrst(i,j) = spval ! GFS does not output time averaged cloud fraction, set acfrst to spval, ncfrst to 1
          ncfrst(i,j) = 1.0
          bgroff(i,j) = spval ! GFS does not have UNDERGROUND RUNOFF
        enddo
      enddo
!     trdlw(i,j)  = 6.0
      ardlw = 1.0 ! GFS incoming sfc longwave has been averaged over 6 hr bucket, set ARDLW to 1

! time averaged incoming sfc longwave
      varname='dlwrf_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,alwin)

! inst incoming sfc longwave
      varname='dlwrf'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,rlwin)
                                                            
! time averaged outgoing sfc longwave
      varname='ulwrf_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,alwout)

! inst outgoing sfc longwave 
      varname='ulwrf'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,radot)

!     where(alwout /= spval) alwout=-alwout ! CLDRAD puts a minus sign before gribbing
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (alwout(i,j) /= spval) alwout(i,j) = -alwout(i,j)
        enddo
      enddo
!     if(debugprint)print*,'sample l',VarName,' = ',1,alwout(isa,jsa)

! time averaged outgoing model top longwave using gfsio
      varname='ulwrf_avetoa'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,alwtoa)
!     if(debugprint)print*,'sample l',VarName,' = ',1,alwtoa(isa,jsa)

! instant outgoing model top longwave
      varname='ulwrf_toa'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,rlwtoa)
!     if(debugprint)print*,'sample l',VarName,' = ',1,rlwtoa(isa,jsa)
      
! GFS incoming sfc longwave has been averaged, set ARDLW to 1
      ardsw=1.0
!     trdsw=6.0

! time averaged incoming sfc shortwave 
      varname='dswrf_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,aswin)
!     if(debugprint)print*,'sample l',VarName,' = ',1,aswin(isa,jsa)

! inst incoming sfc shortwave 
      varname='dswrf'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,rswin)

! inst incoming clear sky sfc shortwave
! FV3 do not output instant incoming clear sky sfc shortwave
      !$omp parallel do private(i,j)
      do j=jsta_2l,jend_2u
        do i=ista_2l,iend_2u
          rswinc(i,j) = spval 
        enddo
      enddo

! time averaged incoming sfc uv-b using getgb
      varname='duvb_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,auvbin)
!     if(debugprint)print*,'sample l',VarName,' = ',1,auvbin(isa,jsa)
       
! time averaged incoming sfc clear sky uv-b using getgb
      varname='cduvb_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,auvbinc)
!     if(debugprint)print*,'sample l',VarName,' = ',1,auvbinc(isa,jsa)
      
! time averaged outgoing sfc shortwave using gfsio
      varname='uswrf_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,aswout)
!     where(aswout /= spval) aswout=-aswout ! CLDRAD puts a minus sign before gribbing 
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (aswout(i,j) /= spval) aswout(i,j) = -aswout(i,j)
        enddo
      enddo
!     if(debugprint)print*,'sample l',VarName,' = ',1,aswout(isa,jsa)

! inst outgoing sfc shortwave using gfsio
      varname='uswrf'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,rswout)

! time averaged model top incoming shortwave
      varname='dswrf_avetoa'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,aswintoa)
!     if(debugprint)print*,'sample l',VarName,' = ',1,aswintoa(isa,jsa)

! time averaged model top outgoing shortwave
      varname='uswrf_avetoa'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,aswtoa)
!     if(debugprint)print*,'sample l',VarName,' = ',1,aswtoa(isa,jsa)

! time averaged surface sensible heat flux, multiplied by -1 because wrf model flux
! has reversed sign convention using gfsio
      varname='shtfl_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sfcshx)
!     where (sfcshx /= spval)sfcshx=-sfcshx
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sfcshx(i,j) /= spval) sfcshx(i,j) = -sfcshx(i,j)
        enddo
      enddo
!     if(debugprint)print*,'sample l',VarName,' = ',1,sfcshx(isa,jsa)

! inst surface sensible heat flux
      varname='shtfl'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,twbs)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (twbs(i,j) /= spval) twbs(i,j) = -twbs(i,j)
        enddo
      enddo

! GFS surface flux has been averaged, set  ASRFC to 1 
      asrfc=1.0  
!      tsrfc=6.0

! time averaged surface latent heat flux, multiplied by -1 because wrf model flux
! has reversed sign vonvention using gfsio
      varname='lhtfl_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sfclhx)
!     where (sfclhx /= spval)sfclhx=-sfclhx
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sfclhx(i,j) /= spval) sfclhx(i,j) = -sfclhx(i,j)
        enddo
      enddo
!     if(debugprint)print*,'sample l',VarName,' = ',1,sfclhx(isa,jsa)

! inst surface latent heat flux
      varname='lhtfl'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,qwbs)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (qwbs(i,j) /= spval) qwbs(i,j) = -qwbs(i,j)
        enddo
      enddo

      if(me==0)print*,'rdaod= ',rdaod
! inst aod550 optical depth
      if(rdaod) then
      varname='aod550'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,aod550)

      varname='du_aod550'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,du_aod550)

      varname='ss_aod550'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ss_aod550)

      varname='su_aod550'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,su_aod550)

      varname='oc_aod550'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,oc_aod550)

      varname='bc_aod550'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,bc_aod550)
      end if

! time averaged ground heat flux using nemsio
      varname='gflux_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,subshx)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) subshx(i,j) = spval
        enddo
      enddo
!     if(debugprint)print*,'sample l',VarName,' = ',1,subshx(isa,jsa)

! inst ground heat flux using nemsio
      varname='gflux'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,grnflx)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) grnflx(i,j) = spval
        enddo
      enddo

! time averaged zonal momentum flux using gfsio
      varname='uflx_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sfcux)
!     if(debugprint)print*,'sample l',VarName,' = ',1,sfcux(isa,jsa)
      
! time averaged meridional momentum flux using nemsio
      varname='vflx_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sfcvx)
!     if(debugprint)print*,'sample l',VarName,' = ',1,sfcvx(isa,jsa)

! dong read in inst surface flux 
! inst zonal momentum flux using gfsio
      varname='uflx'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sfcuxi)
     if(debugprint)print*,'sample l',varname,' = ',1,sfcuxi(isa,jsa)

! inst meridional momentum flux using nemsio
      varname='vflx'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,sfcvxi)
     if(debugprint)print*,'sample l',varname,' = ',1,sfcvxi(isa,jsa)

     
!$omp parallel do private(i,j)
      do j=jsta_2l,jend_2u
        do i=ista_2l,iend_2u
          sfcuvx(i,j) = spval ! GFS does not use total momentum flux
        enddo
      enddo

! time averaged zonal gravity wave stress using nemsio
      varname='u-gwd_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,gtaux)
!     if(debugprint)print*,'sample l',VarName,' = ',1,gtaux(isa,jsa)

! time averaged meridional gravity wave stress using getgb
      varname='v-gwd_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,gtauy)
!     if(debugprint)print*,'sample l',VarName,' = ',1,gtauy(isa,jsa)
                                                     
! time averaged accumulated potential evaporation
      varname='pevpr_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgpotevp)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) avgpotevp(i,j) = spval
        enddo
      enddo
!     if(debugprint)print*,'sample l',VarName,' = ',1,potevp(isa,jsa)

! inst potential evaporation
      varname='pevpr'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,potevp)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) potevp(i,j) = spval
        enddo
      enddo

      do l=1,lm
!$omp parallel do private(i,j)
        do j=jsta_2l,jend_2u
          do i=ista_2l,iend_2u
! GFS does not have temperature tendency due to long wave radiation
            rlwtt(i,j,l)  = spval
! GFS does not have temperature tendency due to short wave radiation
            rswtt(i,j,l)  = spval
! GFS does not have temperature tendency due to latent heating from convection
            tcucn(i,j,l)  = spval
            tcucns(i,j,l) = spval
! GFS does not have temperature tendency due to latent heating from grid scale
            train(i,j,l)  = spval
          enddo
        enddo
      enddo

! set avrain to 1
      avrain=1.0
      avcnvc=1.0
      theat=6.0 ! just in case GFS decides to output T tendency   
      
! 10 m u using nemsio
      varname='ugrd10m'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,u10)

      do j=jsta,jend
        do i=ista,iend
          u10h(i,j)=u10(i,j)
        end do
      end do
!     if(debugprint)print*,'sample l',VarName,' = ',1,u10(isa,jsa)
            
! 10 m v using gfsio
      varname='vgrd10m'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,v10)

      do j=jsta,jend
        do i=ista,iend
          v10h(i,j)=v10(i,j)
        end do
      end do
!     if(debugprint)print*,'sample l',VarName,' = ',1,v10(isa,jsa)
      
! vegetation type, it's in GFS surface file, hopefully will merge into gfsio soon 
      varname='vtype'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,buf)
!     where (buf /= spval)
!      ivgtyp=nint(buf)
!     elsewhere
!      ivgtyp=0 !need to feed reasonable value to crtm
!     end where 
!$omp parallel do private(i,j)
      do j = jsta_2l, jend_2u
        do i=ista,iend
          if (buf(i,j) < spval) then
            ivgtyp(i,j) = nint(buf(i,j))
          else
            ivgtyp(i,j) = 0
          endif
        enddo
      enddo
!     if(debugprint)print*,'sample l',VarName,' = ',1,ivgtyp(isa,jsa)
      
! soil type, it's in GFS surface file, hopefully will merge into gfsio soon
      varname='sotyp'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,buf)
      l=1
!$omp parallel do private(i,j)
      do j = jsta_2l, jend_2u
        do i=ista,iend
          if (buf(i,j) < spval) then
            isltyp(i,j) = nint(buf(i,j))
          else
            isltyp(i,j) = 0 !need to feed reasonable value to crtm
          endif
        enddo
      enddo
!     if(debugprint)print*,'sample l',VarName,' = ',1,isltyp(isa,jsa)
      
      varname='wetness'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,buf)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          smstav(i,j) = buf(i,j)
        enddo
      enddo
      varname='snacc_land'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,buf)
      varname='snacc_ice'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,buf2)
!$omp parallel do private(i,j)
      do j = jsta_2l, jend_2u
        do i=ista,iend
          if(buf(i,j)<spval .and. buf2(i,j)<spval) then
            sndepac(i,j) = buf(i,j) + buf2(i,j)
          else
            sndepac(i,j) = spval
          endif
        enddo
      enddo
!$omp parallel do private(i,j)
      do j=jsta_2l,jend_2u
        do i=ista_2l,iend_2u
!          smstav(i,j) = spval    ! GFS does not have soil moisture availability
!          smstot(i,j) = spval    ! GFS does not have total soil moisture
          sfcevp(i,j) = spval    ! GFS does not have accumulated surface evaporation
          acsnom(i,j) = spval    ! GFS does not have snow melt
!          sst(i,j)    = spval    ! GFS does not have sst????
          thz0(i,j)   = ths(i,j) ! GFS does not have THZ0, use THS to substitute
          qz0(i,j)    = spval    ! GFS does not output humidity at roughness length
          uz0(i,j)    = spval    ! GFS does not output u at roughness length
          vz0(i,j)    = spval    ! GFS does not output humidity at roughness length
        enddo
      enddo
      do l=1,lm
!$omp parallel do private(i,j)
        do j=jsta_2l,jend_2u
          do i=ista_2l,iend_2u
            el_pbl(i,j,l) = spval    ! GFS does not have mixing length
            exch_h(i,j,l) = spval    ! GFS does not output exchange coefficient
          enddo
        enddo
      enddo
!     if(debugprint)print*,'sample l',VarName,' = ',1,thz0(isa,jsa)

! retrieve inst convective cloud top, GFS has cloud top pressure instead of index,
! will need to modify CLDRAD.f to use pressure directly instead of index
!      VarName='pres'
!      VcoordName='convect-cld top' 
!      l=1
!     if(debugprint)print*,'sample l',VarName,' = ',1,ptop(isa,jsa)
      varname='prescnvclt'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ptop)

      
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          htop(i,j) = spval
          if(ptop(i,j) <= 0.0) ptop(i,j) = spval
        enddo
      enddo
      do j=jsta,jend
        do i=ista,iend
          if(ptop(i,j) < spval)then
            do l=1,lm
              if(ptop(i,j) <= pmid(i,j,l))then
                htop(i,j) = l
!                if(i==ii .and. j==jj)print*,'sample ptop,pmid pmid-1,pint= ',   &
!                ptop(i,j),pmid(i,j,l),pmid(i,j,l-1),pint(i,j,l),htop(i,j)
                 exit
              end if
            end do
          end if 
        end do
      end do

! retrieve inst convective cloud bottom, GFS has cloud top pressure instead of index,
! will need to modify CLDRAD.f to use pressure directly instead of index
      varname='prescnvclb'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,pbot)
!     if(debugprint)print*,'sample l',VarName,VcoordName,' = ',1,pbot(isa,jsa)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          hbot(i,j) = spval
          if(pbot(i,j) <= 0.0) pbot(i,j) = spval
        enddo
      enddo
      do j=jsta,jend
        do i=ista,iend
!     if(.not.lb(i,j))print*,'false bitmask for pbot at '
!     +     ,i,j,pbot(i,j)
          if(pbot(i,j) < spval)then
            do l=lm,1,-1
              if(pbot(i,j) >= pmid(i,j,l)) then
                hbot(i,j) = l
!                if(i==ii .and. j==jj)print*,'sample pbot,pmid= ',    &
!                                pbot(i,j),pmid(i,j,l),hbot(i,j)
                exit
              end if
            end do
          end if 
        end do
      end do
      if(debugprint)print*,'sample hbot = ',hbot(isa,jsa)
! retrieve time averaged low cloud top pressure using nemsio
      varname='pres_avelct'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ptopl)
!     if(debugprint)print*,'sample l',VarName,' = ',1,ptopl(isa,jsa)

! retrieve time averaged low cloud bottom pressure using nemsio
      varname='pres_avelcb'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,pbotl)
!     if(debugprint)print*,'sample l',VarName,' = ',1,pbotl(isa,jsa)
     
! retrieve time averaged low cloud top temperature using nemsio
      varname='tmp_avelct'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ttopl)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopl(isa,jsa)

! retrieve time averaged middle cloud top pressure using nemsio
      varname='pres_avemct'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ptopm)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptopm(isa,jsa)
                                                             
! retrieve time averaged middle cloud bottom pressure using  nemsio
      varname='pres_avemcb'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,pbotm)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pbotm(isa,jsa)
      
! retrieve time averaged middle cloud top temperature using nemsio
      varname='tmp_avemct'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ttopm)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopm(isa,jsa)
      
! retrieve time averaged high cloud top pressure using nemsio *********
      varname='pres_avehct'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ptoph)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptoph(isa,jsa)
     
! retrieve time averaged high cloud bottom pressure using  nemsio
      varname='pres_avehcb'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,pboth)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pboth(isa,jsa)

! retrieve time averaged high cloud top temperature using nemsio
      varname='tmp_avehct'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ttoph)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttoph(isa,jsa)
      
! retrieve boundary layer cloud cover using nemsio
      varname='tcdc_avebndcl'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,pblcfr)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,pblcfr(isa,jsa)
!     where (pblcfr /= spval)pblcfr=pblcfr/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j = jsta_2l, jend_2u
        do i=ista,iend
          if (pblcfr(i,j) < spval) pblcfr(i,j) = pblcfr(i,j) * 0.01
        enddo
      enddo
        
! retrieve cloud work function 
      varname='cwork_aveclm'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,cldwork)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,cldwork(isa,jsa)
      
! accumulated total (base+surface) runoff
      varname='watr_acc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,runoff)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) runoff(i,j) = spval
        enddo
      enddo

! total water storage in aquifer
      varname='wa_acc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,twa)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) twa(i,j) = spval
        enddo
      enddo
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,runoff(isa,jsa)

! accumulated evaporation of intercepted water
      varname='ecan_acc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,tecan)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) tecan(i,j) = spval
        enddo
      enddo

! accumulated plant transpiration
      varname='etran_acc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,tetran)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) tetran(i,j) = spval
        enddo
      enddo

! accumulated soil surface evaporation
      varname='edir_acc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,tedir)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) tedir(i,j) = spval
        enddo
      enddo

! retrieve shelter max temperature using nemsio
      varname='t02max'
      if(modelname=='GFS') varname='tmax_max2m'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,maxtshltr)

! retrieve shelter min temperature using nemsio
      varname='t02min'
      if(modelname=='GFS') varname='tmin_min2m'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,mintshltr)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
!     1,mintshltr((ista+iend)/2,(jsta+jend)/2)

! retrieve shelter max RH
      varname='rh02max'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,maxrhshltr)

! retrieve shelter min temperature using nemsio
      varname='rh02min'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,minrhshltr)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
!     1,mintshltr((ista+iend)/2,(jsta+jend)/2)

! retrieve shelter max specific humidity using nemsio
      varname='spfhmax_max2m'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,maxqshltr)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',
!     1,maxqshltr(isa,jsa)

! retrieve shelter min temperature using nemsio
      varname='spfhmin_min2m'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,minqshltr)
 
! retrieve ice thickness using nemsio
      varname='icetk'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,dzice)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,dzice(isa,jsa)

! retrieve wilting point using nemsio
      varname='wilt'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,smcwlt)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) smcwlt(i,j) = spval
        enddo
      enddo
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,smcwlt(isa,jsa)
      
! retrieve sunshine duration using nemsio
      varname='sunsd_acc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,suntime)

! retrieve field capacity using nemsio
      varname='fldcp'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,fieldcapa)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) fieldcapa(i,j) = spval
        enddo
      enddo
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,fieldcapa(isa,jsa)

! retrieve time averaged surface visible beam downward solar flux
      varname='vbdsf_ave'
       call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avisbeamswin)
      vcoordname='sfc'
      l=1

! retrieve time averaged surface visible diffuse downward solar flux
      varname='vddsf_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avisdiffswin)

! retrieve time averaged surface near IR beam downward solar flux
      varname='nbdsf_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,airbeamswin)

! retrieve time averaged surface near IR diffuse downward solar flux
      varname='nddsf_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,airdiffswin)

! retrieve time averaged surface clear sky outgoing LW
      varname='csulf'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,alwoutc)

! retrieve time averaged TOA clear sky outgoing LW
      varname='csulftoa'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,alwtoac)

! retrieve time averaged surface clear sky outgoing SW
      varname='csusf'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,aswoutc)

! retrieve time averaged TOA clear sky outgoing LW
      varname='csusftoa'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,aswtoac)

! retrieve time averaged surface clear sky incoming LW
      varname='csdlf'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,alwinc)

! retrieve time averaged surface clear sky incoming SW
      varname='csdsf'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,aswinc)

! retrieve storm runoff using nemsio
      varname='ssrun_acc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,ssroff)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) ssroff(i,j) = spval
        enddo
      enddo

! retrieve direct soil evaporation
      varname='evbs_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgedir)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) avgedir(i,j) = spval
        enddo
      enddo

! retrieve CANOPY WATER EVAP 
      varname='evcw_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgecan)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) avgecan(i,j) = spval
        enddo
      enddo

! retrieve AVERAGED PRECIP ADVECTED HEAT FLUX
      varname='pah_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,paha)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) paha(i,j) = spval
        enddo
      enddo

! retrieve instantaneous PRECIP ADVECTED HEAT FLUX
      varname='pahi'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,pahi)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) pahi(i,j) = spval
        enddo
      enddo

! retrieve PLANT TRANSPIRATION 
      varname='trans_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgetrans)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) avgetrans(i,j) = spval
        enddo
      enddo

! retrieve snow sublimation
      varname='sbsno_ave'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,avgesnow)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j)==1.0 .and. sice(i,j)==0.) avgesnow(i,j)=spval
        enddo
      enddo

! retrive total soil moisture
      varname='soilm'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,smstot)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) smstot(i,j) = spval
        enddo
      enddo

! retrieve snow phase change heat flux
      varname='snohf'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,snopcx)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          if (sm(i,j) /= 0.0) snopcx(i,j) = spval
        enddo
      enddo

! retrieve pwater
      varname='pwat'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
      spval,varname,pwat)
      
! GFS does not have deep convective cloud top and bottom fields

!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=ista,iend
          htopd(i,j) = spval
          hbotd(i,j) = spval   
          htops(i,j) = spval
          hbots(i,j) = spval 
          cuppt(i,j) = spval 
        enddo
      enddo


      print *, 'gocart_on=',gocart_on
      print *, 'gccpp_on=',gccpp_on
      print *, 'nasa_on=',nasa_on
      if ((gocart_on .or. gccpp_on) .and. d2d_chem) then


! retrieve dust emission fluxes
       do k = 1, nbin_du
       if ( k == 1) varname='duem001'
       if ( k == 2) varname='duem002'
       if ( k == 3) varname='duem003'
       if ( k == 4) varname='duem004'
       if ( k == 5) varname='duem005'

      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)

      duem(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve dust sedimentation fluxes
      do k = 1, nbin_du
       if ( k == 1) varname='dust1sd'
       if ( k == 2) varname='dust2sd'
       if ( k == 3) varname='dust3sd'
       if ( k == 4) varname='dust4sd'
       if ( k == 5) varname='dust5sd'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      dusd(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve dust dry deposition fluxes
      do k = 1, nbin_du
       if ( k == 1) varname='dust1dp'
       if ( k == 2) varname='dust2dp'
       if ( k == 3) varname='dust3dp'
       if ( k == 4) varname='dust4dp'
       if ( k == 5) varname='dust5dp'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      dudp(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve dust wet deposition fluxes
      do k = 1, nbin_du
       if ( k == 1) varname='dust1wtl'
       if ( k == 2) varname='dust2wtl'
       if ( k == 3) varname='dust3wtl'
       if ( k == 4) varname='dust4wtl'
       if ( k == 5) varname='dust5wtl'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      duwt(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve dust scavenging fluxes
      do k = 1, nbin_du
       if ( k == 1) varname='dust1wtc'
       if ( k == 2) varname='dust2wtc'
       if ( k == 3) varname='dust3wtc'
       if ( k == 4) varname='dust4wtc'
       if ( k == 5) varname='dust5wtc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      dusv(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve seasalt emission fluxes
      do k = 1, nbin_ss
       if ( k == 1) varname='ssem001'
       if ( k == 2) varname='ssem002'
       if ( k == 3) varname='ssem003'
       if ( k == 4) varname='ssem004'
       if ( k == 5) varname='ssem005'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      ssem(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve seasalt emission fluxes
      do k = 1, nbin_ss
       if ( k == 1) varname='seas1sd'
       if ( k == 2) varname='seas2sd'
       if ( k == 3) varname='seas3sd'
       if ( k == 4) varname='seas4sd'
       if ( k == 5) varname='seas5sd'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      sssd(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve seasalt dry deposition fluxes
      do k = 1, nbin_ss
       if ( k == 1) varname='seas1dp'
       if ( k == 2) varname='seas2dp'
       if ( k == 3) varname='seas3dp'
       if ( k == 4) varname='seas4dp'
       if ( k == 5) varname='seas5dp'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      ssdp(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve seasalt wet deposition fluxes
      do k = 1, nbin_ss
       if ( k == 1) varname='seas1wtl'
       if ( k == 2) varname='seas2wtl'
       if ( k == 3) varname='seas3wtl'
       if ( k == 4) varname='seas4wtl'
       if ( k == 5) varname='seas5wtl'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      sswt(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve seasalt scavenging fluxes
      do k = 1, nbin_ss
       if ( k == 1) varname='seas1wtc'
       if ( k == 2) varname='seas1wtc'
       if ( k == 3) varname='seas1wtc'
       if ( k == 4) varname='seas1wtc'
       if ( k == 5) varname='seas1wtc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      sssv(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve bc emission fluxes
      do k = 1, nbin_bc
       if ( k == 1) varname='bceman'
       if ( k == 2) varname='bcembb'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      bcem(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve bc sedimentation fluxes
      do k = 1, nbin_bc
       if ( k == 1) varname='bc1sd'
       if ( k == 2) varname='bc2sd'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      bcsd(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve bc dry deposition fluxes
      do k = 1, nbin_bc
       if ( k == 1) varname='bc1dp'
       if ( k == 2) varname='bc2dp'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      bcdp(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve bc large wet deposition fluxes
      do k = 1, nbin_bc
       if ( k == 1) varname='bc1wtl'
       if ( k == 2) varname='bc2wtl'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      bcwt(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve bc convective wet deposition fluxes
      do k = 1, nbin_bc
       if ( k == 1) varname='bc1wtc'
       if ( k == 2) varname='bc2wtc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      bcsv(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve oc emission fluxes
      do k = 1, nbin_oc
       if ( k == 1) varname='oceman'
       if ( k == 2) varname='ocembb'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      ocem(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve oc sedimentation fluxes
      do k = 1, nbin_oc
       if ( k == 1) varname='oc1sd'
       if ( k == 2) varname='oc2sd'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      ocsd(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve oc dry deposition fluxes
      do k = 1, nbin_oc
       if ( k == 1) varname='oc1dp'
       if ( k == 2) varname='oc2dp'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      ocdp(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve oc large wet deposition fluxes
      do k = 1, nbin_oc
       if ( k == 1) varname='oc1wtl'
       if ( k == 2) varname='oc2wtl'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      ocwt(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve oc convective wet deposition fluxes
      do k = 1, nbin_oc
       if ( k == 1) varname='oc1wtc'
       if ( k == 2) varname='oc2wtc'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      ocsv(1:im,jsta_2l:jend_2u,k)=chem_2d(1:im,jsta_2l:jend_2u)
       enddo

! retrieve MIE AOD
       varname='maod'
      call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u,&
      spval,varname,chem_2d)
      maod(1:im,jsta_2l:jend_2u)=chem_2d(1:im,jsta_2l:jend_2u)

       endif ! gocart_on
! done with flux file, close it for now
      status=nf90_close(ncid2d)
!      deallocate(tmp,recname,reclevtyp,reclev)

! pos east
!       call collect_loc(gdlat,dummy)
!       if(me == 0)then
!        latstart = nint(dummy(1,1)*gdsdegr)
!        latlast  = nint(dummy(im,jm)*gdsdegr)
!        print*,'laststart,latlast B bcast= ',latstart,latlast,'gdsdegr=',gdsdegr,&
!          'dummy(1,1)=',dummy(1,1),dummy(im,jm),'gdlat=',gdlat(1,1)
!       end if
!       call mpi_bcast(latstart,1,MPI_INTEGER,0,mpi_comm_comp,irtn)
!       call mpi_bcast(latlast,1,MPI_INTEGER,0,mpi_comm_comp,irtn)
!       write(6,*) 'laststart,latlast,me A calling bcast=',latstart,latlast,me
!       call collect_loc(gdlon,dummy)
!       if(me == 0)then
!        lonstart = nint(dummy(1,1)*gdsdegr)
!        lonlast  = nint(dummy(im,jm)*gdsdegr)
!       end if
!       call mpi_bcast(lonstart,1,MPI_INTEGER,0,mpi_comm_comp,irtn)
!       call mpi_bcast(lonlast, 1,MPI_INTEGER,0,mpi_comm_comp,irtn)

!       write(6,*)'lonstart,lonlast A calling bcast=',lonstart,lonlast
!

! generate look up table for lifted parcel calculations

      thl    = 210.
      plq    = 70000.
      pt_tbl = 10000.          ! this is for 100 hPa added by Moorthi

      CALL table(ptbl,ttbl,pt_tbl,                                     &
                 rdq,rdth,rdp,rdthe,pl,thl,qs0,sqs,sthe,the0)

      CALL tableq(ttblq,rdpq,rdtheq,plq,thl,stheq,the0q)

!     
!     
      IF(me == 0)THEN
        WRITE(6,*)'  SPL (POSTED PRESSURE LEVELS) BELOW: '
        WRITE(6,51) (spl(l),l=1,lsm)
   50   FORMAT(14(f4.1,1x))
   51   FORMAT(8(f8.1,1x))
      ENDIF
!     
!$omp parallel do private(l)
      DO l = 1,lsm
         alsl(l) = log(spl(l))
      END DO
!
!HC WRITE IGDS OUT FOR WEIGHTMAKER TO READ IN AS KGDSIN
      if(me == 0)then
        print*,'writing out igds'
        igdout = 110
!        open(igdout,file='griddef.out',form='unformatted'
!     +  ,status='unknown')
        if(maptype == 1)THEN  ! Lambert conformal
          WRITE(igdout)3
          WRITE(6,*)'igd(1)=',3
          WRITE(igdout)im
          WRITE(igdout)jm
          WRITE(igdout)latstart
          WRITE(igdout)lonstart
          WRITE(igdout)8
          WRITE(igdout)cenlon
          WRITE(igdout)dxval
          WRITE(igdout)dyval
          WRITE(igdout)0
          WRITE(igdout)64
          WRITE(igdout)truelat2
          WRITE(igdout)truelat1
          WRITE(igdout)255
        ELSE IF(maptype  ==  2)THEN  !Polar stereographic
          WRITE(igdout)5
          WRITE(igdout)im
          WRITE(igdout)jm
          WRITE(igdout)latstart
          WRITE(igdout)lonstart
          WRITE(igdout)8
          WRITE(igdout)cenlon
          WRITE(igdout)dxval
          WRITE(igdout)dyval
          WRITE(igdout)0
          WRITE(igdout)64
          WRITE(igdout)truelat2  !Assume projection at +-90
          WRITE(igdout)truelat1
          WRITE(igdout)255
        !  Note: The calculation of the map scale factor at the standard
        !        lat/lon and the PSMAPF
        ! Get map factor at 60 degrees (N or S) for PS projection, which will
        ! be needed to correctly define the DX and DY values in the GRIB GDS
          if (truelat1 < 0.) THEN
            lat = -60.
          else
            lat = 60.
          end if

          CALL msfps(lat,truelat1*0.001,psmapf)

        ELSE IF(maptype == 3) THEN  !Mercator
          WRITE(igdout)1
          WRITE(igdout)im
          WRITE(igdout)jm
          WRITE(igdout)latstart
          WRITE(igdout)lonstart
          WRITE(igdout)8
          WRITE(igdout)latlast
          WRITE(igdout)lonlast
          WRITE(igdout)truelat1
          WRITE(igdout)0
          WRITE(igdout)64
          WRITE(igdout)dxval
          WRITE(igdout)dyval
          WRITE(igdout)255
        ELSE IF(maptype == 0 .OR. maptype == 203)THEN  !A STAGGERED E-GRID
          WRITE(igdout)203
          WRITE(igdout)im
          WRITE(igdout)jm
          WRITE(igdout)latstart
          WRITE(igdout)lonstart
          WRITE(igdout)136
          WRITE(igdout)cenlat
          WRITE(igdout)cenlon
          WRITE(igdout)dxval
          WRITE(igdout)dyval
          WRITE(igdout)64
          WRITE(igdout)0
          WRITE(igdout)0
          WRITE(igdout)0
         ELSE IF(maptype == 207)THEN !Rotated lat-lon grid
           write(flatlon,1001)ifhr
           open(112,file=trim(flatlon),form='formatted', &
             status='unknown')
           write(112,1002)latstart/1000,lonstart/1000,&
             latse/1000,lonse/1000,latnw/1000,lonnw/1000,&
             latlast/1000,lonlast/1000
     1001 format('latlons_corners.txt.f',i3.3)
     1002 format(4(i6,i7,x))
           close(112)
        END IF
      end if
!     
!

      RETURN
      END


      subroutine read_netcdf_3d_para(ncid,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
                 spval,varname,buf,lm)

      use netcdf
      use ctlblk_mod, only : me
      use params_mod, only : small
      implicit none
      include "mpif.h"

      character(len=20),intent(in) :: varname
      real,intent(in)    :: spval
      integer,intent(in) :: ncid,im,jm,lm,jsta_2l,jend_2u,jsta,jend
      integer,intent(in) :: ista_2l,iend_2u,ista,iend
      real,intent(out)   :: buf(ista_2l:iend_2u,jsta_2l:jend_2u,lm)
      integer            :: varid,iret,ii,jj,i,j,l,kk
      integer            :: start(3), count(3), stride(3)
      real,parameter     :: spval_netcdf=9.99e+20
      real               :: fill_value

      iret = nf90_inq_varid(ncid,trim(varname),varid)
      if (iret /= 0) then
        if (me == 0) print*,varname," not found -Assigned missing values"
!$omp parallel do private(i,j,l)
          do l=1,lm
            do j=jsta,jend
              do i=ista,iend
                buf(i,j,l)=spval
              enddo
            enddo
          enddo
      else
        iret = nf90_get_att(ncid,varid,"_FillValue",fill_value)
        if (iret /= 0) fill_value = spval_netcdf
        start = (/ista,jsta,1/)
        ii=iend-ista+1
        jj=jend-jsta+1
        count = (/ii,jj,lm/)
        iret = nf90_get_var(ncid,varid,buf(ista:iend,jsta:jend,1:lm),start=start,count=count)
        if (iret /= 0) then
          print*," iret /=0, Error in reading varid "
        endif
        do l=1,lm
          do j=jsta,jend
            do i=ista,iend
              if(abs(buf(i,j,l)-fill_value)<small)buf(i,j,l)=spval
            end do
          end do
        end do
      endif

      end subroutine read_netcdf_3d_para

      subroutine read_netcdf_2d_para(ncid,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
                 spval,varname,buf)

      use netcdf
      use ctlblk_mod, only : me
      use params_mod, only : small
      implicit none
      include "mpif.h"

      character(len=20),intent(in) :: varname
      real,intent(in)    :: spval
      integer,intent(in) :: ncid,jsta_2l,jend_2u,jsta,jend,ista_2l,iend_2u,ista,iend
      real,intent(out)   :: buf(ista_2l:iend_2u,jsta_2l:jend_2u)
      integer            :: varid,iret,ii,jj,i,j,l,kk
      integer            :: start(2), count(2)
      real,parameter     :: spval_netcdf=9.99e+20
      real               :: fill_value

      iret = nf90_inq_varid(ncid,trim(varname),varid)
      if (iret /= 0) then
        if (me==0) print*,varname," not found -Assigned missing values"
!$omp parallel do private(i,j)
        do j=jsta,jend
          do i=ista,iend
            buf(i,j)=spval
          enddo
        enddo
      else
        iret = nf90_get_att(ncid,varid,"_FillValue",fill_value)
        if (iret /= 0) fill_value = spval_netcdf
        start = (/ista,jsta/)
        ii=iend-ista+1
        jj=jend-jsta+1
        count = (/ii,jj/)
        iret = nf90_get_var(ncid,varid,buf(ista:iend,jsta:jend),start=start,count=count)
        if (iret /= 0) then
          print*," iret /=0, Error in reading varid "
        endif
        do j=jsta,jend
          do i=ista,iend
            if(abs(buf(i,j)-fill_value)<small)buf(i,j)=spval
          end do
        end do
      endif

      end subroutine read_netcdf_2d_para