Changeset 2312 for palm/trunk/SOURCE/microphysics_mod.f90
- Timestamp:
- Jul 14, 2017 8:26:51 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/microphysics_mod.f90 (modified) (43 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/microphysics_mod.f90
r2292 r2312 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 23 ! 22 ! 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 28 ! includes two more prognostic equations for cloud drop concentration (nc) 29 ! and cloud water content (qc). 30 ! - The process of activation is parameterized with a simple Twomey 31 ! activion scheme or with considering solution and curvature 27 ! s1 changed to log_sigma 28 ! 29 ! 2292 2017-06-20 09:51:42Z schwenkel 30 ! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 31 ! includes two more prognostic equations for cloud drop concentration (nc) 32 ! and cloud water content (qc). 33 ! - The process of activation is parameterized with a simple Twomey 34 ! activion scheme or with considering solution and curvature 32 35 ! effects (Khvorostyanov and Curry ,2006). 33 36 ! - The saturation adjustment scheme is replaced by the parameterization 34 37 ! of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128). 35 38 ! - All other microphysical processes of Seifert and Beheng are used. 36 ! Additionally, in those processes the reduction of cloud number concentration 37 ! is considered. 38 ! 39 ! Additionally, in those processes the reduction of cloud number concentration 40 ! is considered. 41 ! 39 42 ! 2233 2017-05-30 18:08:54Z suehring 40 43 ! 41 44 ! 2232 2017-05-30 17:47:52Z suehring 42 45 ! Adjustments to new topography and surface concept 43 ! 46 ! 44 47 ! 2155 2017-02-21 09:57:40Z hoffmann 45 48 ! Bugfix in the calculation of microphysical quantities on ghost points. … … 207 210 limiter_sedimentation, na_init, nc_const, sigma_gc, & 208 211 ventilation_effect 209 212 210 213 211 214 INTERFACE microphysics_control … … 463 466 IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN 464 467 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * & 465 MERGE( 1.0_wp, 0.0_wp, & 468 MERGE( 1.0_wp, 0.0_wp, & 466 469 BTEST( wall_flags_0(k,j,i), 0 ) ) 467 470 ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN 468 471 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * & 469 MERGE( 1.0_wp, 0.0_wp, & 472 MERGE( 1.0_wp, 0.0_wp, & 470 473 BTEST( wall_flags_0(k,j,i), 0 ) ) 471 474 ENDIF … … 478 481 IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN 479 482 nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin * & 480 MERGE( 1.0_wp, 0.0_wp, & 483 MERGE( 1.0_wp, 0.0_wp, & 481 484 BTEST( wall_flags_0(k,j,i), 0 ) ) 482 485 ENDIF … … 521 524 USE particle_attributes, & 522 525 ONLY: molecular_weight_of_solute, molecular_weight_of_water, rho_s, & 523 s1, s2, s3, vanthoff526 log_sigma, vanthoff 524 527 525 528 IMPLICIT NONE … … 578 581 !-- Prescribe parameters for activation 579 582 !-- (see: Bott + Trautmann, 2002, Atm. Res., 64) 580 k_act = 0.7_wp 583 k_act = 0.7_wp 581 584 582 585 IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN 583 586 ! 584 !-- Compute the number of activated Aerosols 587 !-- Compute the number of activated Aerosols 585 588 !-- (see: Twomey, 1959, Pure and applied Geophysics, 43) 586 589 n_act = na_init * sat**k_act … … 592 595 ! 593 596 !-- Compute activation rate after Khairoutdinov and Kogan 594 !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128) 595 sat_max = 1.0_wp / 100.0_wp 596 activ = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN & 597 ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / & 597 !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128) 598 sat_max = 1.0_wp / 100.0_wp 599 activ = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN & 600 ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / & 598 601 dt_micro 599 602 ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk ) THEN 600 603 ! 601 !-- Curvature effect (afactor) with surface tension 604 !-- Curvature effect (afactor) with surface tension 602 605 !-- parameterization by Straka (2009) 603 606 sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp ) … … 609 612 610 613 ! 611 !-- Prescribe power index that describes the soluble fraction 614 !-- Prescribe power index that describes the soluble fraction 612 615 !-- of an aerosol particle (beta) and mean geometric radius of 613 616 !-- dry aerosol spectrum … … 616 619 rd0 = 0.05E-6_wp 617 620 ! 618 !-- Calculate mean geometric supersaturation (s_0) with 621 !-- Calculate mean geometric supersaturation (s_0) with 619 622 !-- parameterization by Khvorostyanov and Curry (2006) 620 623 s_0 = rd0 **(- ( 1.0_wp + beta_act ) ) * & … … 622 625 623 626 ! 624 !-- Calculate number of activated CCN as a function of 627 !-- Calculate number of activated CCN as a function of 625 628 !-- supersaturation and taking Koehler theory into account 626 !-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111) 627 n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( & 628 LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(s1) ) ) )629 !-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111) 630 n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( & 631 LOG( s_0 / sat ) / ( SQRT(2.0_wp) * log_sigma(1) ) ) ) 629 632 activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp ) 630 633 ENDIF … … 644 647 ! Description: 645 648 ! ------------ 646 !> Calculate condensation rate for cloud water content (after Khairoutdinov and 649 !> Calculate condensation rate for cloud water content (after Khairoutdinov and 647 650 !> Kogan, 2000). 648 651 !------------------------------------------------------------------------------! … … 724 727 !-- Mean weight of cloud drops 725 728 IF ( nc(k,j,i) <= 0.0_wp) CYCLE 726 xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin) 729 xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin) 727 730 ! 728 731 !-- Weight averaged diameter of cloud drops: … … 730 733 ! 731 734 !-- Integral diameter of cloud drops 732 nc_0 = nc(k,j,i) * dc 735 nc_0 = nc(k,j,i) * dc 733 736 ! 734 737 !-- Condensation needs only to be calculated in supersaturated regions … … 741 744 cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i) 742 745 cond = MIN( cond, cond_max / dt_micro ) 743 746 744 747 qc(k,j,i) = qc(k,j,i) + cond * dt_micro 745 748 ELSEIF ( sat < 0.0_wp ) THEN … … 891 894 * flag 892 895 IF ( microphysics_morrison ) THEN 893 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp * & 896 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp * & 894 897 autocon / x0 * hyrho(k) * dt_micro * flag ) 895 898 ENDIF … … 1033 1036 ! 1034 1037 !-- Mean weight of cloud drops 1035 xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin) 1038 xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin) 1036 1039 ! 1037 1040 !-- Parameterized turbulence effects on autoconversion (Seifert, … … 1321 1324 1322 1325 USE control_parameters, & 1323 ONLY: call_microphysics_at_all_substeps, & 1326 ONLY: call_microphysics_at_all_substeps, & 1324 1327 intermediate_timestep_count, microphysics_morrison 1325 1328 … … 1366 1369 IF ( microphysics_morrison ) THEN 1367 1370 IF ( qc(k,j,i) > eps_sb .AND. nc(k,j,i) > eps_mr ) THEN 1368 sed_nc(k) = sed_qc_const * & 1371 sed_nc(k) = sed_qc_const * & 1369 1372 ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) * & 1370 1373 ( nc(k,j,i) )**( 1.0_wp / 3.0_wp ) … … 1449 1452 USE surface_mod, & 1450 1453 ONLY : bc_h 1451 1454 1452 1455 IMPLICIT NONE 1453 1456 … … 1455 1458 INTEGER(iwp) :: j !< running index y direction 1456 1459 INTEGER(iwp) :: k !< running index z direction 1457 INTEGER(iwp) :: k_run !< 1460 INTEGER(iwp) :: k_run !< 1458 1461 INTEGER(iwp) :: l !< running index of surface type 1459 1462 INTEGER(iwp) :: m !< running index surface elements … … 1528 1531 ENDDO 1529 1532 ! 1530 !-- Adjust boundary values using surface data type. 1533 !-- Adjust boundary values using surface data type. 1531 1534 !-- Upward-facing 1532 surf_s = bc_h(0)%start_index(j,i) 1535 surf_s = bc_h(0)%start_index(j,i) 1533 1536 surf_e = bc_h(0)%end_index(j,i) 1534 1537 DO m = surf_s, surf_e … … 1539 1542 ! 1540 1543 !-- Downward-facing 1541 surf_s = bc_h(1)%start_index(j,i) 1544 surf_s = bc_h(1)%start_index(j,i) 1542 1545 surf_e = bc_h(1)%end_index(j,i) 1543 1546 DO m = surf_s, surf_e … … 1607 1610 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 1608 1611 ! 1609 !-- Sum up all rain drop number densities which contribute to the flux 1612 !-- Sum up all rain drop number densities which contribute to the flux 1610 1613 !-- through k-1/2 1611 1614 flux = 0.0_wp … … 1730 1733 THEN 1731 1734 ! 1732 !-- Run over all upward-facing surface elements, i.e. non-natural, 1735 !-- Run over all upward-facing surface elements, i.e. non-natural, 1733 1736 !-- natural and urban 1734 1737 DO m = 1, bc_h(0)%ns 1735 i = bc_h(0)%i(m) 1738 i = bc_h(0)%i(m) 1736 1739 j = bc_h(0)%j(m) 1737 1740 k = bc_h(0)%k(m) … … 1961 1964 USE particle_attributes, & 1962 1965 ONLY: molecular_weight_of_solute, molecular_weight_of_water, rho_s, & 1963 s1, s2, s3, vanthoff1966 log_sigma, vanthoff 1964 1967 1965 1968 IMPLICIT NONE … … 2013 2016 !-- Prescribe parameters for activation 2014 2017 !-- (see: Bott + Trautmann, 2002, Atm. Res., 64) 2015 k_act = 0.7_wp 2018 k_act = 0.7_wp 2016 2019 2017 2020 IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN 2018 2021 ! 2019 !-- Compute the number of activated Aerosols 2022 !-- Compute the number of activated Aerosols 2020 2023 !-- (see: Twomey, 1959, Pure and applied Geophysics, 43) 2021 2024 n_act = na_init * sat**k_act … … 2027 2030 ! 2028 2031 !-- Compute activation rate after Khairoutdinov and Kogan 2029 !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128) 2030 sat_max = 0.8_wp / 100.0_wp 2031 activ = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN & 2032 ( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) / & 2032 !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128) 2033 sat_max = 0.8_wp / 100.0_wp 2034 activ = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN & 2035 ( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) / & 2033 2036 dt_micro 2034 2037 … … 2036 2039 ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk ) THEN 2037 2040 ! 2038 !-- Curvature effect (afactor) with surface tension 2041 !-- Curvature effect (afactor) with surface tension 2039 2042 !-- parameterization by Straka (2009) 2040 2043 sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp ) … … 2046 2049 2047 2050 ! 2048 !-- Prescribe power index that describes the soluble fraction 2051 !-- Prescribe power index that describes the soluble fraction 2049 2052 !-- of an aerosol particle (beta) and mean geometric radius of 2050 2053 !-- dry aerosol spectrum … … 2053 2056 rd0 = 0.05E-6_wp 2054 2057 ! 2055 !-- Calculate mean geometric supersaturation (s_0) with 2058 !-- Calculate mean geometric supersaturation (s_0) with 2056 2059 !-- parameterization by Khvorostyanov and Curry (2006) 2057 2060 s_0 = rd0 **(- ( 1.0_wp + beta_act ) ) * & … … 2059 2062 2060 2063 ! 2061 !-- Calculate number of activated CCN as a function of 2064 !-- Calculate number of activated CCN as a function of 2062 2065 !-- supersaturation and taking Koehler theory into account 2063 2066 !-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111) 2064 n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( & 2065 LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(s1) ) ) )2067 n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( & 2068 LOG( s_0 / sat ) / ( SQRT(2.0_wp) * log_sigma(1) ) ) ) 2066 2069 activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp ) 2067 2068 nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init) 2070 2071 nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init) 2069 2072 ENDIF 2070 2073 … … 2076 2079 ! Description: 2077 2080 ! ------------ 2078 !> Calculate condensation rate for cloud water content (after Khairoutdinov and 2081 !> Calculate condensation rate for cloud water content (after Khairoutdinov and 2079 2082 !> Kogan, 2000). 2080 2083 !------------------------------------------------------------------------------! … … 2154 2157 !-- Mean weight of cloud drops 2155 2158 IF ( nc_1d(k) <= 0.0_wp) CYCLE 2156 xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin) 2159 xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin) 2157 2160 ! 2158 2161 !-- Weight averaged diameter of cloud drops: … … 2160 2163 ! 2161 2164 !-- Integral diameter of cloud drops 2162 nc_0 = nc_1d(k) * dc 2165 nc_0 = nc_1d(k) * dc 2163 2166 ! 2164 2167 !-- Condensation needs only to be calculated in supersaturated regions … … 2171 2174 cond_max = q_1d(k) - q_s - qc_1d(k) - qr_1d(k) 2172 2175 cond = MIN( cond, cond_max / dt_micro ) 2173 2176 2174 2177 qc_1d(k) = qc_1d(k) + cond * dt_micro 2175 2178 ELSEIF ( sat < 0.0_wp ) THEN … … 2302 2305 nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag 2303 2306 IF ( microphysics_morrison ) THEN 2304 nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp * & 2307 nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp * & 2305 2308 autocon / x0 * hyrho(k) * dt_micro * flag ) 2306 2309 ENDIF … … 2425 2428 ! 2426 2429 !-- Mean weight of cloud drops 2427 xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin) 2430 xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin) 2428 2431 ! 2429 2432 !-- Parameterized turbulence effects on autoconversion (Seifert, … … 2716 2719 IF ( microphysics_morrison ) THEN 2717 2720 IF ( qc_1d(k) > eps_sb .AND. nc_1d(k) > eps_mr ) THEN 2718 sed_nc(k) = sed_qc_const * & 2721 sed_nc(k) = sed_qc_const * & 2719 2722 ( qc_1d(k) * hyrho(k) )**( 2.0_wp / 3.0_wp ) * & 2720 2723 ( nc_1d(k) )**( 1.0_wp / 3.0_wp ) … … 2790 2793 USE surface_mod, & 2791 2794 ONLY : bc_h 2792 2795 2793 2796 IMPLICIT NONE 2794 2797 … … 2859 2862 ENDDO 2860 2863 ! 2861 !-- Adjust boundary values using surface data type. 2864 !-- Adjust boundary values using surface data type. 2862 2865 !-- Upward facing non-natural 2863 surf_s = bc_h(0)%start_index(j,i) 2866 surf_s = bc_h(0)%start_index(j,i) 2864 2867 surf_e = bc_h(0)%end_index(j,i) 2865 2868 DO m = surf_s, surf_e … … 2870 2873 ! 2871 2874 !-- Downward facing non-natural 2872 surf_s = bc_h(1)%start_index(j,i) 2875 surf_s = bc_h(1)%start_index(j,i) 2873 2876 surf_e = bc_h(1)%end_index(j,i) 2874 2877 DO m = surf_s, surf_e … … 2935 2938 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 2936 2939 ! 2937 !-- Sum up all rain drop number densities which contribute to the flux 2940 !-- Sum up all rain drop number densities which contribute to the flux 2938 2941 !-- through k-1/2 2939 2942 flux = 0.0_wp … … 3044 3047 THEN 3045 3048 3046 surf_s = bc_h(0)%start_index(j,i) 3049 surf_s = bc_h(0)%start_index(j,i) 3047 3050 surf_e = bc_h(0)%end_index(j,i) 3048 3051 DO m = surf_s, surf_e
Note: See TracChangeset
for help on using the changeset viewer.
