Cads for andrew (#616)

NOAA-EMC · Dec 4, 2023 · ea667d9 · ea667d9
1 parent c94bc72
commit ea667d9
Show file tree

Hide file tree

Showing 10 changed files with 2,980 additions and 196 deletions.
diff --git a/src/gsi/cads.f90 b/src/gsi/cads.f90
diff --git a/src/gsi/crtm_interface.f90 b/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
 

diff --git a/src/gsi/gsi_files.cmake b/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

diff --git a/src/gsi/gsimod.F90 b/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