GitHub Issue #540 Modify to assimilate radar reflectivity and convent…

…ional data simultaneously without side effects in EnVar (#543)

This PR removes forcing pure EnVar and adding an option (if_use_w_vr) to
assimilate radar reflectivity and conventional data simultaneously
without side effects in EnVar
(#540). Regression tests for
global 3dvar/4denvar/4dvar are not completed yet, but for other tests,
issues are not found except for "failed the scalability test" and
"exceeded maximum allowable hardware memory limit" on Orion.

Fixes #540

Co-authored-by: Sho Yokota <syokota@Orion-login-1.HPC.MsState.Edu>

src/gsi/gsimod.F90

@@ -24,7 +24,7 @@ module gsimod
 
   use obsmod, only: doradaroneob,oneoblat,oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,&
      radar_no_thinning,ens_hx_dbz_cut,static_gsi_nopcp_dbz,rmesh_dbz,&
-     rmesh_vr,zmesh_dbz,zmesh_vr,if_vterminal, if_model_dbz,if_vrobs_raw,&
+     rmesh_vr,zmesh_dbz,zmesh_vr,if_vterminal, if_model_dbz,if_vrobs_raw,if_use_w_vr,&
      minobrangedbz,maxobrangedbz,maxobrangevr,maxtiltvr,missing_to_nopcp,&
      ntilt_radarfiles,whichradar,&
      minobrangevr,maxtiltdbz,mintiltvr,mintiltdbz,l2rwthin,hurricane_radar 
@@ -765,7 +765,7 @@ module gsimod
        oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,&
        rmesh_vr,zmesh_dbz,zmesh_vr, ntilt_radarfiles, whichradar,&
        radar_no_thinning,ens_hx_dbz_cut,static_gsi_nopcp_dbz,rmesh_dbz,&
-       minobrangevr, maxtiltdbz, mintiltvr,mintiltdbz,if_vterminal,if_vrobs_raw,&
+       minobrangevr, maxtiltdbz, mintiltvr,mintiltdbz,if_vterminal,if_vrobs_raw,if_use_w_vr,&
        if_model_dbz,imp_physics,lupp,netcdf_diag,binary_diag,l_wcp_cwm,aircraft_recon,diag_version,&
        write_fv3_incr,incvars_to_zero,incvars_zero_strat,incvars_efold,diag_version,&
        cao_check,lcalc_gfdl_cfrac,tau_fcst,efsoi_order,lupdqc,lqcoef,cnvw_option,l2rwthin,hurricane_radar,&
@@ -2045,15 +2045,17 @@ subroutine gsimain_initialize
      baldiag_inc =.false.
   end if
 
-! If reflectivity is intended to be assimilated, beta_s0 should be zero.
+! Warning of reflectivity assimilation with static B
   if ( beta_s0 > 0.0_r_kind )then
     ! skipped in case of direct reflectivity DA because it works in Envar and hybrid
     if ( l_use_rw_columntilt .or. l_use_dbz_directDA) then
        do i=1,ndat
           if ( if_model_dbz .and. (index(dtype(i), 'dbz') /= 0) )then
-             write(6,*)'beta_s0 needs to be set to zero in this GSI version, when reflectivity is directly assimilated. &
-                        Static B extended for radar reflectivity assimilation will be included in future version.'
-             call stop2(8888)
+             if (mype==0) then
+                write(6,*)'GSIMOD:  ***WARNING*** static B for reflectivity is regarded as zero in this GSI version &
+                           even though beta_s0 =',beta_s0
+                write(6,*)'Static B extended for radar reflectivity assimilation will be included in future version.'
+             end if
           end if
        end do
     end if

src/gsi/intrw.f90

@@ -96,7 +96,7 @@ subroutine intrw_(rwhead,rval,sval)
 !$$$
   use kinds, only: r_kind,i_kind
   use constants, only: half,one,tiny_r_kind,cg_term,r3600
-  use obsmod, only: lsaveobsens,l_do_adjoint,luse_obsdiag
+  use obsmod, only: lsaveobsens,l_do_adjoint,luse_obsdiag,if_use_w_vr
   use qcmod, only: nlnqc_iter,varqc_iter
   use jfunc, only: jiter
   use gsi_bundlemod, only: gsi_bundle
@@ -128,15 +128,15 @@ subroutine intrw_(rwhead,rval,sval)
   call gsi_bundlegetpointer(sval,'u',su,istatus);ier=istatus+ier
   call gsi_bundlegetpointer(sval,'v',sv,istatus);ier=istatus+ier
   call gsi_bundlegetpointer(sval,'w',sw,istatus)
-  if (istatus==0) then
+  if (if_use_w_vr.and.istatus==0) then
      include_w=.true.
   else
      include_w=.false.
   end if
   call gsi_bundlegetpointer(rval,'u',ru,istatus);ier=istatus+ier
   call gsi_bundlegetpointer(rval,'v',rv,istatus);ier=istatus+ier
   call gsi_bundlegetpointer(rval,'w',rw,istatus)
-  if (istatus==0) then
+  if (if_use_w_vr.and.istatus==0) then
      include_w=.true.
   else
      include_w=.false.

src/gsi/obsmod.F90

@@ -470,7 +470,7 @@ module obsmod
   ! ==== DBZ DA ===
   public :: ntilt_radarfiles
   public :: whichradar
-  public :: vr_dealisingopt, if_vterminal, if_model_dbz, inflate_obserr, if_vrobs_raw, l2rwthin 
+  public :: vr_dealisingopt, if_vterminal, if_model_dbz, inflate_obserr, if_vrobs_raw, if_use_w_vr, l2rwthin
 
   public :: doradaroneob,oneoblat,oneoblon
   public :: oneobddiff,oneobvalue,oneobheight,oneobradid
@@ -617,7 +617,7 @@ module obsmod
 
   logical ::  ta2tb
   logical ::  doradaroneob
-  logical :: vr_dealisingopt, if_vterminal, if_model_dbz, inflate_obserr, if_vrobs_raw, l2rwthin
+  logical :: vr_dealisingopt, if_vterminal, if_model_dbz, inflate_obserr, if_vrobs_raw, if_use_w_vr, l2rwthin
   character(4) :: whichradar,oneobradid
   real(r_kind) :: oneoblat,oneoblon,oneobddiff,oneobvalue,oneobheight
   logical :: radar_no_thinning
@@ -747,6 +747,7 @@ subroutine init_obsmod_dflts
     if_vterminal=.false.
     l2rwthin    =.false.  
     if_vrobs_raw=.false.
+    if_use_w_vr=.true.
     if_model_dbz=.false.
     inflate_obserr=.false.
     whichradar="KKKK"

src/gsi/setuprw.f90

@@ -114,7 +114,7 @@ subroutine setuprw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa
   use obsmod, only: rmiss_single,lobsdiag_forenkf,&
                     lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset,&
                     if_vterminal, ens_hx_dbz_cut, if_model_dbz, &
-                    doradaroneob,oneobddiff,oneobvalue, if_vrobs_raw
+                    doradaroneob,oneobddiff,oneobvalue, if_vrobs_raw, if_use_w_vr
   use obsmod, only: netcdf_diag, binary_diag, dirname,ianldate
   use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, &
        nc_diag_write, nc_diag_data2d
@@ -972,7 +972,7 @@ subroutine check_vars_ (proceed, include_w)
   call gsi_metguess_get ('var::v' , ivar, istatus )
   proceed=proceed.and.ivar>0
   call gsi_metguess_get ('var::w' , ivar, istatus )
-  if (ivar>0) then
+  if (if_use_w_vr.and.ivar>0) then
      include_w=.true.
      if(if_vterminal)then
       if( .not. if_model_dbz ) then

src/gsi/stprw.f90

@@ -83,6 +83,7 @@ subroutine stprw(rwhead,rval,sval,out,sges,nstep)
 !
 !$$$
   use kinds, only: r_kind,i_kind,r_quad
+  use obsmod, only: if_use_w_vr
   use qcmod, only: nlnqc_iter,varqc_iter
   use constants, only: half,one,two,tiny_r_kind,cg_term,zero_quad,r3600
   use gsi_bundlemod, only: gsi_bundle
@@ -124,15 +125,15 @@ subroutine stprw(rwhead,rval,sval,out,sges,nstep)
   call gsi_bundlegetpointer(sval,'u',su,istatus);ier=istatus+ier
   call gsi_bundlegetpointer(sval,'v',sv,istatus);ier=istatus+ier
   call gsi_bundlegetpointer(sval,'w',sw,istatus)
-  if (istatus==0) then
+  if (if_use_w_vr.and.istatus==0) then
      include_w=.true.
   else
      include_w=.false.
   end if
   call gsi_bundlegetpointer(rval,'u',ru,istatus);ier=istatus+ier
   call gsi_bundlegetpointer(rval,'v',rv,istatus);ier=istatus+ier
   call gsi_bundlegetpointer(rval,'w',rw,istatus)
-  if (istatus==0) then
+  if (if_use_w_vr.and.istatus==0) then
      include_w=.true.
   else
      include_w=.false.