INITPOST_GFS_NEMS_MPIIO.f

Go to the documentation of this file.

!         and 2D diag. output (d2d_chem) for GEFS-Aerosols and CCPP-Chem model.
      SUBROUTINE initpost_gfs_nems_mpiio(iostatusAER)


      use vrbls4d, only: dust, salt, suso, soot, waso, pp25, pp10 
      use vrbls3d, only: t, q, uh, vh,wh,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, qqg, 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, ref_10cm
      use vrbls2d, only: f, pd, fis, pblh, ustar, z0, ths, qs, twbs, qwbs, avgcprate,           &
              cprate, avgprec, prec, lspa, sno, si, cldefi, th10, q10, tshltr, pshltr,          &
              tshltr, albase, avgalbedo, avgtcdc, czen, czmean, mxsnal, 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, 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, maxtshltr, mintshltr, maxrhshltr,            &
              minrhshltr, dzice, smcwlt, suntime, fieldcapa, htopd, hbotd, htops, hbots,        &
              cuppt, dusmass, ducmass, dusmass25, ducmass25, aswintoa, &
              maxqshltr, minqshltr, acond, sr, u10h, v10h, &
              avgedir,avgecan,avgetrans,avgesnow,avgprec_cont,avgcprate_cont, &
              avisbeamswin,avisdiffswin,airbeamswin,airdiffswin, &
              alwoutc,alwtoac,aswoutc,aswtoac,alwinc,aswinc,avgpotevp,snoavg, &
              dustcb,bccb,occb,sulfcb,sscb,dustallcb,ssallcb,dustpm,dustpm10,sspm,pp25cb, &
              pp10cb,maod,ti 
      use soil,  only: sldpth, sh2o, smc, stc
      use masks, only: lmv, lmh, htm, vtm, gdlat, gdlon, dx, dy, hbm2, sm, sice
!     use kinds, only: i_llong
!     use nemsio_module, only: nemsio_gfile, nemsio_getfilehead, nemsio_getheadvar, nemsio_close
      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
      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, gocart_on, gccpp_on, pt_tbl, hyb_sigp, filenameflux, filenameaer, &
              isf_surface_physics, d2d_chem
      use gridspec_mod, only: maptype, gridtype, latstart, latlast, lonstart, lonlast, cenlon,  &
              dxval, dyval, truelat2, truelat1, psmapf, cenlat
      use nemsio_module_mpi
      use upp_physics, only: fpvsnew, caldiv, calgradps
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
      implicit none
!
      type(nemsio_gfile) :: nfile,ffile,rfile
!
!     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
      integer,intent(in) :: iostatusaer
      character(len=20)  :: varname, vcoordname
      integer            :: status, fldsize, fldst, recn 
      integer            :: recn_vvel,recn_delz,recn_dpres
      character             startdate*19,sysdepinfo*80,cgar*1
      character             startdate2(19)*4,lprecip_accu*3
! 
!     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 :: reduce_grid = .true. ! set default to true for ops GSM
      CHARACTER*32 label
      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 rinc(5)

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

      character*16,allocatable :: recname(:)
      character*16,allocatable :: reclevtyp(:)
      character*6 :: modelname_nemsio
      integer,     allocatable :: reclev(:), kmsk(:,:)
      real,        allocatable :: glat1d(:), glon1d(:), qstl(:)
      real,        allocatable :: wrk1(:,:), wrk2(:,:)
      real,        allocatable :: p2d(:,:),  t2d(:,:),  q2d(:,:),      &
                                  qs2d(:,:), cw2d(:,:), cfr2d(:,:)
      real(kind=4),allocatable :: vcoord4(:,:,:)
      real, dimension(lm+1)    :: ak5, bk5
      real*8, allocatable      :: pm2d(:,:), pi2d(:,:)
      real,   allocatable      :: tmp(:)
      real                     :: buf(im,jsta_2l:jend_2u)
      integer                  :: lonsperlat(jm/2), numi(jm)

!     real buf(im,jsta_2l:jend_2u),bufsoil(im,nsoil,jsta_2l:jend_2u)   &
!         ,buf3d(im,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
!     REAL,  PARAMETER    :: QMIN = 1.E-15

!      DATA BLANK/'    '/
!
!     integer, parameter    :: npass2=2, npass3=3
!     integer, parameter    :: npass2=20, npass3=30
      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                     :: omg1(im), omg2(im+2)

!***********************************************************************
!     START INIT HERE.
!
      WRITE(6,*)'INITPOST:  ENTER INITPOST_GFS_NEMS_MPIIO'
      WRITE(6,*)'me=',me,'LMV=',size(lmv,1),size(lmv,2),'LMH=',   &
           size(lmh,1),size(lmh,2),'jsta_2l=',jsta_2l,'jend_2u=', &
           jend_2u,'im=',im
!     
      isa = im / 2
      jsa = (jsta+jend) / 2

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

! initialize nemsio using mpi io module
      call nemsio_init()
      call nemsio_open(nfile,trim(filename),'read',mpi_comm_comp,iret=status)
      if ( status /= 0 ) then
        print*,'error opening ',filename, ' Status = ', status ; stop
      endif
      call nemsio_getfilehead(nfile,iret=status,nrec=nrec,idrt=idrt)

! open flux file early yo read imp_physics
      call nemsio_open(ffile,trim(filenameflux),'read',mpi_comm_comp &
       ,iret=status)
      if ( status /= 0 ) then
        print*,'error opening ',filenameflux, ' Status = ', status
      endif

!     
!     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 = 1, im
          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 = 1, im
            htm(i,j,l) = 1.0
            vtm(i,j,l) = 1.0
          end do
        end do
      end do

!     write(*,*)'nrec=',nrec
      allocate(recname(nrec),reclevtyp(nrec),reclev(nrec))
      allocate(glat1d(im*jm),glon1d(im*jm))
      allocate(vcoord4(lm+1,3,2))
! get start date
      call nemsio_getfilehead(nfile,iret=iret                           &
          ,idate=idate(1:7),nfhour=nfhour,recname=recname               &
          ,reclevtyp=reclevtyp,reclev=reclev,lat=glat1d                 &
          ,lon=glon1d,nframe=nframe,vcoord=vcoord4,idrt=maptype         &
          ,modelname=modelname_nemsio)
      if(iret/=0)print*,'error getting idate,nfhour'
      print *,'latstar1=',glat1d(1),glat1d(im*jm)
!
!      modelname_nemsio='FV3GFS'
      print*,'modelname = ',modelname_nemsio
      if(trim(modelname_nemsio)=='FV3GFS')reduce_grid=.false.
      if(reduce_grid)then
!------------------------------
       if (idrt == 4) then
!------------------------------
!       read lonsperlat
        open (201,file='lonsperlat.dat',status='old',form='formatted',     &
                                      action='read',iostat=iret)
        rewind(201)
        read (201,*,iostat=iret) latghf,(lonsperlat(i),i=1,latghf)
        close (201)
        print*,'finished reading lonsperlat' 

        if (jm /= latghf+latghf) then
          write(0,*)' wrong reduced grid - execution skipped'
          stop 777
        endif
        do j=1,jm/2
          numi(j) = lonsperlat(j)
        enddo
        do j=jm/2+1,jm
          numi(j) = lonsperlat(jm+1-j)
        enddo
!------------------------------
       else
!------------------------------
        do j=1,jm
          numi(j) = im
        enddo
!------------------------------
       endif
!------------------------------
      end if

! Specigy grid staggering type
      gridtype = 'A'
      if (me == 0) print *, 'maptype and gridtype is ', &
      maptype,gridtype

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

!$omp parallel do private(i,j,js)
      do j=jsta,jend
        js = (j-1)*im
        do i=1,im
          gdlat(i,j) = glat1d(js+i)
          gdlon(i,j) = glon1d(js+i)
        end do
      end do
!
!      if (hyb_sigp) then
        do l=1,lm+1
         ak5(l) = vcoord4(l,1,1)
         bk5(l) = vcoord4(l,2,1)
        enddo
!      endif

!--Fanglin Yang:  nemsio file created from FV3 does not have vcoord.
      if ( minval(ak5) <0 .or. minval(bk5) <0 ) then
       open (202,file='global_hyblev.txt',status='old',form='formatted',     &
                                       action='read',iostat=iret)
       rewind(202)
       read(202,*)
       do l=1,lm+1
        read (202,*,iostat=iret) ak5(l),bk5(l)                      
       enddo
       close (202)
 
       if (iret == 0  ) then 
         do l=1,lm+1
          vcoord4(l,1,1)=ak5(l)
          vcoord4(l,2,1)=bk5(l)
         enddo
       else
         print *, 'ak5 and bk5 not found, stop !'
         stop
       endif
      endif

      if (me == 0)then
         print *,"ak5",ak5 
         print *,"bk5",bk5 
      endif

!     deallocate(glat1d,glon1d,vcoord4)
      deallocate(glat1d,glon1d)

      print*,'idate = ',(idate(i),i=1,7)
      print*,'idate after broadcast = ',(idate(i),i=1,4)
      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)
      CALL exch(gdlon)
      print *,'after call EXCH,me=',me

!$omp parallel do private(i,j,ip1)
      do j = jsta, jend_m
        do i = 1, im
          ip1 = i + 1
          if (ip1 > im) ip1 = ip1 - im
          dx(i,j) = erad*cos(gdlat(i,j)*dtr) *(gdlon(ip1,j)-gdlon(i,j))*dtr
          dy(i,j) = erad*(gdlat(i,j)-gdlat(i,j+1))*dtr  ! like A*DPH
!     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
      
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          f(i,j) = 1.454441e-4*sin(gdlat(i,j)*dtr)   ! 2*omeg*sin(phi)
        end do
      end do
      
      impf = im
      jmpf = jm
      print*,'impf,jmpf,nframe= ',impf,jmpf,nframe
      
      iyear = idate(1)
      imn   = idate(2) ! ask Jun 
      iday  = idate(3) ! ask Jun
      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 
      
      varname='imp_physics'
      call nemsio_getheadvar(ffile,trim(varname),imp_physics,iret)
      if (iret /= 0) then
        if(me==0)print*,varname, &
        " not found in file-Assigned 99 for Zhao"
        imp_physics=99 !set GFS mp physics to 99 for Zhao scheme
      end if

      if(me==0)print*,'MP_PHYSICS= ',imp_physics

      varname='sf_surface_physi'
      call nemsio_getheadvar(ffile,trim(varname),imp_physics,iret)
      if (iret /= 0) then
        if(me==0)print*,varname, &
        " not found in file-Assigned 2 for NOAH"
        isf_surface_physics=2 !set GFS LSM physics to 2 for NOAH 
      end if

      if(me==0)print*,'SF_SURFACE_PHYSICS= ',isf_surface_physics

! read bucket
      varname='fhzero'
      call nemsio_getheadvar(ffile,trim(varname),fhzero,iret)
      if (iret /= 0) then
        if(me==0)print*,varname, &
        " not found in file-Assign 6 or 12 hours precip bucket"
        tprec   = 6.
        if(ifhr>240)tprec=12.
        tclod   = tprec
        trdlw   = tprec
        trdsw   = tprec
        tsrfc   = tprec
        tmaxmin = tprec
        td3d    = tprec
      else
        tprec=float(fhzero)
        tclod   = tprec
        trdlw   = tprec
        trdsw   = tprec
        tsrfc   = tprec
        tmaxmin = tprec
        td3d    = tprec
      end if


! read meta data to see if precip has zero bucket
!      VarName='lprecip_accu'
!      call nemsio_getheadvar(ffile,trim(VarName),lprecip_accu,iret)
!      if (iret /= 0) then
!        if(me==0)print*,VarName, &
!        " not found in file-Assign non-zero precip bucket"
!        lprecip_accu='no'
!      end if
!      if(lprecip_accu=='yes')tprec=float(ifhr)
      print*,'tprec, tclod, trdlw = ',tprec,tclod,trdlw

      
! Initializes constants for Ferrier microphysics       
      if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95) then
        CALL microinit(imp_physics)
      end if      

! GFS does not need DT to compute accumulated fields, set it to one
!      VarName='DT'
      dt   = 1

      hbm2 = 1.0

      
! start reading nemsio sigma files using parallel read
      fldsize = (jend-jsta+1)*im
      allocate(tmp(fldsize*nrec))
      print*,'allocate tmp successfully'
      tmp = 0.
      call nemsio_denseread(nfile,1,im,jsta,jend,tmp,iret=iret)
      if(iret /=0 ) then
        print*,"fail to read sigma file using mpi io read, stopping"
        stop
      end if
!
!  covert from reduced grid to full grid
!
      if(reduce_grid)then
       print*,'performing reduced grid'
       jtem = jend-jsta+1
       allocate (kmsk(im,jtem))
       kmsk = 0
       do recn=1,nrec
        fldst = (recn-1)*fldsize
        do j=jsta,jend
          js = fldst + (j-jsta)*im
          do i=1,im
            buf(i,j) = tmp(i+js)
          enddo
        enddo
        call gg2rg(im,jtem,numi(jsta),buf(1,jsta))
        call uninterpred(2,kmsk,numi(jsta),im,jtem,buf(1,jsta),tmp(fldst+1))
       enddo
       deallocate(kmsk)
      end if

! Terrain height * G   using nemsio 
      varname='hgt'
      vcoordname = 'sfc'
      l = 1
      call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
      if(recn /=0 ) then
        fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
        do j=jsta,jend
          js = fldst + (j-jsta)*im
          do i=1,im
            fis(i,j) = tmp(i+js)
          enddo
        enddo
      else
        if(me == 0) print*,'fail to read ', varname,vcoordname,l 
      endif

!      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
!      ,l,nrec,fldsize,spval,tmp &
!      ,recname,reclevtyp,reclev,VarName,VcoordName &
!      ,fis)

!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (fis(i,j) /= spval) then
            zint(i,j,lp1) = fis(i,j)
            fis(i,j)      = fis(i,j) * grav
          endif
        enddo
      enddo
      if(debugprint) print*,'sample ',varname,' = ',fis(isa,jsa)

      vcoordname = 'sfc'       ! surface fileds
      l = 1

! Surface pressure  using nemsio
      varname='pres'
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,pint(1,jsta_2l,lp1))

       if(debugprint)print*,'sample surface pressure = ',pint(isa,jsa,lp1)
      
!
!      vertical loop for Layer 3d fields
!      --------------------------------   
      vcoordname = 'mid layer'

      do l=1,lm
        ll = lm-l+1
!                                                     model level T
        !if (me == 0) print*,'start retrieving GFS T using nemsio'
        varname='tmp'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
              t(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          print*,'fail to read ', varname,' at lev ',ll, 'stopping'
          stop
        endif

        if(debugprint)print*,'sample ',ll,varname,' = ',ll,t(isa,jsa,ll)

!                                                     model level q      
        varname='spfh'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
              q(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          print*,'fail to read ', varname,' at lev ',ll, 'stopping'
          stop
        endif

        if(debugprint)print*,'sample ',ll,varname,' = ',ll,q(isa,jsa,ll)

!                                                    model level u      
        varname='ugrd'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
              uh(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          print*,'fail to read ', varname,' at lev ',ll, 'stopping'
          stop
        endif

        if(debugprint)print*,'sample ',ll,varname,' = ',ll,uh(isa,jsa,ll)
      
!                                                     model level v      
        varname='vgrd'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
              vh(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          print*,'fail to read ', varname,' at lev ',ll, 'stopping'
          stop
        endif

        if(debugprint)print*,'sample ',ll,varname,' = ',ll,vh(isa,jsa,ll)
      
! model level pressure      
!        if (.not. hyb_sigp) then
!          VarName='pres'
!          call getrecn(recname,reclevtyp,reclev,nrec,varname,VcoordName,l,recn)
!          if(recn /= 0) then
!            fldst = (recn-1)*fldsize
!!$omp parallel do private(i,j,js)
!            do j=jsta,jend
!              js = fldst + (j-jsta)*im
!              do i=1,im
!                pmid(i,j,ll) = tmp(i+js)
!              enddo
!            enddo
!          else
!            recn_pres=-9999
!            if(me==0)print*,'fail to read ', varname,' at lev ',ll, &
!            'will derive pressure using ak bk later'
!            stop
!          endif

!          if(debugprint)print*,'sample ',ll,VarName,' = ',ll,pmid(isa,jsa,ll)
!        end if      
! GFS is on A grid and does not need PMIDV        

!                                                      dp     
        varname='dpres'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
              dpres(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          recn_dpres = -9999
          if(me==0)print*,'fail to read ', varname,' at lev ',ll, &
            'will derive pressure using ak bk later'
        endif
!                                                      ozone mixing ratio
        varname='o3mr'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
              o3(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          print*,'fail to read ', varname,' at lev ',ll, 'stopping'
          stop
        endif


        if(debugprint)print*,'sample ',ll,varname,' = ',ll,o3(isa,jsa,ll)

! cloud water and ice mixing ratio  for zhao scheme
! need to look up old eta post to derive cloud water/ice from cwm
! Zhao scheme does not produce suspended rain and snow

!!$omp parallel do private(i,j)
        do j = jsta, jend
          do i=1,im
            qqw(i,j,ll) = 0. 
            qqr(i,j,ll) = 0. 
            qqs(i,j,ll) = 0. 
            qqi(i,j,ll) = 0. 
          enddo
        enddo

       if(imp_physics==99 .or. imp_physics==98)then
        varname='clwmr'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
                cwm(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          print*,'fail to read ', varname,' at lev ',ll, 'stopping'
          stop
        endif

        if(debugprint)print*,'sample ',ll,varname,' = ',ll,cwm(isa,jsa,ll)

!$omp parallel do private(i,j)
        do j=jsta,jend
          do i=1,im
            if(t(i,j,ll) < (tfrz-15.) )then ! dividing cloud water from ice
              qqi(i,j,ll) = cwm(i,j,ll)
            else
              qqw(i,j,ll) = cwm(i,j,ll)
            end if
          end do
        end do
       else if(imp_physics==11 .or. imp_physics==8)then ! GFDL or Thompson MP scheme
        varname='clwmr'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
                qqw(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          print*,'fail to read ', varname,' at lev ',ll, 'stopping'
          stop
        endif
        if(debugprint)print*,'sample ',ll,varname,' = ',ll,qqw(isa,jsa,ll)

        varname='icmr'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
                qqi(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          print*,'fail to read ', varname,' at lev ',ll, 'stopping'
          stop
        endif
        if(debugprint)print*,'sample ',ll,varname,' = ',ll,qqi(isa,jsa,ll)
      
        varname='rwmr'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
                qqr(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          print*,'fail to read ', varname,' at lev ',ll, 'stopping'
          stop
        endif
        if(debugprint)print*,'sample ',ll,varname,' = ',ll,qqr(isa,jsa,ll)

        varname ='snmr'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
                qqs(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          print*,'fail to read ', varname,' at lev ',ll, 'stopping'
          stop
        endif
        if(debugprint)print*,'sample ',ll,varname,' = ',ll,qqs(isa,jsa,ll)
        
        varname ='grle'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
                qqg(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          print*,'fail to read ', varname,' at lev ',ll, 'stopping'
          stop
        endif
        if(debugprint)print*,'sample ',ll,varname,' = ',ll,qqg(isa,jsa,ll) 
!define cwm
        do j=jsta,jend
          do i=1,im
            cwm(i,j,ll)=qqg(i,j,ll)+qqs(i,j,ll)+qqr(i,j,ll)+qqi(i,j,ll)+qqw(i,j,ll)
          enddo
        enddo

       end if ! end of reading MP species for diff MP options

!       if (iret /= 0)print*,'Error scattering array';stop

!                                              pressure vertical velocity
       if(trim(modelname_nemsio)=='FV3GFS')then
        recn_vvel = 0 ! do not derive omega 
        varname='dzdt'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
              wh(i,j,ll) = tmp(i+js)
            enddo
          enddo
          if(debugprint)print*,'sample l ',varname,' = ',ll,wh(isa,jsa,ll)
        else 
          if(me==0)print*,'fail to read ', varname,' at lev ',ll 
        end if
       else
        varname='vvel'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
              omga(i,j,ll) = tmp(i+js)
            enddo
          enddo
          if(debugprint)print*,'sample l ',varname,' = ',ll,omga(isa,jsa,ll)
        else
          recn_vvel = -9999 
          if(me==0)print*,'fail to read ', varname,' at lev ',ll, &
            'will derive omega later' 
        endif
       end if

!                                              pressure vertical velocity
        varname='delz'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
! make sure delz is positive.
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
              zint(i,j,ll)=zint(i,j,ll+1)+abs(tmp(i+js))
              if(recn_dpres /= -9999)pmid(i,j,ll)=rgas*dpres(i,j,ll)* &
                      t(i,j,ll)*(q(i,j,ll)*fv+1.0)/grav/abs(tmp(i+js)) 
            enddo
          enddo
          if(debugprint)print*,'sample l ',varname,' = ',ll, &
             zint(isa,jsa,ll)
          if(trim(modelname_nemsio)=='FV3GFS' .and. &
           recn_dpres /= -9999)then
            do j=jsta,jend
              js = fldst + (j-jsta)*im
              do i=1,im
                omga(i,j,ll)=(-1.)*wh(i,j,ll)*dpres(i,j,ll)/abs(tmp(i+js)) 
              end do
            end do
            if(debugprint)print*,'sample l omga for FV3',ll, &
              omga(isa,jsa,ll)
          end if
        else
          recn_delz = -9999
          if(me==0)print*,'fail to read ', varname,' at lev ',ll, &
            'will derive height later'
        endif

! cloud fraction
        varname='cld_amt'
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /= 0) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
              cfr(i,j,ll)=tmp(i+js)
            enddo
          enddo
!          if(debugprint)print*,'sample l ',VarName,' = ',ll, &
!             cfr(isa,jsa,ll)
        endif

        if(imp_physics == 99)then
         allocate(p2d(im,lm),t2d(im,lm),q2d(im,lm),cw2d(im,lm), &
         qs2d(im,lm),cfr2d(im,lm))
         do j=jsta,jend
          do k=1,lm
           do i=1,im
            p2d(i,k)  = pmid(i,j,ll)*0.01
            t2d(i,k)  = t(i,j,ll)
            q2d(i,k)  = q(i,j,ll)
            cw2d(i,k) = cwm(i,j,ll)
            es = min(fpvsnew(t(i,j,ll)),pmid(i,j,ll))
            qs2d(i,k) = eps*es/(pmid(i,j,ll)+epsm1*es)!saturation q for GFS
           enddo
          enddo
          call progcld1 &
!...................................
!  ---  inputs:
             ( p2d,t2d,q2d,qs2d,cw2d,im,lm,0,                         &
!  ---  outputs:
               cfr2d                                                  &
              )
!$omp parallel do private(i,k)
          do k=1,lm
           do i=1,im
            cfr(i,j,k) = cfr2d(i,k)
           enddo
          end do
         end do
         deallocate(p2d,t2d,q2d,qs2d,cw2d,cfr2d)
        end if

! With SHOC NEMS/GSM does output TKE now
        varname='tke'
        recn = 0
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /=0 ) then
          fldst = (recn-1)*fldsize
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
              q2(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
          if(me==0)print*,'fail to read ', varname,' at lev ',ll 
!$omp parallel do private(i,j)
          do j=jsta,jend
            do i=1,im
              q2(i,j,ll) = spval
            end do
          end do
        endif
       if(debugprint)print*,'sample l ',varname,' = ',ll,q2(isa,jsa,ll)

! Read model derived radar ref.
       varname='ref3D'
        recn = 0
        call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
        if(recn /=0 ) then
!$omp parallel do private(i,j,js)
          do j=jsta,jend
            js = fldst + (j-jsta)*im
            do i=1,im
              ref_10cm(i,j,ll) = tmp(i+js)
            enddo
          enddo
        else
!$omp parallel do private(i,j)
          do j=jsta,jend
            do i=1,im
              ref_10cm(i,j,ll) = spval
            end do
          end do
          if(me==0)print*,'fail to read ', varname,' at lev ',ll
        endif
       if(debugprint)print*,'sample l ',varname,' = ',ll,ref_10cm(isa,jsa,ll)


      end do ! do loop for l

! construct interface pressure from model top (which is zero) and dp from top down PDTOP
!     pdtop = spval
!      pt    = 0.
!     pd    = spval           ! GFS does not output PD
      pt=ak5(lp1)

      ii = im/2
      jj = (jsta+jend)/2

!!!!! COMPUTE Z, GFS integrates Z on mid-layer instead
!!! use GFS contants to see if height becomes more aggreable to GFS pressure grib file

      if (recn_dpres == -9999) then
        do l=lm,1,-1
!$omp parallel do private(i,j)
          do j=jsta,jend
            do i=1,im
              pint(i,j,l) = ak5(lm+2-l) + bk5(lm+2-l)*pint(i,j,lp1)
              if(recn_delz == -9999)pmid(i,j,l) = 0.5*(pint(i,j,l)+ &
                         pint(i,j,l+1))  ! for now - Moorthi
            enddo
          enddo
        if (me == 0) print*,'sample pint,pmid' ,ii,jj,l,pint(ii,jj,l),pmid(ii,jj,l)
        enddo
      else
! compute pint using dpres from bot up
!        do l=lm,1,-1
!          do j=jsta,jend
!            do i=1,im
!              pint(i,j,l)   = pint(i,j,l+1) - dpres(i,j,l)
!            enddo
!          enddo
!          if (me == 0) print*,'sample model pint,pmid' ,ii,jj,l &
!           ,pint(ii,jj,l),pmid(ii,jj,l)
!        end do

!Feb 6 2018: per Jun, change to compute pint from top down
        do j=jsta,jend
          do i=1,im
            pint(i,j,1)=ak5(lp1)
          end do
        end do

        do l=2,lp1
          do j=jsta,jend
            do i=1,im
              pint(i,j,l)   = pint(i,j,l-1) + dpres(i,j,l-1)
            enddo
          enddo
          if (me == 0) print*,'sample model pint,pmid' ,ii,jj,l &
          ,pint(ii,jj,l)
        end do
      endif

!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!     Compute omega
!sk05132016

      if (recn_vvel == -9999) then !only compute omega if it's not in model output
        allocate(ps2d(im,jsta_2l:jend_2u),    psx2d(im,jsta_2l:jend_2u),  &
                 psy2d(im,jsta_2l:jend_2u))
        allocate(div3d(im,jsta:jend,lm))

!$omp parallel do private(i,j)
        do j=jsta,jend
          do i=1,im
            ps2d(i,j) = log(pint(i,j,lm+1))
          enddo
        enddo
        call calgradps(ps2d,psx2d,psy2d)

        call caldiv(uh, vh, div3d)

!----------------------------------------------------------------------
        allocate (vcrd(lm+1,2),  d2d(im,lm), u2d(im,lm), v2d(im,lm),    &
                  pi2d(im,lm+1), pm2d(im,lm), omga2d(im,lm))
        omga2d=spval
        idvc    = 2
        idsl    = 2
        nvcoord = 2
        do l=1,lm+1
          vcrd(l,1) = vcoord4(l,1,1)
          vcrd(l,2) = vcoord4(l,2,1)
        enddo

!       jtem = jm / 18 + 1
        jtem = jm / 20 + 1
        do j=jsta,jend
          npass = npass2
!         if (j > jm-jtem+1 .or. j < jtem) npass = npass3
          if (j > jm-jtem+1) then
            npass = npass + nint(0.5*(j-jm+jtem-1))
          elseif (j < jtem)  then
            npass = npass + nint(0.5*(jtem-j))
          endif
!         npass = 0
!$omp parallel do private(i,l,ll)
          do l=1,lm
            ll = lm-l+1
            do i=1,im
              u2d(i,l) = uh(i,j,ll) !flip u & v for calling modstuff
              v2d(i,l) = vh(i,j,ll)
              d2d(i,l) = div3d(i,j,ll)
            end do
          end do

          call modstuff2(im,   im, lm, idvc, idsl, nvcoord,             &
                         vcrd, pint(1,j,lp1), psx2d(1,j), psy2d(1,j),   &
                         d2d,  u2d, v2d, pi2d, pm2d, omga2d, me)
!     if (j ==1 .or. j == jm) &
!     write (*,*)' omga2d=',omga2d(1,:),' j=',j

          if (npass <= 0 ) then
!$omp parallel do private(i,l,ll)
            do l=1,lm
              ll = lm-l+1
              do i=1,im
                omga(i,j,l) = omga2d(i,ll)
!               pmid(i,j,l) = pm2d(i,ll)
!               pint(i,j,l) = pi2d(i,ll+1)
              enddo
            enddo
          else
!$omp parallel do private(i,l,ll,nn,omg1,omg2)
            do l=1,lm
              ll = lm-l+1
              do i=1,im
                omg1(i) = omga2d(i,ll)
              enddo
              do nn=1,npass
                do i=1,im
                  omg2(i+1) = omg1(i)
                enddo
                omg2(1)    = omg2(im+1)
                omg2(im+2) = omg2(2)
                do i=2,im+1
                  omg1(i-1) = third * (omg2(i-1) + omg2(i) + omg2(i+1))
                enddo
              enddo
            
              do i=1,im
                omga(i,j,l) = omg1(i)
              enddo
!     if (j ==1 .or. j == jm) &
!     write (1000+me,*)' omga=',omga(:,j,l),' j=',j,' l=',l

            enddo
          endif

!         Average omega for the last point near the poles - moorthi
          if (j ==1 .or. j == jm) then
            tx1 = 1.0 / im
!     write(*,*)' j=',j,' omgamax=',maxval(omga(:,j,:)),' omgamin=',minval(omga(:,j,:))
!$omp parallel do private(i,l,tx2)
            do l=1,lm
              tx2 = 0.0
              do i=1,im
                tx2 = tx2 + omga(i,j,l)
              enddo
              tx2 = tx2 * tx1
              do i=1,im
                omga(i,j,l) = tx2
              enddo
            enddo
          endif

        enddo                  ! end of j loop

        deallocate (vcrd,d2d,u2d,v2d,pi2d,pm2d,omga2d)
        deallocate (ps2d,psx2d,psy2d,div3d)
      end if
      deallocate (vcoord4)

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

!
      allocate(wrk1(im,jsta:jend),wrk2(im,jsta:jend))
      allocate(fi(im,jsta:jend,2))

      do j=jsta,jend
        do i=1,im
          pd(i,j)         = spval           ! GFS does not output PD
          pint(i,j,1)     = pt
        end do
      end do

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

      if (recn_delz == -9999) then !only compute H if it's not in model output
        do j=jsta,jend
          do i=1,im
            wrk1(i,j)       = log(pmid(i,j,lm))
            wrk2(i,j)       = t(i,j,lm)*(q(i,j,lm)*fv+1.0)
            fi(i,j,1)       = fis(i,j)                      &
                          + wrk2(i,j)*rgas*(alpint(i,j,lp1)-wrk1(i,j))
            zmid(i,j,lm)    = fi(i,j,1) * gravi
          end do
        end do 
        
        DO l=lm,2,-1  ! omit computing model top height because it's infinity
          ll = l - 1
          do j = jsta, jend
            do i = 1, im
              tvll          = t(i,j,ll)*(q(i,j,ll)*fv+1.0)
              pmll          = log(pmid(i,j,ll))

              fi(i,j,2)     = fi(i,j,1) + (0.5*rgas)*(wrk2(i,j)+tvll)   &
                                      * (wrk1(i,j)-pmll)
              zmid(i,j,ll)  = fi(i,j,2) * gravi
!
              fact          = (alpint(i,j,l)-wrk1(i,j)) / (pmll-wrk1(i,j))
              zint(i,j,l)   = zmid(i,j,l) +(zmid(i,j,ll)-zmid(i,j,l))*fact
              fi(i,j,1)     = fi(i,j,2)
              wrk1(i,j)     = pmll
              wrk2(i,j)     = tvll
            ENDDO
          ENDDO

          if (me == 0) print*,'L ZINT= ',l,zint(ii,jj,l),                &
           'alpint=',alpint(ii,jj,l),'pmid=',log(pmid(ii,jj,l)),  &
           'pmid(l-1)=',log(pmid(ii,jj,l-1)),'zmd=',zmid(ii,jj,l), &
           'zmid(l-1)=',zmid(ii,jj,l-1)
        ENDDO
        deallocate(wrk1,wrk2,fi)
      else 
        do l=lm,1,-1
          do j=jsta,jend
            do i=1,im
              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))
            end do
          end do
        end do
      end if

!      do j=jsta,jend
!        do i=1,im
!          pd(i,j)         = spval           ! GFS does not output PD
!          pint(i,j,1)     = PT
!          alpint(i,j,lp1) = log(pint(i,j,lp1))
!          wrk1(i,j)       = log(PMID(I,J,LM))
!          wrk2(i,j)       = T(I,J,LM)*(Q(I,J,LM)*fv+1.0)
!          FI(I,J,1)       = FIS(I,J)                      &
!                          + wrk2(i,j)*rgas*(ALPINT(I,J,Lp1)-wrk1(i,j))
!          ZMID(I,J,LM)    = FI(I,J,1) * gravi
!        end do
!      end do
      
! SECOND, INTEGRATE HEIGHT HYDROSTATICLY, GFS integrate height on mid-layer

!      DO L=LM,2,-1  ! omit computing model top height because it's infinity
!        ll = l - 1
!        do j = jsta, jend
!          do i = 1, im
!            alpint(i,j,l) = log(pint(i,j,l))
!            tvll          = T(I,J,LL)*(Q(I,J,LL)*fv+1.0)
!            pmll          = log(PMID(I,J,LL))

!            FI(I,J,2)     = FI(I,J,1) + (0.5*rgas)*(wrk2(i,j)+tvll)   &
!                                      * (wrk1(i,j)-pmll)
!            ZMID(I,J,LL)  = FI(I,J,2) * gravi
!
!            FACT          = (ALPINT(I,J,L)-wrk1(i,j)) / (pmll-wrk1(i,j))
!            ZINT(I,J,L)   = ZMID(I,J,L) +(ZMID(I,J,LL)-ZMID(I,J,L))*FACT
!            FI(I,J,1)     = FI(I,J,2)
!            wrk1(i,J)     = pmll
!            wrk2(i,j)     = tvll
!          ENDDO
!        ENDDO

!        if (me == 0) print*,'L ZINT= ',l,zint(ii,jj,l),                &
!          'alpint=',ALPINT(ii,jj,l),'pmid=',LOG(PMID(Ii,Jj,L)),  &
!          'pmid(l-1)=',LOG(PMID(Ii,Jj,L-1)),'zmd=',ZMID(Ii,Jj,L), &    
!          'zmid(l-1)=',ZMID(Ii,Jj,L-1)
!      ENDDO
!      deallocate(wrk1,wrk2)


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

! GFS output dust in nemsio (GOCART)
        dustcb=0.0
        dustallcb=0.0
!       DUST = SPVAL
        !VarName='du001'
        varname='dust1'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,dust(1:im,jsta_2l:jend_2u,ll,1))

!        if(debugprint)print*,'sample l ',VarName,' = ',ll,dust(isa,jsa,ll,1)
        end do ! do loop for l      
      
        !VarName='du002'
        varname='dust2'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,dust(1:im,jsta_2l:jend_2u,ll,2))

           dustcb(1:im,jsta_2l:jend_2u)=dustcb(1:im,jsta_2l:jend_2u)+ &
           (dust(1:im,jsta_2l:jend_2u,ll,1)+0.38*dust(1:im,jsta_2l:jend_2u,ll,2))* &
           dpres(1:im,jsta_2l:jend_2u,ll)/grav

!         if(debugprint)print*,'sample l ',VarName,' = ',ll,dust(isa,jsa,ll,2)
        end do ! do loop for l 
      
        !VarName='du003'
        varname='dust3'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,dust(1:im,jsta_2l:jend_2u,ll,3))
!         if(debugprint)print*,'sample l ',VarName,' = ',ll,dust(isa,jsa,ll,3)
        end do ! do loop for l 
      
        !VarName='du004'
        varname='dust4'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,dust(1:im,jsta_2l:jend_2u,ll,4))

!         if(debugprint)print*,'sample l ',VarName,' = ',ll,dust(isa,jsa,ll,4)
        end do ! do loop for l 
      
        !VarName='du005'
        varname='dust5'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,dust(1:im,jsta_2l:jend_2u,ll,5))

           dustallcb(1:im,jsta_2l:jend_2u)=dustallcb(1:im,jsta_2l:jend_2u)+ &
           (dust(1:im,jsta_2l:jend_2u,ll,1)+dust(1:im,jsta_2l:jend_2u,ll,2)+ &
           dust(1:im,jsta_2l:jend_2u,ll,3)+0.74*dust(1:im,jsta_2l:jend_2u,ll,4))* &
           dpres(1:im,jsta_2l:jend_2u,ll)/grav

!         if(debugprint)print*,'sample l ',VarName,' = ',ll,dust(isa,jsa,ll,5)
        end do ! do loop for l 
!
! GFS output sea salt in nemsio (GOCART)
        sscb=0.0
        ssallcb=0.0
!       SALT = SPVAL
        !VarName='ss001'
        varname='seas1'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,salt(1:im,jsta_2l:jend_2u,ll,1))

!         if(debugprint)print*,'sample l ',VarName,' = ',ll,salt(isa,jsa,ll,1)
        end do ! do loop for l

        !VarName='ss002'
        varname='seas2'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,salt(1:im,jsta_2l:jend_2u,ll,2))

!         if(debugprint)print*,'sample l ',VarName,' = ',ll,salt(isa,jsa,ll,2)
        end do ! do loop for l

        !VarName='ss003'
        varname='seas3'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,salt(1:im,jsta_2l:jend_2u,ll,3))

            sscb(1:im,jsta_2l:jend_2u)=sscb(1:im,jsta_2l:jend_2u)+ &
        (salt(1:im,jsta_2l:jend_2u,ll,1)+ &
        salt(1:im,jsta_2l:jend_2u,ll,2)+0.83*salt(1:im,jsta_2l:jend_2u,ll,3))* &
           dpres(1:im,jsta_2l:jend_2u,ll)/grav
     
!         if(debugprint)print*,'sample l ',VarName,' = ',ll,salt(isa,jsa,ll,3)
        end do ! do loop for l

        !VarName='ss004'
        varname='seas4'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,salt(1:im,jsta_2l:jend_2u,ll,4))
!         if(debugprint)print*,'sample l ',VarName,' = ',ll,salt(isa,jsa,ll,4)
        end do ! do loop for l

        !VarName='ss005'
        varname='seas5'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,salt(1:im,jsta_2l:jend_2u,ll,5))

          ssallcb(1:im,jsta_2l:jend_2u)=ssallcb(1:im,jsta_2l:jend_2u)+ &
         (salt(1:im,jsta_2l:jend_2u,ll,1)+salt(1:im,jsta_2l:jend_2u,ll,2)+ &
          salt(1:im,jsta_2l:jend_2u,ll,3)+ &
          salt(1:im,jsta_2l:jend_2u,ll,4))* &
           dpres(1:im,jsta_2l:jend_2u,ll)/grav

!         if(debugprint)print*,'sample l ',VarName,' = ',ll,salt(isa,jsa,ll,5)
        end do ! do loop for l

! GFS output black carbon in nemsio (GOCART)
        bccb=0.0
!       SOOT = SPVAL
        !VarName='bcphobic'
        varname='bc1'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,soot(1:im,jsta_2l:jend_2u,ll,1))

!         if(debugprint)print*,'sample l ',VarName,' = ',ll,soot(isa,jsa,ll,1)
        end do ! do loop for l

        !VarName='bcphilic'
        varname='bc2'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,soot(1:im,jsta_2l:jend_2u,ll,2))

            bccb(1:im,jsta_2l:jend_2u)=bccb(1:im,jsta_2l:jend_2u)+ &
        (soot(1:im,jsta_2l:jend_2u,ll,1)+soot(1:im,jsta_2l:jend_2u,ll,2))* &
           dpres(1:im,jsta_2l:jend_2u,ll)/grav

!         if(debugprint)print*,'sample l ',VarName,' = ',ll,soot(isa,jsa,ll,2)
        end do ! do loop for l

        occb=0.0
! GFS output organic carbon in nemsio (GOCART)
!       WASO = SPVAL
        !VarName='ocphobic'
        varname='oc1'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,waso(1:im,jsta_2l:jend_2u,ll,1))

!         if(debugprint)print*,'sample l ',VarName,' = ',ll,waso(isa,jsa,ll,1)
        end do ! do loop for l

        !VarName='ocphilic'
        varname='oc2'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,waso(1:im,jsta_2l:jend_2u,ll,2))

            occb(1:im,jsta_2l:jend_2u)=occb(1:im,jsta_2l:jend_2u)+ &
        (waso(1:im,jsta_2l:jend_2u,ll,1)+waso(1:im,jsta_2l:jend_2u,ll,2)) * &
           dpres(1:im,jsta_2l:jend_2u,ll)/grav

!         if(debugprint)print*,'sample l ',VarName,' = ',ll,waso(isa,jsa,ll,2)
        end do ! do loop for l

! GFS output sulfate in nemsio (GOCART)
        sulfcb=0.0
!       SUSO = SPVAL
        !VarName='so4'
        varname='sulf'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,suso(1:im,jsta_2l:jend_2u,ll,1))

            sulfcb(1:im,jsta_2l:jend_2u)=sulfcb(1:im,jsta_2l:jend_2u)+ &
        suso(1:im,jsta_2l:jend_2u,ll,1)* &
           dpres(1:im,jsta_2l:jend_2u,ll)/grav

!         if(debugprint)print*,'sample l ',VarName,' = ',ll,suso(isa,jsa,ll,1)
        end do ! do loop for l

! GFS output pp25 in nemsio (GOCART)
        pp25cb=0.0
!       PP25 = SPVAL
        !VarName='so4'
        varname='pp25'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,pp25(1:im,jsta_2l:jend_2u,ll,1))
            pp25cb(1:im,jsta_2l:jend_2u)=pp25cb(1:im,jsta_2l:jend_2u)+ &
        pp25(1:im,jsta_2l:jend_2u,ll,1)* &
           dpres(1:im,jsta_2l:jend_2u,ll)/grav
!         if(debugprint)print*,'sample l ',VarName,' =
!         ',ll,suso(isa,jsa,ll,1)
        end do ! do loop for l
! GFS output pp10 in nemsio (GOCART)
        pp10cb=0.0
!       PP10 = SPVAL
        !VarName='so4'
        varname='pp10'
        vcoordname='mid layer'
        do l=1,lm
          ll=lm-l+1
          call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
          ,l,nrec,fldsize,spval,tmp &
          ,recname,reclevtyp,reclev,varname,vcoordname &
          ,pp10(1:im,jsta_2l:jend_2u,ll,1))
            pp10cb(1:im,jsta_2l:jend_2u)=pp10cb(1:im,jsta_2l:jend_2u)+ &
        pp10(1:im,jsta_2l:jend_2u,ll,1)* &
           dpres(1:im,jsta_2l:jend_2u,ll)/grav
!         if(debugprint)print*,'sample l ',VarName,' =
!         ',ll,suso(isa,jsa,ll,1)
        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,n,tv)
          do j=jsta,jend
            do i=1,im
  
              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

            end do
          end do
        end do
             l=lm
!$omp parallel do private(i,j)
          do j=jsta,jend
            do i=1,im
            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)
            pp25cb(i,j) = max(pp25cb(i,j), 0.0)
            pp10cb(i,j) = max(pp10cb(i,j), 0.0)
!      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 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 
!      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)
            end do
          end do
      endif                     ! endif for gocart_on
!
! done with sigma file, close it for now
      call nemsio_close(nfile,iret=status)
      deallocate(tmp,recname,reclevtyp,reclev)

! open flux file 
! has to move it to beginning to read imp_physics
!      call nemsio_open(ffile,trim(fileNameFlux),'read',mpi_comm_comp &
!       ,iret=status)
!      if ( Status /= 0 ) then
!        print*,'error opening ',fileNameFlux, ' Status = ', Status
!        print*,'skip reading of flux file' 
!      endif
      call nemsio_getfilehead(ffile,iret=status,nrec=nrec)
      print*,'nrec for flux file=',nrec
      allocate(recname(nrec),reclevtyp(nrec),reclev(nrec))
      call nemsio_getfilehead(ffile,iret=iret  &
       ,recname=recname ,reclevtyp=reclevtyp,reclev=reclev)
      if(debugprint)then
       if (me == 0)then
         do i=1,nrec
          print *,'recname,reclevtyp,reclev=',trim(recname(i)),' ', &
          trim(reclevtyp(i)),reclev(i)
         end do
       end if
      end if

! IVEGSRC=1 for IGBP, 0 for USGS, 2 for UMD
      varname='IVEGSRC'
      call nemsio_getheadvar(ffile,trim(varname),ivegsrc,iret)
      if (iret /= 0) then
       print*,varname,' not found in file-use 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

      varname='CU_PHYSICS'
      call nemsio_getheadvar(ffile,trim(varname),icu_physics,iret)
      if (iret /= 0) then
       print*,varname," not found in file-Assigned 4 for SAS as default"
       icu_physics=4
      end if
      if (me == 0) print*,'CU_PHYSICS= ',icu_physics

      varname='dtp'
      call nemsio_getheadvar(ffile,trim(varname),dtp,iret)
      if (iret /= 0) then
       print*,varname," not found in file-Assigned 225. for dtp as default"
       dtp=225.
      end if
      if (me == 0) print*,'dtp= ',dtp

! Chuang: zhour is when GFS empties bucket last so using this
! to compute buket will result in changing bucket with forecast time.
! set default bucket for now

!     call nemsio_getheadvar(ffile,'zhour',zhour,iret=iret)
!     if(iret == 0) then
!        tprec   = 1.0*ifhr-zhour
!        tclod   = tprec
!        trdlw   = tprec
!        trdsw   = tprec
!        tsrfc   = tprec
!        tmaxmin = tprec
!        td3d    = tprec
!        print*,'tprec from flux file header= ',tprec
!     else
!        print*,'Error reading accumulation bucket from flux file', &
!            'header - will try to read from env variable FHZER'
!        CALL GETENV('FHZER',ENVAR)
!        read(ENVAR, '(I2)')idum
!        tprec   = idum*1.0
!        tclod   = tprec
!        trdlw   = tprec
!        trdsw   = tprec
!        tsrfc   = tprec
!        tmaxmin = tprec
!        td3d    = tprec
!        print*,'TPREC from FHZER= ',tprec
!     end if


! start reading nemsio flux files using parallel read
      fldsize = (jend-jsta+1)*im
      allocate(tmp(fldsize*nrec))
      print*,'allocate tmp successfully'
      tmp=0.
      call nemsio_denseread(ffile,1,im,jsta,jend,tmp,iret=iret)
! can't stop because anl does not have flux file
!      if(iret/=0)then
!        print*,"fail to read flux file using mpi io read, stopping"
!        stop
!      end if

      vcoordname='sfc'       ! surface fileds
      varname='land'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                  &
                          ,l,nrec,fldsize,spval,tmp                      &
                          ,recname,reclevtyp,reclev,varname,vcoordname,sm)
      if(debugprint)print*,'sample ',varname,' =',sm(im/2,(jsta+jend)/2)

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

! sea ice mask 

      varname='icec'
      vcoordname = 'sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                    &
                          ,l,nrec,fldsize,spval,tmp                        &
                          ,recname,reclevtyp,reclev,varname,vcoordname,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=1,im
          if (sm(i,j) /= spval .and. sm(i,j) == 0.0) sice(i,j) = 0.0
        enddo
      enddo


! PBL height using nemsio
      varname='hpbl'
      vcoordname = 'sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,pblh)
     if(debugprint)print*,'sample ',varname,' = ',pblh(isa,jsa)

! frictional velocity using nemsio
      varname='fricv'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,ustar) 
!     if(debugprint)print*,'sample ',VarName,' = ',ustar(isa,jsa)

! roughness length using getgb
      varname='sfcr'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,z0)
!     if(debugprint)print*,'sample ',VarName,' = ',z0(isa,jsa)

! sfc exchange coeff
      varname='sfexc'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,sfcexc)

! aerodynamic conductance
      varname='acond'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,acond)

! surface potential T  using getgb
      varname='tmp'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,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=1,im
          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) 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)

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

! convective precip in m per physics time step using getgb
! read 6 hour bucket
      varname='cpratb_ave'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,avgcprate)
!     where(avgcprate /= spval)avgcprate=avgcprate*dtq2/1000. ! convert to m
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (avgcprate(i,j) /= spval) avgcprate(i,j) = avgcprate(i,j) * (dtq2*0.001)
!wm          cprate(i,j) = avgcprate(i,j)
        enddo
      enddo
! read continuous bucket
      varname='cprat_ave'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,avgcprate_cont)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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(isa,jsa)
      
!      print*,'maxval CPRATE: ', maxval(CPRATE)

! time averaged bucketed precip rate 
      varname='prateb_ave'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,avgprec)
!     where(avgprec /= spval)avgprec=avgprec*dtq2/1000. ! convert to m
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (avgprec(i,j) /= spval) avgprec(i,j) = avgprec(i,j) * (dtq2*0.001)
        enddo
      enddo

!     if(debugprint)print*,'sample ',VarName,' = ',avgprec(isa,jsa)

! time averaged continuous precip rate in m per physics time step using getgb
      varname='prate_ave'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,avgprec_cont)
!     where(avgprec /= spval)avgprec=avgprec*dtq2/1000. ! convert to m
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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(isa,jsa)

! precip rate in m per physics time step 
      varname='tprcp'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,prec)
!$omp parallel do private(i,j)
! unit of prec and cprate in post is supposed to be m per physics time step
! it will be converted back to kg/m^2/s by multiplying by density in SURFCE
      do j=jsta,jend
        do i=1,im
          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'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,cprate)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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
      if(debugprint)print*,'sample ',varname,' = ',cprate(isa,jsa)
      
! GFS does not have accumulated total, gridscale, and convective precip, will use inst precip to derive in SURFCE.f


! inst snow water eqivalent using nemsio
      varname='weasd'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,sno)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,snoavg)
! snow cover is multipled by 100 in SURFCE before writing it out
      do j=jsta,jend
        do i=1,im
          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'
!     VcoordName='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,si)
!     where(si /= spval)si=si*1000. ! convert to mm
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='tmp'
      vcoordname='2 m above gnd'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,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=1,im
          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='spfh'
      vcoordname='2 m above gnd'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,qshltr)
!     if(debugprint)print*,'sample ',VarName,' = ',qshltr(isa,jsa)
      
! mid day avg albedo in fraction using nemsio
      varname='albdo_ave'
      vcoordname='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,avgalbedo)
!     where(avgalbedo /= spval)avgalbedo=avgalbedo/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (avgalbedo(i,j) /= spval) avgalbedo(i,j) = avgalbedo(i,j) * 0.01
        enddo
      enddo
!     if(debugprint)print*,'sample ',VarName,' = ',avgalbedo(isa,jsa)
     
! time averaged column cloud fractionusing nemsio
      varname='tcdc_ave'
      vcoordname='atmos col'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,avgtcdc)
!     where(avgtcdc /= spval)avgtcdc=avgtcdc/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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=1,im
          czen(i,j)   = spval
          czmean(i,j) = spval      
        enddo
      enddo

! maximum snow albedo in fraction using nemsio
      varname='mxsalb'
      vcoordname='sfc'
!     l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,mxsnal)
!     where(mxsnal /= spval)mxsnal=mxsnal/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (mxsnal(i,j) /= spval) mxsnal(i,j) = mxsnal(i,j) * 0.01
        enddo
      enddo
!     if(debugprint)print*,'sample ',VarName,' = ',mxsnal(isa,jsa)
     
!$omp parallel do private(i,j)
      do j=jsta_2l,jend_2u
        do i=1,im
          radot(i,j) = spval ! GFS does not have inst surface outgoing longwave
        enddo
      enddo

! 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=1,im
          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=1,im
          cfrach(i,j) = spval
          cfracl(i,j) = spval
          cfracm(i,j) = spval
        enddo
      enddo

! ave high cloud fraction using nemsio
      varname='tcdc_ave'
      vcoordname='high cld lay'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,avgcfrach)
!     where(avgcfrach /= spval)avgcfrach=avgcfrach/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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_ave'
      vcoordname='low cld lay'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,avgcfracl)
!     where(avgcfracl /= spval)avgcfracl=avgcfracl/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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_ave'
      vcoordname='mid cld lay'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,avgcfracm)
!     where(avgcfracm /= spval)avgcfracm=avgcfracm/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='tcdc'
      vcoordname='convect-cld laye'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,cnvcfr)
!     where(cnvcfr /= spval)cnvcfr=cnvcfr/100. ! convert to fraction
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,buf)
!     where(buf /= spval)islope=nint(buf) 
!$omp parallel do private(i,j)
      do j = jsta_2l, jend_2u
        do i=1,im
          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 using nemsio
      varname='cnwat'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,cmc)
!     where(cmc /= spval)cmc=cmc/1000. ! convert from kg*m^2 to m
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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=1,im
          grnflx(i,j) = spval ! GFS does not have inst ground heat flux
        enddo
      enddo

! frozen precip fraction using nemsio
      varname='cpofp'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,sr)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='ti'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,ti)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (sice(i,j) == spval .or. sice(i,j) == 0.) ti(i,j)=spval  
        enddo
      enddo

! vegetation fraction in fraction. using nemsio
      varname='veg'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,vegfrc)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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=1,im
          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
 
! liquid volumetric soil mpisture in fraction using nemsio
      varname='soill'
      vcoordname='0-10 cm down'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,sh2o(1,jsta_2l,1))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='soill'
      vcoordname='10-40 cm down'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,sh2o(1,jsta_2l,2))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='soill'
      vcoordname='40-100 cm down'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,sh2o(1,jsta_2l,3))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='soill'
      vcoordname='100-200 cm down'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,sh2o(1,jsta_2l,4))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='soilw'
      vcoordname='0-10 cm down'
      l=1
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,smc(1,jsta_2l,1))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='soilw'
      vcoordname='10-40 cm down'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,smc(1,jsta_2l,2))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='soilw'
      vcoordname='40-100 cm down'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,smc(1,jsta_2l,3))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='soilw'
      vcoordname='100-200 cm down'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,smc(1,jsta_2l,4))
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (sm(i,j) /= 0.0) smc(i,j,4) = spval
        enddo
      enddo
!     if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,4)

! soil temperature using nemsio
      varname='tmp'
      vcoordname='0-10 cm down'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,stc(1,jsta_2l,1))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='tmp'
      vcoordname='10-40 cm down'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,stc(1,jsta_2l,2))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='tmp'
      vcoordname='40-100 cm down'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,stc(1,jsta_2l,3))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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='tmp'
      vcoordname='100-200 cm down'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,stc(1,jsta_2l,4))
!     mask open water areas, combine with sea ice tmp
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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)

!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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
          ssroff(i,j) = spval ! GFS does not have storm runoff
          bgroff(i,j) = spval ! GFS does not have UNDERGROUND RUNOFF
          rlwin(i,j)  = spval  ! GFS does not have inst incoming sfc longwave
          rlwtoa(i,j) = spval ! GFS does not have inst model top outgoing longwave
        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 using nemsio
      varname='dlwrf_ave'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,alwin)

! inst incoming sfc longwave using nemsio
      varname='dlwrf'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,rlwin)
                                                            
! time averaged outgoing sfc longwave using gfsio
      varname='ulwrf_ave'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,alwout)
! inst outgoing sfc longwave using nemsio
      varname='ulwrf'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,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=1,im
          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_ave'
      vcoordname='nom. top' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,alwtoa)
!     if(debugprint)print*,'sample l',VarName,' = ',1,alwtoa(isa,jsa)
      
!$omp parallel do private(i,j)
      do j=jsta_2l,jend_2u
        do i=1,im
          rswin(i,j)  = spval  ! GFS does not have inst incoming sfc shortwave
          rswinc(i,j) = spval  ! GFS does not have inst incoming clear sky sfc shortwave 
          rswout(i,j) = spval  ! GFS does not have inst outgoing sfc shortwave
        enddo
      enddo
           
! GFS incoming sfc longwave has been averaged, set ARDLW to 1
      ardsw=1.0
!     trdsw=6.0

! time averaged incoming sfc shortwave using gfsio
      varname='dswrf_ave'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,aswin)
!     if(debugprint)print*,'sample l',VarName,' = ',1,aswin(isa,jsa)

! inst incoming sfc shortwave using nemsio
      varname='dswrf'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,rswin)

! time averaged incoming sfc uv-b using getgb
      varname='duvb_ave'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,auvbin)
!     if(debugprint)print*,'sample l',VarName,' = ',1,auvbin(isa,jsa)
       
! time averaged incoming sfc clear sky uv-b using getgb
      varname='cduvb_ave'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,auvbinc)
!     if(debugprint)print*,'sample l',VarName,' = ',1,auvbinc(isa,jsa)
      
! time averaged outgoing sfc shortwave using gfsio
      varname='uswrf_ave'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,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=1,im
          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'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,rswout)

! time averaged model top incoming shortwave
      varname='dswrf_ave'
      vcoordname='nom. top'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,aswintoa)
      
!     if(debugprint)print*,'sample l',VarName,' = ',1,aswintoa(isa,jsa)      

! time averaged model top outgoing shortwave
      varname='uswrf_ave'
      vcoordname='nom. top' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,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'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,sfcshx)
!     where (sfcshx /= spval)sfcshx=-sfcshx
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,twbs)
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,sfclhx)
!     where (sfclhx /= spval)sfclhx=-sfclhx
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,qwbs)
!     where (sfclhx /= spval)sfclhx=-sfclhx
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (qwbs(i,j) /= spval) qwbs(i,j) = -qwbs(i,j)
        enddo
      enddo

! time averaged ground heat flux using nemsio
      varname='gflux_ave'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,subshx)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,grnflx)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,sfcux)
!     if(debugprint)print*,'sample l',VarName,' = ',1,sfcux(isa,jsa)
      
! time averaged meridional momentum flux using nemsio
      varname='vflx_ave'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,sfcvx)
!     if(debugprint)print*,'sample l',VarName,' = ',1,sfcvx(isa,jsa)
     
!$omp parallel do private(i,j)
      do j=jsta_2l,jend_2u
        do i=1,im
!          snopcx(i,j)  =spval ! GFS does not have snow phase change heat flux
          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'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,gtaux)

!     if(debugprint)print*,'sample l',VarName,' = ',1,gtaux(isa,jsa)

! time averaged meridional gravity wave stress using getgb
      varname='v-gwd_ave'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,gtauy)
!     if(debugprint)print*,'sample l',VarName,' = ',1,gtauy(isa,jsa)
                                                     
! time averaged accumulated potential evaporation
      varname='pevpr_ave'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,avgpotevp)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,potevp)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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=1,im
! 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='ugrd'
      vcoordname='10 m above gnd' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,u10)

      do j=jsta,jend
        do i=1,im
          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='vgrd'
      vcoordname='10 m above gnd' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,v10)

      do j=jsta,jend
        do i=1,im
          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='vgtyp'
!Use for fv3 model output
      varname='vtype'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,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=1,im
          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'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,buf)
!     where (buf /= spval)
!      isltyp=nint(buf)
!     elsewhere
!      isltyp=0 !need to feed reasonable value to crtm
!     end where 
!$omp parallel do private(i,j)
      do j = jsta_2l, jend_2u
        do i=1,im
          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)
      
!$omp parallel do private(i,j)
      do j=jsta_2l,jend_2u
        do i=1,im
          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
          acsnow(i,j) = spval    ! GFS does not have averaged accumulated snow
          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=1,im
            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
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,ptop)
!     if(debugprint)print*,'sample l',VarName,' = ',1,ptop(isa,jsa)
      
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          htop(i,j) = spval
          if(ptop(i,j) <= 0.0) ptop(i,j) = spval
        enddo
      enddo
      do j=jsta,jend
        do i=1,im
          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='pres'
      vcoordname='convect-cld bot' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,pbot)
!     if(debugprint)print*,'sample l',VarName,VcoordName,' = ',1,pbot(isa,jsa)
      
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          hbot(i,j) = spval
          if(pbot(i,j) <= 0.0) pbot(i,j) = spval
        enddo
      enddo
      do j=jsta,jend
        do i=1,im
!     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

! retrieve time averaged low cloud top pressure using nemsio
      varname='pres_ave'
      vcoordname='low cld top' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,ptopl)
!     if(debugprint)print*,'sample l',VarName,' = ',1,ptopl(isa,jsa)

! retrieve time averaged low cloud bottom pressure using nemsio
      varname='pres_ave'
      vcoordname='low cld bot' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,pbotl)
!     if(debugprint)print*,'sample l',VarName,' = ',1,pbotl(isa,jsa)
     
! retrieve time averaged low cloud top temperature using nemsio
      varname='tmp_ave'
      vcoordname='low cld top' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,ttopl)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopl(isa,jsa)

! retrieve time averaged middle cloud top pressure using nemsio
      varname='pres_ave'
      vcoordname='mid cld top' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,ptopm)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptopm(isa,jsa)
                                                             
! retrieve time averaged middle cloud bottom pressure using  nemsio
      varname='pres_ave'
      vcoordname='mid cld bot' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,pbotm)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pbotm(isa,jsa)
      
! retrieve time averaged middle cloud top temperature using nemsio
      varname='tmp_ave'
      vcoordname='mid cld top' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,ttopm)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopm(isa,jsa)
      
! retrieve time averaged high cloud top pressure using nemsio *********
      varname='pres_ave'
      vcoordname='high cld top' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,ptoph)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptoph(isa,jsa)
     
! retrieve time averaged high cloud bottom pressure using  nemsio
      varname='pres_ave'
      vcoordname='high cld bot' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,pboth)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pboth(isa,jsa)

! retrieve time averaged high cloud top temperature using nemsio
      varname='tmp_ave'
      vcoordname='high cld top' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,ttoph)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttoph(isa,jsa)
      
! retrieve boundary layer cloud cover using nemsio
      varname='tcdc_ave'
      vcoordname='bndary-layer cld' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,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=1,im
          if (pblcfr(i,j) < spval) pblcfr(i,j) = pblcfr(i,j) * 0.01
        enddo
      enddo
        
! retrieve cloud work function using nemsio
      varname='cwork_ave'
      vcoordname='atmos col' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,cldwork)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,cldwork(isa,jsa)
      
! retrieve water runoff using nemsio
      varname='watr_acc'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,runoff)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (sm(i,j) /= 0.0) runoff(i,j) = spval
        enddo
      enddo
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,runoff(isa,jsa)
      
! retrieve shelter max temperature using nemsio
      varname='tmax_max'
      vcoordname='2 m above gnd' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,maxtshltr)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,maxtshltr(isa,jsa)

! retrieve shelter min temperature using nemsio
      varname='tmin_min'
      vcoordname='2 m above gnd' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,mintshltr)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
!     1,mintshltr(im/2,(jsta+jend)/2)
 
!$omp parallel do private(i,j)
      do j=jsta_2l,jend_2u
        do i=1,im
          maxrhshltr(i,j) = spval
          minrhshltr(i,j) = spval
        enddo
      enddo
      
! retrieve ice thickness using nemsio
      varname='icetk'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,dzice)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,dzice(isa,jsa)

! retrieve wilting point using nemsio
      varname='wilt'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,smcwlt)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,suntime)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,suntime(isa,jsa)

! retrieve field capacity using nemsio
      varname='fldcp'
      vcoordname='sfc' 
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u                &
                          ,l,nrec,fldsize,spval,tmp                    &
                          ,recname,reclevtyp,reclev,varname,vcoordname &
                          ,fieldcapa)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          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'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,avisbeamswin)

! retrieve time averaged surface visible diffuse downward solar flux
      varname='vddsf_ave'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,avisdiffswin)

! retrieve time averaged surface near IR beam downward solar flux
      varname='nbdsf_ave'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,airbeamswin)

! retrieve time averaged surface near IR diffuse downward solar flux
      varname='nddsf_ave'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,airdiffswin)

! retrieve time averaged surface clear sky outgoing LW
      varname='csulf'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,alwoutc)

! retrieve time averaged TOA clear sky outgoing LW
      varname='csulf'
      vcoordname='nom. top'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,alwtoac)

! retrieve time averaged surface clear sky outgoing SW
      varname='csusf'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,aswoutc)

! retrieve time averaged TOA clear sky outgoing LW
      varname='csusf'
      vcoordname='nom. top'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,aswtoac)

! retrieve time averaged surface clear sky incoming LW
      varname='csdlf'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,alwinc)

! retrieve time averaged surface clear sky incoming SW
      varname='csdsf'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,aswinc)

! retrieve shelter max specific humidity using nemsio
      varname='spfhmax_max'
      vcoordname='2 m above gnd'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,maxqshltr)
!     if(debugprint)print*,'sample l',VcoordName,VarName,' = ',
!     1,maxtshltr(isa,jsa)

! retrieve shelter min temperature using nemsio
      varname='spfhmin_min'
      vcoordname='2 m above gnd'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,minqshltr)

! retrieve storm runoff using nemsio
      varname='ssrun_acc'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,ssroff)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (sm(i,j) /= 0.0) ssroff(i,j) = spval
        enddo
      enddo

! retrieve direct soil evaporation
      varname='evbs_ave'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,avgedir)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (sm(i,j) /= 0.0) avgedir(i,j) = spval
        enddo
      enddo

! retrieve CANOPY WATER EVAP 
      varname='evcw_ave'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,avgecan)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (sm(i,j) /= 0.0) avgecan(i,j) = spval
        enddo
      enddo

! retrieve PLANT TRANSPIRATION 
      varname='trans_ave'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,avgetrans)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (sm(i,j) /= 0.0) avgetrans(i,j) = spval
        enddo
      enddo

! retrieve snow sublimation
      varname='sbsno_ave'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,avgesnow)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (sm(i,j)==1.0 .and. sice(i,j)==0.) avgesnow(i,j)=spval
        enddo
      enddo

! retrive total soil moisture
      varname='soilm'
      vcoordname='0-200 cm down'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,smstot)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (sm(i,j) /= 0.0) smstot(i,j) = spval
        enddo
      enddo

! retrieve snow phase change heat flux
      varname='snohf'
      vcoordname='sfc'
      l=1
      call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
      ,l,nrec,fldsize,spval,tmp &
      ,recname,reclevtyp,reclev,varname,vcoordname &
      ,snopcx)
!     mask water areas
!$omp parallel do private(i,j)
      do j=jsta,jend
        do i=1,im
          if (sm(i,j) /= 0.0) snopcx(i,j) = spval
        enddo
      enddo
      
! GFS does not have deep convective cloud top and bottom fields

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

      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'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,duem(1,jsta_2l,k))
!     if(debugprint)print*,'sample ',VarName,' = ',duem(isa,jsa,k)
      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='dsut5sd'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,dusd(1,jsta_2l,k))
!      if(debugprint)print*,'sample ',VarName,' = ',dusd(isa,jsa,k)
      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'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,dudp(1,jsta_2l,k))
        print *,'dudp,ck=',maxval(dudp(1:im,jsta:jend,k)), &
                 minval(dudp(1:im,jsta:jend,k))
!      if(debugprint)print*,'sample ',VarName,' = ',dudp(isa,jsa,k)
      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'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,duwt(1,jsta_2l,k))
      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'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,dusv(1,jsta_2l,k))
      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'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,ssem(1,jsta_2l,k))
      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'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,sssd(1,jsta_2l,k))
      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'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,ssdp(1,jsta_2l,k))
      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'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,sswt(1,jsta_2l,k))
      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'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,sssv(1,jsta_2l,k))
      enddo

! retrieve bc emission fluxes
      do k = 1, nbin_bc
       if ( k == 1) varname='bceman'
       if ( k == 2) varname='bcembb'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,bcem(1,jsta_2l,k))
      enddo

! retrieve bc sedimentation fluxes
      do k = 1, nbin_bc
       if ( k == 1) varname='bc1sd'
       if ( k == 2) varname='bc2sd'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,bcsd(1,jsta_2l,k))
      enddo

! retrieve bc dry deposition fluxes
      do k = 1, nbin_bc
       if ( k == 1) varname='bc1dp'
       if ( k == 2) varname='bc2dp'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,bcdp(1,jsta_2l,k))
      enddo

! retrieve bc large wet deposition fluxes
      do k = 1, nbin_bc
       if ( k == 1) varname='bc1wtl'
       if ( k == 2) varname='bc2wtl'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,bcwt(1,jsta_2l,k))
      enddo

! retrieve bc convective wet deposition fluxes
      do k = 1, nbin_bc
       if ( k == 1) varname='bc1wtc'
       if ( k == 2) varname='bc2wtc'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,bcsv(1,jsta_2l,k))
      enddo

! retrieve oc emission fluxes
      do k = 1, nbin_oc
       if ( k == 1) varname='oceman'
       if ( k == 2) varname='ocembb'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,ocem(1,jsta_2l,k))
      enddo

! retrieve oc sedimentation fluxes
      do k = 1, nbin_oc
       if ( k == 1) varname='oc1sd'
       if ( k == 2) varname='oc2sd'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,ocsd(1,jsta_2l,k))
      enddo

! retrieve oc dry deposition fluxes
      do k = 1, nbin_oc
       if ( k == 1) varname='oc1dp'
       if ( k == 2) varname='oc2dp'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,ocdp(1,jsta_2l,k))
      enddo

! retrieve oc large wet deposition fluxes
      do k = 1, nbin_oc
       if ( k == 1) varname='oc1wtl'
       if ( k == 2) varname='oc2wtl'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,ocwt(1,jsta_2l,k))
      enddo

! retrieve oc convective wet deposition fluxes
      do k = 1, nbin_oc
       if ( k == 1) varname='oc1wtc'
       if ( k == 2) varname='oc2wtc'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,ocsv(1,jsta_2l,k))
      enddo

! retrieve MIE AOD
       varname='maod'
       vcoordname='sfc'
       l=1
       call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u               &
                           ,l,nrec,fldsize,spval,tmp                   &
                           ,recname,reclevtyp,reclev,varname,vcoordname&
                           ,maod(1,jsta_2l))

       endif
! done with flux file, close it for now
      call nemsio_close(ffile,iret=status)
      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
        END IF
      end if
!     
!

      RETURN
      END
      subroutine rg2gg(im,jm,numi,a)
!
      implicit none
!
      integer,intent(in):: im,jm,numi(jm)
      real,intent(inout):: a(im,jm)
      integer j,ir,ig
      real r,t(im)
      do j=1,jm
        r = real(numi(j))/real(im)
        do ig=1,im
          ir    = mod(nint((ig-1)*r),numi(j)) + 1
          t(ig) = a(ir,j)
        enddo
        do ig=1,im
          a(ig,j) = t(ig)
        enddo
      enddo
      end subroutine rg2gg
      subroutine gg2rg(im,jm,numi,a)
!
      implicit none
!
      integer,intent(in):: im,jm,numi(jm)
      real,intent(inout):: a(im,jm)
      integer j,ir,ig
      real r,t(im)
      do j=1,jm
        r = real(numi(j))/real(im)
        do ir=1,numi(j)
          ig    = nint((ir-1)/r) + 1
          t(ir) = a(ig,j)
        enddo
        do ir=1,numi(j)
          a(ir,j) = t(ir)
        enddo
      enddo
      end subroutine gg2rg

      subroutine uninterpred(iord,kmsk,lonsperlat,lonr,latr,fi,f)
!!
      implicit none
!!
      integer, intent(in)  :: iord, lonr, latr
      integer, intent(in)  :: kmsk(lonr,latr), lonsperlat(latr)
      real,    intent(in)  :: fi(lonr,latr)
      real,    intent(out) :: f(lonr,latr)
      integer j,lons
!!
!!$omp parallel do private(j,lons)
      do j=1,latr
        lons = lonsperlat(j)
        if(lons /= lonr) then
          call intlon(iord,1,lons,lonr,kmsk(1,j),fi(1,j),f(1,j))
        else
          f(:,j) = fi(:,j)
        endif
      enddo
      end subroutine
      subroutine intlon(iord,imsk,m1,m2,k1,f1,f2)
      implicit none
      integer,intent(in) :: iord,imsk,m1,m2
      integer,intent(in) :: k1(m1)
      real,   intent(in) :: f1(m1)
      real,   intent(out):: f2(m2)
      integer i2,in,il,ir
      real  r,x1
      r = real(m1)/real(m2)
      do i2=1,m2
         x1 = (i2-1)*r
         il = int(x1)+1
         ir = mod(il,m1)+1
          if(iord == 2 .and. (imsk == 0 .or. k1(il) == k1(ir))) then
            f2(i2) = f1(il)*(il-x1) + f1(ir)*(x1-il+1)
          else
            in = mod(nint(x1),m1) + 1
            f2(i2) = f1(in)
          endif
      enddo
      end subroutine intlon