!> @file bulk_cloud_model_mod.f90 !------------------------------------------------------------------------------! ! This file is part of the PALM model system. ! ! PALM is free software: you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation, either version 3 of the License, or (at your option) any later ! version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 1997-2018 Leibniz Universitaet Hannover !------------------------------------------------------------------------------! ! ! Current revisions: ! ------------------ ! ! ! Former revisions: ! ----------------- ! $Id: bulk_cloud_model_mod.f90 3371 2018-10-18 13:40:12Z knoop $ ! Modularization of all bulk cloud physics code components ! ! ! unused variables removed ! ! 3026 2018-05-22 10:30:53Z schwenkel ! Changed the name specific humidity to mixing ratio, since we are computing ! mixing ratios. ! ! 2718 2018-01-02 08:49:38Z maronga ! Corrected "Former revisions" section ! ! 2701 2017-12-15 15:40:50Z suehring ! Changes from last commit documented ! ! 2698 2017-12-14 18:46:24Z suehring ! Bugfix in get_topography_top_index ! ! 2696 2017-12-14 17:12:51Z kanani ! Change in file header (GPL part) ! ! 2608 2017-11-13 14:04:26Z schwenkel ! Calculation of supersaturation in external module (diagnostic_quantities_mod). ! Change: correct calculation of saturation specific humidity to saturation ! mixing ratio (the factor of 0.378 vanishes). ! ! 2522 2017-10-05 14:20:37Z schwenkel ! Minor bugfix ! ! 2375 2017-08-29 14:10:28Z schwenkel ! Improved aerosol initilization and some minor bugfixes ! for droplet sedimenation ! ! 2318 2017-07-20 17:27:44Z suehring ! Get topography top index via Function call ! ! 2317 2017-07-20 17:27:19Z suehring ! s1 changed to log_sigma ! ! 2292 2017-06-20 09:51:42Z schwenkel ! Implementation of new microphysic scheme: cloud_scheme = 'morrison' ! includes two more prognostic equations for cloud drop concentration (nc) ! and cloud water content (qc). ! - The process of activation is parameterized with a simple Twomey ! activion scheme or with considering solution and curvature ! effects (Khvorostyanov and Curry ,2006). ! - The saturation adjustment scheme is replaced by the parameterization ! of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128). ! - All other microphysical processes of Seifert and Beheng are used. ! Additionally, in those processes the reduction of cloud number concentration ! is considered. ! ! 2233 2017-05-30 18:08:54Z suehring ! ! 2232 2017-05-30 17:47:52Z suehring ! Adjustments to new topography and surface concept ! ! 2155 2017-02-21 09:57:40Z hoffmann ! Bugfix in the calculation of microphysical quantities on ghost points. ! ! 2031 2016-10-21 15:11:58Z knoop ! renamed variable rho to rho_ocean ! ! 2000 2016-08-20 18:09:15Z knoop ! Forced header and separation lines into 80 columns ! ! 1850 2016-04-08 13:29:27Z maronga ! Module renamed ! Adapted for modularization of microphysics. ! ! 1845 2016-04-08 08:29:13Z raasch ! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface ! point (instead of one point above surface) ! ! 1831 2016-04-07 13:15:51Z hoffmann ! turbulence renamed collision_turbulence, ! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also ! avaialble for microphysics_kessler. ! ! 1822 2016-04-07 07:49:42Z hoffmann ! Unused variables removed. ! Kessler scheme integrated. ! ! 1691 2015-10-26 16:17:44Z maronga ! Added new routine calc_precipitation_amount. The routine now allows to account ! for precipitation due to sedimenation of cloud (fog) droplets ! ! 1682 2015-10-07 23:56:08Z knoop ! Code annotations made doxygen readable ! ! 1646 2015-09-02 16:00:10Z hoffmann ! Bugfix: Wrong computation of d_mean. ! ! 1361 2014-04-16 15:17:48Z hoffmann ! Bugfix in sedimentation_rain: Index corrected. ! Vectorized version of adjust_cloud added. ! Little reformatting of the code. ! ! 1353 2014-04-08 15:21:23Z heinze ! REAL constants provided with KIND-attribute ! ! 1346 2014-03-27 13:18:20Z heinze ! Bugfix: REAL constants provided with KIND-attribute especially in call of ! intrinsic function like MAX, MIN, SIGN ! ! 1334 2014-03-25 12:21:40Z heinze ! Bugfix: REAL constants provided with KIND-attribute ! ! 1322 2014-03-20 16:38:49Z raasch ! REAL constants defined as wp-kind ! ! 1320 2014-03-20 08:40:49Z raasch ! ONLY-attribute added to USE-statements, ! kind-parameters added to all INTEGER and REAL declaration statements, ! kinds are defined in new module kinds, ! comment fields (!:) to be used for variable explanations added to ! all variable declaration statements ! ! 1241 2013-10-30 11:36:58Z heinze ! hyp and rho_ocean have to be calculated at each time step if data from external ! file LSF_DATA are used ! ! 1115 2013-03-26 18:16:16Z hoffmann ! microphyical tendencies are calculated in bcm_actions in an optimized ! way; unrealistic values are prevented; bugfix in evaporation; some reformatting ! ! 1106 2013-03-04 05:31:38Z raasch ! small changes in code formatting ! ! 1092 2013-02-02 11:24:22Z raasch ! unused variables removed ! file put under GPL ! ! 1065 2012-11-22 17:42:36Z hoffmann ! Sedimentation process implemented according to Stevens and Seifert (2008). ! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens ! and Stevens, 2010). ! ! 1053 2012-11-13 17:11:03Z hoffmann ! initial revision ! ! Description: ! ------------ !> Calculate bulk cloud microphysics. !------------------------------------------------------------------------------! MODULE bulk_cloud_model_mod USE arrays_3d, & #if defined (__nopointer) ONLY: ddzu, diss, dzu, dzw, hyp, hyrho, & nc, nc_p, nr, nr_p, & precipitation_amount, prr, pt, d_exner, pt_init, q, ql, & qc, qc_p, qr, qr_p, & exner, zu, tnc_m, tnr_m, tqc_m, tqr_m #else ONLY: ddzu, diss, dzu, dzw, hyp, hyrho, & nc, nc_1, nc_2, nc_3, nc_p, nr, nr_1, nr_2, nr_3, nr_p, & precipitation_amount, prr, pt, d_exner, pt_init, q, ql, ql_1, & qc, qc_1, qc_2, qc_3, qc_p, qr, qr_1, qr_2, qr_3, qr_p, & exner, zu, tnc_m, tnr_m, tqc_m, tqr_m #endif USE averaging, & ONLY: nc_av, nr_av, prr_av, qc_av, ql_av, qr_av USE basic_constants_and_equations_mod, & ONLY: c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, molecular_weight_of_solute,& molecular_weight_of_water, pi, rho_l, rho_s, r_d, r_v, vanthoff,& exner_function, exner_function_invers, ideal_gas_law_rho, & ideal_gas_law_rho_pt, barometric_formula, rd_d_rv USE control_parameters, & ONLY: dt_3d, dt_do2d_xy, intermediate_timestep_count, & intermediate_timestep_count_max, large_scale_forcing, & lsf_surf, pt_surface, rho_surface, surface_pressure, & time_do2d_xy, message_string USE cpulog, & ONLY: cpu_log, log_point_s USE grid_variables, & ONLY: dx, dy USE indices, & ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt, & wall_flags_0 USE kinds USE statistics, & ONLY: weight_pres, weight_substep USE surface_mod, & ONLY : bc_h, get_topography_top_index_ji, surf_bulk_cloud_model, & surf_microphysics_morrison, surf_microphysics_seifert IMPLICIT NONE CHARACTER (LEN=20) :: aerosol_bulk = 'nacl' !< namelist parameter CHARACTER (LEN=20) :: cloud_scheme = 'saturation_adjust' !< namelist parameter LOGICAL :: aerosol_nacl =.TRUE. !< nacl aerosol for bulk scheme LOGICAL :: aerosol_c3h4o4 =.FALSE. !< malonic acid aerosol for bulk scheme LOGICAL :: aerosol_nh4no3 =.FALSE. !< malonic acid aerosol for bulk scheme LOGICAL :: bulk_cloud_model = .FALSE. !< namelist parameter LOGICAL :: cloud_water_sedimentation = .FALSE. !< cloud water sedimentation LOGICAL :: curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory LOGICAL :: limiter_sedimentation = .TRUE. !< sedimentation limiter LOGICAL :: collision_turbulence = .FALSE. !< turbulence effects LOGICAL :: ventilation_effect = .TRUE. !< ventilation effect LOGICAL :: call_microphysics_at_all_substeps = .FALSE. !< namelist parameter LOGICAL :: microphysics_sat_adjust = .FALSE. !< use saturation adjust bulk scheme? LOGICAL :: microphysics_kessler = .FALSE. !< use kessler bulk scheme? LOGICAL :: microphysics_morrison = .FALSE. !< use 2-moment Morrison (add. prog. eq. for nc and qc) LOGICAL :: microphysics_seifert = .FALSE. !< use 2-moment Seifert and Beheng scheme LOGICAL :: precipitation = .FALSE. !< namelist parameter REAL(wp) :: precipitation_amount_interval = 9999999.9_wp !< namelist parameter REAL(wp) :: a_1 = 8.69E-4_wp !< coef. in turb. parametrization (cm-2 s3) REAL(wp) :: a_2 = -7.38E-5_wp !< coef. in turb. parametrization (cm-2 s3) REAL(wp) :: a_3 = -1.40E-2_wp !< coef. in turb. parametrization REAL(wp) :: a_term = 9.65_wp !< coef. for terminal velocity (m s-1) REAL(wp) :: a_vent = 0.78_wp !< coef. for ventilation effect REAL(wp) :: b_1 = 11.45E-6_wp !< coef. in turb. parametrization (m) REAL(wp) :: b_2 = 9.68E-6_wp !< coef. in turb. parametrization (m) REAL(wp) :: b_3 = 0.62_wp !< coef. in turb. parametrization REAL(wp) :: b_term = 9.8_wp !< coef. for terminal velocity (m s-1) REAL(wp) :: b_vent = 0.308_wp !< coef. for ventilation effect REAL(wp) :: beta_cc = 3.09E-4_wp !< coef. in turb. parametrization (cm-2 s3) REAL(wp) :: c_1 = 4.82E-6_wp !< coef. in turb. parametrization (m) REAL(wp) :: c_2 = 4.8E-6_wp !< coef. in turb. parametrization (m) REAL(wp) :: c_3 = 0.76_wp !< coef. in turb. parametrization REAL(wp) :: c_const = 0.93_wp !< const. in Taylor-microscale Reynolds number REAL(wp) :: c_evap = 0.7_wp !< constant in evaporation REAL(wp) :: c_term = 600.0_wp !< coef. for terminal velocity (m-1) REAL(wp) :: diff_coeff_l = 0.23E-4_wp !< diffusivity of water vapor (m2 s-1) REAL(wp) :: eps_sb = 1.0E-10_wp !< threshold in two-moments scheme REAL(wp) :: eps_mr = 0.0_wp !< threshold for morrison scheme REAL(wp) :: k_cc = 9.44E09_wp !< const. cloud-cloud kernel (m3 kg-2 s-1) REAL(wp) :: k_cr0 = 4.33_wp !< const. cloud-rain kernel (m3 kg-1 s-1) REAL(wp) :: k_rr = 7.12_wp !< const. rain-rain kernel (m3 kg-1 s-1) REAL(wp) :: k_br = 1000.0_wp !< const. in breakup parametrization (m-1) REAL(wp) :: k_st = 1.2E8_wp !< const. in drizzle parametrization (m-1 s-1) REAL(wp) :: kin_vis_air = 1.4086E-5_wp !< kin. viscosity of air (m2 s-1) REAL(wp) :: prec_time_const = 0.001_wp !< coef. in Kessler scheme (s-1) REAL(wp) :: ql_crit = 0.0005_wp !< coef. in Kessler scheme (kg kg-1) REAL(wp) :: schmidt_p_1d3=0.8921121_wp !< Schmidt number**0.33333, 0.71**0.33333 REAL(wp) :: sigma_gc = 1.3_wp !< geometric standard deviation cloud droplets REAL(wp) :: thermal_conductivity_l = 2.43E-2_wp !< therm. cond. air (J m-1 s-1 K-1) REAL(wp) :: w_precipitation = 9.65_wp !< maximum terminal velocity (m s-1) REAL(wp) :: x0 = 2.6E-10_wp !< separating drop mass (kg) REAL(wp) :: xamin = 5.24E-19_wp !< average aerosol mass (kg) (~ 0.05µm) REAL(wp) :: xcmin = 4.18E-15_wp !< minimum cloud drop size (kg) (~ 1µm) REAL(wp) :: xrmin = 2.6E-10_wp !< minimum rain drop size (kg) REAL(wp) :: xrmax = 5.0E-6_wp !< maximum rain drop site (kg) REAL(wp) :: c_sedimentation = 2.0_wp !< Courant number of sedimentation process REAL(wp) :: dpirho_l !< 6.0 / ( pi * rho_l ) REAL(wp) :: dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius REAL(wp) :: dt_micro !< microphysics time step REAL(wp) :: sigma_bulk = 2.0_wp !< width of aerosol spectrum REAL(wp) :: na_init = 100.0E6_wp !< Total particle/aerosol concentration (cm-3) REAL(wp) :: nc_const = 70.0E6_wp !< cloud droplet concentration REAL(wp) :: dt_precipitation = 100.0_wp !< timestep precipitation (s) REAL(wp) :: sed_qc_const !< const. for sedimentation of cloud water REAL(wp) :: pirho_l !< pi * rho_l / 6.0; REAL(wp) :: e_s !< saturation water vapor pressure REAL(wp) :: q_s !< saturation mixing ratio REAL(wp) :: sat !< supersaturation REAL(wp) :: t_l !< actual temperature SAVE PRIVATE PUBLIC bcm_parin, & bcm_check_parameters, & bcm_check_data_output, & bcm_check_data_output_pr, & bcm_header, & bcm_init_arrays, & bcm_init, & bcm_3d_data_averaging, & bcm_data_output_2d, & bcm_data_output_3d, & bcm_swap_timelevel, & bcm_rrd_global, & bcm_rrd_local, & bcm_wrd_global, & bcm_wrd_local, & bcm_actions, & calc_liquid_water_content PUBLIC call_microphysics_at_all_substeps, & cloud_water_sedimentation, & bulk_cloud_model, & cloud_scheme, & collision_turbulence, & dt_precipitation, & microphysics_morrison, & microphysics_sat_adjust, & microphysics_seifert, & na_init, & nc_const, & precipitation, & sigma_gc INTERFACE bcm_parin MODULE PROCEDURE bcm_parin END INTERFACE bcm_parin INTERFACE bcm_check_parameters MODULE PROCEDURE bcm_check_parameters END INTERFACE bcm_check_parameters INTERFACE bcm_check_data_output MODULE PROCEDURE bcm_check_data_output END INTERFACE bcm_check_data_output INTERFACE bcm_check_data_output_pr MODULE PROCEDURE bcm_check_data_output_pr END INTERFACE bcm_check_data_output_pr INTERFACE bcm_header MODULE PROCEDURE bcm_header END INTERFACE bcm_header INTERFACE bcm_init_arrays MODULE PROCEDURE bcm_init_arrays END INTERFACE bcm_init_arrays INTERFACE bcm_init MODULE PROCEDURE bcm_init END INTERFACE bcm_init INTERFACE bcm_3d_data_averaging MODULE PROCEDURE bcm_3d_data_averaging END INTERFACE bcm_3d_data_averaging INTERFACE bcm_data_output_2d MODULE PROCEDURE bcm_data_output_2d END INTERFACE bcm_data_output_2d INTERFACE bcm_data_output_3d MODULE PROCEDURE bcm_data_output_3d END INTERFACE bcm_data_output_3d INTERFACE bcm_swap_timelevel MODULE PROCEDURE bcm_swap_timelevel END INTERFACE bcm_swap_timelevel INTERFACE bcm_rrd_global MODULE PROCEDURE bcm_rrd_global END INTERFACE bcm_rrd_global INTERFACE bcm_rrd_local MODULE PROCEDURE bcm_rrd_local END INTERFACE bcm_rrd_local INTERFACE bcm_wrd_global MODULE PROCEDURE bcm_wrd_global END INTERFACE bcm_wrd_global INTERFACE bcm_wrd_local MODULE PROCEDURE bcm_wrd_local END INTERFACE bcm_wrd_local INTERFACE bcm_actions MODULE PROCEDURE bcm_actions MODULE PROCEDURE bcm_actions_ij END INTERFACE bcm_actions INTERFACE adjust_cloud MODULE PROCEDURE adjust_cloud MODULE PROCEDURE adjust_cloud_ij END INTERFACE adjust_cloud INTERFACE activation MODULE PROCEDURE activation MODULE PROCEDURE activation_ij END INTERFACE activation INTERFACE condensation MODULE PROCEDURE condensation MODULE PROCEDURE condensation_ij END INTERFACE condensation INTERFACE autoconversion MODULE PROCEDURE autoconversion MODULE PROCEDURE autoconversion_ij END INTERFACE autoconversion INTERFACE autoconversion_kessler MODULE PROCEDURE autoconversion_kessler MODULE PROCEDURE autoconversion_kessler_ij END INTERFACE autoconversion_kessler INTERFACE accretion MODULE PROCEDURE accretion MODULE PROCEDURE accretion_ij END INTERFACE accretion INTERFACE selfcollection_breakup MODULE PROCEDURE selfcollection_breakup MODULE PROCEDURE selfcollection_breakup_ij END INTERFACE selfcollection_breakup INTERFACE evaporation_rain MODULE PROCEDURE evaporation_rain MODULE PROCEDURE evaporation_rain_ij END INTERFACE evaporation_rain INTERFACE sedimentation_cloud MODULE PROCEDURE sedimentation_cloud MODULE PROCEDURE sedimentation_cloud_ij END INTERFACE sedimentation_cloud INTERFACE sedimentation_rain MODULE PROCEDURE sedimentation_rain MODULE PROCEDURE sedimentation_rain_ij END INTERFACE sedimentation_rain INTERFACE calc_precipitation_amount MODULE PROCEDURE calc_precipitation_amount MODULE PROCEDURE calc_precipitation_amount_ij END INTERFACE calc_precipitation_amount INTERFACE supersaturation MODULE PROCEDURE supersaturation END INTERFACE supersaturation INTERFACE calc_liquid_water_content MODULE PROCEDURE calc_liquid_water_content END INTERFACE calc_liquid_water_content CONTAINS !------------------------------------------------------------------------------! ! Description: ! ------------ !> Parin for &bulk_cloud_parameters for the microphysics module !------------------------------------------------------------------------------! SUBROUTINE bcm_parin IMPLICIT NONE CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file NAMELIST /bulk_cloud_parameters/ & aerosol_bulk, & c_sedimentation, & call_microphysics_at_all_substeps, & bulk_cloud_model, & cloud_scheme, & cloud_water_sedimentation, & collision_turbulence, & curvature_solution_effects_bulk, & dry_aerosol_radius, & limiter_sedimentation, & na_init, & nc_const, & precipitation, & precipitation_amount_interval, & sigma_bulk, & ventilation_effect line = ' ' ! !-- Try to find microphysics module package REWIND ( 11 ) line = ' ' DO WHILE ( INDEX( line, '&bulk_cloud_parameters' ) == 0 ) READ ( 11, '(A)', END=10 ) line ENDDO BACKSPACE ( 11 ) ! !-- Read user-defined namelist READ ( 11, bulk_cloud_parameters ) ! !-- Set flag that indicates that the microphysics module is switched on !bulk_cloud_model = .TRUE. 10 CONTINUE END SUBROUTINE bcm_parin !------------------------------------------------------------------------------! ! Description: ! ------------ !> Check parameters routine for microphysics module !------------------------------------------------------------------------------! SUBROUTINE bcm_check_parameters IMPLICIT NONE ! !-- Check cloud scheme IF ( cloud_scheme == 'saturation_adjust' ) THEN microphysics_sat_adjust = .TRUE. microphysics_seifert = .FALSE. microphysics_kessler = .FALSE. precipitation = .FALSE. ELSEIF ( cloud_scheme == 'seifert_beheng' ) THEN microphysics_sat_adjust = .FALSE. microphysics_seifert = .TRUE. microphysics_kessler = .FALSE. microphysics_morrison = .FALSE. precipitation = .TRUE. ELSEIF ( cloud_scheme == 'kessler' ) THEN microphysics_sat_adjust = .FALSE. microphysics_seifert = .FALSE. microphysics_kessler = .TRUE. microphysics_morrison = .FALSE. precipitation = .TRUE. ELSEIF ( cloud_scheme == 'morrison' ) THEN microphysics_sat_adjust = .FALSE. microphysics_seifert = .TRUE. microphysics_kessler = .FALSE. microphysics_morrison = .TRUE. precipitation = .TRUE. ELSE message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // & TRIM( cloud_scheme ) // '"' CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 ) ENDIF ! !-- Set the default value for the integration interval of precipitation amount IF ( microphysics_seifert .OR. microphysics_kessler ) THEN IF ( precipitation_amount_interval == 9999999.9_wp ) THEN precipitation_amount_interval = dt_do2d_xy ELSE IF ( precipitation_amount_interval > dt_do2d_xy ) THEN WRITE( message_string, * ) 'precipitation_amount_interval = ', & precipitation_amount_interval, ' must not be larger than ', & 'dt_do2d_xy = ', dt_do2d_xy CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 ) ENDIF ENDIF ENDIF ! TODO: find better sollution for circular dependency problem surf_bulk_cloud_model = bulk_cloud_model surf_microphysics_morrison = microphysics_morrison surf_microphysics_seifert = microphysics_seifert ! !-- Check aerosol IF ( aerosol_bulk == 'nacl' ) THEN aerosol_nacl = .TRUE. aerosol_c3h4o4 = .FALSE. aerosol_nh4no3 = .FALSE. ELSEIF ( aerosol_bulk == 'c3h4o4' ) THEN aerosol_nacl = .FALSE. aerosol_c3h4o4 = .TRUE. aerosol_nh4no3 = .FALSE. ELSEIF ( aerosol_bulk == 'nh4no3' ) THEN aerosol_nacl = .FALSE. aerosol_c3h4o4 = .FALSE. aerosol_nh4no3 = .TRUE. ELSE message_string = 'unknown aerosol = "' // TRIM( aerosol_bulk ) // '"' CALL message( 'check_parameters', 'PA0469', 1, 2, 0, 6, 0 ) ENDIF END SUBROUTINE bcm_check_parameters !------------------------------------------------------------------------------! ! Description: ! ------------ !> Check data output for microphysics module !------------------------------------------------------------------------------! SUBROUTINE bcm_check_data_output( var, unit ) IMPLICIT NONE CHARACTER (LEN=*) :: unit !< CHARACTER (LEN=*) :: var !< SELECT CASE ( TRIM( var ) ) CASE ( 'nc' ) IF ( .NOT. microphysics_morrison ) THEN message_string = 'output of "' // TRIM( var ) // '" ' // & 'requires ' // & 'cloud_scheme = "morrison"' CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 ) ENDIF unit = '1/m3' CASE ( 'nr' ) IF ( .NOT. microphysics_seifert ) THEN message_string = 'output of "' // TRIM( var ) // '" ' // & 'requires ' // & 'cloud_scheme = "seifert_beheng"' CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 ) ENDIF unit = '1/m3' CASE ( 'prr' ) IF ( microphysics_sat_adjust ) THEN message_string = 'output of "' // TRIM( var ) // '" ' // & 'is not available for ' // & 'cloud_scheme = "saturation_adjust"' CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 ) ENDIF unit = 'kg/kg m/s' CASE ( 'qc' ) unit = 'kg/kg' CASE ( 'qr' ) IF ( .NOT. microphysics_seifert ) THEN message_string = 'output of "' // TRIM( var ) // '" ' // & 'requires ' // & 'cloud_scheme = "seifert_beheng"' CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 ) ENDIF unit = 'kg/kg' CASE ( 'pra*' ) IF ( .NOT. microphysics_kessler .AND. & .NOT. microphysics_seifert ) THEN message_string = 'output of "' // TRIM( var ) // '" ' // & 'requires ' // & 'cloud_scheme = "kessler" or "seifert_beheng"' CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 ) ENDIF ! TODO: find sollution (maybe connected to flow_statistics redesign?) ! IF ( j == 1 ) THEN ! message_string = 'temporal averaging of precipitation ' // & ! 'amount "' // TRIM( var ) // '" is not possible' ! CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 ) ! ENDIF unit = 'mm' CASE ( 'prr*' ) IF ( .NOT. microphysics_kessler .AND. & .NOT. microphysics_seifert ) THEN message_string = 'output of "' // TRIM( var ) // '"' // & ' requires' // & ' cloud_scheme = "kessler" or "seifert_beheng"' CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 ) ENDIF unit = 'mm/s' CASE DEFAULT unit = 'illegal' END SELECT END SUBROUTINE bcm_check_data_output !------------------------------------------------------------------------------! ! Description: ! ------------ !> Check data output of profiles for microphysics module !------------------------------------------------------------------------------! SUBROUTINE bcm_check_data_output_pr( variable, var_count, unit, dopr_unit ) USE arrays_3d, & ONLY: zu USE control_parameters, & ONLY: data_output_pr USE profil_parameter, & ONLY: dopr_index USE statistics, & ONLY: hom, statistic_regions, pr_palm IMPLICIT NONE CHARACTER (LEN=*) :: unit !< CHARACTER (LEN=*) :: variable !< CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit INTEGER(iwp) :: var_count !< INTEGER(iwp) :: pr_index !< SELECT CASE ( TRIM( variable ) ) ! TODO: make index generic: pr_index = pr_palm+1 CASE ( 'nc' ) IF ( .NOT. microphysics_morrison ) THEN message_string = 'data_output_pr = ' // & TRIM( data_output_pr(var_count) ) // & ' is not implemented for' // & ' cloud_scheme /= morrison' CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 ) ENDIF pr_index = 89 dopr_index(var_count) = pr_index dopr_unit = '1/m3' unit = dopr_unit hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 ) CASE ( 'nr' ) IF ( .NOT. microphysics_seifert ) THEN message_string = 'data_output_pr = ' // & TRIM( data_output_pr(var_count) ) // & ' is not implemented for' // & ' cloud_scheme /= seifert_beheng' CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 ) ENDIF pr_index = 73 dopr_index(var_count) = pr_index dopr_unit = '1/m3' unit = dopr_unit hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 ) CASE ( 'prr' ) IF ( microphysics_sat_adjust ) THEN message_string = 'data_output_pr = ' // & TRIM( data_output_pr(var_count) ) // & ' is not available for' // & ' cloud_scheme = saturation_adjust' CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 ) ENDIF pr_index = 76 dopr_index(var_count) = pr_index dopr_unit = 'kg/kg m/s' unit = dopr_unit hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 ) CASE ( 'qc' ) pr_index = 75 dopr_index(var_count) = pr_index dopr_unit = 'kg/kg' unit = dopr_unit hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 ) CASE ( 'qr' ) IF ( .NOT. microphysics_seifert ) THEN message_string = 'data_output_pr = ' // & TRIM( data_output_pr(var_count) ) // & ' is not implemented for' // & ' cloud_scheme /= seifert_beheng' CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 ) ENDIF pr_index = 744 dopr_index(var_count) = pr_index dopr_unit = 'kg/kg' unit = dopr_unit hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 ) CASE DEFAULT unit = 'illegal' END SELECT END SUBROUTINE bcm_check_data_output_pr !------------------------------------------------------------------------------! ! Description: ! ------------ !> Allocate microphysics module arrays and define pointers !------------------------------------------------------------------------------! SUBROUTINE bcm_init_arrays USE indices, & ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< ! !-- Liquid water content #if defined( __nopointer ) ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) #else ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) #endif ! !-- 3D-cloud water content IF ( .NOT. microphysics_morrison ) THEN #if defined( __nopointer ) ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) #else ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) #endif ENDIF ! !-- Precipitation amount and rate (only needed if output is switched) ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg) ) ! !-- 3d-precipitation rate ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) IF ( microphysics_morrison ) THEN ! !-- 3D-cloud drop water content, cloud drop concentration arrays #if defined( __nopointer ) ALLOCATE( nc(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & nc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & tnc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & tqc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) #else ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) #endif ENDIF IF ( microphysics_seifert ) THEN ! !-- 3D-rain water content, rain drop concentration arrays #if defined( __nopointer ) ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) #else ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) #endif ENDIF #if ! defined( __nopointer ) ! !-- Initial assignment of the pointers ql => ql_1 IF ( .NOT. microphysics_morrison ) THEN qc => qc_1 ENDIF IF ( microphysics_morrison ) THEN qc => qc_1; qc_p => qc_2; tqc_m => qc_3 nc => nc_1; nc_p => nc_2; tnc_m => nc_3 ENDIF IF ( microphysics_seifert ) THEN qr => qr_1; qr_p => qr_2; tqr_m => qr_3 nr => nr_1; nr_p => nr_2; tnr_m => nr_3 ENDIF #endif END SUBROUTINE bcm_init_arrays !------------------------------------------------------------------------------! ! Description: ! ------------ !> Initialization of the microphysics module !------------------------------------------------------------------------------! SUBROUTINE bcm_init !( dots_label, dots_unit, dots_num, dots_max ) IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< ! INTEGER(iwp) :: dots_num ! INTEGER(iwp) :: dots_max ! CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit ! CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label CALL location_message( 'initializing microphysics module', .FALSE. ) IF ( bulk_cloud_model ) THEN ! dots_label(dots_num+1) = 'some_var' ! dots_unit(dots_num+1) = 'm/s' ! ! dots_num_palm = dots_num ! dots_num = dots_num + 1 ! ! ! !-- Stuff for the module ! ! abs_velocity = 0.0_wp ! !-- Initialize the remaining quantities IF ( microphysics_morrison ) THEN DO i = nxlg, nxrg DO j = nysg, nyng qc(:,j,i) = 0.0_wp nc(:,j,i) = 0.0_wp ENDDO ENDDO ENDIF IF ( microphysics_seifert ) THEN DO i = nxlg, nxrg DO j = nysg, nyng qr(:,j,i) = 0.0_wp nr(:,j,i) = 0.0_wp ENDDO ENDDO ENDIF ! !-- Liquid water content and precipitation amount !-- are zero at beginning of the simulation IF ( bulk_cloud_model ) THEN ql = 0.0_wp ! TODO ??? qc = 0.0_wp precipitation_amount = 0.0_wp ENDIF ! !-- Initialize old and new time levels. IF ( microphysics_morrison ) THEN tqc_m = 0.0_wp tnc_m = 0.0_wp qc_p = qc nc_p = nc ENDIF IF ( microphysics_seifert ) THEN tqr_m = 0.0_wp tnr_m = 0.0_wp qr_p = qr nr_p = nr ENDIF ! !-- constant for the sedimentation of cloud water (2-moment cloud physics) sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ) & )**( 2.0_wp / 3.0_wp ) * & EXP( 5.0_wp * LOG( sigma_gc )**2 ) ! !-- Calculate timestep according to precipitation IF ( microphysics_seifert ) THEN dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) / & w_precipitation ENDIF ! !-- Set constants for certain aerosol type IF ( microphysics_morrison ) THEN IF ( aerosol_nacl ) THEN molecular_weight_of_solute = 0.05844_wp rho_s = 2165.0_wp vanthoff = 2.0_wp ELSEIF ( aerosol_c3h4o4 ) THEN molecular_weight_of_solute = 0.10406_wp rho_s = 1600.0_wp vanthoff = 1.37_wp ELSEIF ( aerosol_nh4no3 ) THEN molecular_weight_of_solute = 0.08004_wp rho_s = 1720.0_wp vanthoff = 2.31_wp ENDIF ENDIF ! !-- Pre-calculate frequently calculated fractions of pi and rho_l pirho_l = pi * rho_l / 6.0_wp dpirho_l = 1.0_wp / pirho_l CALL location_message( 'finished', .TRUE. ) ELSE CALL location_message( 'skipped', .TRUE. ) ENDIF END SUBROUTINE bcm_init !------------------------------------------------------------------------------! ! Description: ! ------------ !> Swapping of timelevels !------------------------------------------------------------------------------! SUBROUTINE bcm_swap_timelevel ( mod_count ) IMPLICIT NONE INTEGER, INTENT(IN) :: mod_count IF ( bulk_cloud_model ) THEN #if defined( __nopointer ) IF ( microphysics_morrison ) THEN qc = qc_p nc = nc_p ENDIF IF ( microphysics_seifert ) THEN qr = qr_p nr = nr_p ENDIF #else SELECT CASE ( mod_count ) CASE ( 0 ) IF ( microphysics_morrison ) THEN qc => qc_1; qc_p => qc_2 nc => nc_1; nc_p => nc_2 ENDIF IF ( microphysics_seifert ) THEN qr => qr_1; qr_p => qr_2 nr => nr_1; nr_p => nr_2 ENDIF CASE ( 1 ) IF ( microphysics_morrison ) THEN qc => qc_2; qc_p => qc_1 nc => nc_2; nc_p => nc_1 ENDIF IF ( microphysics_seifert ) THEN qr => qr_2; qr_p => qr_1 nr => nr_2; nr_p => nr_1 ENDIF END SELECT #endif ENDIF END SUBROUTINE bcm_swap_timelevel !------------------------------------------------------------------------------! ! Description: ! ------------ !> Header output for microphysics module !------------------------------------------------------------------------------! SUBROUTINE bcm_header ( io ) IMPLICIT NONE INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file ! !-- Write microphysics module header WRITE ( io, 1 ) WRITE ( io, 2 ) WRITE ( io, 3 ) IF ( microphysics_kessler ) THEN WRITE ( io, 4 ) 'Kessler-Scheme' ENDIF IF ( microphysics_seifert ) THEN WRITE ( io, 4 ) 'Seifert-Beheng-Scheme' IF ( cloud_water_sedimentation ) WRITE ( io, 5 ) IF ( collision_turbulence ) WRITE ( io, 6 ) IF ( ventilation_effect ) WRITE ( io, 7 ) IF ( limiter_sedimentation ) WRITE ( io, 8 ) ENDIF WRITE ( io, 20 ) WRITE ( io, 21 ) surface_pressure WRITE ( io, 22 ) r_d WRITE ( io, 23 ) rho_surface WRITE ( io, 24 ) c_p WRITE ( io, 25 ) l_v IF ( microphysics_seifert ) THEN WRITE ( io, 26 ) 1.0E-6_wp * nc_const WRITE ( io, 27 ) c_sedimentation ENDIF 1 FORMAT ( //' Bulk cloud microphysics module information:'/ & ' ------------------------------------------'/ ) 2 FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/ & ' total water content is used.' ) 3 FORMAT ( '--> Condensation is parameterized via 0% - or 100% scheme.' ) 4 FORMAT ( '--> Precipitation parameterization via ', A ) 5 FORMAT ( '--> Cloud water sedimentation parameterization via Stokes law' ) 6 FORMAT ( '--> Turbulence effects on precipitation process' ) 7 FORMAT ( '--> Ventilation effects on evaporation of rain drops' ) 8 FORMAT ( '--> Slope limiter used for sedimentation process' ) 20 FORMAT ( '--> Essential parameters:' ) 21 FORMAT ( ' Surface pressure : p_0 = ', F7.2, ' hPa') 22 FORMAT ( ' Gas constant : R = ', F5.1, ' J/(kg K)') 23 FORMAT ( ' Density of air : rho_0 = ', F6.3, ' kg/m**3') 24 FORMAT ( ' Specific heat cap. : c_p = ', F6.1, ' J/(kg K)') 25 FORMAT ( ' Vapourization heat : L_v = ', E9.2, ' J/kg') 26 FORMAT ( ' Droplet density : N_c = ', F6.1, ' 1/cm**3' ) 27 FORMAT ( ' Sedimentation Courant number : C_s = ', F4.1 ) END SUBROUTINE bcm_header !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine for averaging 3D data !------------------------------------------------------------------------------! SUBROUTINE bcm_3d_data_averaging( mode, variable ) USE control_parameters, & ONLY: average_count_3d IMPLICIT NONE CHARACTER (LEN=*) :: mode !< CHARACTER (LEN=*) :: variable !< INTEGER(iwp) :: i !< local index INTEGER(iwp) :: j !< local index INTEGER(iwp) :: k !< local index IF ( mode == 'allocate' ) THEN SELECT CASE ( TRIM( variable ) ) CASE ( 'nc' ) IF ( .NOT. ALLOCATED( nc_av ) ) THEN ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF nc_av = 0.0_wp CASE ( 'nr' ) IF ( .NOT. ALLOCATED( nr_av ) ) THEN ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF nr_av = 0.0_wp CASE ( 'prr' ) IF ( .NOT. ALLOCATED( prr_av ) ) THEN ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF prr_av = 0.0_wp CASE ( 'qc' ) IF ( .NOT. ALLOCATED( qc_av ) ) THEN ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF qc_av = 0.0_wp CASE ( 'ql' ) IF ( .NOT. ALLOCATED( ql_av ) ) THEN ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF ql_av = 0.0_wp CASE ( 'qr' ) IF ( .NOT. ALLOCATED( qr_av ) ) THEN ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF qr_av = 0.0_wp CASE DEFAULT CONTINUE END SELECT ELSEIF ( mode == 'sum' ) THEN SELECT CASE ( TRIM( variable ) ) CASE ( 'nc' ) IF ( ALLOCATED( nc_av ) ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 nc_av(k,j,i) = nc_av(k,j,i) + nc(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'nr' ) IF ( ALLOCATED( nr_av ) ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 nr_av(k,j,i) = nr_av(k,j,i) + nr(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'prr' ) IF ( ALLOCATED( prr_av ) ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 prr_av(k,j,i) = prr_av(k,j,i) + prr(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'qc' ) IF ( ALLOCATED( qc_av ) ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 qc_av(k,j,i) = qc_av(k,j,i) + qc(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'ql' ) IF ( ALLOCATED( ql_av ) ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 ql_av(k,j,i) = ql_av(k,j,i) + ql(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'qr' ) IF ( ALLOCATED( qr_av ) ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 qr_av(k,j,i) = qr_av(k,j,i) + qr(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE DEFAULT CONTINUE END SELECT ELSEIF ( mode == 'average' ) THEN SELECT CASE ( TRIM( variable ) ) CASE ( 'nc' ) IF ( ALLOCATED( nc_av ) ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 nc_av(k,j,i) = nc_av(k,j,i) / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO ENDIF CASE ( 'nr' ) IF ( ALLOCATED( nr_av ) ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 nr_av(k,j,i) = nr_av(k,j,i) / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO ENDIF CASE ( 'prr' ) IF ( ALLOCATED( prr_av ) ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 prr_av(k,j,i) = prr_av(k,j,i) / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO ENDIF CASE ( 'qc' ) IF ( ALLOCATED( qc_av ) ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 qc_av(k,j,i) = qc_av(k,j,i) / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO ENDIF CASE ( 'ql' ) IF ( ALLOCATED( ql_av ) ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 ql_av(k,j,i) = ql_av(k,j,i) / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO ENDIF CASE ( 'qr' ) IF ( ALLOCATED( qr_av ) ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 qr_av(k,j,i) = qr_av(k,j,i) / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO ENDIF CASE DEFAULT CONTINUE END SELECT ENDIF END SUBROUTINE bcm_3d_data_averaging !------------------------------------------------------------------------------! ! Description: ! ------------ !> Define 2D output variables. !------------------------------------------------------------------------------! SUBROUTINE bcm_data_output_2d( av, variable, found, grid, mode, local_pf, & two_d, nzb_do, nzt_do ) IMPLICIT NONE CHARACTER (LEN=*), INTENT(INOUT) :: grid !< name of vertical grid CHARACTER (LEN=*), INTENT(IN) :: mode !< either 'xy', 'xz' or 'yz' CHARACTER (LEN=*), INTENT(IN) :: variable !< name of variable INTEGER(iwp), INTENT(IN) :: av !< flag for (non-)average output INTEGER(iwp), INTENT(IN) :: nzb_do !< vertical output index (bottom) INTEGER(iwp), INTENT(IN) :: nzt_do !< vertical output index (top) INTEGER(iwp) :: flag_nr !< number of masking flag INTEGER(iwp) :: i !< loop index along x-direction INTEGER(iwp) :: j !< loop index along y-direction INTEGER(iwp) :: k !< loop index along z-direction LOGICAL, INTENT(INOUT) :: found !< flag if output variable is found LOGICAL, INTENT(INOUT) :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections) LOGICAL :: resorted !< flag if output is already resorted REAL(wp), PARAMETER :: fill_value = -999.0_wp !< value for the _FillValue attribute REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) :: local_pf !< local !< array to which output data is resorted to REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable found = .TRUE. resorted = .FALSE. ! !-- Set masking flag for topography for not resorted arrays flag_nr = 0 ! 0 = scalar, 1 = u, 2 = v, 3 = w SELECT CASE ( TRIM( variable ) ) CASE ( 'nc_xy', 'nc_xz', 'nc_yz' ) IF ( av == 0 ) THEN to_be_resorted => nc ELSE IF ( .NOT. ALLOCATED( nc_av ) ) THEN ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) nc_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => nc_av ENDIF IF ( mode == 'xy' ) grid = 'zu' CASE ( 'nr_xy', 'nr_xz', 'nr_yz' ) IF ( av == 0 ) THEN to_be_resorted => nr ELSE IF ( .NOT. ALLOCATED( nr_av ) ) THEN ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) nr_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => nr_av ENDIF IF ( mode == 'xy' ) grid = 'zu' CASE ( 'pra*_xy' ) ! 2d-array / integral quantity => no av ! CALL exchange_horiz_2d( precipitation_amount ) DO i = nxl, nxr DO j = nys, nyn local_pf(i,j,nzb+1) = precipitation_amount(j,i) ENDDO ENDDO precipitation_amount = 0.0_wp ! reset for next integ. interval resorted = .TRUE. two_d = .TRUE. IF ( mode == 'xy' ) grid = 'zu1' CASE ( 'prr_xy', 'prr_xz', 'prr_yz' ) IF ( av == 0 ) THEN ! CALL exchange_horiz( prr, nbgp ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 local_pf(i,j,k) = prr(k,j,i) * hyrho(nzb+1) ENDDO ENDDO ENDDO ELSE IF ( .NOT. ALLOCATED( prr_av ) ) THEN ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) prr_av = REAL( fill_value, KIND = wp ) ENDIF ! CALL exchange_horiz( prr_av, nbgp ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 local_pf(i,j,k) = prr_av(k,j,i) * hyrho(nzb+1) ENDDO ENDDO ENDDO ENDIF resorted = .TRUE. IF ( mode == 'xy' ) grid = 'zu' CASE ( 'qc_xy', 'qc_xz', 'qc_yz' ) IF ( av == 0 ) THEN to_be_resorted => qc ELSE IF ( .NOT. ALLOCATED( qc_av ) ) THEN ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) qc_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => qc_av ENDIF IF ( mode == 'xy' ) grid = 'zu' CASE ( 'ql_xy', 'ql_xz', 'ql_yz' ) IF ( av == 0 ) THEN to_be_resorted => ql ELSE IF ( .NOT. ALLOCATED( ql_av ) ) THEN ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ql_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => ql_av ENDIF IF ( mode == 'xy' ) grid = 'zu' CASE ( 'qr_xy', 'qr_xz', 'qr_yz' ) IF ( av == 0 ) THEN to_be_resorted => qr ELSE IF ( .NOT. ALLOCATED( qr_av ) ) THEN ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) qr_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => qr_av ENDIF IF ( mode == 'xy' ) grid = 'zu' CASE DEFAULT found = .FALSE. grid = 'none' END SELECT IF ( found .AND. .NOT. resorted ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = MERGE( & to_be_resorted(k,j,i), & REAL( fill_value, KIND = wp ), & BTEST( wall_flags_0(k,j,i), flag_nr ) & ) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE bcm_data_output_2d !------------------------------------------------------------------------------! ! Description: ! ------------ !> Define 3D output variables. !------------------------------------------------------------------------------! SUBROUTINE bcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) IMPLICIT NONE CHARACTER (LEN=*), INTENT(IN) :: variable !< name of variable INTEGER(iwp), INTENT(IN) :: av !< flag for (non-)average output INTEGER(iwp), INTENT(IN) :: nzb_do !< lower limit of the data output (usually 0) INTEGER(iwp), INTENT(IN) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d) INTEGER(iwp) :: flag_nr !< number of masking flag INTEGER(iwp) :: i !< loop index along x-direction INTEGER(iwp) :: j !< loop index along y-direction INTEGER(iwp) :: k !< loop index along z-direction LOGICAL, INTENT(INOUT) :: found !< flag if output variable is found LOGICAL :: resorted !< flag if output is already resorted REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) :: local_pf !< local !< array to which output data is resorted to REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable found = .TRUE. resorted = .FALSE. ! !-- Set masking flag for topography for not resorted arrays flag_nr = 0 ! 0 = scalar, 1 = u, 2 = v, 3 = w SELECT CASE ( TRIM( variable ) ) CASE ( 'nc' ) IF ( av == 0 ) THEN to_be_resorted => nc ELSE IF ( .NOT. ALLOCATED( nc_av ) ) THEN ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) nc_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => nc_av ENDIF CASE ( 'nr' ) IF ( av == 0 ) THEN to_be_resorted => nr ELSE IF ( .NOT. ALLOCATED( nr_av ) ) THEN ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) nr_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => nr_av ENDIF CASE ( 'prr' ) IF ( av == 0 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = prr(k,j,i) ENDDO ENDDO ENDDO ELSE IF ( .NOT. ALLOCATED( prr_av ) ) THEN ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) prr_av = REAL( fill_value, KIND = wp ) ENDIF DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = prr_av(k,j,i) ENDDO ENDDO ENDDO ENDIF resorted = .TRUE. CASE ( 'qc' ) IF ( av == 0 ) THEN to_be_resorted => qc ELSE IF ( .NOT. ALLOCATED( qc_av ) ) THEN ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) qc_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => qc_av ENDIF CASE ( 'ql' ) IF ( av == 0 ) THEN to_be_resorted => ql ELSE IF ( .NOT. ALLOCATED( ql_av ) ) THEN ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ql_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => ql_av ENDIF CASE ( 'qr' ) IF ( av == 0 ) THEN to_be_resorted => qr ELSE IF ( .NOT. ALLOCATED( qr_av ) ) THEN ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) qr_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => qr_av ENDIF CASE DEFAULT found = .FALSE. END SELECT IF ( found .AND. .NOT. resorted ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = MERGE( & to_be_resorted(k,j,i), & REAL( fill_value, KIND = wp ), & BTEST( wall_flags_0(k,j,i), flag_nr ) & ) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE bcm_data_output_3d !------------------------------------------------------------------------------! ! Description: ! ------------ !> This routine reads the respective restart data for the microphysics module. !------------------------------------------------------------------------------! SUBROUTINE bcm_rrd_global( found ) USE control_parameters, & ONLY: length, restart_string IMPLICIT NONE LOGICAL, INTENT(OUT) :: found found = .TRUE. SELECT CASE ( restart_string(1:length) ) CASE ( 'c_sedimentation' ) READ ( 13 ) c_sedimentation CASE ( 'bulk_cloud_model' ) READ ( 13 ) bulk_cloud_model CASE ( 'cloud_scheme' ) READ ( 13 ) cloud_scheme CASE ( 'cloud_water_sedimentation' ) READ ( 13 ) cloud_water_sedimentation CASE ( 'collision_turbulence' ) READ ( 13 ) collision_turbulence CASE ( 'limiter_sedimentation' ) READ ( 13 ) limiter_sedimentation CASE ( 'nc_const' ) READ ( 13 ) nc_const CASE ( 'precipitation' ) READ ( 13 ) precipitation CASE ( 'ventilation_effect' ) READ ( 13 ) ventilation_effect ! CASE ( 'global_paramter' ) ! READ ( 13 ) global_parameter ! CASE ( 'global_array' ) ! IF ( .NOT. ALLOCATED( global_array ) ) ALLOCATE( global_array(1:10) ) ! READ ( 13 ) global_array CASE DEFAULT found = .FALSE. END SELECT END SUBROUTINE bcm_rrd_global !------------------------------------------------------------------------------! ! Description: ! ------------ !> This routine reads the respective restart data for the microphysics module. !------------------------------------------------------------------------------! SUBROUTINE bcm_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, & nxr_on_file, nynf, nync, nyn_on_file, nysf, & nysc, nys_on_file, tmp_2d, tmp_3d, found ) USE control_parameters USE indices USE pegrid IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: k !< INTEGER(iwp) :: nxlc !< INTEGER(iwp) :: nxlf !< INTEGER(iwp) :: nxl_on_file !< INTEGER(iwp) :: nxrc !< INTEGER(iwp) :: nxrf !< INTEGER(iwp) :: nxr_on_file !< INTEGER(iwp) :: nync !< INTEGER(iwp) :: nynf !< INTEGER(iwp) :: nyn_on_file !< INTEGER(iwp) :: nysc !< INTEGER(iwp) :: nysf !< INTEGER(iwp) :: nys_on_file !< LOGICAL, INTENT(OUT) :: found REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d !< REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !< ! !-- Here the reading of user-defined restart data follows: !-- Sample for user-defined output found = .TRUE. SELECT CASE ( restart_string(1:length) ) CASE ( 'nc' ) IF ( k == 1 ) READ ( 13 ) tmp_3d nc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'nc_av' ) IF ( .NOT. ALLOCATED( nc_av ) ) THEN ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF IF ( k == 1 ) READ ( 13 ) tmp_3d nc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'nr' ) IF ( k == 1 ) READ ( 13 ) tmp_3d nr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'nr_av' ) IF ( .NOT. ALLOCATED( nr_av ) ) THEN ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF IF ( k == 1 ) READ ( 13 ) tmp_3d nr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'precipitation_amount' ) IF ( k == 1 ) READ ( 13 ) tmp_2d precipitation_amount(nysc-nbgp:nync+nbgp, & nxlc-nbgp:nxrc+nbgp) = & tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'prr' ) IF ( .NOT. ALLOCATED( prr ) ) THEN ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF IF ( k == 1 ) READ ( 13 ) tmp_3d prr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'prr_av' ) IF ( .NOT. ALLOCATED( prr_av ) ) THEN ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF IF ( k == 1 ) READ ( 13 ) tmp_3d prr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'qc' ) IF ( k == 1 ) READ ( 13 ) tmp_3d qc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'qc_av' ) IF ( .NOT. ALLOCATED( qc_av ) ) THEN ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF IF ( k == 1 ) READ ( 13 ) tmp_3d qc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'ql' ) IF ( k == 1 ) READ ( 13 ) tmp_3d ql(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'ql_av' ) IF ( .NOT. ALLOCATED( ql_av ) ) THEN ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF IF ( k == 1 ) READ ( 13 ) tmp_3d ql_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'qr' ) IF ( k == 1 ) READ ( 13 ) tmp_3d qr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'qr_av' ) IF ( .NOT. ALLOCATED( qr_av ) ) THEN ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF IF ( k == 1 ) READ ( 13 ) tmp_3d qr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) ! CASE DEFAULT found = .FALSE. END SELECT END SUBROUTINE bcm_rrd_local !------------------------------------------------------------------------------! ! Description: ! ------------ !> This routine writes the respective restart data for the microphysics module. !------------------------------------------------------------------------------! SUBROUTINE bcm_wrd_global IMPLICIT NONE CALL wrd_write_string( 'c_sedimentation' ) WRITE ( 14 ) c_sedimentation CALL wrd_write_string( 'bulk_cloud_model' ) WRITE ( 14 ) bulk_cloud_model CALL wrd_write_string( 'cloud_scheme' ) WRITE ( 14 ) cloud_scheme CALL wrd_write_string( 'cloud_water_sedimentation' ) WRITE ( 14 ) cloud_water_sedimentation CALL wrd_write_string( 'collision_turbulence' ) WRITE ( 14 ) collision_turbulence CALL wrd_write_string( 'limiter_sedimentation' ) WRITE ( 14 ) limiter_sedimentation CALL wrd_write_string( 'nc_const' ) WRITE ( 14 ) nc_const CALL wrd_write_string( 'precipitation' ) WRITE ( 14 ) precipitation CALL wrd_write_string( 'ventilation_effect' ) WRITE ( 14 ) ventilation_effect ! needs preceeding allocation if array ! CALL wrd_write_string( 'global_parameter' ) ! WRITE ( 14 ) global_parameter ! IF ( ALLOCATED( inflow_damping_factor ) ) THEN ! CALL wrd_write_string( 'inflow_damping_factor' ) ! WRITE ( 14 ) inflow_damping_factor ! ENDIF END SUBROUTINE bcm_wrd_global !------------------------------------------------------------------------------! ! Description: ! ------------ !> This routine writes the respective restart data for the microphysics module. !------------------------------------------------------------------------------! SUBROUTINE bcm_wrd_local IMPLICIT NONE IF ( ALLOCATED( prr ) ) THEN CALL wrd_write_string( 'prr' ) WRITE ( 14 ) prr ENDIF IF ( ALLOCATED( prr_av ) ) THEN CALL wrd_write_string( 'prr_av' ) WRITE ( 14 ) prr_av ENDIF IF ( ALLOCATED( precipitation_amount ) ) THEN CALL wrd_write_string( 'precipitation_amount' ) WRITE ( 14 ) precipitation_amount ENDIF CALL wrd_write_string( 'ql' ) WRITE ( 14 ) ql IF ( ALLOCATED( ql_av ) ) THEN CALL wrd_write_string( 'ql_av' ) WRITE ( 14 ) ql_av ENDIF CALL wrd_write_string( 'qc' ) WRITE ( 14 ) qc IF ( ALLOCATED( qc_av ) ) THEN CALL wrd_write_string( 'qc_av' ) WRITE ( 14 ) qc_av ENDIF IF ( microphysics_morrison ) THEN CALL wrd_write_string( 'nc' ) WRITE ( 14 ) nc IF ( ALLOCATED( nc_av ) ) THEN CALL wrd_write_string( 'nc_av' ) WRITE ( 14 ) nc_av ENDIF ENDIF IF ( microphysics_seifert ) THEN CALL wrd_write_string( 'nr' ) WRITE ( 14 ) nr IF ( ALLOCATED( nr_av ) ) THEN CALL wrd_write_string( 'nr_av' ) WRITE ( 14 ) nr_av ENDIF CALL wrd_write_string( 'qr' ) WRITE ( 14 ) qr IF ( ALLOCATED( qr_av ) ) THEN CALL wrd_write_string( 'qr_av' ) WRITE ( 14 ) qr_av ENDIF ENDIF END SUBROUTINE bcm_wrd_local !------------------------------------------------------------------------------! ! Description: ! ------------ !> Control of microphysics for all grid points !------------------------------------------------------------------------------! SUBROUTINE bcm_actions IMPLICIT NONE IF ( large_scale_forcing .AND. lsf_surf ) THEN ! !-- Calculate vertical profile of the hydrostatic pressure (hyp) hyp = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp) d_exner = exner_function_invers(hyp) exner = 1.0_wp / exner_function_invers(hyp) hyrho = ideal_gas_law_rho_pt(hyp, pt_init) ! !-- Compute reference density rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp)) ENDIF ! !-- Compute length of time step IF ( call_microphysics_at_all_substeps ) THEN dt_micro = dt_3d * weight_pres(intermediate_timestep_count) ELSE dt_micro = dt_3d ENDIF ! !-- Reset precipitation rate IF ( intermediate_timestep_count == 1 ) prr = 0.0_wp ! !-- Compute cloud physics IF ( microphysics_kessler ) THEN CALL autoconversion_kessler IF ( cloud_water_sedimentation ) CALL sedimentation_cloud ELSEIF ( microphysics_seifert ) THEN CALL adjust_cloud IF ( microphysics_morrison ) CALL activation IF ( microphysics_morrison ) CALL condensation CALL autoconversion CALL accretion CALL selfcollection_breakup CALL evaporation_rain CALL sedimentation_rain IF ( cloud_water_sedimentation ) CALL sedimentation_cloud ENDIF CALL calc_precipitation_amount END SUBROUTINE bcm_actions !------------------------------------------------------------------------------! ! Description: ! ------------ !> Adjust number of raindrops to avoid nonlinear effects in sedimentation and !> evaporation of rain drops due to too small or too big weights !> of rain drops (Stevens and Seifert, 2008). !------------------------------------------------------------------------------! SUBROUTINE adjust_cloud IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: flag !< flag to indicate first grid level above CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( qr(k,j,i) <= eps_sb ) THEN qr(k,j,i) = 0.0_wp nr(k,j,i) = 0.0_wp ELSE IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag ENDIF ENDIF IF ( microphysics_morrison ) THEN IF ( qc(k,j,i) <= eps_sb ) THEN qc(k,j,i) = 0.0_wp nc(k,j,i) = 0.0_wp ELSE IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin * flag ENDIF ENDIF ENDIF ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' ) END SUBROUTINE adjust_cloud !------------------------------------------------------------------------------! ! Description: ! ------------ !> Calculate number of activated condensation nucleii after simple activation !> scheme of Twomey, 1959. !------------------------------------------------------------------------------! SUBROUTINE activation IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: activ !< REAL(wp) :: afactor !< REAL(wp) :: beta_act !< REAL(wp) :: bfactor !< REAL(wp) :: k_act !< REAL(wp) :: n_act !< REAL(wp) :: n_ccn !< REAL(wp) :: s_0 !< REAL(wp) :: sat_max !< REAL(wp) :: sigma !< REAL(wp) :: sigma_act !< REAL(wp) :: flag !< flag to indicate first grid level above CALL cpu_log( log_point_s(65), 'activation', 'start' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation located !-- in diagnostic_quantities_mod CALL supersaturation ( i, j, k ) ! !-- Prescribe parameters for activation !-- (see: Bott + Trautmann, 2002, Atm. Res., 64) k_act = 0.7_wp activ = 0.0_wp IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN ! !-- Compute the number of activated Aerosols !-- (see: Twomey, 1959, Pure and applied Geophysics, 43) n_act = na_init * sat**k_act ! !-- Compute the number of cloud droplets !-- (see: Morrison + Grabowski, 2007, JAS, 64) ! activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro ! !-- Compute activation rate after Khairoutdinov and Kogan !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128) sat_max = 1.0_wp / 100.0_wp activ = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN & ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / & dt_micro ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN ! !-- Curvature effect (afactor) with surface tension !-- parameterization by Straka (2009) sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp ) afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l ) ! !-- Solute effect (bfactor) bfactor = vanthoff * molecular_weight_of_water * & rho_s / ( molecular_weight_of_solute * rho_l ) ! !-- Prescribe power index that describes the soluble fraction !-- of an aerosol particle (beta) !-- (see: Morrison + Grabowski, 2007, JAS, 64) beta_act = 0.5_wp sigma_act = sigma_bulk**( 1.0_wp + beta_act ) ! !-- Calculate mean geometric supersaturation (s_0) with !-- parameterization by Khvorostyanov and Curry (2006) s_0 = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) * & ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp ! !-- Calculate number of activated CCN as a function of !-- supersaturation and taking Koehler theory into account !-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111) n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( & LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) ) activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp ) ENDIF nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init) ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(65), 'activation', 'stop' ) END SUBROUTINE activation !------------------------------------------------------------------------------! ! Description: ! ------------ !> Calculate condensation rate for cloud water content (after Khairoutdinov and !> Kogan, 2000). !------------------------------------------------------------------------------! SUBROUTINE condensation IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: cond !< REAL(wp) :: cond_max !< REAL(wp) :: dc !< REAL(wp) :: evap !< REAL(wp) :: g_fac !< REAL(wp) :: nc_0 !< REAL(wp) :: temp !< REAL(wp) :: xc !< REAL(wp) :: flag !< flag to indicate first grid level above CALL cpu_log( log_point_s(66), 'condensation', 'start' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation located !-- in diagnostic_quantities_mod CALL supersaturation ( i, j, k ) ! !-- Actual temperature: temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) ) g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * & l_v / ( thermal_conductivity_l * temp ) & + r_v * temp / ( diff_coeff_l * e_s ) & ) ! !-- 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 ! !-- 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) cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i) cond = MIN( cond, cond_max / dt_micro ) 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 = MAX( evap, -qc(k,j,i) / dt_micro ) qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag ENDIF IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin ENDIF ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(66), 'condensation', 'stop' ) END SUBROUTINE condensation !------------------------------------------------------------------------------! ! Description: ! ------------ !> Autoconversion rate (Seifert and Beheng, 2006). !------------------------------------------------------------------------------! SUBROUTINE autoconversion IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: alpha_cc !< REAL(wp) :: autocon !< REAL(wp) :: dissipation !< REAL(wp) :: flag !< flag to mask topography grid points REAL(wp) :: k_au !< REAL(wp) :: l_mix !< REAL(wp) :: nc_auto !< REAL(wp) :: nu_c !< REAL(wp) :: phi_au !< REAL(wp) :: r_cc !< REAL(wp) :: rc !< REAL(wp) :: re_lambda !< REAL(wp) :: sigma_cc !< REAL(wp) :: tau_cloud !< REAL(wp) :: xc !< CALL cpu_log( log_point_s(55), 'autoconversion', 'start' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( microphysics_morrison ) THEN nc_auto = nc(k,j,i) ELSE nc_auto = nc_const ENDIF IF ( qc(k,j,i) > eps_sb .AND. nc_auto > eps_mr ) THEN k_au = k_cc / ( 20.0_wp * x0 ) ! !-- Intern time scale of coagulation (Seifert and Beheng, 2006): !-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )) tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + & qc(k,j,i) ), 0.0_wp ) ! !-- Universal function for autoconversion process !-- (Seifert and Beheng, 2006): phi_au = 600.0_wp * tau_cloud**0.68_wp * & ( 1.0_wp - tau_cloud**0.68_wp )**3 ! !-- Shape parameter of gamma distribution (Geoffroy et al., 2010): !-- (Use constant nu_c = 1.0_wp instead?) nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp ) ! !-- Mean weight of cloud droplets: xc = MAX( hyrho(k) * qc(k,j,i) / nc_auto, xcmin) ! !-- Parameterized turbulence effects on autoconversion (Seifert, !-- Nuijens and Stevens, 2010) IF ( collision_turbulence ) THEN ! !-- Weight averaged radius of cloud droplets: rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp ) alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c ) r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c ) sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c ) ! !-- Mixing length (neglecting distance to ground and !-- stratification) l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp ) ! !-- Limit dissipation rate according to Seifert, Nuijens and !-- Stevens (2010) dissipation = MIN( 0.06_wp, diss(k,j,i) ) ! !-- Compute Taylor-microscale Reynolds number: re_lambda = 6.0_wp / 11.0_wp * & ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) * & SQRT( 15.0_wp / kin_vis_air ) * & dissipation**( 1.0_wp / 6.0_wp ) ! !-- The factor of 1.0E4 is needed to convert the dissipation !-- rate from m2 s-3 to cm2 s-3. k_au = k_au * ( 1.0_wp + & dissipation * 1.0E4_wp * & ( re_lambda * 1.0E-3_wp )**0.25_wp * & ( alpha_cc * EXP( -1.0_wp * ( ( rc - & r_cc ) / & sigma_cc )**2 & ) + beta_cc & ) & ) ENDIF ! !-- Autoconversion rate (Seifert and Beheng, 2006): autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / & ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 * & ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * & rho_surface autocon = MIN( autocon, qc(k,j,i) / dt_micro ) qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro & * flag IF ( microphysics_morrison ) THEN nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp * & autocon / x0 * hyrho(k) * dt_micro * flag ) ENDIF ENDIF ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' ) END SUBROUTINE autoconversion !------------------------------------------------------------------------------! ! Description: ! ------------ !> Autoconversion process (Kessler, 1969). !------------------------------------------------------------------------------! SUBROUTINE autoconversion_kessler IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< INTEGER(iwp) :: k_wall !< topgraphy top index REAL(wp) :: dqdt_precip !< REAL(wp) :: flag !< flag to mask topography grid points DO i = nxlg, nxrg DO j = nysg, nyng ! !-- Determine vertical index of topography top k_wall = get_topography_top_index_ji( j, i, 's' ) DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( qc(k,j,i) > ql_crit ) THEN dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit ) ELSE dqdt_precip = 0.0_wp ENDIF qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro * flag pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * & d_exner(k) * flag ! !-- Compute the rain rate (stored on surface grid point) prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag ENDDO ENDDO ENDDO END SUBROUTINE autoconversion_kessler !------------------------------------------------------------------------------! ! Description: ! ------------ !> Accretion rate (Seifert and Beheng, 2006). !------------------------------------------------------------------------------! SUBROUTINE accretion IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: accr !< REAL(wp) :: flag !< flag to mask topography grid points REAL(wp) :: k_cr !< REAL(wp) :: nc_accr !< REAL(wp) :: phi_ac !< REAL(wp) :: tau_cloud !< REAL(wp) :: xc !< CALL cpu_log( log_point_s(56), 'accretion', 'start' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( microphysics_morrison ) THEN nc_accr = nc(k,j,i) ELSE nc_accr = nc_const ENDIF IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) & .AND. ( nc_accr > eps_mr ) ) THEN ! !-- Intern time scale of coagulation (Seifert and Beheng, 2006): tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ) ! !-- Universal function for accretion process (Seifert and !-- Beheng, 2001): phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4 ! !-- Mean weight of cloud drops xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin) ! !-- Parameterized turbulence effects on autoconversion (Seifert, !-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to !-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3. IF ( collision_turbulence ) THEN k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * & MIN( 600.0_wp, & diss(k,j,i) * 1.0E4_wp )**0.25_wp & ) ELSE k_cr = k_cr0 ENDIF ! !-- Accretion rate (Seifert and Beheng, 2006): accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * & SQRT( rho_surface * hyrho(k) ) accr = MIN( accr, qc(k,j,i) / dt_micro ) qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag IF ( microphysics_morrison ) THEN nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), & accr / xc * hyrho(k) * dt_micro * flag) ENDIF ENDIF ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(56), 'accretion', 'stop' ) END SUBROUTINE accretion !------------------------------------------------------------------------------! ! Description: ! ------------ !> Collisional breakup rate (Seifert, 2008). !------------------------------------------------------------------------------! SUBROUTINE selfcollection_breakup IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: breakup !< REAL(wp) :: dr !< REAL(wp) :: flag !< flag to mask topography grid points REAL(wp) :: phi_br !< REAL(wp) :: selfcoll !< CALL cpu_log( log_point_s(57), 'selfcollection', 'start' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( qr(k,j,i) > eps_sb ) THEN ! !-- Selfcollection rate (Seifert and Beheng, 2001): selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * & SQRT( hyrho(k) * rho_surface ) ! !-- Weight averaged diameter of rain drops: dr = ( hyrho(k) * qr(k,j,i) / & nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp ) ! !-- Collisional breakup rate (Seifert, 2008): IF ( dr >= 0.3E-3_wp ) THEN phi_br = k_br * ( dr - 1.1E-3_wp ) breakup = selfcoll * ( phi_br + 1.0_wp ) ELSE breakup = 0.0_wp ENDIF selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro ) nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag ENDIF ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' ) END SUBROUTINE selfcollection_breakup !------------------------------------------------------------------------------! ! Description: ! ------------ !> Evaporation of precipitable water. Condensation is neglected for !> precipitable water. !------------------------------------------------------------------------------! SUBROUTINE evaporation_rain IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: dr !< REAL(wp) :: evap !< REAL(wp) :: evap_nr !< REAL(wp) :: f_vent !< REAL(wp) :: flag !< flag to mask topography grid points REAL(wp) :: g_evap !< REAL(wp) :: lambda_r !< REAL(wp) :: mu_r !< REAL(wp) :: mu_r_2 !< REAL(wp) :: mu_r_5d2 !< REAL(wp) :: nr_0 !< REAL(wp) :: temp !< REAL(wp) :: xr !< CALL cpu_log( log_point_s(58), 'evaporation', 'start' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( qr(k,j,i) > eps_sb ) THEN ! !-- Call calculation of supersaturation located !-- in diagnostic_quantities_mod CALL supersaturation ( i, j, k ) ! !-- Evaporation needs only to be calculated in subsaturated regions IF ( sat < 0.0_wp ) THEN ! !-- Actual temperature: temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) ) g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * & l_v / ( thermal_conductivity_l * temp ) & + r_v * temp / ( diff_coeff_l * e_s ) & ) ! !-- Mean weight of rain drops xr = hyrho(k) * qr(k,j,i) / nr(k,j,i) ! !-- Weight averaged diameter of rain drops: dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp ) ! !-- Compute ventilation factor and intercept parameter !-- (Seifert and Beheng, 2006; Seifert, 2008): IF ( ventilation_effect ) THEN ! !-- Shape parameter of gamma distribution (Milbrandt and Yau, !-- 2005; Stevens and Seifert, 2008): mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * & ( dr - 1.4E-3_wp ) ) ) ! !-- Slope parameter of gamma distribution (Seifert, 2008): lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & ( mu_r + 1.0_wp ) & )**( 1.0_wp / 3.0_wp ) / dr mu_r_2 = mu_r + 2.0_wp mu_r_5d2 = mu_r + 2.5_wp f_vent = a_vent * gamm( mu_r_2 ) * & lambda_r**( -mu_r_2 ) + b_vent * & schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *& gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) * & ( 1.0_wp - & 0.5_wp * ( b_term / a_term ) * & ( lambda_r / ( c_term + lambda_r ) & )**mu_r_5d2 - & 0.125_wp * ( b_term / a_term )**2 * & ( lambda_r / ( 2.0_wp * c_term + lambda_r ) & )**mu_r_5d2 - & 0.0625_wp * ( b_term / a_term )**3 * & ( lambda_r / ( 3.0_wp * c_term + lambda_r ) & )**mu_r_5d2 - & 0.0390625_wp * ( b_term / a_term )**4 * & ( lambda_r / ( 4.0_wp * c_term + lambda_r ) & )**mu_r_5d2 & ) nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / & gamm( mu_r + 1.0_wp ) ELSE f_vent = 1.0_wp nr_0 = nr(k,j,i) * dr ENDIF ! !-- Evaporation rate of rain water content (Seifert and !-- Beheng, 2006): evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / & hyrho(k) evap = MAX( evap, -qr(k,j,i) / dt_micro ) evap_nr = MAX( c_evap * evap / xr * hyrho(k), & -nr(k,j,i) / dt_micro ) qr(k,j,i) = qr(k,j,i) + evap * dt_micro * flag nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag ENDIF ENDIF ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(58), 'evaporation', 'stop' ) END SUBROUTINE evaporation_rain !------------------------------------------------------------------------------! ! Description: ! ------------ !> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR). !------------------------------------------------------------------------------! SUBROUTINE sedimentation_cloud IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: flag !< flag to mask topography grid points REAL(wp) :: nc_sedi !< REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !< REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !< CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' ) sed_qc(nzt+1) = 0.0_wp sed_nc(nzt+1) = 0.0_wp DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzt, nzb+1, -1 ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( microphysics_morrison ) THEN nc_sedi = nc(k,j,i) ELSE nc_sedi = nc_const ENDIF ! !-- Sedimentation fluxes for number concentration are only calculated !-- for cloud_scheme = 'morrison' IF ( microphysics_morrison ) THEN IF ( qc(k,j,i) > eps_sb .AND. nc(k,j,i) > eps_mr ) THEN sed_nc(k) = sed_qc_const * & ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) * & ( nc(k,j,i) )**( 1.0_wp / 3.0_wp ) ELSE sed_nc(k) = 0.0_wp ENDIF sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * & nc(k,j,i) / dt_micro + sed_nc(k+1) & ) * flag nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) * & ddzu(k+1) / hyrho(k) * dt_micro * flag ENDIF IF ( qc(k,j,i) > eps_sb .AND. nc_sedi > eps_mr ) THEN sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) * & ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * & flag ELSE sed_qc(k) = 0.0_wp ENDIF sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / & dt_micro + sed_qc(k+1) & ) * flag q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * & ddzu(k+1) / hyrho(k) * dt_micro * flag qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * & ddzu(k+1) / hyrho(k) * dt_micro * flag pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * & ddzu(k+1) / hyrho(k) * lv_d_cp * & d_exner(k) * dt_micro * flag ! !-- Compute the precipitation rate due to cloud (fog) droplets IF ( call_microphysics_at_all_substeps ) THEN prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) & * weight_substep(intermediate_timestep_count) & * flag ELSE prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag ENDIF ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' ) END SUBROUTINE sedimentation_cloud !------------------------------------------------------------------------------! ! Description: ! ------------ !> Computation of sedimentation flux. Implementation according to Stevens !> and Seifert (2008). Code is based on UCLA-LES. !------------------------------------------------------------------------------! SUBROUTINE sedimentation_rain IMPLICIT NONE INTEGER(iwp) :: i !< running index x direction INTEGER(iwp) :: j !< running index y direction INTEGER(iwp) :: k !< running index z direction INTEGER(iwp) :: k_run !< INTEGER(iwp) :: m !< running index surface elements INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint REAL(wp) :: c_run !< REAL(wp) :: d_max !< REAL(wp) :: d_mean !< REAL(wp) :: d_min !< REAL(wp) :: dr !< REAL(wp) :: flux !< REAL(wp) :: flag !< flag to mask topography grid points REAL(wp) :: lambda_r !< REAL(wp) :: mu_r !< REAL(wp) :: z_run !< REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !< REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !< REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !< REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !< REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !< REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !< REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !< REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !< CALL cpu_log( log_point_s(60), 'sed_rain', 'start' ) ! !-- Compute velocities DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( qr(k,j,i) > eps_sb ) THEN ! !-- Weight averaged diameter of rain drops: dr = ( hyrho(k) * qr(k,j,i) / & nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp ) ! !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; !-- Stevens and Seifert, 2008): mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * & ( dr - 1.4E-3_wp ) ) ) ! !-- Slope parameter of gamma distribution (Seifert, 2008): lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, & a_term - b_term * ( 1.0_wp + & c_term / & lambda_r )**( -1.0_wp * & ( mu_r + 1.0_wp ) ) & ) & ) * flag w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, & a_term - b_term * ( 1.0_wp + & c_term / & lambda_r )**( -1.0_wp * & ( mu_r + 4.0_wp ) ) & ) & ) * flag ELSE w_nr(k) = 0.0_wp w_qr(k) = 0.0_wp ENDIF ENDDO ! !-- Adjust boundary values using surface data type. !-- Upward-facing surf_s = bc_h(0)%start_index(j,i) surf_e = bc_h(0)%end_index(j,i) DO m = surf_s, surf_e k = bc_h(0)%k(m) w_nr(k-1) = w_nr(k) w_qr(k-1) = w_qr(k) ENDDO ! !-- Downward-facing surf_s = bc_h(1)%start_index(j,i) surf_e = bc_h(1)%end_index(j,i) DO m = surf_s, surf_e k = bc_h(1)%k(m) w_nr(k+1) = w_nr(k) w_qr(k+1) = w_qr(k) ENDDO ! !-- Model top boundary value w_nr(nzt+1) = 0.0_wp w_qr(nzt+1) = 0.0_wp ! !-- Compute Courant number DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) c_nr(k) = 0.25_wp * ( w_nr(k-1) + & 2.0_wp * w_nr(k) + w_nr(k+1) ) * & dt_micro * ddzu(k) * flag c_qr(k) = 0.25_wp * ( w_qr(k-1) + & 2.0_wp * w_qr(k) + w_qr(k+1) ) * & dt_micro * ddzu(k) * flag ENDDO ! !-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977): IF ( limiter_sedimentation ) THEN DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) ) d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i) qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, & 2.0_wp * d_max, & ABS( d_mean ) ) & * flag d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) ) d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i) nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, & 2.0_wp * d_max, & ABS( d_mean ) ) ENDDO ELSE nr_slope = 0.0_wp qr_slope = 0.0_wp ENDIF sed_nr(nzt+1) = 0.0_wp sed_qr(nzt+1) = 0.0_wp ! !-- Compute sedimentation flux DO k = nzt, nzb+1, -1 ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ! !-- Sum up all rain drop number densities which contribute to the flux !-- through k-1/2 flux = 0.0_wp z_run = 0.0_wp ! height above z(k) k_run = k c_run = MIN( 1.0_wp, c_nr(k) ) DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt ) flux = flux + hyrho(k_run) * & ( nr(k_run,j,i) + nr_slope(k_run) * & ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run) & * flag z_run = z_run + dzu(k_run) * flag k_run = k_run + 1 * flag c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) & * flag ENDDO ! !-- It is not allowed to sediment more rain drop number density than !-- available flux = MIN( flux, & hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * & dt_micro & ) sed_nr(k) = flux / dt_micro * flag nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * & ddzu(k+1) / hyrho(k) * dt_micro * flag ! !-- Sum up all rain water content which contributes to the flux !-- through k-1/2 flux = 0.0_wp z_run = 0.0_wp ! height above z(k) k_run = k c_run = MIN( 1.0_wp, c_qr(k) ) DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt ) flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + & qr_slope(k_run) * ( 1.0_wp - c_run ) * & 0.5_wp ) * c_run * dzu(k_run) * flag z_run = z_run + dzu(k_run) * flag k_run = k_run + 1 * flag c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) & * flag ENDDO ! !-- It is not allowed to sediment more rain water content than !-- available flux = MIN( flux, & hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * & dt_micro & ) sed_qr(k) = flux / dt_micro * flag qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * & ddzu(k+1) / hyrho(k) * dt_micro * flag q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * & ddzu(k+1) / hyrho(k) * dt_micro * flag pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * & ddzu(k+1) / hyrho(k) * lv_d_cp * & d_exner(k) * dt_micro * flag ! !-- Compute the rain rate IF ( call_microphysics_at_all_substeps ) THEN prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) & * weight_substep(intermediate_timestep_count) & * flag ELSE prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag ENDIF ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' ) END SUBROUTINE sedimentation_rain !------------------------------------------------------------------------------! ! Description: ! ------------ !> Computation of the precipitation amount due to gravitational settling of !> rain and cloud (fog) droplets !------------------------------------------------------------------------------! SUBROUTINE calc_precipitation_amount IMPLICIT NONE INTEGER(iwp) :: i !< running index x direction INTEGER(iwp) :: j !< running index y direction INTEGER(iwp) :: k !< running index y direction INTEGER(iwp) :: m !< running index surface elements IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.& ( .NOT. call_microphysics_at_all_substeps .OR. & intermediate_timestep_count == intermediate_timestep_count_max ) ) & THEN ! !-- Run over all upward-facing surface elements, i.e. non-natural, !-- natural and urban DO m = 1, bc_h(0)%ns i = bc_h(0)%i(m) j = bc_h(0)%j(m) k = bc_h(0)%k(m) precipitation_amount(j,i) = precipitation_amount(j,i) + & prr(k,j,i) * hyrho(k) * dt_3d ENDDO ENDIF END SUBROUTINE calc_precipitation_amount !------------------------------------------------------------------------------! ! Description: ! ------------ !> Control of microphysics for grid points i,j !------------------------------------------------------------------------------! SUBROUTINE bcm_actions_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< IF ( large_scale_forcing .AND. lsf_surf ) THEN ! !-- Calculate vertical profile of the hydrostatic pressure (hyp) hyp = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp) d_exner = exner_function_invers(hyp) exner = 1.0_wp / exner_function_invers(hyp) hyrho = ideal_gas_law_rho_pt(hyp, pt_init) ! !-- Compute reference density rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp)) ENDIF ! !-- Compute length of time step IF ( call_microphysics_at_all_substeps ) THEN dt_micro = dt_3d * weight_pres(intermediate_timestep_count) ELSE dt_micro = dt_3d ENDIF ! !-- Reset precipitation rate IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp ! !-- Compute cloud physics IF( microphysics_kessler ) THEN CALL autoconversion_kessler( i,j ) IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j ) ELSEIF ( microphysics_seifert ) THEN CALL adjust_cloud( i,j ) IF ( microphysics_morrison ) CALL activation( i,j ) IF ( microphysics_morrison ) CALL condensation( i,j ) CALL autoconversion( i,j ) CALL accretion( i,j ) CALL selfcollection_breakup( i,j ) CALL evaporation_rain( i,j ) CALL sedimentation_rain( i,j ) IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j ) ENDIF CALL calc_precipitation_amount( i,j ) END SUBROUTINE bcm_actions_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Adjust number of raindrops to avoid nonlinear effects in !> sedimentation and evaporation of rain drops due to too small or !> too big weights of rain drops (Stevens and Seifert, 2008). !> The same procedure is applied to cloud droplets if they are determined !> prognostically. Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE adjust_cloud_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: flag !< flag to indicate first grid level above surface DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( qr(k,j,i) <= eps_sb ) THEN qr(k,j,i) = 0.0_wp nr(k,j,i) = 0.0_wp ELSE ! !-- Adjust number of raindrops to avoid nonlinear effects in !-- sedimentation and evaporation of rain drops due to too small or !-- too big weights of rain drops (Stevens and Seifert, 2008). IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag ENDIF ENDIF IF ( microphysics_morrison ) THEN IF ( qc(k,j,i) <= eps_sb ) THEN qc(k,j,i) = 0.0_wp nc(k,j,i) = 0.0_wp ELSE IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN nc(k,j,i) = qc(k,j,i) * hyrho(k) / xamin * flag ENDIF ENDIF ENDIF ENDDO END SUBROUTINE adjust_cloud_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Calculate number of activated condensation nucleii after simple activation !> scheme of Twomey, 1959. !------------------------------------------------------------------------------! SUBROUTINE activation_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: activ !< REAL(wp) :: afactor !< REAL(wp) :: alpha !< REAL(wp) :: beta_act !< REAL(wp) :: bfactor !< REAL(wp) :: e_s !< REAL(wp) :: flag !< flag to indicate first grid level above surface REAL(wp) :: k_act !< REAL(wp) :: n_act !< REAL(wp) :: n_ccn !< REAL(wp) :: q_s !< REAL(wp) :: s_0 !< REAL(wp) :: sat !< REAL(wp) :: sat_max !< REAL(wp) :: sigma !< REAL(wp) :: sigma_act !< REAL(wp) :: t_int !< REAL(wp) :: t_l !< DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ! !-- Actual liquid water temperature: t_l = exner(k) * pt(k,j,i) ! !-- Calculate actual temperature t_int = pt(k,j,i) * exner_function( hyp(k) ) ! !-- Saturation vapor pressure at t_l: e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / & ( t_l - 35.86_wp ) & ) ! !-- Computation of saturation mixing ratio: q_s = rd_d_rv * e_s / ( hyp(k) - e_s ) alpha = rd_d_rv * lv_d_rd * lv_d_cp / ( t_l * t_l ) q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / & ( 1.0_wp + alpha * q_s ) !-- Supersaturation: sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp ! !-- Prescribe parameters for activation !-- (see: Bott + Trautmann, 2002, Atm. Res., 64) k_act = 0.7_wp activ = 0.0_wp IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN ! !-- Compute the number of activated Aerosols !-- (see: Twomey, 1959, Pure and applied Geophysics, 43) n_act = na_init * sat**k_act ! !-- Compute the number of cloud droplets !-- (see: Morrison + Grabowski, 2007, JAS, 64) ! activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro ! !-- Compute activation rate after Khairoutdinov and Kogan !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128) sat_max = 0.8_wp / 100.0_wp activ = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN & ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / & dt_micro nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init) ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk ) THEN ! !-- Curvature effect (afactor) with surface tension !-- parameterization by Straka (2009) sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp ) afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int ) ! !-- Solute effect (bfactor) bfactor = vanthoff * molecular_weight_of_water * & rho_s / ( molecular_weight_of_solute * rho_l ) ! !-- Prescribe power index that describes the soluble fraction !-- of an aerosol particle (beta). !-- (see: Morrison + Grabowski, 2007, JAS, 64) beta_act = 0.5_wp sigma_act = sigma_bulk**( 1.0_wp + beta_act ) ! !-- Calculate mean geometric supersaturation (s_0) with !-- parameterization by Khvorostyanov and Curry (2006) s_0 = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) * & ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp ! !-- Calculate number of activated CCN as a function of !-- supersaturation and taking Koehler theory into account !-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111) n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( & LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) ) activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp ) nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init) ENDIF ENDDO END SUBROUTINE activation_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Calculate condensation rate for cloud water content (after Khairoutdinov and !> Kogan, 2000). !------------------------------------------------------------------------------! SUBROUTINE condensation_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: alpha !< REAL(wp) :: cond !< REAL(wp) :: cond_max !< REAL(wp) :: dc !< REAL(wp) :: e_s !< REAL(wp) :: evap !< REAL(wp) :: flag !< flag to indicate first grid level above surface REAL(wp) :: g_fac !< REAL(wp) :: nc_0 !< REAL(wp) :: q_s !< REAL(wp) :: sat !< REAL(wp) :: t_l !< REAL(wp) :: temp !< REAL(wp) :: xc !< DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ! !-- Actual liquid water temperature: t_l = exner(k) * pt(k,j,i) ! !-- Saturation vapor pressure at t_l: e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / & ( t_l - 35.86_wp ) & ) ! !-- Computation of saturation mixing ratio: q_s = rd_d_rv * e_s / ( hyp(k) - e_s ) alpha = rd_d_rv * lv_d_rd * lv_d_cp / ( t_l * t_l ) q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / & ( 1.0_wp + alpha * q_s ) !-- Supersaturation: sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp ! !-- Actual temperature: temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) ) g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * & l_v / ( thermal_conductivity_l * temp ) & + r_v * temp / ( diff_coeff_l * e_s ) & ) ! !-- 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 ! !-- 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) cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i) cond = MIN( cond, cond_max / dt_micro ) 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 = MAX( evap, -qc(k,j,i) / dt_micro ) qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag ENDIF ENDDO END SUBROUTINE condensation_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE autoconversion_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: alpha_cc !< REAL(wp) :: autocon !< REAL(wp) :: dissipation !< REAL(wp) :: flag !< flag to indicate first grid level above surface REAL(wp) :: k_au !< REAL(wp) :: l_mix !< REAL(wp) :: nc_auto !< REAL(wp) :: nu_c !< REAL(wp) :: phi_au !< REAL(wp) :: r_cc !< REAL(wp) :: rc !< REAL(wp) :: re_lambda !< REAL(wp) :: sigma_cc !< REAL(wp) :: tau_cloud !< REAL(wp) :: xc !< DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) nc_auto = MERGE ( nc(k,j,i), nc_const, microphysics_morrison ) IF ( qc(k,j,i) > eps_sb .AND. nc_auto > eps_mr ) THEN k_au = k_cc / ( 20.0_wp * x0 ) ! !-- Intern time scale of coagulation (Seifert and Beheng, 2006): !-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )) tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ), & 0.0_wp ) ! !-- Universal function for autoconversion process !-- (Seifert and Beheng, 2006): phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3 ! !-- Shape parameter of gamma distribution (Geoffroy et al., 2010): !-- (Use constant nu_c = 1.0_wp instead?) nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp ) ! !-- Mean weight of cloud droplets: xc = hyrho(k) * qc(k,j,i) / nc_auto ! !-- Parameterized turbulence effects on autoconversion (Seifert, !-- Nuijens and Stevens, 2010) IF ( collision_turbulence ) THEN ! !-- Weight averaged radius of cloud droplets: rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp ) alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c ) r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c ) sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c ) ! !-- Mixing length (neglecting distance to ground and stratification) l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp ) ! !-- Limit dissipation rate according to Seifert, Nuijens and !-- Stevens (2010) dissipation = MIN( 0.06_wp, diss(k,j,i) ) ! !-- Compute Taylor-microscale Reynolds number: re_lambda = 6.0_wp / 11.0_wp * & ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) * & SQRT( 15.0_wp / kin_vis_air ) * & dissipation**( 1.0_wp / 6.0_wp ) ! !-- The factor of 1.0E4 is needed to convert the dissipation rate !-- from m2 s-3 to cm2 s-3. k_au = k_au * ( 1.0_wp + & dissipation * 1.0E4_wp * & ( re_lambda * 1.0E-3_wp )**0.25_wp * & ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / & sigma_cc )**2 & ) + beta_cc & ) & ) ENDIF ! !-- Autoconversion rate (Seifert and Beheng, 2006): autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / & ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 * & ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * & rho_surface autocon = MIN( autocon, qc(k,j,i) / dt_micro ) qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag IF ( microphysics_morrison ) THEN nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp * & autocon / x0 * hyrho(k) * dt_micro * flag ) ENDIF ENDIF ENDDO END SUBROUTINE autoconversion_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Autoconversion process (Kessler, 1969). !------------------------------------------------------------------------------! SUBROUTINE autoconversion_kessler_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< INTEGER(iwp) :: k_wall !< topography top index REAL(wp) :: dqdt_precip !< REAL(wp) :: flag !< flag to indicate first grid level above surface ! !-- Determine vertical index of topography top k_wall = get_topography_top_index_ji( j, i, 's' ) DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( qc(k,j,i) > ql_crit ) THEN dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit ) ELSE dqdt_precip = 0.0_wp ENDIF qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro * flag pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k) & * flag ! !-- Compute the rain rate (stored on surface grid point) prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag ENDDO END SUBROUTINE autoconversion_kessler_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE accretion_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: accr !< REAL(wp) :: flag !< flag to indicate first grid level above surface REAL(wp) :: k_cr !< REAL(wp) :: nc_accr !< REAL(wp) :: phi_ac !< REAL(wp) :: tau_cloud !< REAL(wp) :: xc !< DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) nc_accr = MERGE ( nc(k,j,i), nc_const, microphysics_morrison ) IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) .AND. & ( nc_accr > eps_mr ) ) THEN ! !-- Intern time scale of coagulation (Seifert and Beheng, 2006): tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ) ! !-- Universal function for accretion process !-- (Seifert and Beheng, 2001): phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4 ! !-- Mean weight of cloud drops xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin) ! !-- Parameterized turbulence effects on autoconversion (Seifert, !-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to !-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3. IF ( collision_turbulence ) THEN k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * & MIN( 600.0_wp, & diss(k,j,i) * 1.0E4_wp )**0.25_wp & ) ELSE k_cr = k_cr0 ENDIF ! !-- Accretion rate (Seifert and Beheng, 2006): accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * & SQRT( rho_surface * hyrho(k) ) accr = MIN( accr, qc(k,j,i) / dt_micro ) qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag IF ( microphysics_morrison ) THEN nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc * & hyrho(k) * dt_micro * flag & ) ENDIF ENDIF ENDDO END SUBROUTINE accretion_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Collisional breakup rate (Seifert, 2008). Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE selfcollection_breakup_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: breakup !< REAL(wp) :: dr !< REAL(wp) :: flag !< flag to indicate first grid level above surface REAL(wp) :: phi_br !< REAL(wp) :: selfcoll !< DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( qr(k,j,i) > eps_sb ) THEN ! !-- Selfcollection rate (Seifert and Beheng, 2001): selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * SQRT( hyrho(k) * rho_surface ) ! !-- Weight averaged diameter of rain drops: dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp ) ! !-- Collisional breakup rate (Seifert, 2008): IF ( dr >= 0.3E-3_wp ) THEN phi_br = k_br * ( dr - 1.1E-3_wp ) breakup = selfcoll * ( phi_br + 1.0_wp ) ELSE breakup = 0.0_wp ENDIF selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro ) nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag ENDIF ENDDO END SUBROUTINE selfcollection_breakup_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Evaporation of precipitable water. Condensation is neglected for !> precipitable water. Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE evaporation_rain_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: alpha !< REAL(wp) :: dr !< REAL(wp) :: e_s !< REAL(wp) :: evap !< REAL(wp) :: evap_nr !< REAL(wp) :: f_vent !< REAL(wp) :: flag !< flag to indicate first grid level above surface REAL(wp) :: g_evap !< REAL(wp) :: lambda_r !< REAL(wp) :: mu_r !< REAL(wp) :: mu_r_2 !< REAL(wp) :: mu_r_5d2 !< REAL(wp) :: nr_0 !< REAL(wp) :: q_s !< REAL(wp) :: sat !< REAL(wp) :: t_l !< REAL(wp) :: temp !< REAL(wp) :: xr !< DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( qr(k,j,i) > eps_sb ) THEN ! !-- Actual liquid water temperature: t_l = exner(k) * pt(k,j,i) ! !-- Saturation vapor pressure at t_l: e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / & ( t_l - 35.86_wp ) & ) ! !-- Computation of saturation mixing ratio: q_s = rd_d_rv * e_s / ( hyp(k) - e_s ) alpha = rd_d_rv * lv_d_rd * lv_d_cp / ( t_l * t_l ) q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / ( 1.0_wp + alpha * q_s ) ! !-- Supersaturation: sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp ! !-- Evaporation needs only to be calculated in subsaturated regions IF ( sat < 0.0_wp ) THEN ! !-- Actual temperature: temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) ) g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / & ( thermal_conductivity_l * temp ) + & r_v * temp / ( diff_coeff_l * e_s ) & ) ! !-- Mean weight of rain drops xr = hyrho(k) * qr(k,j,i) / nr(k,j,i) ! !-- Weight averaged diameter of rain drops: dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp ) ! !-- Compute ventilation factor and intercept parameter !-- (Seifert and Beheng, 2006; Seifert, 2008): IF ( ventilation_effect ) THEN ! !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; !-- Stevens and Seifert, 2008): mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) ) ! !-- Slope parameter of gamma distribution (Seifert, 2008): lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & ( mu_r + 1.0_wp ) & )**( 1.0_wp / 3.0_wp ) / dr mu_r_2 = mu_r + 2.0_wp mu_r_5d2 = mu_r + 2.5_wp f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + & b_vent * schmidt_p_1d3 * & SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * & lambda_r**( -mu_r_5d2 ) * & ( 1.0_wp - & 0.5_wp * ( b_term / a_term ) * & ( lambda_r / ( c_term + lambda_r ) & )**mu_r_5d2 - & 0.125_wp * ( b_term / a_term )**2 * & ( lambda_r / ( 2.0_wp * c_term + lambda_r ) & )**mu_r_5d2 - & 0.0625_wp * ( b_term / a_term )**3 * & ( lambda_r / ( 3.0_wp * c_term + lambda_r ) & )**mu_r_5d2 - & 0.0390625_wp * ( b_term / a_term )**4 * & ( lambda_r / ( 4.0_wp * c_term + lambda_r ) & )**mu_r_5d2 & ) nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / & gamm( mu_r + 1.0_wp ) ELSE f_vent = 1.0_wp nr_0 = nr(k,j,i) * dr ENDIF ! !-- Evaporation rate of rain water content (Seifert and Beheng, 2006): evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k) evap = MAX( evap, -qr(k,j,i) / dt_micro ) evap_nr = MAX( c_evap * evap / xr * hyrho(k), & -nr(k,j,i) / dt_micro ) qr(k,j,i) = qr(k,j,i) + evap * dt_micro * flag nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag ENDIF ENDIF ENDDO END SUBROUTINE evaporation_rain_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR). !> Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE sedimentation_cloud_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: flag !< flag to indicate first grid level above surface REAL(wp) :: nc_sedi !< REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !< REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !< sed_qc(nzt+1) = 0.0_wp sed_nc(nzt+1) = 0.0_wp DO k = nzt, nzb+1, -1 ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) nc_sedi = MERGE( nc(k,j,i), nc_const, microphysics_morrison ) ! !-- Sedimentation fluxes for number concentration are only calculated !-- for cloud_scheme = 'morrison' IF ( microphysics_morrison ) THEN IF ( qc(k,j,i) > eps_sb .AND. nc(k,j,i) > eps_mr ) THEN sed_nc(k) = sed_qc_const * & ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) * & ( nc(k,j,i) )**( 1.0_wp / 3.0_wp ) ELSE sed_nc(k) = 0.0_wp ENDIF sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * & nc(k,j,i) / dt_micro + sed_nc(k+1) & ) * flag nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) * & ddzu(k+1) / hyrho(k) * dt_micro * flag ENDIF IF ( qc(k,j,i) > eps_sb .AND. nc_sedi > eps_mr ) THEN sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) * & ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * flag ELSE sed_qc(k) = 0.0_wp ENDIF sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / & dt_micro + sed_qc(k+1) & ) * flag q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & hyrho(k) * dt_micro * flag qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & hyrho(k) * dt_micro * flag pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & hyrho(k) * lv_d_cp * d_exner(k) * dt_micro & * flag ! !-- Compute the precipitation rate of cloud (fog) droplets IF ( call_microphysics_at_all_substeps ) THEN prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * & weight_substep(intermediate_timestep_count) * flag ELSE prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag ENDIF ENDDO END SUBROUTINE sedimentation_cloud_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Computation of sedimentation flux. Implementation according to Stevens !> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE sedimentation_rain_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< running index x direction INTEGER(iwp) :: j !< running index y direction INTEGER(iwp) :: k !< running index z direction INTEGER(iwp) :: k_run !< INTEGER(iwp) :: m !< running index surface elements INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint REAL(wp) :: c_run !< REAL(wp) :: d_max !< REAL(wp) :: d_mean !< REAL(wp) :: d_min !< REAL(wp) :: dr !< REAL(wp) :: flux !< REAL(wp) :: flag !< flag to indicate first grid level above surface REAL(wp) :: lambda_r !< REAL(wp) :: mu_r !< REAL(wp) :: z_run !< REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !< REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !< REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !< REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !< REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !< REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !< REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !< REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !< ! !-- Compute velocities DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) IF ( qr(k,j,i) > eps_sb ) THEN ! !-- Weight averaged diameter of rain drops: dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp ) ! !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; !-- Stevens and Seifert, 2008): mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) ) ! !-- Slope parameter of gamma distribution (Seifert, 2008): lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, & a_term - b_term * ( 1.0_wp + & c_term / lambda_r )**( -1.0_wp * & ( mu_r + 1.0_wp ) ) & ) & ) * flag w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, & a_term - b_term * ( 1.0_wp + & c_term / lambda_r )**( -1.0_wp * & ( mu_r + 4.0_wp ) ) & ) & ) * flag ELSE w_nr(k) = 0.0_wp w_qr(k) = 0.0_wp ENDIF ENDDO ! !-- Adjust boundary values using surface data type. !-- Upward facing non-natural surf_s = bc_h(0)%start_index(j,i) surf_e = bc_h(0)%end_index(j,i) DO m = surf_s, surf_e k = bc_h(0)%k(m) w_nr(k-1) = w_nr(k) w_qr(k-1) = w_qr(k) ENDDO ! !-- Downward facing non-natural surf_s = bc_h(1)%start_index(j,i) surf_e = bc_h(1)%end_index(j,i) DO m = surf_s, surf_e k = bc_h(1)%k(m) w_nr(k+1) = w_nr(k) w_qr(k+1) = w_qr(k) ENDDO ! !-- Neumann boundary condition at model top w_nr(nzt+1) = 0.0_wp w_qr(nzt+1) = 0.0_wp ! !-- Compute Courant number DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * & dt_micro * ddzu(k) * flag c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * & dt_micro * ddzu(k) * flag ENDDO ! !-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977): IF ( limiter_sedimentation ) THEN DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) ) d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i) qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, & 2.0_wp * d_max, & ABS( d_mean ) ) * flag d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) ) d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i) nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, & 2.0_wp * d_max, & ABS( d_mean ) ) * flag ENDDO ELSE nr_slope = 0.0_wp qr_slope = 0.0_wp ENDIF sed_nr(nzt+1) = 0.0_wp sed_qr(nzt+1) = 0.0_wp ! !-- Compute sedimentation flux DO k = nzt, nzb+1, -1 ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ! !-- Sum up all rain drop number densities which contribute to the flux !-- through k-1/2 flux = 0.0_wp z_run = 0.0_wp ! height above z(k) k_run = k c_run = MIN( 1.0_wp, c_nr(k) ) DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt ) flux = flux + hyrho(k_run) * & ( nr(k_run,j,i) + nr_slope(k_run) * ( 1.0_wp - c_run ) * & 0.5_wp ) * c_run * dzu(k_run) * flag z_run = z_run + dzu(k_run) * flag k_run = k_run + 1 * flag c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag ENDDO ! !-- It is not allowed to sediment more rain drop number density than !-- available flux = MIN( flux, & hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * dt_micro ) sed_nr(k) = flux / dt_micro * flag nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / & hyrho(k) * dt_micro * flag ! !-- Sum up all rain water content which contributes to the flux !-- through k-1/2 flux = 0.0_wp z_run = 0.0_wp ! height above z(k) k_run = k c_run = MIN( 1.0_wp, c_qr(k) ) DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt ) flux = flux + hyrho(k_run) * & ( qr(k_run,j,i) + qr_slope(k_run) * ( 1.0_wp - c_run ) * & 0.5_wp ) * c_run * dzu(k_run) * flag z_run = z_run + dzu(k_run) * flag k_run = k_run + 1 * flag c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag ENDDO ! !-- It is not allowed to sediment more rain water content than available flux = MIN( flux, & hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * dt_micro ) sed_qr(k) = flux / dt_micro * flag qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & hyrho(k) * dt_micro * flag q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & hyrho(k) * dt_micro * flag pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & hyrho(k) * lv_d_cp * d_exner(k) * dt_micro & * flag ! !-- Compute the rain rate IF ( call_microphysics_at_all_substeps ) THEN prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) & * weight_substep(intermediate_timestep_count) * flag ELSE prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag ENDIF ENDDO END SUBROUTINE sedimentation_rain_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> This subroutine computes the precipitation amount due to gravitational !> settling of rain and cloud (fog) droplets !------------------------------------------------------------------------------! SUBROUTINE calc_precipitation_amount_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< running index x direction INTEGER(iwp) :: j !< running index y direction INTEGER(iwp) :: k !< running index z direction INTEGER(iwp) :: m !< running index surface elements INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.& ( .NOT. call_microphysics_at_all_substeps .OR. & intermediate_timestep_count == intermediate_timestep_count_max ) ) & THEN surf_s = bc_h(0)%start_index(j,i) surf_e = bc_h(0)%end_index(j,i) DO m = surf_s, surf_e k = bc_h(0)%k(m) precipitation_amount(j,i) = precipitation_amount(j,i) + & prr(k,j,i) * hyrho(k) * dt_3d ENDDO ENDIF END SUBROUTINE calc_precipitation_amount_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Computation of the diagnostic supersaturation sat, actual liquid water !< temperature t_l and saturation water vapor mixing ratio q_s !------------------------------------------------------------------------------! SUBROUTINE supersaturation ( i,j,k ) IMPLICIT NONE INTEGER(iwp) :: i !< running index INTEGER(iwp) :: j !< running index INTEGER(iwp) :: k !< running index REAL(wp) :: alpha !< correction factor ! !-- Actual liquid water temperature: t_l = exner(k) * pt(k,j,i) ! !-- Calculate water vapor saturation pressure e_s = magnus( t_l ) ! !-- Computation of saturation mixing ratio: q_s = rd_d_rv * e_s / ( hyp(k) - e_s ) ! !-- Correction factor alpha = rd_d_rv * lv_d_rd * lv_d_cp / ( t_l * t_l ) ! !-- Correction of the approximated value !-- (see: Cuijpers + Duynkerke, 1993, JAS, 23) q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / ( 1.0_wp + alpha * q_s ) ! !-- Supersaturation: !-- Not in case of microphysics_kessler or microphysics_sat_adjust !-- since qr is unallocated IF ( .NOT. microphysics_kessler .AND. & .NOT. microphysics_sat_adjust ) THEN sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp ENDIF END SUBROUTINE supersaturation !------------------------------------------------------------------------------! ! Description: ! ------------ !> Calculation of the liquid water content (0%-or-100%-scheme). This scheme is !> used by the one and the two moment cloud physics scheme. Using the two moment !> scheme, this calculation results in the cloud water content. !------------------------------------------------------------------------------! SUBROUTINE calc_liquid_water_content IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: flag !< flag to indicate first grid level above surface DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation located !-- in diagnostic_quantities_mod CALL supersaturation( i, j, k ) ! !-- Compute the liquid water content IF ( microphysics_seifert .AND. .NOT. microphysics_morrison ) & THEN IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp ) THEN qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) ) * flag ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag ELSE IF ( q(k,j,i) < qr(k,j,i) ) q(k,j,i) = qr(k,j,i) qc(k,j,i) = 0.0_wp ql(k,j,i) = qr(k,j,i) * flag ENDIF ELSEIF ( microphysics_morrison ) THEN ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag ELSE IF ( ( q(k,j,i) - q_s ) > 0.0_wp ) THEN qc(k,j,i) = ( q(k,j,i) - q_s ) * flag ql(k,j,i) = qc(k,j,i) * flag ELSE qc(k,j,i) = 0.0_wp ql(k,j,i) = 0.0_wp ENDIF ENDIF ENDDO ENDDO ENDDO END SUBROUTINE calc_liquid_water_content !------------------------------------------------------------------------------! ! Description: ! ------------ !> This function computes the gamma function (Press et al., 1992). !> The gamma function is needed for the calculation of the evaporation !> of rain drops. !------------------------------------------------------------------------------! FUNCTION gamm( xx ) IMPLICIT NONE INTEGER(iwp) :: j !< REAL(wp) :: gamm !< REAL(wp) :: ser !< REAL(wp) :: tmp !< REAL(wp) :: x_gamm !< REAL(wp) :: xx !< REAL(wp) :: y_gamm !< REAL(wp), PARAMETER :: stp = 2.5066282746310005_wp !< REAL(wp), PARAMETER :: cof(6) = (/ 76.18009172947146_wp, & -86.50532032941677_wp, & 24.01409824083091_wp, & -1.231739572450155_wp, & 0.1208650973866179E-2_wp, & -0.5395239384953E-5_wp /) !< x_gamm = xx y_gamm = x_gamm tmp = x_gamm + 5.5_wp tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp ser = 1.000000000190015_wp DO j = 1, 6 y_gamm = y_gamm + 1.0_wp ser = ser + cof( j ) / y_gamm ENDDO ! !-- Until this point the algorithm computes the logarithm of the gamma !-- function. Hence, the exponential function is used. ! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) ) gamm = EXP( tmp ) * stp * ser / x_gamm RETURN END FUNCTION gamm END MODULE bulk_cloud_model_mod