make some Thompson tuning parameters externally accessible

Author: grantfirl

physics/MP/Thompson/module_mp_thompson.F90

@@ -92,10 +92,17 @@ module module_mp_thompson
 !.. droplet concentration and nu_c is also variable depending on local
 !.. droplet number concentration.
    !real(wp), parameter :: Nt_c = 100.e6
-   real(wp), parameter :: Nt_c_o = 50.e6
-   real(wp), parameter :: Nt_c_l = 150.e6
    real(wp), parameter, private :: Nt_c_max = 1999.e6
 
+   ! Tuning parameters
+   real(wp)            :: Nt_c_l = 150.e6   ! Cloud number concentration over land (set in thompson_init)
+   real(wp)            :: Nt_c_o = 50.e6    ! Cloud number concentration over ocean (set in thompson_init)
+   real(wp)            :: av_i
+   real(wp)            :: xnc_max = 1000.e3
+   real(wp)            :: ssati_min = 0.15
+   real(wp)            :: Nt_i_max = 4999.e3_dp
+   real(wp)            :: rr_min = 1000.0
+   
 !..Declaration of constants for assumed CCN/IN aerosols when none in
 !.. the input data.  Look inside the init routine for modifications
 !.. due to surface land-sea points or vegetation characteristics.
@@ -456,7 +463,7 @@ subroutine thompson_init(is_aerosol_aware_in,       &
          logical:: micro_init
          real(wp) :: stime, etime
          logical, parameter :: precomputed_tables = .FALSE.
-
+         
 ! Set module derived constants
          am_r = PI*rho_w/6.0
          am_g = PI*rho_g/6.0
@@ -472,6 +479,9 @@ subroutine thompson_init(is_aerosol_aware_in,       &
 
          lsub = lvap0 + lfus
          olfus = 1./lfus
+         
+        ! Transition value of coefficient matching at crossover from cloud ice to snow
+         av_i = av_s * D0s ** (bv_s - bv_i)
 
 ! Set module variable is_aerosol_aware/merra2_aerosol_aware
          is_aerosol_aware = is_aerosol_aware_in
@@ -2017,7 +2027,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d,     &
       real(wp) :: Ef_ra, Ef_sa, Ef_ga
       real(wp) :: dtsave, odts, odt, odzq, hgt_agl, SR
       real(wp) :: xslw1, ygra1, zans1, eva_factor
-      real(wp) av_i
       integer :: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq
       integer, dimension(5) :: ksed1
       integer :: nir, nis, nig, nii, nic, niin
@@ -2042,8 +2051,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d,     &
       odt = 1./dt
       odts = 1./dtsave
       iexfrq = 1
-! Transition value of coefficient matching at crossover from cloud ice to snow
-      av_i = av_s * D0s ** (bv_s - bv_i)
 
 !+---+-----------------------------------------------------------------+
 !> - Initialize Source/sink terms.  First 2 chars: "pr" represents source/sink of
@@ -2269,15 +2276,15 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d,     &
             ni(k) = max(R2, ni1d(k)*rho(k))
             if (ni(k).le. R2) then
                lami = cie(2)/5.E-6
-               ni(k) = min(4999.e3_dp, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
+               ni(k) = min(Nt_i_max, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
             endif
             L_qi(k) = .true.
             lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
             ilami = 1./lami
             xDi = (bm_i + mu_i + 1.) * ilami
             if (xDi.lt. 5.E-6) then
                lami = cie(2)/5.E-6
-               ni(k) = min(4999.e3_dp, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
+               ni(k) = min(Nt_i_max, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
             elseif (xDi.gt. 300.E-6) then
                lami = cie(2)/300.E-6
                ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i
@@ -2933,13 +2940,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d,     &
 
 !>  - Deposition nucleation of dust/mineral from DeMott et al (2010)
 !! we may need to relax the temperature and ssati constraints.
-               if ( (ssati(k).ge. 0.15) .or. (ssatw(k).gt. eps &
+               if ( (ssati(k).ge. ssati_min) .or. (ssatw(k).gt. eps &
                                     .and. temp(k).lt.253.15) ) then
                   if (dustyIce .AND. (is_aerosol_aware .or. merra2_aerosol_aware)) then
                      xnc = iceDeMott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k))
                      xnc = xnc*(1.0 + 50.*rand3)
                   else
-                     xnc = min(1000.E3, TNO*EXP(ATO*(T_0-temp(k))))
+                     xnc = min(xnc_max, TNO*EXP(ATO*(T_0-temp(k))))
                   endif
                   xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave
                   pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
@@ -2949,7 +2956,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d,     &
 
 !>  - Freezing of aqueous aerosols based on Koop et al (2001, Nature)
                xni = smo0(k)+ni(k) + (pni_rfz(k)+pni_wfz(k)+pni_inu(k))*dtsave
-               if ((is_aerosol_aware .or. merra2_aerosol_aware) .AND. homogIce .AND. (xni.le.4999.E3)    & 
+               if ((is_aerosol_aware .or. merra2_aerosol_aware) .AND. homogIce .AND. (xni.le.Nt_i_max)    & 
                               .AND.(temp(k).lt.238).AND.(ssati(k).ge.0.4) ) then
                   xnc = iceKoop(temp(k),qv(k),qvs(k),nwfa(k), dtsave)
                   pni_iha(k) = xnc*odts
@@ -3282,7 +3289,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d,     &
             xDi = (bm_i + mu_i + 1.) * ilami
             if (xDi.lt. 5.E-6) then
                lami = cie(2)/5.E-6
-               xni = min(4999.e3_dp, cig(1)*oig2*xri/am_i*lami**bm_i)
+               xni = min(Nt_i_max, cig(1)*oig2*xri/am_i*lami**bm_i)
                niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
             elseif (xDi.gt. 300.E-6) then 
                lami = cie(2)/300.E-6
@@ -3293,8 +3300,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d,     &
             niten(k) = -ni1d(k)*odts
          endif
          xni=max(0.,(ni1d(k) + niten(k)*dtsave)*rho(k))
-         if (xni.gt.4999.E3) &
-                niten(k) = (4999.E3-ni1d(k)*rho(k))*odts*orho
+         if (xni.gt.Nt_i_max) &
+                niten(k) = (Nt_i_max-ni1d(k)*rho(k))*odts*orho
 
 !>  - Rain tendency
          qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) &
@@ -3977,7 +3984,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d,     &
                   pfll1(k) = pfll1(k) + sed_r(k)*DT*onstep(1)
                enddo
 
-               if (rr(kts).gt.R1*1000.) then
+               if (rr(kts).gt.R1*rr_min) then
                   pptrain = pptrain + sed_r(kts)*DT*onstep(1)
                endif 
             enddo
@@ -4072,7 +4079,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d,     &
                pfil1(k) = pfil1(k) + sed_i(k)*DT*onstep(2)
             enddo
 
-            if (ri(kts).gt.R1*1000.) then
+            if (ri(kts).gt.R1*rr_min) then
                pptice = pptice + sed_i(kts)*DT*onstep(2)
             endif 
          enddo
@@ -4102,7 +4109,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d,     &
                pfil1(k) = pfil1(k) + sed_s(k)*DT*onstep(3)
             enddo
 
-            if (rs(kts).gt.R1*1000.) then
+            if (rs(kts).gt.R1*rr_min) then
                pptsnow = pptsnow + sed_s(kts)*DT*onstep(3)
             endif 
          enddo
@@ -4133,7 +4140,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d,     &
                   pfil1(k) = pfil1(k) + sed_g(k)*DT*onstep(4)
                enddo
 
-               if (rg(kts).gt.R1*1000.) then
+               if (rg(kts).gt.R1*rr_min) then
                   pptgraul = pptgraul + sed_g(kts)*DT*onstep(4)
                endif
             enddo
@@ -4261,7 +4268,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d,     &
                lami = cie(2)/300.E-6
             endif
             ni1d(k) = min(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i,           &
-                           4999.e3_dp/rho(k))
+                           Nt_i_max/rho(k))
          endif
          qr1d(k) = qr1d(k) + qrten(k)*DT
          nr1d(k) = max(R2/rho(k), nr1d(k) + nrten(k)*DT)

physics/MP/Thompson/mp_thompson.F90

@@ -11,7 +11,7 @@ module mp_thompson
       use machine, only : kind_phys
 
       use module_mp_thompson, only : thompson_init, mp_gt_driver, thompson_finalize, calc_effectRad
-      use module_mp_thompson, only : naIN0, naIN1, naCCN0, naCCN1, eps, Nt_c_l, Nt_c_o
+      use module_mp_thompson, only : naIN0, naIN1, naCCN0, naCCN1, eps
       use module_mp_thompson, only : re_qc_min, re_qc_max, re_qi_min, re_qi_max, re_qs_min, re_qs_max
 
       use module_mp_thompson_make_number_concentrations, only: make_IceNumber, make_DropletNumber, make_RainNumber
@@ -34,6 +34,9 @@ subroutine mp_thompson_init(ncol, nlev, con_pi, con_t0c, con_rv,     &
                                   con_cp, con_rgas, con_boltz, con_amd,    &
                                   con_amw, con_avgd, con_hvap, con_hfus,   &
                                   con_g, con_rd, con_eps,                  &
+                                  con_Nt_c_l, con_Nt_c_o, con_av_i,        &
+                                  con_xnc_max, con_ssati_min, con_Nt_i_max,&
+                                  con_rr_min,                              &
                                   restart, imp_physics,                    &
                                   imp_physics_thompson, convert_dry_rho,   &
                                   spechum, qc, qr, qi, qs, qg, ni, nr,     &
@@ -45,6 +48,7 @@ subroutine mp_thompson_init(ncol, nlev, con_pi, con_t0c, con_rv,     &
                                   is_initialized, errmsg, errflg)
          use module_mp_thompson, only : PI, T_0, Rv, R, RoverRv, Cp
          use module_mp_thompson, only : R_uni, k_b, M_w, M_a, N_avo, lvap0, lfus
+         use module_mp_thompson, only : nt_c_l, nt_c_o, av_i, xnc_max, ssati_min, Nt_i_max, rr_min
 
          implicit none
 
@@ -54,6 +58,8 @@ subroutine mp_thompson_init(ncol, nlev, con_pi, con_t0c, con_rv,     &
          real(kind_phys),           intent(in   ) :: con_pi, con_t0c, con_rv, con_cp, con_rgas, &
                                                      con_boltz, con_amd, con_amw, con_avgd,     &
                                                      con_hvap, con_hfus, con_g, con_rd, con_eps
+         real(kind_phys), optional, intent(in   ) :: con_Nt_c_l, con_Nt_c_o, con_av_i, con_xnc_max, &
+                                                     con_ssati_min, con_Nt_i_max, con_rr_min
          logical,                   intent(in   ) :: restart
          logical,                   intent(inout) :: is_initialized
          integer,                   intent(in   ) :: imp_physics
@@ -126,6 +132,13 @@ subroutine mp_thompson_init(ncol, nlev, con_pi, con_t0c, con_rv,     &
          lvap0 = con_hvap
          lfus = con_hfus
 
+         if (present(con_nt_c_l)) nt_c_l = con_nt_c_l
+         if (present(con_nt_c_o)) nt_c_o = con_nt_c_o
+         if (present(con_xnc_max)) xnc_max = con_xnc_max
+         if (present(con_ssati_min)) ssati_min = con_ssati_min
+         if (present(con_Nt_i_max)) Nt_i_max = con_Nt_i_max
+         if (present(con_rr_min)) rr_min = con_rr_min
+         
          ! Consistency checks
          if (imp_physics/=imp_physics_thompson) then
             write(errmsg,'(*(a))') "Logic error: namelist choice of microphysics is different from Thompson MP"
@@ -153,7 +166,10 @@ subroutine mp_thompson_init(ncol, nlev, con_pi, con_t0c, con_rv,     &
                             mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, &
                             threads=threads, errmsg=errmsg, errflg=errflg)
          if (errflg /= 0) return
-
+         
+         ! override the value of av_i with the constant value
+         if (present(con_av_i)) av_i = con_av_i
+         
          ! For restart runs, the init is done here
          if (restart) then
            is_initialized = .true.

physics/MP/Thompson/mp_thompson.meta

@@ -135,6 +135,69 @@
   type = real
   kind = kind_phys
   intent = in
+[con_Nt_c_l]
+  standard_name = prescribed_number_concentration_of_cloud_liquid_water_particles_in_air_over_land
+  long_name = volumetric number concentration of cloud liquid water droplets in air over land points
+  units = m-3
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = True
+[con_Nt_c_o]
+  standard_name = prescribed_number_concentration_of_cloud_liquid_water_particles_in_air_over_ocean
+  long_name = volumetric number concentration of cloud liquid water droplets in air over ocean points
+  units = m-3
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = True
+[con_av_i]
+  standard_name = cloud_ice_to_snow_tuning_parameter
+  long_name = transition value of coefficient matching at crossover from cloud ice to snow
+  units = 1
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = True
+[con_xncmax]
+  standard_name = maximum_mass_number_concentration_of_cloud_liquid_water_particles_in_air_used_in_deposition_nucleation
+  long_name = maximum mass number concentration of cloud liquid water particles in air
+  units = kg-1
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = True
+[con_ssati_min]
+  standard_name = minimum_threshold_supersaturation_over_ice_for_deposition_nucleation
+  long_name = minimum supersaturation over ice threshold for deposition nucleation
+  units = fraction
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = True
+[con_Nt_i_max]
+  standard_name = maximum_threshold_number_concentration_of_cloud_ice_water_crystals_in_air
+  long_name = maximum threshold number concentration of cloud ice water crystals in air
+  units = m-3
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = True
+[con_rr_min]
+  standard_name = multiplicative_tuning_parameter_for_microphysical_sedimentation_minimum_threshold
+  long_name = multiplicative tuning parameter for microphysical sedimentation minimum threshold
+  units = 1
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = True
 [restart]
   standard_name = flag_for_restart
   long_name = flag for restart (warmstart) or coldstart