Updates for soil moisture and soil temperature analysis (#675)

Co-authored-by: jswhit2 <Jeffrey.S.Whitaker@noaa.gov>
Co-authored-by: guoqing.ge <Guoqing.Ge@noaa.gov>

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/read_prepbufr.f90

@@ -1669,13 +1669,15 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
               pmq(k)=nint(qcmark(8,k))
            end do
 
-!          181, 183, 187, and 188 are the screen-level obs over land
-           global_2m_land = ( (kx==181 .or. kx==183 .or. kx==188 .or. kx==188 ) .and. hofx_2m_sfcfile )
+!          187, 181, and 183 are the screen-level obs over land
+!          note: don't need the hofx_2m_sfcfile if set usage in convinfo, and qm updated in the input file
+           global_2m_land = ( (kx==187 .or. kx==181 .or. kx==183) .and. hofx_2m_sfcfile )
 
 !          If temperature ob, extract information regarding virtual
 !          versus sensible temperature
            if(tob) then
               ! use tvirtual if tsensible flag not set, and not in either 2Dregional or global_2m DA mode
+              ! for now, keeping 2m obs as sensible, for global system.
               if ( (.not. tsensible)  .and. .not. (twodvar_regional .or. global_2m_land) ) then
 
                  do k=1,levs
@@ -1979,18 +1981,17 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
                      pqm(k)=2 ! otherwise, type 183 will be discarded.
                      qm=2
                      tqm(k)=2
-                     if (kx==187) obserr(3,k)=2.2
-                     if (kx==181) obserr(3,k)=1.5
-                     if (kx==183) obserr(3,k)=2.6
+                     if (kx==187) obserr(3,k)=2.0_r_double
+                     if (kx==181) obserr(3,k)=2.0_r_double
+                     if (kx==183) obserr(3,k)=2.0_r_double
                 endif
                 if (qob .and. qm == 9 ) then 
                      qm = 2
                      ! qob err specified as fraction of qsat, multiplied by 10.
-                     if (kx==187) obserr(2,k)=1.0
-                     if (kx==181) obserr(2,k)=1.0
-                     if (kx==183) obserr(2,k)=1.0
+                     if (kx==187) obserr(2,k)=1.0_r_double
+                     if (kx==181) obserr(2,k)=1.0_r_double
+                     if (kx==183) obserr(2,k)=1.0_r_double
                 endif
-
               endif
 !             Set usage variable              
               usage = zero

src/gsi/setupq.f90

@@ -376,7 +376,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
   do k=1,nobs
      ikx=nint(data(ikxx,k))
      itype=ictype(ikx)
-     landsfctype =( itype==181 .or. itype==183 .or. itype==187 .or. itype==188  )
+     landsfctype =( itype==181 .or. itype==183 .or. itype==187 )
      do l=k+1,nobs
         if (twodvar_regional .or. (hofx_2m_sfcfile .and. landsfctype) ) then
            duplogic=data(ilat,k) == data(ilat,l) .and.  &

src/gsi/setupt.f90

@@ -448,7 +448,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
   do k=1,nobs
      ikx=nint(data(ikxx,k))
      itype=ictype(ikx)
-     landsfctype =( itype==181 .or. itype==183 .or. itype==187 .or. itype==188 )
+     landsfctype =( itype==181 .or. itype==183 .or. itype==187 )
      do l=k+1,nobs
         if (twodvar_regional .or. (hofx_2m_sfcfile .and. landsfctype) ) then
            duplogic=data(ilat,k) == data(ilat,l) .and.  &
@@ -483,9 +483,6 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
 
 ! Run a buddy-check
 ! Note: buddy check crashes for hofx_2m_sfcfile option.
-! Ccurrent params have buddy radius of 108 km, max diff of 8 K.
-! The gross error check removes O-F > 7., so this is probably removing
-! most obs that fail the buddy check already
   if (twodvar_regional .and. buddycheck_t) call buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse)
 
 ! If requested, save select data for output to diagnostic file