Merge pull request #702 from PUMASDevelopment/katetc/rainbows_pr

cam6_3_159: Diagnostic rainbows and GPU fixes in PUMAS microphysics
Fixes #683 Add Diagnostic Rainbow Calculation
Partially addresses #1007 Broken PUMAS GPU code and GPU regression test by bringing in an updated PUMAS external tag with fixed GPU code. However, the needed ccs_config tag ( ccs_config_cesm0.0.99 ) caused a failure in the single column tests right now. So the updates to the ccs_config external will be left for the next GPU fixing tag.

RAINBOWS!!! I'M NOT KIDDING!!!!111!
Also a new PUMAS tag that fixes GPU OpenACC calls. B4b

Does not change answers.

Externals_CAM.cfg

@@ -32,7 +32,7 @@ required = True
 local_path = src/physics/pumas
 protocol = git
 repo_url = https://github.com/ESCOMP/PUMAS
-tag = pumas_cam-release_v1.35
+tag = pumas_cam-release_v1.36
 required = True
 
 [pumas-frozen]

cime_config/testdefs/testlist_cam.xml

@@ -1353,6 +1353,7 @@
   <test compset="QPC5" grid="T42_T42_mg17" name="SCT_D_Ln7" testmods="cam/scm_prep">
     <machines>
       <machine name="derecho" compiler="intel" category="aux_cam"/>
+      <machine name="derecho" compiler="intel" category="prealpha"/>
     </machines>
     <options>
       <option name="comment" >001 tsc</option>
@@ -1520,6 +1521,15 @@
       <option name="wallclock">01:00:00</option>
       <option name="comment">Test PUMAS ML emulated warm rain option</option>
     </options>
+  </test>
+    <test compset="F2000dev" grid="ne30pg3_ne30pg3_mg17" name="SMS_D_Ln9" testmods="cam/outfrq9s_pumas_rainbows" supported="false">
+    <machines>
+      <machine name="derecho" compiler="intel" category="aux_pumas"/>
+    </machines>
+    <options>
+      <option name="wallclock">01:00:00</option>
+      <option name="comment">Test PUMAS diagnostics rainbows output</option>
+    </options>
   </test>
   <test compset="F2000dev" grid="f09_f09_mg17" name="ERP_Ln9" testmods="cam/outfrq9s_mg3" supported="false">
     <machines>

cime_config/testdefs/testmods_dirs/cam/outfrq9s_pumas_rainbows/shell_commands

@@ -0,0 +1,3 @@
+./xmlchange ROF_NCPL=\$ATM_NCPL
+./xmlchange GLC_NCPL=\$ATM_NCPL
+./xmlchange CAM_CONFIG_OPTS=' -microphys pumas' --append

cime_config/testdefs/testmods_dirs/cam/outfrq9s_pumas_rainbows/user_nl_cam

@@ -0,0 +1,5 @@
+mfilt=1,1,1,1,1,1
+ndens=1,1,1,1,1,1
+nhtfrq=9,9,9,9,9,9
+inithist='ENDOFRUN'
+fincl1 = 'RBFRAC','RBFREQ','rbSZA'

cime_config/testdefs/testmods_dirs/cam/outfrq9s_pumas_rainbows/user_nl_clm

@@ -0,0 +1,26 @@
+!----------------------------------------------------------------------------------
+! Users should add all user specific namelist changes below in the form of 
+! namelist_var = new_namelist_value 
+!
+! Include namelist variables for drv_flds_in ONLY if -megan and/or -drydep options
+! are set in the CLM_NAMELIST_OPTS env variable.
+!
+! EXCEPTIONS: 
+! Set use_cndv           by the compset you use and the CLM_BLDNML_OPTS -dynamic_vegetation setting
+! Set use_vichydro       by the compset you use and the CLM_BLDNML_OPTS -vichydro           setting
+! Set use_cn             by the compset you use and CLM_BLDNML_OPTS -bgc  setting
+! Set use_crop           by the compset you use and CLM_BLDNML_OPTS -crop setting
+! Set spinup_state       by the CLM_BLDNML_OPTS -bgc_spinup      setting
+! Set irrigate           by the CLM_BLDNML_OPTS -irrig           setting
+! Set dtime              with L_NCPL                             option
+! Set fatmlndfrc         with LND_DOMAIN_PATH/LND_DOMAIN_FILE    options
+! Set finidat            with RUN_REFCASE/RUN_REFDATE/RUN_REFTOD options for hybrid or branch cases
+!                        (includes $inst_string for multi-ensemble cases)
+! Set glc_grid           with CISM_GRID                          option
+! Set glc_smb            with GLC_SMB                            option
+! Set maxpatch_glcmec    with GLC_NEC                            option
+! Set glc_do_dynglacier  with GLC_TWO_WAY_COUPLING               env variable
+!----------------------------------------------------------------------------------
+hist_nhtfrq = 9
+hist_mfilt  = 1
+hist_ndens  = 1

doc/ChangeLog

@@ -1,3 +1,68 @@
+===============================================================
+
+Tag name: cam6_3_159
+Originator(s): katetc, andrewgettelman, sjsprecious
+Date: April 26, 2024
+One-line Summary: Diagnostic rainbows and new PUMAS external with fixed GPU directives
+Github PR URL: https://github.com/ESCOMP/CAM/pull/702
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+ - Add diagnostic rainbow calculation: https://github.com/ESCOMP/CAM/issues/683
+ - Partially addresses Broken PUMAS GPU code and GPU regression test: https://github.com/ESCOMP/CAM/issues/1007
+
+Describe any changes made to build system: N/A
+
+Describe any changes made to the namelist: N/A
+
+List any changes to the defaults for the boundary datasets: N/A
+
+Describe any substantial timing or memory changes: N/A
+
+Code reviewed by: sjsprecious, andrewgettelman, cacraigucar, nusbaume
+
+List all files eliminated: N/A
+
+List all files added and what they do: N/A
+
+List all existing files that have been modified, and describe the changes:
+M	Externals_CAM.cgf
+	- Point to new PUMAS tag
+
+M	cime_config/testdefs/testlist_cam.xml
+	- Add new test for rainbows output to aux_pumas suite
+	- Add SCT test to prealpha test suite to ensure it is not broken by the next ccs_config tags
+
+M	cime_config/testdefs/testmods_dir/cam/outfrq9s_pumas_rainbows/shell_comands
+	cime_config/testdefs/testmods_dir/cam/outfrq9s_pumas_rainbows/user_nl_cam
+	cime_config/testdefs/testmods_dir/cam/outfrq9s_pumas_rainbows/user_nl_clm
+	- Added test to the aux_pumas suite to make sure rainbows functionality is maintained
+
+M	src/physics/cam/micro_pumas_cam.F90
+	src/physics/cam_dev/micro_pumas_cam.F90
+	- Diagnostic rainbows parameterization added identically in both versions of the file.
+
+If there were any failures reported from running test_driver.sh on any test
+platform, and checkin with these failures has been OK'd by the gatekeeper,
+then copy the lines from the td.*.status files for the failed tests to the
+appropriate machine below.  All failed tests must be justified.
+
+derecho/intel/aux_cam: all BFB, except:
+  ERP_Ln9.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq9s (Overall: FAIL) details:
+    - pre-existing failure
+
+  SMS_Ld1.f09_f09_mg17.FCHIST_GC.derecho_intel.cam-outfrq1d (Overall: DIFF) details:
+    FAIL SMS_Ld1.f09_f09_mg17.FCHIST_GC.derecho_intel.cam-outfrq1d BASELINE /glade/campaign/cesm/community/amwg/cam_baselines/cam6_3_158: DIFF
+  SMS_Lh12.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq3h (Overall: DIFF) details:
+    FAIL SMS_Lh12.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq3h BASELINE /glade/campaign/cesm/community/amwg/cam_baselines/cam6_3_158: DIFF
+    - Unexpected baseline comparison failures. Documented in ESCOMP/cam issue #1018
+
+izumi/nag/aux_cam: all B4B, except:
+  DAE.f45_f45_mg37.FHS94.izumi_nag.cam-dae (Overall: FAIL) details:
+     - pre-existing failure
+
+izumi/gnu/aux_cam: all BFB
+
+Summarize any changes to answers: bit-for-bit unchanged except GEOS-Chem and HEMCO tests described in issue #1018
 
 ===============================================================
 
@@ -236,7 +301,7 @@ Summarize any changes to answers: bit-for-bit unchanged
 
 Tag name: cam6_3_155
 Originator(s): katec,vlarson,bstephens82,huebleruwm,zarzycki,JulioTBacmeister
-Date: April 11, 2024
+Date: April 12, 2024
 One-line Summary: Update CLUBB and SILHS to new UWM external
 Github PR URL: https://github.com/ESCOMP/CAM/pull/960
 

src/physics/cam/micro_pumas_cam.F90

@@ -12,8 +12,11 @@ module micro_pumas_cam
 use physconst,      only: gravit, rair, tmelt, cpair, rh2o, rhoh2o, &
                           latvap, latice, mwh2o
 use phys_control,   only: phys_getopts, use_hetfrz_classnuc
-
-
+use shr_const_mod,  only: pi => shr_const_pi
+use time_manager,   only: get_curr_date, get_curr_calday
+use phys_grid,      only: get_rlat_all_p, get_rlon_all_p
+use orbit,          only: zenith
+
 use physics_types,  only: physics_state, physics_ptend, &
                           physics_ptend_init, physics_state_copy, &
                           physics_update, physics_state_dealloc, &
@@ -197,6 +200,8 @@ module micro_pumas_cam
    ast_idx = -1,            &
    cld_idx = -1,            &
    concld_idx = -1,         &
+   prec_dp_idx = -1,        &
+   prec_sh_idx = -1,        &
    qsatfac_idx = -1
 
 ! Pbuf fields needed for subcol_SILHS
@@ -1019,6 +1024,10 @@ subroutine micro_pumas_cam_init(pbuf2d)
 
    end if
 
+   call addfld ('RBFRAC',  horiz_only, 'A', 'Fraction',  'Fraction of sky covered by a potential rainbow' )
+   call addfld ('RBFREQ',  horiz_only, 'A', 'Frequency',  'Potential rainbow frequency' )
+   call addfld( 'rbSZA', horiz_only, 'I', 'degrees', 'solar zenith angle' )
+
    ! History variables for CAM5 microphysics
    call addfld ('MPDT',       (/ 'lev' /), 'A', 'W/kg',     'Heating tendency - Morrison microphysics'                )
    call addfld ('MPDQ',       (/ 'lev' /), 'A', 'kg/kg/s',  'Q tendency - Morrison microphysics'                      )
@@ -1274,6 +1283,8 @@ subroutine micro_pumas_cam_init(pbuf2d)
    ast_idx      = pbuf_get_index('AST')
    cld_idx      = pbuf_get_index('CLD')
    concld_idx   = pbuf_get_index('CONCLD')
+   prec_dp_idx  = pbuf_get_index('PREC_DP')
+   prec_sh_idx  = pbuf_get_index('PREC_SH')
 
    naai_idx     = pbuf_get_index('NAAI')
    naai_hom_idx = pbuf_get_index('NAAI_HOM')
@@ -1603,7 +1614,10 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf)
 
    real(r8), pointer :: cld(:,:)          ! Total cloud fraction
    real(r8), pointer :: concld(:,:)       ! Convective cloud fraction
-   real(r8), pointer :: iciwpst(:,:)      ! Stratiform in-cloud ice water path for radiation
+   real(r8), pointer :: prec_dp(:)      ! Deep Convective precip
+   real(r8), pointer :: prec_sh(:)      ! Shallow Convective precip    
+
+   real(r8), pointer :: iciwpst(:,:) ! Stratiform in-cloud ice water path for radiation
    real(r8), pointer :: iclwpst(:,:)      ! Stratiform in-cloud liquid water path for radiation
    real(r8), pointer :: cldfsnow(:,:)     ! Cloud fraction for liquid+snow
    real(r8), pointer :: icswp(:,:)        ! In-cloud snow water path
@@ -1834,6 +1848,34 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf)
    real(r8), parameter :: mucon  = 5.3_r8            ! Convective size distribution shape parameter
    real(r8), parameter :: deicon = 50._r8            ! Convective ice effective diameter (meters)
 
+   ! Rainbows: solar zenith angle (SZA)
+   real(r8) :: zen_angle(state%psetcols)                   ! Daytime solar zenith angles (radians)
+   real(r8) :: rlats(state%psetcols), rlons(state%psetcols)   ! chunk latitudes and longitudes (radains)
+   real(r8) :: sza(state%psetcols)                         ! solar zenith angles (degrees)
+   real(r8), parameter :: rad2deg = 180._r8/pi                ! radians to degrees conversion factor
+   real(r8) :: calday  !current calendar day
+
+   real(r8) :: precc(state%psetcols)                ! convective precip rate
+
+! Rainbow frequency and fraction for output
+
+   real(r8) :: rbfreq(state%psetcols)
+   real(r8) :: rbfrac(state%psetcols)
+
+!Rainbows: parameters
+
+   real(r8), parameter :: rb_rmin  =1.e-6_r8                    ! Strat Rain threshold (mixing ratio)
+   real(r8), parameter :: rb_rcmin = 5._r8/(86400._r8*1000._r8) ! Conv Rain Threshold (mm/d--> m/s)
+   real(r8), parameter :: rb_pmin =85000._r8                       ! Minimum pressure for surface layer
+   real(r8), parameter :: deg2rad = pi/180._r8            ! Conversion factor
+   integer :: top_idx  !Index for top level below rb_pmin
+   real(r8) :: convmx
+   real(r8) :: cldmx
+   real(r8) :: frlow
+   real(r8) :: cldtot
+   real(r8) :: rmax
+   logical :: rval
+
    !-------------------------------------------------------------------------------
 
    lchnk = state%lchnk
@@ -1871,6 +1913,29 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf)
    call pbuf_get_field(pbuf, ast_idx,         ast,     start=(/1,1,itim_old/), kount=(/psetcols,pver,1/), &
         col_type=col_type, copy_if_needed=use_subcol_microp)
 
+   ! Get convective precip
+   if (prec_dp_idx > 0) then
+      call pbuf_get_field(pbuf, prec_dp_idx, prec_dp, col_type=col_type, copy_if_needed=use_subcol_microp)
+   else
+      nullify(prec_dp)
+   end if
+   if (prec_sh_idx > 0) then
+      call pbuf_get_field(pbuf, prec_sh_idx, prec_sh, col_type=col_type, copy_if_needed=use_subcol_microp)
+   else
+      nullify(prec_sh)
+   end if
+
+! Merge Precipitation rates (multi-process)
+   if (associated(prec_dp) .and. associated(prec_sh)) then
+      precc(:ncol) = prec_dp(:ncol)  + prec_sh(:ncol)
+   else if (associated(prec_dp)) then
+      precc(:ncol) = prec_dp(:ncol)
+   else if (associated(prec_sh)) then
+      precc(:ncol) = prec_sh(:ncol)
+   else
+      precc(:ncol) = 0._r8
+   end if
+
    if (.not. do_cldice) then
       ! If we are NOT prognosing ice and snow tendencies, then get them from the Pbuf
       call pbuf_get_field(pbuf, tnd_qsnow_idx,   tnd_qsnow,   col_type=col_type, copy_if_needed=use_subcol_microp)
@@ -2043,6 +2108,26 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf)
    call pbuf_get_field(pbuf, evpsnow_st_idx,  evpsnow_st_grid)
    call pbuf_get_field(pbuf, am_evp_st_idx,   am_evp_st_grid)
 
+   !-----------------------------------------------------------------------    
+   !        ... Calculate cosine of zenith angle
+   !            then cast back to angle (radians)
+   !-----------------------------------------------------------------------    
+   zen_angle(:) = 0.0_r8
+   rlats(:) = 0.0_r8
+   rlons(:) = 0.0_r8
+   calday = get_curr_calday()
+   call get_rlat_all_p( lchnk, ncol, rlats )
+   call get_rlon_all_p( lchnk, ncol, rlons )
+   call zenith( calday, rlats, rlons, zen_angle, ncol )
+   where (zen_angle(:) <= 1.0_r8 .and. zen_angle(:) >= -1.0_r8)
+      zen_angle(:) = acos( zen_angle(:) )
+   elsewhere
+      zen_angle(:) = 0.0_r8
+   end where
+
+   sza(:) = zen_angle(:) * rad2deg
+   call outfld( 'rbSZA',   sza,    ncol, lchnk )
+
    !-------------------------------------------------------------------------------------
    ! Microphysics assumes 'liquid stratus frac = ice stratus frac
    !                      = max( liquid stratus frac, ice stratus frac )'.
@@ -2134,6 +2219,10 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf)
       state_loc_numgraup(:ncol,:) = 0._r8
    end if
 
+   ! Zero out diagnostic rainbow arrays
+   rbfreq = 0._r8
+   rbfrac = 0._r8   
+
    ! Zero out values above top_lev before passing into _tend for some pbuf variables that are inputs
    naai(:ncol,:top_lev-1) = 0._r8
    npccn(:ncol,:top_lev-1) = 0._r8
@@ -3123,6 +3212,63 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf)
 
    racau_grid = min(racau_grid, 1.e10_r8)
 
+!-----------------------------------------------------------------------
+! Diagnostic Rainbow Calculation. Seriously.
+!-----------------------------------------------------------------------
+
+! Rainbows currently calculated on the grid, not subcolumn specific
+   do i = 1, ngrdcol
+
+      top_idx = pver
+      convmx = 0._r8
+      frlow  = 0._r8
+      cldmx  = 0._r8
+      cldtot  = maxval(ast(i,top_lev:))
+
+! Find levels in surface layer
+      do k = top_lev, pver
+         if (state%pmid(i,k) > rb_pmin) then     
+            top_idx = min(k,top_idx)
+         end if 
+      end do 
+
+!For all fractional precip calculated below, use maximum in surface layer.
+!For convective precip, base on convective cloud area
+      convmx = maxval(concld(i,top_idx:))
+!For stratiform precip, base on precip fraction
+      cldmx= maxval(freqr(i,top_idx:))
+! Combine and use maximum of strat or conv fraction
+      frlow= max(cldmx,convmx)
+
+!max precip
+      rmax=maxval(qrout_grid(i,top_idx:))
+
+!     Stratiform precip mixing ratio OR some convective precip
+!      (rval = true  if any sig precip)
+
+      rval = ((precc(i) > rb_rcmin) .or. (rmax > rb_rmin))
+
+!Now can find conditions for a rainbow:
+! Maximum cloud cover (CLDTOT) < 0.5
+! 48 < SZA < 90
+! freqr (below rb_pmin) > 0.25
+! Some rain (liquid > 1.e-6 kg/kg, convective precip > 1.e-7 m/s 
+
+      if ((cldtot < 0.5_r8) .and. (sza(i) > 48._r8) .and. (sza(i) < 90._r8) .and. rval) then  
+
+!Rainbow 'probability' (area) derived from solid angle theory 
+!as the fraction of the hemisphere for a spherical cap with angle phi=sza-48.
+! This is only valid between 48 < sza < 90 (controlled for above).
+
+          rbfrac(i) =  max(0._r8,(1._r8-COS((sza(i)-48._r8)*deg2rad))/2._r8) * frlow
+          rbfreq(i) =  1.0_r8 
+      end if 
+
+   end do                    ! end column loop for rainbows
+
+   call outfld('RBFRAC',   rbfrac,      psetcols, lchnk, avg_subcol_field=use_subcol_microp)
+   call outfld('RBFREQ',   rbfreq,      psetcols, lchnk, avg_subcol_field=use_subcol_microp)   
+
    ! --------------------- !
    ! History Output Fields !
    ! --------------------- !