Merge remote-tracking branch 'upstream/develop' into optimize6

src/enkf/gridinfo_gfs.f90

@@ -67,6 +67,7 @@ module gridinfo
 character(len=max_varname_length),public, dimension(13) :: vars3d_supported = (/'u   ', 'v   ', 'tv  ', 'q   ', 'oz  ', 'cw  ', 'tsen', 'prse', &
                                                                                 'ql  ', 'qi  ', 'qr  ', 'qs  ', 'qg  '/) 
 character(len=max_varname_length),public, dimension(13)  :: vars2d_supported = (/'ps ', 'pst', 'sst', 't2m', 'q2m', 'st1', 'st2', 'st3', 'st4', 'sl1', 'sl2', 'sl3', 'sl4' /)
+character(len=max_varname_length),public, dimension(8)  :: vars2d_landonly = (/'st1', 'st2', 'st3', 'st4', 'sl1', 'sl2', 'sl3', 'sl4' /)
 ! supported variable names in anavinfo
 contains
 

src/enkf/gridio_gfs.f90

@@ -926,6 +926,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,
         endif  ! use_full_hydro
      enddo
   else if (use_gfs_ncio) then
+     clip=tiny_r_kind
      call read_vardata(dset, 'ugrd', ug3d,errcode=iret)
      if (iret /= 0) then
         print *,'error reading ugrd'
@@ -1203,36 +1204,36 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,
         call copytogrdin(ug,grdin(:,levels(n3d) + soilt4_ind,nb,ne))
      endif
      if (slc1_ind > 0) then
-        call read_vardata(dset_sfc, 'slc1', values_2d, errcode=iret)
+        call read_vardata(dset_sfc, 'soill1', values_2d, errcode=iret)
         if (iret /= 0) then
-                print *,'error reading slc1'
+                print *,'error reading soill1'
                 call stop2(22)
         endif
         ug = reshape(values_2d,(/nlons*nlats/))
         call copytogrdin(ug,grdin(:,levels(n3d) + slc1_ind,nb,ne))
      endif
      if (slc2_ind > 0) then
-        call read_vardata(dset_sfc, 'slc2', values_2d, errcode=iret)
+        call read_vardata(dset_sfc, 'soill2', values_2d, errcode=iret)
         if (iret /= 0) then
-                print *,'error reading slc2'
+                print *,'error reading soill2'
                 call stop2(22)
         endif
         ug = reshape(values_2d,(/nlons*nlats/))
         call copytogrdin(ug,grdin(:,levels(n3d) + slc2_ind,nb,ne))
      endif
      if (slc3_ind > 0) then
-        call read_vardata(dset_sfc, 'slc3', values_2d, errcode=iret)
+        call read_vardata(dset_sfc, 'soill3', values_2d, errcode=iret)
         if (iret /= 0) then
-                print *,'error reading slc3'
+                print *,'error reading soill3'
                 call stop2(22)
         endif
         ug = reshape(values_2d,(/nlons*nlats/))
         call copytogrdin(ug,grdin(:,levels(n3d) + slc3_ind,nb,ne))
      endif
      if (slc4_ind > 0) then
-        call read_vardata(dset_sfc, 'slc4', values_2d, errcode=iret)
+        call read_vardata(dset_sfc, 'soill4', values_2d, errcode=iret)
         if (iret /= 0) then
-                print *,'error reading slc2'
+                print *,'error reading soill4'
                 call stop2(22)
         endif
         ug = reshape(values_2d,(/nlons*nlats/))
@@ -2122,6 +2123,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n
   character(nemsio_charkind) :: field
   character(len=nf90_max_name) :: time_units
   logical :: hasfield
+  character(len=max_varname_length), dimension(n3d) :: no_vars3d
 
   real(r_kind) kap,kapr,kap1,clip
   real(r_single) compress_err
@@ -2143,10 +2145,12 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n
 
   call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file)
 
-  if (write_sfc_file .and. nproc==0 ) then
+  if (write_sfc_file ) then
         ! adding the sfc increments requires adjusting several other variables. This is done is a separate
         ! program.
-        write(6,*)'gridio/writegriddata: not coded to write sfc analysis, use separate add_incr program instead'
+        if (nproc == 0) write(6,*)'gridio/writegriddata: not coded to write sfc analysis, will write increment for sfc fields'
+        no_vars3d=''
+        call  writeincrement(nanal1,nanal2,no_vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag)
   endif
 
   nocompress = .true.
@@ -3584,7 +3588,7 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,
   integer :: ql_ind, qi_ind, qr_ind, qs_ind, qg_ind
 
   ! netcdf things
-  integer(i_kind) :: dimids3(3), ncstart(3), nccount(3)
+  integer(i_kind) :: dimids3(3), ncstart(3), nccount(3), dimids2(2) 
   integer(i_kind) :: ncid_out, lon_dimid, lat_dimid, lev_dimid, ilev_dimid
   integer(i_kind) :: lonvarid, latvarid, levvarid, pfullvarid, ilevvarid, &
                      hyaivarid, hybivarid, uvarid, vvarid, delpvarid, delzvarid, &
@@ -3615,14 +3619,14 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,
 
   call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file)
 
-  if ( write_atm_file) then
   use_full_hydro = .false.
   clip = tiny_r_kind
   read(datestring,*) iadateout
 
   ncstart = (/1, 1, 1/)
   nccount = (/nlons, nlats, nlevs/)
 
+  if ( write_atm_file) then
   ne = 0
   ensmemloop: do nanal=nanal1,nanal2
   ne = ne + 1
@@ -3978,20 +3982,21 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,
      ! create dimensions based on analysis resolution, not guess
      call nccheck_incr(nf90_def_dim(ncid_out, "longitude", nlons, lon_dimid))
      call nccheck_incr(nf90_def_dim(ncid_out, "latitude", nlats, lat_dimid))
+     dimids2 = (/ lon_dimid, lat_dimid /)
      ! create variables
      call nccheck_incr(nf90_def_var(ncid_out, "longitude", nf90_real, (/lon_dimid/), lonvarid))
      call nccheck_incr(nf90_def_var(ncid_out, "latitude", nf90_real, (/lat_dimid/), latvarid))
-     call nccheck_incr(nf90_def_var(ncid_out, "tmp2m_inc", nf90_real, dimids3(1:2), tmp2mvarid))
-     call nccheck_incr(nf90_def_var(ncid_out, "spfh2m_inc", nf90_real, dimids3(1:2), spfh2mvarid))
-     call nccheck_incr(nf90_def_var(ncid_out, "soilt1_inc", nf90_real, dimids3(1:2), soilt1varid))
-     call nccheck_incr(nf90_def_var(ncid_out, "soilt2_inc", nf90_real, dimids3(1:2), soilt2varid))
-     call nccheck_incr(nf90_def_var(ncid_out, "soilt3_inc", nf90_real, dimids3(1:2), soilt3varid))
-     call nccheck_incr(nf90_def_var(ncid_out, "soilt4_inc", nf90_real, dimids3(1:2), soilt4varid))
-     call nccheck_incr(nf90_def_var(ncid_out, "slc1_inc", nf90_real, dimids3(1:2), slc1varid))
-     call nccheck_incr(nf90_def_var(ncid_out, "slc2_inc", nf90_real, dimids3(1:2), slc2varid))
-     call nccheck_incr(nf90_def_var(ncid_out, "slc3_inc", nf90_real, dimids3(1:2), slc3varid))
-     call nccheck_incr(nf90_def_var(ncid_out, "slc4_inc", nf90_real, dimids3(1:2), slc4varid))
-     call nccheck_incr(nf90_def_var(ncid_out, "soilsnow_mask", nf90_int, dimids3(1:2), maskvarid))
+     call nccheck_incr(nf90_def_var(ncid_out, "tmp2m_inc", nf90_real, dimids2, tmp2mvarid))
+     call nccheck_incr(nf90_def_var(ncid_out, "spfh2m_inc", nf90_real, dimids2, spfh2mvarid))
+     call nccheck_incr(nf90_def_var(ncid_out, "soilt1_inc", nf90_real, dimids2, soilt1varid))
+     call nccheck_incr(nf90_def_var(ncid_out, "soilt2_inc", nf90_real, dimids2, soilt2varid))
+     call nccheck_incr(nf90_def_var(ncid_out, "soilt3_inc", nf90_real, dimids2, soilt3varid))
+     call nccheck_incr(nf90_def_var(ncid_out, "soilt4_inc", nf90_real, dimids2, soilt4varid))
+     call nccheck_incr(nf90_def_var(ncid_out, "slc1_inc", nf90_real, dimids2, slc1varid))
+     call nccheck_incr(nf90_def_var(ncid_out, "slc2_inc", nf90_real, dimids2, slc2varid))
+     call nccheck_incr(nf90_def_var(ncid_out, "slc3_inc", nf90_real, dimids2, slc3varid))
+     call nccheck_incr(nf90_def_var(ncid_out, "slc4_inc", nf90_real, dimids2, slc4varid))
+     call nccheck_incr(nf90_def_var(ncid_out, "soilsnow_mask", nf90_int, dimids2, maskvarid))
      ! place global attributes to serial calc_increment output
      call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "source", "GSI EnKF"))
      call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "comment", &
@@ -4036,7 +4041,7 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,
      ! note: same logic/threshold used in global_cycle to produce
      ! mask on model grid.
 
-     call read_vardata(dsfg, 'slc1', values_2d, errcode=iret)
+     call read_vardata(dsfg, 'soill1', values_2d, errcode=iret)
 
      mask = 0
      do j=1,nlats

src/enkf/inflation.f90

@@ -77,7 +77,8 @@ module inflation
 use constants, only: one, zero, rad2deg, deg2rad
 use covlocal, only: latval, taper
 use controlvec, only: ncdim, cvars3d, cvars2d, nc3d, nc2d, clevels
-use gridinfo, only: latsgrd, logp, npts, nlevs_pres
+! note: vars2d_landonly currently only defined for gridio_gfs, but smoothing only coded for gfs.
+use gridinfo, only: latsgrd, logp, npts, nlevs_pres, vars2d_landonly
 use loadbal, only: indxproc, numptsperproc, npts_max, anal_chunk, anal_chunk_prior
 use smooth_mod, only: smooth
 
@@ -101,9 +102,10 @@ subroutine inflate_ens()
 real(r_single),dimension(ndiag) :: sumcoslat,suma,suma2,sumi,sumf,sumitot,sumatot, &
      sumcoslattot,suma2tot,sumftot
 real(r_single) fnanalsml,coslat
-integer(i_kind) i,nn,iunit,ierr,nb,nnlvl,ps_ind
+integer(i_kind) i,nn,iunit,ierr,nb,nnlvl,ps_ind, this_ind, ind
+integer(i_kind), dimension(8) :: soil_index
 character(len=500) filename
-real(r_single), allocatable, dimension(:,:) :: tmp_chunk2,covinfglobal
+real(r_single), allocatable, dimension(:,:) :: tmp_chunk2,covinfglobal,store_presmooth
 real(r_single) r
 
 fnanalsml = one/(real(nanals-1,r_single))
@@ -231,7 +233,33 @@ subroutine inflate_ens()
    do nn=1,ncdim
      call mpi_allreduce(mpi_in_place,covinfglobal(1,nn),npts,mpi_real4,mpi_sum,mpi_comm_world,ierr)
    enddo
+   ! do not apply smoothing to soil temp. or soil moisture (not globally defined)
+
+   ind = 0
+   do i = 1,8
+       this_ind = getindex(cvars2d, vars2d_landonly(i))
+       if (this_ind>0) then
+            ind=ind+1
+            soil_index(ind)=this_ind
+       endif
+   enddo
+
+   if (ind>0) then
+        allocate(store_presmooth(npts,ind))
+        do i = 1, ind
+           store_presmooth(:,i) = covinfglobal(:,clevels(nc3d)+soil_index(i))
+        enddo
+   endif
+
    call smooth(covinfglobal)
+
+   if (ind>0) then
+        do i = 1, ind
+           covinfglobal(:,clevels(nc3d) + soil_index(i))  = store_presmooth(:,i)
+        enddo
+        deallocate(store_presmooth)
+   endif
+
    where (covinfglobal < covinflatemin) covinfglobal = covinflatemin
    where (covinfglobal > covinflatemax) covinfglobal = covinflatemax
    do i=1,numptsperproc(nproc+1)

src/enkf/readconvobs.f90

@@ -24,7 +24,6 @@ module readconvobs
 !                       reflectivity and radial velocity assimilation. POC: xuguang.wang@ou.edu
 !   2017-12-13  shlyaeva - added netcdf diag read/write capability
 !   2019-03-21  CAPS(C. Tong) - added direct reflectivity DA capability
-!   2022-03-23  draper - added option to not scale qobs by forecast qsat.
 !
 ! attributes:
 !   language: f95

src/enkf/statevec.f90

@@ -136,7 +136,7 @@ subroutine init_statevec()
 do i = 1, ns2d
   if (getindex(vars2d_supported, svars2d(i))<0) then
     if (nproc .eq. 0) then
-      print *,'Error: 2D variable ', svars2d(i), ' is not supported in current version.'
+      print *,'Error: state 2D variable ', svars2d(i), ' is not supported in current version.'
       print *,'Supported variables: ', vars2d_supported
     endif
     call stop2(502)
@@ -145,7 +145,7 @@ subroutine init_statevec()
 do i = 1, ns3d
   if (getindex(vars3d_supported, svars3d(i))<0) then
     if (nproc .eq. 0) then 
-       print *,'Error: 3D variable ', svars3d(i), ' is not supported in current version.'
+       print *,'Error: state 3D variable ', svars3d(i), ' is not supported in current version.'
        print *,'Supported variables: ', vars3d_supported
     endif
     call stop2(502)

src/gsi/gsi_rfv3io_mod.f90

@@ -25,6 +25,7 @@ module gsi_rfv3io_mod
 !   2022-08-10  Wang    - add IO for regional FV3-SMOKE (RRFS-SMOKE) model 
 !   2023-07-30  Zhao    - add IO for the analysis of the significant wave height
 !                         (SWH, aka howv in GSI) in fv3-lam based DA (eg., RRFS-3DRTMA)
+!   2024-01-24  X.Zhang - bug fix for reading the soil temp and mois from the wram start file 
 !
 ! subroutines included:
 !   sub gsi_rfv3io_get_grid_specs
@@ -118,6 +119,11 @@ module gsi_rfv3io_mod
       vgsiname = [character(len=max_varname_length) :: &
         'u','v','w','tsen','delp','q','oz','ql','qi','qr','qs','qg','qnr','dbz','fed','ps','delzinc', &
         aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk',aeronames_smoke_fv3]
+
+  integer(i_kind) ,parameter:: nnonnegtracer=7
+  character(len=max_varname_length), dimension(nnonnegtracer), parameter :: &
+    vnames_nonnegativetracers = [character(len=max_varname_length) :: &
+      "sphum","o3mr","liq_wat","ice_wat","rainwat","snowwat","graupel"]
   character(len=max_varname_length),dimension(:),allocatable:: name_metvars2d
   character(len=max_varname_length),dimension(:),allocatable:: name_metvars3d
   character(len=max_varname_length),dimension(:),allocatable:: name_chemvars3d
@@ -2043,6 +2049,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
     integer(i_kind),allocatable :: gfile_loc_layout(:)
     character(len=180)  :: filename_layout
 
+! for sfc 2d vaiable exist or not
+    logical, dimension(n2d) :: sfc_var_exist
+
     sfcdata= fv3filenamegin%sfcdata
     dynvars= fv3filenamegin%dynvars
 
@@ -2054,6 +2063,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
 !-- initialisation of the array for howv
     sfcn2d(:,:,k_howv) = zero
 
+!-- initialisation of the array for sfc_var_exist 
+    sfc_var_exist = .false.
+
     if(mype==mype_2d ) then
        allocate(sfc_fulldomain(nx,ny))
 
@@ -2104,30 +2116,43 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
           iret=nf90_inquire_variable(gfile_loc,i,name,len)
           if( trim(name)=='f10m'.or.trim(name)=='F10M' ) then
              k=k_f10m
+             sfc_var_exist(k) = .true.
           else if( trim(name)=='stype'.or.trim(name)=='STYPE' ) then
              k=k_stype
+             sfc_var_exist(k) = .true.
           else if( trim(name)=='vfrac'.or.trim(name)=='VFRAC' ) then
              k=k_vfrac
+             sfc_var_exist(k) = .true.
           else if( trim(name)=='vtype'.or.trim(name)=='VTYPE' ) then
              k=k_vtype
+             sfc_var_exist(k) = .true.
           else if( trim(name)=='zorl'.or.trim(name)=='ZORL' ) then
              k=k_zorl
+             sfc_var_exist(k) = .true.
           else if( trim(name)=='tsea'.or.trim(name)=='TSEA' ) then
              k=k_tsea
+             sfc_var_exist(k) = .true.
           else if( trim(name)=='sheleg'.or.trim(name)=='SHELEG' ) then
              k=k_snwdph
-          else if( trim(name)=='stc'.or.trim(name)=='STC' ) then
+             sfc_var_exist(k) = .true.
+          else if( trim(name)=='stc'.or.trim(name)=='tslb' ) then
              k=k_stc 
-          else if( trim(name)=='smc'.or.trim(name)=='SMC' ) then
+             sfc_var_exist(k) = .true.
+          else if( trim(name)=='smc'.or.trim(name)=='smois' ) then
              k=k_smc
+             sfc_var_exist(k) = .true.
           else if( trim(name)=='SLMSK'.or.trim(name)=='slmsk' ) then
              k=k_slmsk
+             sfc_var_exist(k) = .true.
           else if( trim(name)=='T2M'.or.trim(name)=='t2m' ) then
              k=k_t2m
+             sfc_var_exist(k) = .true.
           else if( trim(name)=='Q2M'.or.trim(name)=='q2m' ) then
              k=k_q2m
+             sfc_var_exist(k) = .true.
           else if( trim(name)=='HOWV'.or.trim(name)=='howv' ) then
              k=k_howv
+             sfc_var_exist(k) = .true.
           else
              cycle 
           endif
@@ -2263,29 +2288,32 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
 !-- broadcast the updated i_howv_3dda to all tasks (!!!!)
     call mpi_bcast(i_howv_3dda, 1, mpi_itype, mype_2d, mpi_comm_world, iret_bcast)
 
+!-- broadcast the updated sfc_var_exist to all tasks (!!!!)
+    call mpi_bcast(sfc_var_exist, n2d, mpi_itype, mype_2d, mpi_comm_world, iret_bcast)
+
 !!!!!!! scatter !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     call mpi_scatterv(work,ijns2d,displss2d,mpi_rtype,&
       sfcn2d,ijns2d(mm1),mpi_rtype,mype_2d,mpi_comm_world,ierror)
 
     deallocate ( work )
 
-    fact10(:,:,it)=sfcn2d(:,:,k_f10m)
-    soil_type(:,:,it)=sfcn2d(:,:,k_stype)
-    veg_frac(:,:,it)=sfcn2d(:,:,k_vfrac)
-    veg_type(:,:,it)=sfcn2d(:,:,k_vtype)
-    sfc_rough(:,:,it)=sfcn2d(:,:,k_zorl)
-    sfct(:,:,it)=sfcn2d(:,:,k_tsea)
-    sno(:,:,it)=sfcn2d(:,:,k_snwdph)
-    soil_temp(:,:,it)=sfcn2d(:,:,k_stc)
-    soil_moi(:,:,it)=sfcn2d(:,:,k_smc)
-    ges_z(:,:)=sfcn2d(:,:,k_orog)/grav
-    isli(:,:,it)=nint(sfcn2d(:,:,k_slmsk))
+    if ( sfc_var_exist(k_f10m)  ) fact10(:,:,it)=sfcn2d(:,:,k_f10m)
+    if ( sfc_var_exist(k_stype) ) soil_type(:,:,it)=sfcn2d(:,:,k_stype)
+    if ( sfc_var_exist(k_vfrac) ) veg_frac(:,:,it)=sfcn2d(:,:,k_vfrac)
+    if ( sfc_var_exist(k_vtype) ) veg_type(:,:,it)=sfcn2d(:,:,k_vtype)
+    if ( sfc_var_exist(k_zorl)  ) sfc_rough(:,:,it)=sfcn2d(:,:,k_zorl)
+    if ( sfc_var_exist(k_tsea)  ) sfct(:,:,it)=sfcn2d(:,:,k_tsea)
+    if ( sfc_var_exist(k_snwdph)) sno(:,:,it)=sfcn2d(:,:,k_snwdph)
+    if ( sfc_var_exist(k_stc)   ) soil_temp(:,:,it)=sfcn2d(:,:,k_stc)
+    if ( sfc_var_exist(k_smc)   ) soil_moi(:,:,it)=sfcn2d(:,:,k_smc)
+    if ( sfc_var_exist(k_orog)  ) ges_z(:,:)=sfcn2d(:,:,k_orog)/grav
+    if ( sfc_var_exist(k_slmsk) ) isli(:,:,it)=nint(sfcn2d(:,:,k_slmsk))
     if(i_use_2mq4b > 0 .and. i_use_2mt4b > 0 ) then
-       ges_t2m(:,:)=sfcn2d(:,:,k_t2m)
-       ges_q2m(:,:)=sfcn2d(:,:,k_q2m)
+       if ( sfc_var_exist(k_t2m) ) ges_t2m(:,:)=sfcn2d(:,:,k_t2m)
+       if ( sfc_var_exist(k_q2m)   ) ges_q2m(:,:)=sfcn2d(:,:,k_q2m)
     endif
     if ( i_howv_3dda == 1 ) then
-       ges_howv(:,:)=sfcn2d(:,:,k_howv)
+       if ( sfc_var_exist(k_howv)  ) ges_howv(:,:)=sfcn2d(:,:,k_howv)
     endif
     deallocate (sfcn2d,a)
     return
@@ -4687,6 +4715,9 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file
                 call fv3_ll_to_h(work_a(:,:),work_b(:,:),nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag)
              endif
           endif
+          if (ifindstrloc(vnames_nonnegativetracers,trim(varname))> 0) then
+              where (work_b <0.0_r_kind)  work_b=0.0_r_kind
+          endif
           if(fv3_io_layout_y > 1) then
              do nio=0,fv3_io_layout_y-1
                 countloc=(/nxcase,ny_layout_len(nio),1/)
@@ -4847,6 +4878,9 @@ subroutine gsi_fv3ncdf_write_v1(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3f
       else
         call fv3_ll_to_h(work_a(:,:),work_b(:,:),nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag)
       endif
+      if (ifindstrloc(vnames_nonnegativetracers,trim(varname))> 0) then
+           where (work_b <0.0_r_kind)  work_b=0.0_r_kind
+      endif
       call check( nf90_put_var(gfile_loc,VarId,work_b,start=startloc,count=countloc) )
     enddo  !ilevtot
     call check(nf90_close(gfile_loc))