Merge pull request #167 from GEOS-ESM/bugfix/rtodling/divo_tropoprs

Fix to Div-Vor calculation - and consequences

CHANGELOG.md

@@ -11,7 +11,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 
+<<<<<<< HEAD
+- bugfix: update calculation of div/vor
+- merge of 5.30.3 with 5.40.0 and 5.29.5-p7
+=======
 - add offline RO bufr handling program
+>>>>>>> refs/remotes/origin/bugfix/rtodling/divo_tropoprs
 
 ### Fixed
 

GEOSaana_GridComp/GSI_GridComp/GSI_GridCompMod.F90

@@ -1,7 +1,5 @@
 #include "MAPL_ErrLog.h" 
-#ifdef _REAL8_
-#define _GMAO_FVGSI_
-#endif
+!#define _GMAO_FVGSI_
 !#define PRINT_STATES
 !#define SFCverbose
 !#define UPAverbose
@@ -2194,6 +2192,7 @@ subroutine GSI_GridCompComputeVorDiv_(lit)
    use compact_diffs, only: uv2vordiv
    use xhat_vordivmod, only: xhat_vordiv_calc2
 #endif /* _GMAO_FVGSI_ */
+   use mpeu_util, only: die
 
    implicit none
 
@@ -2288,13 +2287,18 @@ subroutine GSI_GridCompComputeVorDiv_(lit)
    VERIFY_(STATUS)
 
 #else /* _GMAO_FVGSI_ */
+   if(IamRoot) print *,trim(Iam),': Using GSI-based div/vor procedure'
    if(.not.cdiff_created()) call create_cdiff_coefs()
    if(.not.cdiff_initialized()) call inisph(rearth,rlats(2),wgtlats(2),nlon,nlat-2)
    ier=0
+   call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'div', ges_div_nnn, istatus );ier=ier+istatus
+   call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'vor', ges_vor_nnn, istatus );ier=ier+istatus
    call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'u', ges_u_it, istatus );ier=ier+istatus
    call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'v', ges_v_it, istatus );ier=ier+istatus
    if(ier==0) then
       call xhat_vordiv_calc2 (ges_u_it,ges_v_it,ges_vor_nnn,ges_div_nnn)
+   else
+      call die(Iam,': unable to calculate vor/div',99)
    endif
 !!   call destroy_cdiff_coefs
 #endif /* _GMAO_FVGSI_ */

GEOSaana_GridComp/GSI_GridComp/analyzer

@@ -649,7 +649,8 @@ sub init {
    if ( $ENV{CRTM_COEFFS} ) {
       Assignfn( "$CRTM_COEFFS","CRTM_Coeffs");
    } else {
-      Assignfn( "$fvInput/$myetc/fix_ncep20221018/REL-2.4.0-jcsda/CRTM_Coeffs/Little_Endian","CRTM_Coeffs");
+#     Assignfn( "$fvInput/$myetc/fix_ncep20221018/REL-2.4.0-jcsda/CRTM_Coeffs/Little_Endian","CRTM_Coeffs");
+      Assignfn( "$fvInput/$myetc/r21c_ncep20221018/Little_Endian","CRTM_Coeffs");
    }
    Assignfn( "$fvhome/run/prepobs_errtable.global","errtable");
 

GEOSaana_GridComp/GSI_GridComp/bkgvar_rewgt.f90

@@ -376,7 +376,7 @@ subroutine bkgvar_rewgt(sfvar,vpvar,tvar,psvar,mype)
   nsmth=8
   call smooth2d(sfvar,vpvar,tvar,psvar,nsig,nsmth,mype)
 
-  if (bkgv_write) call write_bkgvars_grid(sfvar,vpvar,tvar,psvar,mype)
+  if (bkgv_write) call write_bkgvars_grid(sfvar,vpvar,tvar,psvar,'bkgvar',mype)
   if(mype==0) write(6,*) 'bkgvar_rewgt: Flow-dependent feature on: nt=',nfldsig, ' minus nt= 1'
 
   return

GEOSaana_GridComp/GSI_GridComp/gsimod.F90

@@ -19,6 +19,7 @@ module gsimod
      dtbduv_on,time_window_max,offtime_data,init_directories,oberror_tune,ext_sonde, &
      blacklst,init_obsmod_vars,lobsdiagsave,lobskeep,lobserver,hilbert_curve,&
      lread_obs_save,lread_obs_skip,time_window_rad,lgpsbnd_revint
+  use obsmod, only: saberTbot,saberTtop
   use gsi_dbzOper, only: diag_radardbz
 
   use obsmod, only: doradaroneob,oneoblat,oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,&
@@ -642,7 +643,7 @@ module gsimod
        minobrangevr, maxtiltdbz, mintiltvr,mintiltdbz,if_vterminal,if_vrobs_raw,&
        if_model_dbz,imp_physics,lupp,netcdf_diag,binary_diag,l_wcp_cwm,diag_version,&
        cao_check,lcalc_gfdl_cfrac,tau_fcst,evfsoi_order,lupdqc,lqcoef,lgpsbnd_revint,&
-       wrtgeovals
+       wrtgeovals,saberTbot,saberTtop
 
 ! GRIDOPTS (grid setup variables,including regional specific variables):
 !     jcap     - spectral resolution

GEOSaana_GridComp/GSI_GridComp/obsmod.F90

@@ -389,6 +389,8 @@ module obsmod
 !                        includes cwm for both swcp and lwcp or not
 !   def lgpsbnd_revint - namelist logical determining whether or not to use
 !                        revised integration of bending angle profiles
+!   def saberTbot      - pressure level below which to ignore SABER-T
+!   def saberTtop      - pressure level above which to ignore SABER-T
 !
 ! attributes:
 !   langauge: f90
@@ -445,6 +447,7 @@ module obsmod
   public :: use_limit,lrun_subdirs
   public :: l_foreaft_thin,luse_obsdiag
   public :: lgpsbnd_revint
+  public :: saberTbot, saberTtop
 
   ! ==== DBZ DA ===
   public :: ntilt_radarfiles
@@ -557,6 +560,7 @@ module obsmod
   real(r_kind) ::static_gsi_nopcp_dbz
   real(r_kind) ::rmesh_dbz,zmesh_dbz
   real(r_kind) ::rmesh_vr,zmesh_vr
+  real(r_kind) ::saberTbot,saberTtop
 
   logical :: debugmode
   real(r_kind) :: minobrangevr,maxobrangevr,mintiltvr,maxtiltvr
@@ -666,6 +670,10 @@ subroutine init_obsmod_dflts
     maxtiltvr=20.0_r_kind
     missing_to_nopcp=.false.
 
+!   Default for SABER use set to ignore all
+    saberTbot =    0.0_r_kind ! mb
+    saberTtop = 1000.0_r_kind ! mb
+
 !   Set logical flag
     perturb_obs = .false.   ! .true. = perturb observations
     oberror_tune = .false.   ! .true. = tune oberror

GEOSaana_GridComp/GSI_GridComp/prewgt.f90

@@ -438,7 +438,7 @@ subroutine prewgt(mype)
   if (bkgv_flowdep) then
       call bkgvar_rewgt(sfvar,vpvar,tvar,psvar,mype)
   else
-      if (bkgv_write) call write_bkgvars_grid(sfvar,vpvar,tvar,psvar,mype)
+      if (bkgv_write) call write_bkgvars_grid(sfvar,vpvar,tvar,psvar,'prewgt',mype)
   endif
 
 ! vertical length scales

GEOSaana_GridComp/GSI_GridComp/setupt.f90

@@ -33,6 +33,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
        lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset
   use obsmod, only: dplat
   use obsmod, only: wrtgeovals
+  use obsmod, only: saberTbot,saberTtop
   use m_obsNode, only: obsNode
   use m_tNode, only: tNode
   use m_tNode, only: tNode_appendto
@@ -853,6 +854,14 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
         endif
      endif
 
+! Select SABER-T observations within given slab
+     if(itype ==294 ) then
+        if(prest>saberTtop .or. prest<saberTbot) then
+           error=zero
+           ratio_errors=zero
+        endif
+     endif
+
 ! Compute innovation
      if(i_use_2mt4b>0 .and. sfctype) then
         ddiff = tob-tges2m

GEOSaana_GridComp/GSI_GridComp/tpause.f90

@@ -61,6 +61,8 @@ subroutine tpause(mype,method)
   real(r_kind),parameter:: r40=40.0_r_kind
   real(r_kind),parameter:: r1e5=1.0e5_r_kind
 
+  logical, parameter :: tpause_debug = .true.
+
 ! Declare local variables
   logical t_method,pvoz_capable
 
@@ -82,6 +84,7 @@ subroutine tpause(mype,method)
   real(r_kind),dimension(:,:,:),pointer :: ges_tv_nt=>NULL()
   real(r_kind),dimension(:,:,:),pointer :: ges_oz_nt=>NULL()
   real(r_kind),dimension(:,:,:),pointer :: ges_vor_nt=>NULL()
+  real(r_kind),dimension(:,:,:),pointer :: ges_div_nt=>NULL()
 
 !================================================================================
 ! Set local constants
@@ -207,6 +210,11 @@ subroutine tpause(mype,method)
         end do
      end do
 
+     if (tpause_debug) then
+        call gsi_bundlegetpointer (gsi_metguess_bundle(nt),'div',ges_div_nt,istatus)
+        call write_bkgvars_grid(ges_div_nt,ges_vor_nt,ges_oz_nt,tropprs,'tpause',mype)
+     endif
+
 ! End of tropopause location method blocks     
   endif
 

GEOSaana_GridComp/GSI_GridComp/write_bkgvars_grid.f90

@@ -1,4 +1,4 @@
-subroutine write_bkgvars_grid(a,b,c,d,mype)
+subroutine write_bkgvars_grid(a,b,c,d,fname,mype)
 !$$$  subroutine documentation block
 !
 ! subprogram:    write_bkgvars_grid
@@ -28,22 +28,26 @@ subroutine write_bkgvars_grid(a,b,c,d,mype)
 !$$$
   use kinds, only: r_kind,i_kind,r_single
   use gridmod, only: nlat,nlon,nsig,lat2,lon2
+  use guess_grids, only: get_ref_gesprs
+  use mpeu_util, only: get_lun => luavail
   implicit none
 
   integer(i_kind)                       ,intent(in   ) :: mype
+  character(len=*)                      ,intent(in   ) :: fname
 
   real(r_kind),dimension(lat2,lon2,nsig),intent(in   ) :: a,b,c
   real(r_kind),dimension(lat2,lon2)     ,intent(in   ) :: d
 
-  character(255):: grdfile
 
   real(r_kind),dimension(nlat,nlon,nsig):: ag,bg,cg
   real(r_kind),dimension(nlat,nlon):: dg
 
   real(r_single),dimension(nlon,nlat,nsig):: a4,b4,c4
   real(r_single),dimension(nlon,nlat):: d4
 
-  integer(i_kind) ncfggg,iret,i,j,k
+  character(len=80) :: grdfile
+  integer(i_kind) iret,i,j,k,lu
+  real(r_kind),dimension(nsig+1)::prs
 
 ! gather stuff to processor 0
   do k=1,nsig
@@ -72,15 +76,38 @@ subroutine write_bkgvars_grid(a,b,c,d,mype)
      end do
 
 ! Create byte-addressable binary file for grads
-     grdfile='bkgvar_rewgt.grd'
-     ncfggg=len_trim(grdfile)
-     call baopenwt(22,grdfile(1:ncfggg),iret)
+     grdfile = trim(fname)//'.grd' 
+     call baopenwt(22,trim(grdfile),iret)
      call wryte(22,4*nlat*nlon*nsig,a4)
      call wryte(22,4*nlat*nlon*nsig,b4)
      call wryte(22,4*nlat*nlon*nsig,c4)
      call wryte(22,4*nlat*nlon,d4)
      call baclose(22,iret)
-  end if
+
+     call get_ref_gesprs(prs)
+
+!    Now create corresponding grads table file
+     lu=get_lun()
+     open(lu,file=trim(fname)//'.ctl',form='formatted')
+     write(lu,'(2a)') 'DSET  ^', trim(grdfile)
+     write(lu,'(2a)') 'TITLE ', 'gsi berror variances'
+     write(lu,'(a,2x,e13.6)') 'UNDEF', 1.E+15 ! any other preference for this?
+     write(lu,'(a,2x,i4,2x,a,2x,f5.1,2x,f9.6)') 'XDEF',nlon, 'LINEAR',   0.0, 360./nlon
+     write(lu,'(a,2x,i4,2x,a,2x,f5.1,2x,f9.6)') 'YDEF',nlat, 'LINEAR', -90.0, 180./(nlat-1.)
+     write(lu,'(a,2x,i4,2x,a,1x,f8.3)')    'ZDEF',nsig, 'LEVELS', prs(1) ! prs is in mb
+     do k=2,nsig  ! grads tends not to like a long line of pressures, thus split
+        write(lu,'(1x,f8.3)')   prs(k) ! prs is in mb
+     enddo
+     write(lu,'(a,2x,i4,2x,a)')   'TDEF', 1, 'LINEAR 12:00Z04JUL1776 6hr' ! any date suffices
+     write(lu,'(a,2x,i4)')        'VARS', 4
+     write(lu,'(a,1x,2(i4,1x),a)') 'a',nsig,0,'a'
+     write(lu,'(a,1x,2(i4,1x),a)') 'b',nsig,0,'b'
+     write(lu,'(a,1x,2(i4,1x),a)') 'c',nsig,0,'c'
+     write(lu,'(a,1x,2(i4,1x),a)') 'd',   1,0,'d'
+     write(lu,'(a)') 'ENDVARS'
+     close(lu)
+
+  end if ! mype=0
 
   return
 end subroutine write_bkgvars_grid