Merge branch 'NOAA-EMC:develop' into feature/ThompsonMP_Azadeh

modulefiles/gsi_jet.intel.lua

@@ -1,7 +1,7 @@
 help([[
 ]])
 
-prepend_path("MODULEPATH", "/mnt/lfs4/HFIP/hfv3gfs/role.epic/spack-stack/spack-stack-1.6.0/envs/gsi-addon-dev/install/modulefiles/Core")
+prepend_path("MODULEPATH", "/mnt/lfs4/HFIP/hfv3gfs/role.epic/spack-stack/spack-stack-1.6.0/envs/gsi-addon-dev-rocky8/install/modulefiles/Core")
 
 local python_ver=os.getenv("python_ver") or "3.11.6"
 local stack_intel_ver=os.getenv("stack_intel_ver") or "2021.5.0"

modulefiles/gsi_noaacloud.intel.lua

@@ -0,0 +1,25 @@
+help([[
+]])
+
+prepend_path("MODULEPATH", "/contrib/spack-stack/spack-stack-1.6.0/envs/gsi-addon-env/install/modulefiles/Core")
+
+local python_ver=os.getenv("python_ver") or "3.10.13"
+local stack_intel_ver=os.getenv("stack_intel_ver") or "2021.3.0"
+local stack_impi_ver=os.getenv("stack_impi_ver") or "2021.3.0"
+local cmake_ver=os.getenv("cmake_ver") or "3.23.1"
+local prod_util_ver=os.getenv("prod_util_ver") or "2.1.1"
+
+load(pathJoin("stack-intel", stack_intel_ver))
+load(pathJoin("stack-intel-oneapi-mpi", stack_impi_ver))
+load(pathJoin("python", python_ver))
+load(pathJoin("cmake", cmake_ver))
+
+load("gsi_common")
+load(pathJoin("prod_util", prod_util_ver))
+
+pushenv("CFLAGS", "-xHOST")
+pushenv("FFLAGS", "-xHOST")
+
+pushenv("GSI_BINARY_SOURCE_DIR", "/contrib/Wei.Huang/data/hack-orion/fix/gsi/20240208")
+
+whatis("Description: GSI environment on NOAA Cloud with Intel Compilers")

regression/global_4denvar.sh

@@ -128,6 +128,8 @@ errtable=$fixgsi/prepobs_errtable.global
 aeroinfo=$fixgsi/global_aeroinfo.txt
 atmsbeaminfo=$fixgsi/atms_beamwidth.txt
 cloudyinfo=$fixgsi/cloudy_radiance_info.txt
+cris_clddet=$fixgsi/CRIS_CLDDET.NL
+iasi_clddet=$fixgsi/IASI_CLDDET.NL
 
 emiscoef_IRwater=$fixcrtm/Nalli.IRwater.EmisCoeff.bin
 emiscoef_IRice=$fixcrtm/NPOESS.IRice.EmisCoeff.bin
@@ -169,6 +171,8 @@ $ncp $errtable ./errtable
 $ncp $aeroinfo ./aeroinfo
 $ncp $atmsbeaminfo ./atms_beamwidth.txt
 $ncp $cloudyinfo   ./cloudy_radiance_info.txt
+$ncp $cris_clddet ./CRIS_CLDDET.NL
+$ncp $iasi_clddet ./IASI_CLDDET.NL
 
 $ncp $bufrtable ./prepobs_prep.bufrtable
 $ncp $bftab_sst ./bftab_sstphr

regression/regression_namelists.sh

@@ -68,7 +68,7 @@ export gsi_namelist="
    dfact=0.75,dfact1=3.0,noiqc=.true.,oberrflg=.false.,c_varqc=0.04,
    use_poq7=.true.,qc_noirjaco3_pole=.true.,vqc=.false.,nvqc=.true.,hub_norm=.true.,
    aircraft_t_bc=.true.,biaspredt=1.0e5,upd_aircraft=.true.,cleanup_tail=.true.,
-   tcp_width=70.0,tcp_ermax=7.35,
+   tcp_width=70.0,tcp_ermax=7.35,cris_cads=.true.,iasi_cads=.true.,
    $OBSQC
  /
  &OBS_INPUT

src/gsi/correlated_obsmod.F90

@@ -961,14 +961,18 @@ subroutine upd_varch_
             enddo
             nchanl1=jc
 
-            if(nchanl1==0) call die(myname_,' improperly set GSI_BundleErrorCov')
             if(.not.amiset_(GSI_BundleErrorCov(itbl))) then 
-               if (iamroot_) write(6,*) 'WARNING: Error Covariance not set for ',trim(idnames(itbl))
+               if (iamroot_) write(6,*) trim(myname_), ' WARNING: Error Covariance not set for ',trim(idnames(itbl))
                cycle read_tab
             endif
 
             nch_active=GSI_BundleErrorCov(itbl)%nch_active
-            if(nch_active<0) return
+            if(nch_active<0) then
+               if (iamroot_) write(6,*) trim(myname_), ' WARNING: No active channels for ',trim(idnames(itbl))
+               return
+            endif
+
+            if(nchanl1==0) call die(myname_,' improperly set GSI_BundleErrorCov')            
 
             if(GMAO_ObsErrorCov)then
                do jj=1,nch_active

src/gsi/gsi_rfv3io_mod.f90

@@ -61,7 +61,7 @@ module gsi_rfv3io_mod
   use rapidrefresh_cldsurf_mod, only: i_use_2mq4b,i_use_2mt4b
   use chemmod, only: naero_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3,laeroana_fv3cmaq
   use chemmod, only: naero_smoke_fv3,aeronames_smoke_fv3,laeroana_fv3smoke  
-  use rapidrefresh_cldsurf_mod, only: i_howv_3dda
+  use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda
 
   implicit none
   public type_fv3regfilenameg
@@ -147,7 +147,7 @@ module gsi_rfv3io_mod
   public :: mype_u,mype_v,mype_t,mype_q,mype_p,mype_oz,mype_ql
   public :: mype_qi,mype_qr,mype_qs,mype_qg,mype_qnr,mype_w
   public :: k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc
-  public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv
+  public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv,k_gust
   public :: ijns,ijns2d,displss,displss2d,ijnz,displsz_g
   public :: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv
   public :: fv3lam_io_tracerchemvars3d_nouv,fv3lam_io_tracersmokevars3d_nouv
@@ -158,7 +158,7 @@ module gsi_rfv3io_mod
   integer(i_kind) mype_qi,mype_qr,mype_qs,mype_qg,mype_qnr,mype_w
 
   integer(i_kind) k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc
-  integer(i_kind) k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv
+  integer(i_kind) k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv,k_gust
   parameter(                   &  
     k_f10m =1,                  &   !fact10
     k_stype=2,                  &   !soil_type
@@ -174,7 +174,8 @@ module gsi_rfv3io_mod
     k_q2m  =12,                 & ! 2 m Q
     k_orog =13,                 & !terrain
     k_howv =14,                 &   ! significant wave height (aka howv in GSI)
-    n2d=14                   )
+    k_gust =15,                 &   ! wind gust (aka gust in GSI)
+    n2d=15                   )
   logical :: grid_reverse_flag
   character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_dynmetvars3d_nouv 
                                     ! copy of cvars3d excluding uv 3-d fields   
@@ -996,6 +997,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
     real(r_kind),dimension(:,:),pointer::ges_t2m=>NULL()
     real(r_kind),dimension(:,:),pointer::ges_q2m=>NULL()
     real(r_kind),dimension(:,:),pointer::ges_howv=>NULL()
+    real(r_kind),dimension(:,:),pointer::ges_gust=>NULL()
 
     real(r_kind),dimension(:,:,:),pointer::ges_ql=>NULL()
     real(r_kind),dimension(:,:,:),pointer::ges_qi=>NULL()
@@ -1274,6 +1276,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
           else if(trim(vartem)=='t2m') then
           else if(trim(vartem)=='q2m') then
           else if(trim(vartem)=='howv') then
+          else if(trim(vartem)=='gust') then
           else 
             write(6,*)'the metvarname2 ',trim(vartem),' has not been considered yet, stop'
             call stop2(333)
@@ -1294,7 +1297,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
         vartem=trim(name_metvars2d(i))
         if(.not.( (trim(vartem)=='ps'.and.fv3sar_bg_opt==0).or.(trim(vartem)=="z") &
                   .or.(trim(vartem)=="t2m").or.(trim(vartem)=="q2m")               &
-                  .or.(trim(vartem)=="howv")))  then ! z is treated separately
+                  .or.(trim(vartem)=="howv").or.(trim(vartem)=="gust")))  then ! z is treated separately
           if (ifindstrloc(vardynvars,trim(vartem)) > 0) then
             jdynvar=jdynvar+1
             fv3lam_io_dynmetvars2d_nouv(jdynvar)=trim(vartem)
@@ -1557,6 +1560,12 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
             if (ier/=0) call die(trim(myname),'cannot get pointers for howv, ier=',ier)
          endif
 
+!---     wind gust (gust)
+         if ( i_gust_3dda == 1 ) then
+            call GSI_BundleGetPointer(GSI_MetGuess_Bundle(it),'gust',ges_gust,istatus ); ier=ier+istatus
+            if (ier/=0) call die(trim(myname),'cannot get pointers for gust, ier=',ier)
+         endif
+
          if(mype == 0 ) then
            call check(nf90_open(fv3filenamegin(it)%dynvars,nf90_nowrite,loc_id))
            call check(nf90_inquire(loc_id,formatNum=ncfmt))
@@ -1739,7 +1748,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
          endif
 
 
-         call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m,ges_howv)
+         call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m,   &
+                                 ges_howv,ges_gust)
 
          if(i_use_2mq4b > 0 .and. i_use_2mt4b > 0 ) then
 ! Convert 2m guess mixing ratio to specific humidity
@@ -1975,7 +1985,8 @@ end subroutine gsi_bundlegetpointer_fv3lam_tracerchem_nouv
 
 end subroutine read_fv3_netcdf_guess
 
-subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
+subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,          &
+                              ges_howv,ges_gust)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    gsi_fv3ncdf2d_read       
@@ -2022,6 +2033,7 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
     real(r_kind),               intent(in),dimension(:,:),pointer::ges_t2m
     real(r_kind),               intent(in),dimension(:,:),pointer::ges_q2m
     real(r_kind),               intent(in),dimension(:,:),pointer::ges_howv
+    real(r_kind),               intent(in),dimension(:,:),pointer::ges_gust
     type (type_fv3regfilenameg),intent(in) :: fv3filenamegin
 
     character(len=max_varname_length) :: name
@@ -2036,8 +2048,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
     integer(i_kind) kk,n,ns,j,ii,jj,mm1
       character(len=:),allocatable :: sfcdata   !='fv3_sfcdata'
       character(len=:),allocatable :: dynvars   !='fv3_dynvars'
-! for checking the existence of howv in firstguess file
-    integer(i_kind) id_howv
+! for checking the existence of howv/gust in firstguess file
+    integer(i_kind) id_howv, id_gust
     integer(i_kind) iret_bcast
 
 ! for io_layout > 1
@@ -2057,8 +2069,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
     allocate(work(itotsub*n2d))
     allocate( sfcn2d(lat2,lon2,n2d))
 
-!-- initialisation of the array for howv
+!-- initialisation of the array for howv/gust
     sfcn2d(:,:,k_howv) = zero
+    sfcn2d(:,:,k_gust) = zero
 
 !-- initialisation of the array for sfc_var_exist 
     sfc_var_exist = .false.
@@ -2107,6 +2120,21 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
              trim(sfcdata), ', iret, varid = ',iret, id_howv,' (on pe: ', mype,').'
          end if
        end if
+!---   check the existence of wind gust (gust) in 2D FV3-LAM firstguess file
+!      (similar as done above for howv)
+       if ( i_gust_3dda == 1 ) then
+         iret = nf90_inq_varid(gfile_loc,'gust',id_gust)
+         if ( iret /= nf90_noerr ) then
+           iret = nf90_inq_varid(gfile_loc,'GUST',id_gust) ! double check with name in uppercase
+         end if
+         if ( iret /= nf90_noerr ) then
+           i_gust_3dda = 0                ! gust does not exist in firstguess, then stop GSI run.
+           call die('gsi_fv3ncdf2d_read','Warning: CANNOT find gust in firstguess, aborting..., iret = ', iret)
+         else
+           write(6,'(1x,A,1x,A,1x,A,1x,I4,1x,I4,1x,A,1x,I4.4,A)') 'gsi_fv3ncdf2d_read:: Found gust in firstguess ',  &
+             trim(sfcdata), ', iret, varid = ',iret, id_gust,' (on pe: ', mype,').'
+         end if
+       end if
 
    !!!!!!!!!!!! read in 2d variables !!!!!!!!!!!!!!!!!!!!!!!!!!
        do i=ndimensions+1,nvariables
@@ -2150,6 +2178,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
           else if( trim(name)=='HOWV'.or.trim(name)=='howv' ) then
              k=k_howv
              sfc_var_exist(k) = .true.
+          else if( trim(name)=='GUST'.or.trim(name)=='gust' ) then
+             k=k_gust
+             sfc_var_exist(k) = .true.
           else
              cycle 
           endif
@@ -2283,8 +2314,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
        if(allocated(sfc_fulldomain)) deallocate (sfc_fulldomain)
     endif  ! mype
 
-!-- broadcast the updated i_howv_3dda to all tasks (!!!!)
+!-- broadcast the updated i_howv_3dda, i_gust_3dda to all tasks (!!!!)
     call mpi_bcast(i_howv_3dda, 1, mpi_itype, mype_2d, mpi_comm_world, iret_bcast)
+    call mpi_bcast(i_gust_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)
@@ -2313,6 +2345,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
     if ( i_howv_3dda == 1 ) then
        if ( sfc_var_exist(k_howv)  ) ges_howv(:,:)=sfcn2d(:,:,k_howv)
     endif
+    if ( i_gust_3dda == 1 ) then
+       if ( sfc_var_exist(k_gust)  ) ges_gust(:,:)=sfcn2d(:,:,k_gust)
+    endif
     deallocate (sfcn2d,a)
     return
 end subroutine gsi_fv3ncdf2d_read
@@ -3628,6 +3663,7 @@ subroutine wrfv3_netcdf(fv3filenamegin)
     real(r_kind),pointer,dimension(:,:  ):: ges_t2m =>NULL()
     real(r_kind),pointer,dimension(:,:  ):: ges_q2m  =>NULL()
     real(r_kind),pointer,dimension(:,:  ):: ges_howv =>NULL()
+    real(r_kind),pointer,dimension(:,:  ):: ges_gust =>NULL()
 
     integer(i_kind) i,k
 
@@ -3750,6 +3786,9 @@ subroutine wrfv3_netcdf(fv3filenamegin)
     if ( i_howv_3dda == 1 ) then
        call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'howv',ges_howv,istatus); ier=ier+istatus
     endif
+    if ( i_gust_3dda == 1 ) then
+       call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'gust',ges_gust,istatus); ier=ier+istatus
+    endif
     if (ier/=0) call die('wrfv3_netcdf','cannot get pointers for fv3 met-fields, ier =',ier)
 
     if (laeroana_fv3cmaq) then
@@ -3964,6 +4003,10 @@ subroutine wrfv3_netcdf(fv3filenamegin)
     if ( i_howv_3dda == 1 ) then
       call gsi_fv3ncdf_write_sfc(fv3filenamegin,'howv',ges_howv,add_saved)
     endif
+!-- output analysis of gust
+    if ( i_gust_3dda == 1 ) then
+      call gsi_fv3ncdf_write_sfc(fv3filenamegin,'gust',ges_gust,add_saved)
+    endif
 
     if(allocated(g_prsi)) deallocate(g_prsi)
 

src/gsi/gsimod.F90

@@ -184,7 +184,7 @@ module gsimod
                             cld_bld_coverage,cld_clr_coverage,&
                             i_cloud_q_innovation,i_ens_mean,DTsTmax,&
                             i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check, &
-                            corp_howv, hwllp_howv
+                            corp_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust
   use gsi_metguess_mod, only: gsi_metguess_init,gsi_metguess_final
   use gsi_chemguess_mod, only: gsi_chemguess_init,gsi_chemguess_final
   use tcv_mod, only: init_tcps_errvals,tcp_refps,tcp_width,tcp_ermin,tcp_ermax
@@ -1604,6 +1604,9 @@ module gsimod
 !                           = 0.42 meters (default)
 !      hwllp_howv    - real, background error de-correlation length scale of howv 
 !                           = 170,000.0 meters (default 170 km)
+!      corp_gust     - real, static background error of gust (stddev error)
+!      hwllp_gust    - real, background error de-correlation length scale of gust 
+!      oerr_gust     - real, observation error of gust
 !
   namelist/rapidrefresh_cldsurf/dfi_radar_latent_heat_time_period, &
                                 metar_impact_radius,metar_impact_radius_lowcloud, &
@@ -1625,7 +1628,7 @@ module gsimod
                                 cld_bld_coverage,cld_clr_coverage,&
                                 i_cloud_q_innovation,i_ens_mean,DTsTmax, &
                                 i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check, &
-                                corp_howv, hwllp_howv
+                                corp_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust
 
 ! chem(options for gsi chem analysis) :
 !     berror_chem       - .true. when background  for chemical species that require

src/gsi/m_berror_stats_reg.f90

@@ -313,6 +313,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt
       use mpeu_util,only: getindex
       use radiance_mod, only: icloud_cv,n_clouds_fwd,cloud_names_fwd
       use chemmod, only: berror_fv3_cmaq_regional,berror_fv3_sd_regional
+      use rapidrefresh_cldsurf_mod, only: corp_gust, hwllp_gust, l_rtma3d
 
       implicit none
 
@@ -825,11 +826,35 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt
 !       end if
      else if (n==nrf2_gust) then
         do i=1,mlat
-           corp(i,n)=three
+           corp(i,n)=three    ! background error stddev of wind gust = 3 m/s (default: legacy code from 2DRTMA) 
         end do
         do i=0,mlat+1
-           hwllp(i,n)=hwll(i,1,nrf3_q)
+           hwllp(i,n)=hwll(i,1,nrf3_q)  ! de-correlation length of bkgd error of gust is
+                                        ! same as the value of q at bottom level (default: legacy code from 2DRTMA)
+                                        ! for other DA apps, it is recommended to change it 
+                                        ! by setting hwllp_gust in GSI namelist.
         end do
+        if ( l_rtma3d ) then  ! For 3drtma only: allowing to change the stddev and 
+                              ! de-correlation length of bkgd error of gust:
+                              !   corp_gust : set in namelist(if <=0, using default value above (3.0)
+                              !   hwllp_gust: set in namelist(if <=0, using default value above (value of q)
+           if ( corp_gust  .gt. 0.0_r_kind ) then
+              corp(1:mlat,   n) = corp_gust
+              if (mype==0) write(6,'(1x,A,A,I5.5,A,F8.3)') &
+                           myname_,"@pe=",mype," (3drtma) set b_error stddev of gust = ",corp_gust
+           else
+              if (mype==0) write(6,'(1x,A,A,I5.5,A,F8.3)') &
+                           myname_,"@pe=",mype," (3drtma) set b_error stddev of gust (default) = ",three
+           end if
+           if ( hwllp_gust .gt. 0.0_r_kind ) then
+              hwllp(0:mlat+1,n) = hwllp_gust
+              if (mype==0) write(6,'(1x,A,A,I5.5,A,F12.3)') &
+                           myname_,"@pe=",mype," (3drtma) set b_error de-corr length of gust = ",hwllp_gust
+           else
+              if (mype==0) write(6,'(1x,A,A,I5.5,A)') &
+                           myname_,"@pe=",mype," (3drtma) set b_error de-corr length of gust is same as length of q."
+           end if
+        end if
      else if (n==nrf2_vis) then
         do i=1,mlat
            corp(i,n)=3.0_r_kind