Changeset 4532 for palm/trunk/SOURCE/bulk_cloud_model_mod.f90

Timestamp:

May 14, 2020 1:41:35 PM (5 years ago)

Author:

schwenkel

Message:

Calculate mean droplet radius assuming gamma distibution for condensation

File:

• palm/trunk/SOURCE/bulk_cloud_model_mod.f90 (modified) (9 diffs)

palm/trunk/SOURCE/bulk_cloud_model_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Calculate mean droplet radius assuming gamma distibution for condensation
+! 
+! 4521 2020-05-06 11:39:49Z schwenkel
 ! Rename variable
 ! 
@@ -4317 +4320 @@
                    CALL supersaturation_ice ( i, j, k )
                    nucle = 0.0_wp
-                   !IF ( zu(k) >= 1500.0_wp ) CYCLE
                    IF ( sat_ice >= 0.05_wp  .OR.  ql(k,j,i) >= 0.001E-3_wp  ) THEN
 !
 !--                   Calculate ice nucleation
                       nucle = MAX( ( in_init - ni(k,j,i) ) / dt_micro, 0.0_wp )
-                      !qi(k,j,i) = qi(k,j,i) + nucle * dt_micro * ximin
                       ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), in_init)
                    ENDIF
@@ -4386 +4387 @@
              CALL supersaturation_ice ( i, j, k )
              nucle = 0.0_wp
-             !IF ( zu(k) >= 1500.0_wp ) CYCLE
              IF ( sat_ice >= 0.05_wp  .OR.  ql(k,j,i) >= 0.001E-3_wp ) THEN
 !
 !--             Calculate ice nucleation
                 nucle = MAX( ( in_init - ni(k,j,i) ) / dt_micro, 0.0_wp )
-                !qi(k,j,i) = qi(k,j,i) + nucle * dt_micro * ximin
                 ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), in_init)
              ENDIF
@@ -4431 +4430 @@
        INTEGER(iwp) ::  k                 !< loop index
 
+       REAL(wp)     ::  alpha_rc = 1.0_wp !< Tuning parameter see (Seifert and Stevens, 2010)
        REAL(wp)     ::  cond              !< condensation rate
        REAL(wp)     ::  cond_max          !< maximum condensation rate
-       REAL(wp)     ::  dc                !< weight avageraed diameter
        REAL(wp)     ::  evap              !< evaporation rate
        REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
        REAL(wp)     ::  g_fac             !< factor 1 / Fk + Fd
-       REAL(wp)     ::  nc_0              !< integral diameter
+       REAL(wp)     ::  nu = 1.0_wp       !< Shape parameter of gernerlized gama distribution
+       REAL(wp)     ::  rc                !< mean cloud droplets radius assuming gamma distribution
        REAL(wp)     ::  temp              !< actual temperature
-       REAL(wp)     ::  xc                !< mean mass droplets
 
        CALL cpu_log( log_point_s(66), 'condensation', 'start' )
@@ -4467 +4466 @@
 !
 !--             Mean weight of cloud drops
-                IF ( nc(k,j,i) <= 0.0_wp) CYCLE
-                xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
-!
-!--             Weight averaged diameter of cloud drops:
-                dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
-!
-!--             Integral diameter of cloud drops
-                nc_0 = nc(k,j,i) * dc
+                IF ( nc(k,j,i) <= 0.0_wp ) CYCLE
+!
+!--             Calculating mean radius of cloud droplets assuming gamma distribution with shape
+!--             parameter nu=1 (Seifert and Beheng, 2006). Tuning factor alpha_rc (intorduced
+!--             in Seifert and Stevens, 2010 ) is switched off. Minimum radius is set to 1µm
+!--             following (Seifert and Beheng, 2006, Kogan and Khairoutdinov, 2000,
+!--             Seifert and Stevens, 2010)
+                rc = MAX( alpha_rc  * gamma(nu + 1.33_wp) / gamma(nu + 1.0_wp) *                   &
+                        ( 3.0_wp * qc(k,j,i) /                                                     &
+                                   ( 4.0_wp * pi * rho_l  * ( nu + 2.0_wp ) * nc(k,j,i) )          &
+                        )**0.33_wp, 1.0E-6_wp )
 !
 !--             Condensation needs only to be calculated in supersaturated regions
                 IF ( sat > 0.0_wp )  THEN
 !
-!--                Condensation rate of cloud water content
-!--                after KK scheme.
-!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
-                   cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
+!--                Condensation rate of cloud water content after KK scheme.
+!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128, Morrison and Grabowski,
+!--                2007, and Seifert and Stevens, 2010)
+                   cond      = 4.0_wp * pi * nc(k,j,i) * g_fac * sat * rc / hyrho(k)
                    IF ( microphysics_seifert ) THEN
                       cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
@@ -4492 +4494 @@
                    qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
                 ELSEIF ( sat < 0.0_wp ) THEN
-                   evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
+                   evap      = 4.0_wp * pi * nc(k,j,i) * g_fac * sat * rc / hyrho(k)
                    evap      = MAX( evap, -qc(k,j,i) / dt_micro )
 
@@ -4522 +4524 @@
        INTEGER(iwp) ::  k                 !< loop index
 
+       REAL(wp)     ::  alpha_rc = 1.0_wp !< Tuning parameter see (Seifert and Stevens, 2010)
        REAL(wp)     ::  cond              !< condensation rate
        REAL(wp)     ::  cond_max          !< maximum condensation rate
-       REAL(wp)     ::  dc                !< weight avageraed diameter
        REAL(wp)     ::  evap              !< evaporation rate
        REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
        REAL(wp)     ::  g_fac             !< factor 1 / Fk + Fd
-       REAL(wp)     ::  nc_0              !< integral diameter
+       REAL(wp)     ::  nu = 1.0_wp       !< Shape parameter of gernerlized gama distribution
+       REAL(wp)     ::  rc                !< mean cloud droplets radius assuming gamma distribution
        REAL(wp)     ::  temp              !< actual temperature
-       REAL(wp)     ::  xc                !< mean mass droplets
 
        DO  k = nzb+1, nzt
@@ -4555 +4557 @@
 !--       Mean weight of cloud drops
           IF ( nc(k,j,i) <= 0.0_wp) CYCLE
-          xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
-!
-!--       Weight averaged diameter of cloud drops:
-          dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
-!
-!--       Integral diameter of cloud drops
-          nc_0 = nc(k,j,i) * dc
+!
+!--       Calculating mean radius of cloud droplets assuming gamma distribution with shape
+!--       parameter nu=1 (Seifert and Beheng, 2006). Tuning factor alpha_rc (intorduced
+!--       in Seifert and Stevens, 2010 ) is switched off. Minimum radius is set to 1µm following
+!--       (Seifert and Beheng, 2006, Kogan and Khairoutdinov, 2000, Seifert and Stevens, 2010)
+          rc = MAX( alpha_rc  * gamma(nu + 1.33_wp) / gamma(nu + 1.0_wp) *                         &
+               ( 3.0_wp * qc(k,j,i) / ( 4.0_wp * pi * rho_l  * ( nu + 2.0_wp ) * nc(k,j,i) )       &
+               )**0.33_wp, 1.0E-6_wp )
 !
 !--       Condensation needs only to be calculated in supersaturated regions
           IF ( sat > 0.0_wp )  THEN
 !
-!--          Condensation rate of cloud water content
-!--          after KK scheme.
-!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
-             cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
+!--          Condensation rate of cloud water content after KK scheme.
+!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128, Morrison and Grabowski, 2007
+!--          and Seifert and Stevens, 2010)
+             cond      = 4.0_wp * pi * nc(k,j,i) * g_fac * sat * rc / hyrho(k)
              IF ( microphysics_seifert ) THEN
                 cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
@@ -4579 +4582 @@
              qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
           ELSEIF ( sat < 0.0_wp ) THEN
-             evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
+             evap      = 4.0_wp * pi * nc(k,j,i) * g_fac * sat * rc / hyrho(k)
              evap      = MAX( evap, -qc(k,j,i) / dt_micro )