Merge remote-tracking branch 'upstream/develop' into feature/hafsv2_b…

…aseline

src/gsi/constants.f90

@@ -31,7 +31,6 @@ module constants
 !   2016-02-15 Johnson, Y. Wang, X. Wang - define additional constant values for
 !                                          radar DA, POC: [email protected]
 !   2019-09-25  X.Su     - put stndrd_atmos_ps constant values
-!   2022-03-01  X.Lu & X.Wang  - increased max_varname_length for HAFS dual ens. POC: [email protected]
 !
 ! Subroutines Included:
 !   sub init_constants_derived - compute derived constants

src/gsi/crtm_interface.f90

@@ -977,7 +977,7 @@ end subroutine destroy_crtm
 subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
                    h,q,qs,clw_guess,ciw_guess,rain_guess,snow_guess,prsl,prsi, &
                    trop5,tzbgr,dtsavg,sfc_speed,&
-                   tsim,emissivity,ptau5,ts, &
+                   tsim,emissivity,chan_level,ptau5,ts, &
                    emissivity_k,temp,wmix,jacobian,error_status,tsim_clr,tcc, & 
                    tcwv,hwp_ratio,stability,layer_od,jacobian_aero)  
 !$$$  subprogram documentation block
@@ -1097,6 +1097,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
   real(r_kind)                          ,intent(  out) :: sfc_speed,dtsavg
   real(r_kind),dimension(nchanl+nreal)  ,intent(in   ) :: data_s
   real(r_kind),dimension(nchanl)        ,intent(  out) :: tsim,emissivity,ts,emissivity_k
+  real(r_kind),dimension(nchanl)        ,intent(  out) :: chan_level 
   character(10)                         ,intent(in   ) :: obstype
   integer(i_kind)                       ,intent(  out) :: error_status
   real(r_kind),dimension(nsig,nchanl)   ,intent(  out) :: temp,ptau5,wmix
@@ -1150,6 +1151,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
   real(r_kind):: sno00,sno01,sno10,sno11,secant_term
   real(r_kind):: hwp_total,theta_700,theta_sfc,hs
   real(r_kind):: dlon,dlat,dxx,dyy,yy,zz,garea
+  real(r_kind):: radiance, radiance_overcast, radiance_ratio
   real(r_kind),dimension(0:3):: wgtavg
   real(r_kind),dimension(nsig,nchanl):: omix
   real(r_kind),dimension(nsig,nchanl,n_aerosols_jac):: jaero
@@ -2217,8 +2219,10 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
        end do
     end if
 
+    chan_level = zero
+
 !$omp parallel do  schedule(dynamic,1) private(i) &
-!$omp private(total_od,k,kk,m,term,ii,cwj)
+!$omp private(total_od,k,kk,m,term,ii,cwj,radiance,radiance_overcast,radiance_ratio)
     do i=1,nchanl
 !   Zero jacobian and transmittance arrays
        do k=1,nsig
@@ -2228,6 +2232,16 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
          wmix(k,i)=zero
        end do
 
+       radiance=rtsolution(i,1)%radiance
+       do k=msig, 1, -1
+          radiance_overcast = rtsolution(i,1)%upwelling_overcast_radiance(k)
+          radiance_ratio = abs(radiance_overcast/radiance)
+          if (radiance_ratio < 0.99_r_kind) then
+             chan_level(i) = atmosphere(1)%pressure(k) / r10
+             exit
+          endif
+       enddo
+
 !  Simulated brightness temperatures
        tsim(i)=rtsolution(i,1)%brightness_temperature
 

src/gsi/gsi_files.cmake

@@ -101,6 +101,7 @@ bkgvar_rewgt.f90
 blacklist.f90
 blendmod.f90
 buddycheck_mod.f90
+cads.f90
 calc_fov_conical.f90
 calc_fov_crosstrk.f90
 calctends.f90

src/gsi/gsi_rfv3io_mod.f90

@@ -559,7 +559,7 @@ subroutine gsi_rfv3io_get_ens_grid_specs(grid_spec,ierr)
   integer(i_kind) :: nio,nylen
   integer(i_kind),allocatable :: gfile_loc_layout(:)
   character(len=180)  :: filename_layout
-  integer(i_kind) imiddle,jmiddle
+  integer(i_kind) imiddle,jmiddle,grid_ens_type_fv3_regional
 
 
     iret=nf90_open(trim(grid_spec),nf90_nowrite,gfile_grid_spec)
@@ -635,33 +635,19 @@ subroutine gsi_rfv3io_get_ens_grid_specs(grid_spec,ierr)
     if(mype==0)write(6,*),'ny_layout_bens=',ny_layout_bens
     if(mype==0)write(6,*),'ny_layout_eens=',ny_layout_eens
 
-!
-!  need to decide the grid orientation of the FV regional model   
-!
-!   grid_type_fv3_regional = 0 : decide grid orientation based on
-!                                grid_lat/grid_lon
-!                            1 : input is E-W N-S grid
-!                            2 : input is W-E S-N grid
-!
-    if(grid_type_fv3_regional == 0) then
-        imiddle=nxens/2
-        jmiddle=nyens/2
-        if( (grid_latt(imiddle,1) < grid_latt(imiddle,nyens)) .and. &
-            (grid_lont(1,jmiddle) < grid_lont(nxens,jmiddle)) ) then
-            grid_type_fv3_regional = 2
-        else
-            grid_type_fv3_regional = 1
-        endif
+    imiddle=nxens/2
+    jmiddle=nyens/2
+    if( (grid_latt(imiddle,1) < grid_latt(imiddle,nyens)) .and. &
+        (grid_lont(1,jmiddle) < grid_lont(nxens,jmiddle)) ) then
+        grid_ens_type_fv3_regional = 2
+    else
+        grid_ens_type_fv3_regional = 1
     endif
 ! check the grid type
-    if( grid_type_fv3_regional == 1 ) then
-       if(mype==0) write(6,*) 'FV3 regional input grid is  E-W N-S grid'
-       grid_reverse_flag=.true.    ! grid is revered comparing to usual map view
-    else if(grid_type_fv3_regional == 2) then
-       if(mype==0) write(6,*) 'FV3 regional input grid is  W-E S-N grid'
-       grid_reverse_flag=.false.   ! grid orientated just like we see on map view
+    if( grid_type_fv3_regional == grid_ens_type_fv3_regional ) then
+       if(mype==0) write(6,*) 'Ensemble has the same orientation as the control, Cool!'
     else
-       write(6,*) 'Error: FV3 regional input grid is unknown grid'
+       write(6,*) 'Warning! Ensemble has a different orientation as the control. This case needs further tests, Abort!'
        call stop2(678)
     endif
 !

src/gsi/gsimod.F90

@@ -93,8 +93,14 @@ module gsimod
       erradar_inflate,tdrerr_inflate,use_poq7,qc_satwnds,&
       init_qcvars,vadfile,noiqc,c_varqc,gps_jacqc,qc_noirjaco3,qc_noirjaco3_pole,&
       buddycheck_t,buddydiag_save,njqc,vqc,nvqc,hub_norm,vadwnd_l2rw_qc, &
-      pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres,cao_check
+      pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres,cao_check, &
+      cris_cads, iasi_cads, airs_cads
   use qcmod, only: troflg,lat_c,nrand
+  use cads, only: M__Sensor,N__Num_Bands,N__GradChkInterval,N__Band_Size,N__Bands,N__Window_Width, &
+      N__Window_Bounds,R__BT_Threshold,R__Grad_Threshold,R__Window_Grad_Threshold, L__Do_Quick_Exit, &
+      L__Do_CrossBand, N__BandToUse,L__Do_Imager_Cloud_Detection, N__Num_Imager_Chans, &
+      N__Num_Imager_Clusters,N__Imager_Chans,R__Stddev_Threshold,R__Coverage_Threshold, &
+      R__FG_Departure_Threshold, CADS_Setup_Cloud
   use pcpinfo, only: npredp,diag_pcp,dtphys,deltim,init_pcp
   use jfunc, only: iout_iter,iguess,miter,factqmin,factqmax,superfact,limitqobs, &
      factql,factqi,factqr,factqs,factqg, &  
@@ -507,6 +513,9 @@ module gsimod
 !                           2. fv3_regional =      .true.  
 !                           3. fv3_cmaq_regional = .true. 
 !                           4. berror_fv3_cmaq_regional = .true. 
+!  09-02-2022 Jung      Added namelist entries to call a new IR cloud detection routine
+!                       the original cloud detection routine is the default.  To use the new 
+!                       cloud detection routine, set the flags to .true.
 !  09-15-2022 yokota  - add scale/variable/time-dependent localization
 !  2023-07-30 Zhao    - added namelist options for analysis of significant wave height
 !                       (aka howv in GSI code): corp_howv, hwllp_howv
@@ -1051,6 +1060,13 @@ module gsimod
 !                                                wind observations.
 
 !     vad_near_analtime - assimilate newvadwnd obs around analysis time only
+!
+!     Flags to use the new IR cloud detection routine.  Flag must be set to true to use the new routine.  The default
+!     (no flag or .false.) will use the default.
+!     airs_cads: use the clod and aerosool detection software for the AIRS instrument
+!     cris_cads: use the cloud and aerosol detection software for CrIS instruments
+!     iasi_cads: use the cloud and aerosol detection software for IASI instruments
+!     
 
   namelist/obsqc/dfact,dfact1,erradar_inflate,tdrerr_inflate,oberrflg,&
        vadfile,noiqc,c_varqc,blacklst,use_poq7,hilbert_curve,tcp_refps,tcp_width,&
@@ -1061,7 +1077,7 @@ module gsimod
        q_doe_a_136,q_doe_a_137,q_doe_b_136,q_doe_b_137, &
        t_doe_a_136,t_doe_a_137,t_doe_b_136,t_doe_b_137, &
        uv_doe_a_236,uv_doe_a_237,uv_doe_a_213,uv_doe_b_236,uv_doe_b_237,uv_doe_b_213, &
-       vad_near_analtime
+       vad_near_analtime,airs_cads,cris_cads,iasi_cads
 
 ! OBS_INPUT (controls input data):
 !      dmesh(max(dthin))- thinning mesh for each group
@@ -1663,6 +1679,40 @@ module gsimod
 !     fac_tsl  - index to apply thermal skin layer or not: 0 = no; 1 = yes.
    namelist/nst/nst_gsi,nstinfo,zsea1,zsea2,fac_dtl,fac_tsl
 
+!  Initialize the Cloud and Aerosol Detection Software (CADS) 
+!
+!     M__Sensor                     Unique ID for sensor
+!     N__Num_Bands                  Number of channel bands
+!     N__Band_Size(:)               Number of channels in each band
+!     N__Bands(:,:)                 Channel lists
+!     N__Window_Width(:)            Smoothing filter window widths per band
+!     N__Window_Bounds(:,:)         Channels in the spectral window gradient check
+!     N__GradChkInterval(:)         Window width used in gradient calculation
+!     R__BT_Threshold(:)            BT threshold for cloud contamination
+!     R__Grad_Threshold(:)          Gradient threshold for cloud contamination
+!     R__Window_Grad_Threshold(:)   Threshold for window gradient check in QE
+!     L__Do_Quick_Exit              On/off switch for the Quick Exit scenario
+!     L__Do_CrossBand               On/off switch for the cross-band method
+!     N__BandToUse(:)               Band number assignment for each channel
+!     L__Do_Imager_Cloud_Detection  On/off switch for the imager cloud detection
+!     N__Num_Imager_Chans           No. of imager channels
+!     N__Num_Imager_Clusters        No. of clusters to be expected
+!     N__Imager_Chans(:)            List of imager channels
+!     R__Stddev_Threshold(:)        St. Dev. threshold, one for each imager channel
+!     R__Coverage_Threshold         Threshold for fractional coverage of a cluster
+!     R__FG_Departure_Threshold     Threshold for imager FG departure
+
+   NAMELIST / Cloud_Detect_Coeffs / M__Sensor, N__Num_Bands,            &
+           N__Band_Size, N__Bands, N__Window_Width, N__Window_Bounds,   &
+           N__GradChkInterval, R__BT_Threshold, R__Grad_Threshold,      &
+           R__Window_Grad_Threshold, L__Do_Quick_Exit,                  &
+           L__Do_CrossBand, N__BandToUse,                               &
+           L__Do_Imager_Cloud_Detection, N__Num_Imager_Chans,           &
+           N__Num_Imager_Clusters, N__Imager_Chans,                     &
+           R__Stddev_Threshold, R__Coverage_Threshold,                  &
+           R__FG_Departure_Threshold
+
+
 !EOC
 
 !---------------------------------------------------------------------------
@@ -1749,6 +1799,7 @@ subroutine gsimain_initialize
   call set_fgrid2agrid
   call gsi_nstcoupler_init_nml
   call init_radaruse_directDA
+  call CADS_Setup_Cloud
 
  if(mype==0) write(6,*)' at 0 in gsimod, use_gfs_stratosphere,nems_nmmb_regional = ', &
                        use_gfs_stratosphere,nems_nmmb_regional

src/gsi/hybrid_ensemble_isotropic.F90

@@ -1393,10 +1393,6 @@ subroutine load_ensemble
                 endif
 
 !     regional_ensemble_option = 5: ensembles are fv3 regional.
-                if (l_both_fv3sar_gfs_ens) then ! first read in gfs ensembles for regional
-                   call get_gefs_for_regional
-!clthink do we need mpi_bar here?
-                endif
                 call fv3_regional_enspert%get_fv3_regional_ensperts(en_perts,nelen,ps_bar)
 
 
@@ -4042,7 +4038,7 @@ subroutine hybens_grid_setup
   logical,allocatable::vector(:)
   real(r_kind) eps,r_e
   real(r_kind) rlon_a(nlon),rlat_a(nlat),rlon_e(nlon),rlat_e(nlat)
-  character(:),allocatable:: fv3_spec_grid_filename
+  character(:),allocatable:: fv3_ens_spec_grid_filename
   integer :: ierr
 
   nord_e2a=4       !   soon, move this to hybrid_ensemble_parameters
@@ -4130,9 +4126,9 @@ subroutine hybens_grid_setup
   else
      if(dual_res) then
         call get_region_dx_dy_ens(region_dx_ens,region_dy_ens)
-        if(regional_ensemble_option) then
-           fv3_spec_grid_filename="fv3_ens_grid_spec"
-           call gsi_rfv3io_get_ens_grid_specs(fv3_spec_grid_filename,ierr)
+        if(regional_ensemble_option == 5) then
+           fv3_ens_spec_grid_filename="fv3_ens_grid_spec"
+           call gsi_rfv3io_get_ens_grid_specs(fv3_ens_spec_grid_filename,ierr)
         endif
      else
         region_dx_ens=region_dx