Progress on substep_compare.

Author: quantheory

build_micro.csh

@@ -46,7 +46,7 @@ set FILES="${DIR}/wv_sat_methods.F90 ${DIR}/micro_mg_utils.F90 ${DIR}/micro_mg2_
 #are separated by spaces - the code still compiles if commas are used instead,
 #but it crashes when you try to load it into python. 
 
-set PUB_FNS='wv_sat_methods_init wv_sat_svp_water wv_sat_svp_ice wv_sat_svp_to_qsat wv_sat_qsat_water wv_sat_qsat_ice micro_mg_init micro_mg_tend calc_precip_frac micro_mg_utils_init ice_deposition_sublimation kk2000_liq_autoconversion ice_autoconversion immersion_freezing contact_freezing snow_self_aggregation accrete_cloud_water_snow secondary_ice_production accrete_rain_snow heterogeneous_rain_freezing accrete_cloud_water_rain self_collection_rain accrete_cloud_ice_snow evaporate_sublimate_precip bergeron_process_snow sedimentation avg_diameter rain_fall_speeds'
+set PUB_FNS='wv_sat_methods_init wv_sat_svp_water wv_sat_svp_ice wv_sat_svp_to_qsat wv_sat_qsat_water wv_sat_qsat_ice micro_mg_init micro_mg_tend calc_precip_frac micro_mg_utils_init ice_deposition_sublimation kk2000_liq_autoconversion ice_autoconversion immersion_freezing contact_freezing snow_self_aggregation accrete_cloud_water_snow secondary_ice_production accrete_rain_snow heterogeneous_rain_freezing accrete_cloud_water_rain self_collection_rain accrete_cloud_ice_snow evaporate_rain sublimate_snow bergeron_process_snow sedimentation avg_diameter rain_fall_speeds'
 
 #CREATE .pyf FILE WITH INTERFACES TO PYTHON:
 #=================

micro_mg2_0.F90

@@ -347,7 +347,8 @@ end subroutine micro_mg_init
 !microphysics routine for each timestep goes here...
 
 subroutine micro_mg_tend ( &
-     mgncol,             nlev,               deltatin,           &
+     mgncol,             nlev,               evap_col_steps,     &
+     evap_steps,         col_steps,          deltatin,           &
      t,                            q,                            &
      qcn,                          qin,                          &
      ncn,                          nin,                          &
@@ -442,7 +443,8 @@ subroutine micro_mg_tend ( &
        accrete_cloud_water_rain, &
        self_collection_rain, &
        accrete_cloud_ice_snow, &
-       evaporate_sublimate_precip, &
+       evaporate_rain, &
+       sublimate_snow, &
        bergeron_process_snow, &
        sedimentation
 
@@ -452,6 +454,9 @@ subroutine micro_mg_tend ( &
   ! input arguments
   integer,  intent(in) :: mgncol         ! number of microphysics columns
   integer,  intent(in) :: nlev           ! number of layers
+  integer,  intent(in) :: evap_col_steps ! Evap/collection substeps per MG2 steps
+  integer,  intent(in) :: evap_steps     ! Rain evaporation substeps per evap/col step
+  integer,  intent(in) :: col_steps      ! Self-collection substeps per evap/col step
   real(r8), intent(in) :: deltatin       ! time step (s)
   real(r8), intent(in) :: t(mgncol,nlev)        ! input temperature (K)
   real(r8), intent(in) :: q(mgncol,nlev)        ! input h20 vapor mixing ratio (kg/kg)
@@ -842,13 +847,17 @@ subroutine micro_mg_tend ( &
   ! of sedimentation to MG2 timestep, respectively.
   real(r8) :: sed_deltat, sed_time, sed_step_ratio
 
+  real(r8) :: dumt(mgncol), dumq(mgncol), dumrho(mgncol), dumqvl(mgncol)
+  real(r8) :: dumlamr(mgncol), dumn0r(mgncol), dumarn(mgncol), dumr2(mgncol), dumnr2(mgncol)
+  real(r8) :: dumnragg(mgncol), dumnragg2(mgncol), dumpre(mgncol)
+
   ! Processes that can be disabled.
   logical :: do_sed_loc, do_inst_loc
 
   ! loop array variables
   ! "i" and "k" are column/level iterators for internal (MG) variables
-  ! "n" is used for other looping (currently just sedimentation)
-  integer i, k, n
+  ! "m" and "n" are used for other looping
+  integer i, k, m, n
 
   ! number of sub-steps for loops over "n" (for sedimentation)
   integer nstep
@@ -1534,8 +1543,6 @@ subroutine micro_mg_tend ( &
      call accrete_cloud_water_rain(microp_uniform, qric(:,k), qcic(:,k), &
           ncic(:,k), relvar(:,k), accre_enhan(:,k), pra(:,k), npra(:,k))
 
-     call self_collection_rain(rho(:,k), qric(:,k), nric(:,k), nragg(:,k))
-
      if (do_cldice) then
         call accrete_cloud_ice_snow(t(:,k), rho(:,k), asn(:,k), qiic(:,k), niic(:,k), &
              qsic(:,k), lams(:,k), n0s(:,k), prai(:,k), nprai(:,k))
@@ -1544,11 +1551,10 @@ subroutine micro_mg_tend ( &
         nprai(:,k) = 0._r8
      end if
 
-     call evaporate_sublimate_precip(t(:,k), rho(:,k), &
+     call sublimate_snow(t(:,k), rho(:,k), &
           dv(:,k), mu(:,k), sc(:,k), q(:,k), qvl(:,k), qvi(:,k), &
-          lcldm(:,k), precip_frac(:,k), arn(:,k), asn(:,k), qcic(:,k), qiic(:,k), &
-          qric(:,k), qsic(:,k), lamr(:,k), n0r(:,k), lams(:,k), n0s(:,k), &
-          pre(:,k), prds(:,k))
+          lcldm(:,k), precip_frac(:,k), asn(:,k), qcic(:,k), qiic(:,k), &
+          qsic(:,k), lams(:,k), n0s(:,k), prds(:,k))
 
      call bergeron_process_snow(t(:,k), rho(:,k), dv(:,k), mu(:,k), sc(:,k), &
           qvl(:,k), qvi(:,k), asn(:,k), qcic(:,k), qsic(:,k), lams(:,k), n0s(:,k), &
@@ -1649,6 +1655,62 @@ subroutine micro_mg_tend ( &
 
      end do
 
+     dumr(:,k) = qric(:,k)
+     dumnr(:,k) = nric(:,k)
+     dumt = t(:,k)
+     dumq = q(:,k)
+     nragg(:,k) = 0._r8
+     pre(:,k) = 0._r8
+     do m = 1, evap_col_steps
+        dumrho = p(:,k)/(r*dumt)
+
+        ! Need separate temporary nr for the two loops because they have to jointly
+        ! conserve rain number.
+        dumnr2 = dumnr(:,k)
+        dumnragg2 = 0._r8
+        do n = 1, col_steps
+           call self_collection_rain(dumrho, dumr(:,k), dumnr2, dumnragg)
+           dumnragg = max(dumnragg, -dumnr2*col_steps*evap_col_steps/deltat)
+           nragg(:,k) = nragg(:,k) + dumnragg / col_steps
+           dumnragg2 = dumnragg2 + dumnragg / col_steps
+           dumnr2 = max(dumnr2 + dumnragg*deltat/(col_steps*evap_col_steps), 0._r8)
+        end do
+
+        dumnr2 = dumnr(:,k)
+        do n = 1, evap_steps
+           dumrho = p(:,k)/(r*dumt)
+           dumarn = ar*(rhosu/dumrho)**0.54_r8
+           do i = 1, mgncol
+              call qsat_water(dumt(i), p(i,k), dum, dumqvl(i))
+           end do
+           call size_dist_param_basic(mg_rain_props, dumr(:,k), dumnr2, &
+                dumlamr, dumn0r)
+           call evaporate_rain(dumt, dumrho, p(:,k), dumq, dumqvl, &
+                lcldm(:,k), precip_frac(:,k), dumarn, qcic(:,k), qiic(:,k), &
+                dumr(:,k), dumlamr, dumn0r, dumpre)
+           dumpre = max(dumpre, -dumr(:,k)*evap_steps*evap_col_steps/deltat)
+           pre(:,k) = pre(:,k) + dumpre / evap_steps
+           dumnr2 = max(dumnr2 + (dumnr2/dumr(:,k))*dumpre*deltat/(evap_steps*evap_col_steps), 0._r8)
+           dumr(:,k) = max(dumr(:,k) + dumpre*deltat/(evap_steps*evap_col_steps), 0._r8)
+           ! Not that dumq and dumt are grid-mean, not in-precip quantities.
+           dumq = dumq - precip_frac(:,k)*dumpre*deltat/(evap_steps*evap_col_steps)
+           dumt = dumt + xxlv*precip_frac(:,k)*dumpre*deltat/(evap_steps*evap_col_steps)
+        end do
+
+        ! Joint conservation of rain number; adjust self-collection only for
+        ! consistency with existing checks, and since evaporation should not be
+        ! limited by number directly anyway.
+        where (-dumnragg2*deltat/evap_col_steps > dumnr2)
+           dumnragg = dumnragg2
+           dumnragg2 = -dumnr2 * evap_col_steps / deltat
+           nragg(:,k) = nragg(:,k) + dumnragg2 - dumnragg
+           ! Note: might be better to use size_dist_param_basic to set min here.
+           dumnr(:,k) = 0._r8
+        elsewhere
+           dumnr(:,k) = dumnr2 + dumnragg2*deltat/evap_col_steps
+        end where
+     end do
+
      do i=1,mgncol
 
         !===================================================================

micro_mg_utils.F90

@@ -66,7 +66,8 @@ module micro_mg_utils
      accrete_cloud_water_rain, &
      self_collection_rain, &
      accrete_cloud_ice_snow, &
-     evaporate_sublimate_precip, &
+     evaporate_rain, &
+     sublimate_snow, &
      bergeron_process_snow, &
      sedimentation
 
@@ -1242,46 +1243,39 @@ end subroutine accrete_cloud_ice_snow
 ! in-cloud condensation/deposition of rain and snow is neglected
 ! except for transfer of cloud water to snow through bergeron process
 
-elemental subroutine evaporate_sublimate_precip(t, rho, dv, mu, sc, q, qvl, qvi, &
-     lcldm, precip_frac, arn, asn, qcic, qiic, qric, qsic, lamr, n0r, lams, n0s, &
-     pre, prds)
+elemental subroutine evaporate_rain(t, rho, p, q, qvl, &
+     lcldm, precip_frac, arn, qcic, qiic, qric, lamr, n0r, pre)
 
   real(r8), intent(in) :: t    ! temperature
   real(r8), intent(in) :: rho  ! air density
-  real(r8), intent(in) :: dv   ! water vapor diffusivity
-  real(r8), intent(in) :: mu   ! viscosity
-  real(r8), intent(in) :: sc   ! schmidt number
+  real(r8), intent(in) :: p    ! pressure
   real(r8), intent(in) :: q    ! humidity
   real(r8), intent(in) :: qvl  ! saturation humidity (water)
-  real(r8), intent(in) :: qvi  ! saturation humidity (ice)
   real(r8), intent(in) :: lcldm  ! liquid cloud fraction
   real(r8), intent(in) :: precip_frac ! precipitation fraction (maximum overlap)
 
   ! fallspeed parameters
   real(r8), intent(in) :: arn  ! rain
-  real(r8), intent(in) :: asn  ! snow
 
   ! In-cloud MMRs
   real(r8), intent(in) :: qcic ! cloud liquid
   real(r8), intent(in) :: qiic ! cloud ice
   real(r8), intent(in) :: qric ! rain
-  real(r8), intent(in) :: qsic ! snow
 
   ! Size parameters
   ! rain
   real(r8), intent(in) :: lamr
   real(r8), intent(in) :: n0r
-  ! snow
-  real(r8), intent(in) :: lams
-  real(r8), intent(in) :: n0s
 
   ! Output tendencies
   real(r8), intent(out) :: pre
-  real(r8), intent(out) :: prds
-
+ 
   real(r8) :: qclr   ! water vapor mixing ratio in clear air
   real(r8) :: ab     ! correction to account for latent heat
   real(r8) :: eps    ! 1/ sat relaxation timescale
+  real(r8) :: dv   ! water vapor diffusivity
+  real(r8) :: mu   ! viscosity
+  real(r8) :: sc   ! schmidt number
 
   real(r8) :: dum
 
@@ -1304,6 +1298,10 @@ elemental subroutine evaporate_sublimate_precip(t, rho, dv, mu, sc, q, qvl, qvi,
      ! evaporation of rain
      if (qric >= qsmall) then
 
+        dv = 8.794E-5_r8 * t**1.81_r8 / p
+        mu = 1.496E-6_r8 * t**1.5_r8 / (t + 120._r8)
+        sc = mu/(rho*dv)
+
         ab = calc_ab(t, qvl, xxlv)
         eps = 2._r8*pi*n0r*rho*Dv* &
              (f1r/(lamr*lamr)+ &
@@ -1321,6 +1319,70 @@ elemental subroutine evaporate_sublimate_precip(t, rho, dv, mu, sc, q, qvl, qvi,
         pre = 0._r8
      end if
 
+  else
+     pre = 0._r8
+  end if
+
+end subroutine evaporate_rain
+
+! calculate evaporation/sublimation of rain and snow
+!===================================================================
+! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
+! in-cloud condensation/deposition of rain and snow is neglected
+! except for transfer of cloud water to snow through bergeron process
+
+elemental subroutine sublimate_snow(t, rho, dv, mu, sc, q, qvl, qvi, &
+     lcldm, precip_frac, asn, qcic, qiic, qsic, lams, n0s, prds)
+
+  real(r8), intent(in) :: t    ! temperature
+  real(r8), intent(in) :: rho  ! air density
+  real(r8), intent(in) :: dv   ! water vapor diffusivity
+  real(r8), intent(in) :: mu   ! viscosity
+  real(r8), intent(in) :: sc   ! schmidt number
+  real(r8), intent(in) :: q    ! humidity
+  real(r8), intent(in) :: qvl  ! saturation humidity (water)
+  real(r8), intent(in) :: qvi  ! saturation humidity (ice)
+  real(r8), intent(in) :: lcldm  ! liquid cloud fraction
+  real(r8), intent(in) :: precip_frac ! precipitation fraction (maximum overlap)
+
+  ! fallspeed parameters
+  real(r8), intent(in) :: asn  ! snow
+
+  ! In-cloud MMRs
+  real(r8), intent(in) :: qcic ! cloud liquid
+  real(r8), intent(in) :: qiic ! cloud ice
+  real(r8), intent(in) :: qsic ! snow
+
+  ! Size parameters
+  ! snow
+  real(r8), intent(in) :: lams
+  real(r8), intent(in) :: n0s
+
+  ! Output tendencies
+  real(r8), intent(out) :: prds
+
+  real(r8) :: qclr   ! water vapor mixing ratio in clear air
+  real(r8) :: ab     ! correction to account for latent heat
+  real(r8) :: eps    ! 1/ sat relaxation timescale
+
+  real(r8) :: dum
+
+  ! set temporary cloud fraction to zero if cloud water + ice is very small
+  ! this will ensure that evaporation/sublimation of precip occurs over
+  ! entire grid cell, since min cloud fraction is specified otherwise
+  if (qcic+qiic < 1.e-6_r8) then
+     dum = 0._r8
+  else
+     dum = lcldm
+  end if
+
+  ! only calculate if there is some precip fraction > cloud fraction
+
+  if (precip_frac > dum) then
+
+     ! calculate q for out-of-cloud region
+     qclr=(q-dum*qvl)/(1._r8-dum)
+
      ! sublimation of snow
      if (qsic >= qsmall) then
         ab = calc_ab(t, qvi, xxls)
@@ -1340,10 +1402,9 @@ elemental subroutine evaporate_sublimate_precip(t, rho, dv, mu, sc, q, qvl, qvi,
 
   else
      prds = 0._r8
-     pre = 0._r8
   end if
 
-end subroutine evaporate_sublimate_precip
+end subroutine sublimate_snow
 
 ! bergeron process - evaporation of droplets and deposition onto snow
 !===================================================================