Recompute pgam, lamc, degrau, dei for radiation from updated state in pumas_post_main

Mirrors CAM micro_pumas_cam post-call recomputation, which discards the
per-substep values PUMAS returns. Match dei standard name with RRTMGP and
the CAM-SIMA registry (singular crystal). Run pumas_post_main before
apply_heating_rate (air density uses the pre-heating temperature) and add
the low stratiform cloud fraction optics limiter to the suite.

Author: Haipeng Lin

schemes/pumas/micro_pumas_ccpp_dimensions_post.meta

@@ -464,7 +464,7 @@
   kind = pumas_r8
   intent = in
 [deffi]
-  standard_name = effective_diameter_of_stratiform_cloud_ice_crystals_for_radiation
+  standard_name = effective_diameter_of_stratiform_cloud_ice_crystal_for_radiation
   long_name = effective diameter of stratiform cloud ice particles for radiation
   units = um
   dimensions = (horizontal_loop_extent, vertical_layer_dimension)

schemes/pumas/micro_pumas_ccpp_optics_limiter.meta

@@ -53,7 +53,7 @@
   kind = kind_phys
   intent = inout
 [deffi]
-  standard_name = effective_diameter_of_stratiform_cloud_ice_particle_for_radiation
+  standard_name = effective_diameter_of_stratiform_cloud_ice_crystal_for_radiation
   long_name = effective diameter of stratiform cloud ice particles for radiation
   units = um
   dimensions = (horizontal_loop_extent, vertical_layer_dimension)

schemes/pumas/pumas_post_main.F90

@@ -9,29 +9,53 @@ module pumas_post_main
 
   !> \section arg_table_pumas_post_main_run Argument Table
   !! \htmlinclude pumas_post_main_run.html
-  subroutine pumas_post_main_run(ncol, cldice, numice, strat_cldfrc, effi, errmsg, errcode)
+  subroutine pumas_post_main_run(ncol, nlev, trop_cloud_top_lev, rair, pmid, temp, &
+                cldliq, numliq, cldice, numice, graupel, numgraup, strat_cldfrc,   &
+                effi, dei, pgam, lamc, degrau, errmsg, errcode)
     use ccpp_kinds,        only: kind_phys
     use pumas_kinds,       only: pumas_r8=>kind_r8
-    use micro_pumas_utils, only: size_dist_param_basic, mg_ice_props, mincld, qsmall
+    use micro_pumas_utils, only: size_dist_param_basic, size_dist_param_liq, &
+                                 mg_ice_props, mg_liq_props, avg_diameter,   &
+                                 mincld, qsmall, rhog, rhows, rhoi
 
     integer,                 intent(in)  :: ncol
+    integer,                 intent(in)  :: nlev
+    integer,                 intent(in)  :: trop_cloud_top_lev !Index of the top model level for which
+                                                               !tropospheric cloud physics is run (index)
+    real(kind_phys),                 intent(in)  :: rair        ! gas constant of dry air (J kg-1 K-1)
+    real(kind_phys), dimension(:,:), intent(in)  :: pmid        ! layer midpoint pressure (Pa)
+    real(kind_phys), dimension(:,:), intent(in)  :: temp        ! air temperature, before microphysics heating is applied (K)
+    real(kind_phys), dimension(:,:), intent(in)  :: cldliq      ! updated cloud liquid water mixing ratio (kg/kg)
+    real(kind_phys), dimension(:,:), intent(in)  :: numliq      ! updated cloud droplet number concentration (kg-1)
     real(kind_phys), dimension(:,:), intent(in)  :: cldice      ! updated cloud ice mixing ratio (kg/kg)
     real(kind_phys), dimension(:,:), intent(in)  :: numice      ! updated cloud ice number concentration (kg-1)
+    real(kind_phys), dimension(:,:), intent(in)  :: graupel     ! updated graupel mixing ratio (kg/kg)
+    real(kind_phys), dimension(:,:), intent(in)  :: numgraup    ! updated graupel number concentration (kg-1)
     real(kind_phys), dimension(:,:), intent(in)  :: strat_cldfrc! total stratiform cloud area fraction (= ast)
     real(kind_phys), dimension(:,:), intent(out) :: effi       ! ice effective radius (micron)
+    real(kind_phys), dimension(:,:), intent(out) :: dei        ! ice effective diameter for radiation (micron)
+    real(kind_phys), dimension(:,:), intent(out) :: pgam       ! droplet size distribution shape parameter for radiation (1)
+    real(kind_phys), dimension(:,:), intent(out) :: lamc       ! droplet size distribution slope for radiation (1)
+    real(kind_phys), dimension(:,:), intent(out) :: degrau     ! graupel effective diameter for radiation (m)
     character(len=512),          intent(out) :: errmsg
     integer,                 intent(out) :: errcode
 
-    integer        :: k, nlev
-    real(pumas_r8) :: icimrst(ncol, size(cldice,2))  ! in-cloud (grid-mean) ice mixing ratio
-    real(pumas_r8) :: niic(ncol, size(cldice,2))     ! in-cloud (grid-mean) ice number conc
-    real(pumas_r8) :: rei(ncol, size(cldice,2))      ! ice slope param, then effective radius
+    integer        :: k
+    real(pumas_r8) :: icimrst(ncol, nlev)  ! in-cloud (grid-mean) ice mixing ratio
+    real(pumas_r8) :: niic(ncol, nlev)     ! in-cloud (grid-mean) ice number conc
+    real(pumas_r8) :: rei(ncol, nlev)      ! ice slope param, then effective radius
+    real(pumas_r8) :: icwmrst(ncol, nlev)  ! in-cloud (grid-mean) liquid mixing ratio
+    real(pumas_r8) :: ncic(ncol, nlev)     ! in-cloud (grid-mean) droplet number conc
+    real(pumas_r8) :: rho(ncol, nlev)      ! air density
+    real(pumas_r8) :: mu(ncol, nlev)       ! droplet size distribution shape parameter
+    real(pumas_r8) :: lambdac(ncol, nlev)  ! droplet size distribution slope
+    real(pumas_r8) :: qg(ncol, nlev)       ! grid-mean graupel mixing ratio
+    real(pumas_r8) :: ng(ncol, nlev)       ! grid-mean graupel number conc
+    real(pumas_r8) :: dgout2(ncol, nlev)   ! mean graupel particle diameter
 
     errmsg = ' '
     errcode = 0
 
-    nlev = size(cldice, 2)
-
     ! Ice effective radius is recomputed here from the post-microphysics
     ! grid-mean in-cloud ice, NOT from the raw per-substep value PUMAS returns.
     ! PUMAS uses the total stratiform fraction for ice cloud (icecldf = ast).
@@ -54,6 +78,59 @@ subroutine pumas_post_main_run(ncol, cldice, numice, strat_cldfrc, effi, errmsg,
 
     effi(:ncol,:) = real(rei(:,:), kind_phys)
 
+    ! Ice effective diameter for radiation follows from the recomputed
+    ! effective radius. Mirrors CAM micro_pumas_cam.F90:3080 (dei).
+    dei(:ncol,:) = real(rei(:,:) * rhoi/rhows * 2._pumas_r8, kind_phys)
+
+    ! Droplet size distribution parameters for radiation are likewise recomputed
+    ! from the post-microphysics grid-mean in-cloud liquid. Air density uses the
+    ! temperature BEFORE the microphysics heating is applied (in CAM the input
+    ! state, not the substepped local state, "to preserve answers").
+    ! PUMAS uses the total stratiform fraction for liquid cloud (liqcldf = ast).
+    ! Mirrors CAM micro_pumas_cam.F90:2919 (rho), :2621 (icwmrst), :2957 (ncic),
+    ! :2952-2978 (mu, lambdac).
+    rho(:,:)     = real(pmid(:ncol,:), pumas_r8) / &
+                (real(rair, pumas_r8) * real(temp(:ncol,:), pumas_r8))
+    icwmrst(:,:) = min(real(cldliq(:ncol,:), pumas_r8) / &
+                max(mincld, real(strat_cldfrc(:ncol,:), pumas_r8)), 0.005_pumas_r8)
+    ncic(:,:)   =    real(numliq(:ncol,:), pumas_r8) / &
+                max(mincld, real(strat_cldfrc(:ncol,:), pumas_r8))
+
+    mu(:,:)      = 0._pumas_r8
+    lambdac(:,:) = 0._pumas_r8
+    do k = trop_cloud_top_lev, nlev
+      call size_dist_param_liq(mg_liq_props, icwmrst(:,k), ncic(:,k), rho(:,k), &
+                               mu(:,k), lambdac(:,k), ncol)
+    end do
+
+    ! size_dist_param_liq flags no-cloud points with mu = -100; CAM resets
+    ! those to zero (lambdac is already zero there).
+    where (icwmrst(:,trop_cloud_top_lev:) < qsmall)
+      mu(:,trop_cloud_top_lev:) = 0._pumas_r8
+    end where
+
+    pgam(:ncol,:) = real(mu(:,:), kind_phys)
+    lamc(:ncol,:) = real(lambdac(:,:), kind_phys)
+
+    ! Graupel effective diameter for radiation, from the post-microphysics
+    ! grid-mean graupel. Mirrors CAM micro_pumas_cam.F90:3043-3056 (degrau).
+    qg(:,:) = real(graupel(:ncol,:), pumas_r8)
+    ng(:,:) = real(numgraup(:ncol,:), pumas_r8)
+
+    dgout2(:,:) = 0._pumas_r8
+    where (qg(:,trop_cloud_top_lev:) >= 1.e-7_pumas_r8)
+      dgout2(:,trop_cloud_top_lev:) = avg_diameter( &
+           qg(:,trop_cloud_top_lev:), &
+           ng(:,trop_cloud_top_lev:) * rho(:,trop_cloud_top_lev:), &
+           rho(:,trop_cloud_top_lev:), rhog)
+    end where
+
+    degrau(:ncol,:) = 0._kind_phys
+    where (qg(:,trop_cloud_top_lev:) >= 1.e-7_pumas_r8)
+      degrau(:ncol,trop_cloud_top_lev:) = real(dgout2(:,trop_cloud_top_lev:) * &
+           3._pumas_r8 * rhog/rhows, kind_phys)
+    end where
+
   end subroutine pumas_post_main_run
 
 end module pumas_post_main

schemes/pumas/pumas_post_main.meta

@@ -13,6 +13,62 @@
   dimensions = ()
   type = integer
   intent = in
+[nlev]
+  standard_name = vertical_layer_dimension
+  long_name = vertical layer dimension
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[trop_cloud_top_lev]
+  standard_name = vertical_layer_index_of_troposphere_cloud_physics_top
+  long_name = vertical layer index of troposphere cloud physics top
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[rair]
+  standard_name = gas_constant_of_dry_air
+  long_name = ideal gas constant for dry air
+  units = J kg-1 K-1
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[pmid]
+  standard_name = air_pressure
+  long_name = air pressure
+  units = Pa
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[temp]
+  standard_name = air_temperature
+  long_name = air temperature
+  units = K
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[cldliq]
+  standard_name = cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water
+  long_name = cloud liquid water mixing ratio wrt moist air and condensed water
+  units = kg kg-1
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+  advected = true
+[numliq]
+  standard_name = mass_number_concentration_of_cloud_liquid_water_droplets_in_moist_air_and_condensed_water
+  long_name = mass number concentration of cloud liquid wrt moist air and condensed water
+  units = kg-1
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+  advected = true
 [cldice]
   standard_name = cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water
   long_name = cloud ice mixing ratio wrt moist air and condensed water
@@ -31,6 +87,24 @@
   kind = kind_phys
   intent = in
   advected = true
+[graupel]
+  standard_name = graupel_water_mixing_ratio_wrt_moist_air_and_condensed_water
+  long_name = graupel mixing ratio wrt moist air and condensed water
+  units = kg kg-1
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+  advected = true
+[numgraup]
+  standard_name = mass_number_concentration_of_graupel_particles_in_moist_air_and_condensed_water
+  long_name = mass number concentration of graupel wrt moist air and condensed water
+  units = kg-1
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+  advected = true
 [strat_cldfrc]
   standard_name = stratiform_cloud_area_fraction
   long_name = stratiform cloud area fraction
@@ -47,6 +121,38 @@
   type = real
   kind = kind_phys
   intent = out
+[dei]
+  standard_name = effective_diameter_of_stratiform_cloud_ice_crystal_for_radiation
+  long_name = effective diameter of stratiform cloud ice particles for radiation
+  units = um
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = out
+[pgam]
+  standard_name = size_distribution_shape_parameter_for_microphysics
+  long_name = cloud particle size distribution shape (gamma) parameter
+  units = 1
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = out
+[lamc]
+  standard_name = slope_of_droplet_distribution_for_optics
+  long_name = cloud particle size distribution slope (lambda) parameter
+  units = 1
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = out
+[degrau]
+  standard_name = effective_diameter_of_stratiform_cloud_graupel_particle_for_radiation
+  long_name = effective diameter of stratiform graupel particles for radiation
+  units = m
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = out
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP

test/test_suites/suite_pumas.xml

@@ -21,10 +21,13 @@
     <!-- Apply tendencies -->
     <scheme>apply_constituent_tendencies</scheme>
     <scheme>clamp_number_concentrations</scheme>
+    <!-- Recompute size distribution quantities for radiation from the updated
+         constituents; needs the pre-heating air temperature, so it must run
+         before apply_heating_rate -->
+    <scheme>pumas_post_main</scheme>
+    <scheme>micro_pumas_ccpp_optics_limiter</scheme>
     <scheme>apply_heating_rate</scheme>
     <!--<scheme>qneg</scheme>-->
     <scheme>geopotential_temp</scheme>
-    <!-- Post PUMAS for REI - needs updated numice, numliq -->
-    <scheme>pumas_post_main</scheme>
   </group>
 </suite>