!> @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-2020 Leibniz Universitaet Hannover !--------------------------------------------------------------------------------------------------! ! ! Current revisions: ! ------------------ ! ! ! Former revisions: ! ----------------- ! $Id: bulk_cloud_model_mod.f90 4577 2020-06-25 09:53:58Z suehring $ ! further re-formatting concerning Fortran parameter variables ! ! 4542 2020-05-19 15:45:12Z raasch ! file re-formatted to follow the PALM coding standard ! ! 4535 2020-05-15 12:07:23Z raasch ! bugfix for restart data format query ! ! 4533 2020-05-14 14:46:46Z schwenkel ! Reformat intrinsic function ! ! 4532 2020-05-14 13:41:35Z schwenkel ! Calculate mean droplet radius assuming gamma distibution for condensation ! ! 4521 2020-05-06 11:39:49Z schwenkel ! Rename variable ! ! 4517 2020-05-03 14:29:30Z raasch ! added restart with MPI-IO for reading local arrays ! ! 4506 2020-04-21 10:57:45Z schwenkel ! Use correct magnus formula for liquid water temperature ! ! 4502 2020-04-17 16:14:16Z schwenkel ! Implementation of ice microphysics ! ! 4495 2020-04-13 20:11:20Z raasch ! restart data handling with MPI-IO added ! ! 4457 2020-03-11 14:20:43Z raasch ! use statement for exchange horiz added ! ! 4418 2020-02-21 09:41:13Z raasch ! bugfix for raindrop number adjustment ! ! 4370 2020-01-10 14:00:44Z raasch ! vector directives added to force vectorization on Intel19 compiler ! ! 4360 2020-01-07 11:25:50Z suehring ! Introduction of wall_flags_total_0, which currently sets bits based on static ! topography information used in wall_flags_static_0 ! ! 4329 2019-12-10 15:46:36Z motisi ! Renamed wall_flags_0 to wall_flags_static_0 ! ! 4289 2019-11-05 14:33:41Z knoop ! Removed parameters precipitation and precipitation_amount_interval from namelist ! ! 4268 2019-10-17 11:29:38Z schwenkel ! Introducing bcm_boundary_conditions ! ! 4182 2019-08-22 15:20:23Z scharf ! Corrected "Former revisions" section ! ! 4168 2019-08-16 13:50:17Z suehring ! Replace function get_topography_top_index by topo_top_ind ! ! 4110 2019-07-22 17:05:21Z suehring ! Pass integer flag array as well as boundary flags to WS scalar advection ! routine ! ! 4109 2019-07-22 17:00:34Z suehring ! Added microphyics scheme 'morrision_no_rain' ! ! 3931 2019-04-24 16:34:28Z schwenkel ! Added bcm_exchange_horiz which is called after non_transport_physics ! ! 3885 2019-04-11 11:29:34Z kanani ! Changes related to global restructuring of location messages and introduction ! of additional debug messages ! ! 3874 2019-04-08 16:53:48Z knoop ! Implemented non_transport_physics module interfaces ! ! 3870 2019-04-08 13:44:34Z knoop ! Moving prognostic equations of bcm into bulk_cloud_model_mod ! ! 3869 2019-04-08 11:54:20Z knoop ! moving the furniture around ;-) ! ! 3786 2019-03-06 16:58:03Z raasch ! unsed variables removed ! ! 3767 2019-02-27 08:18:02Z raasch ! unused variable for file index removed from rrd-subroutines parameter list ! ! 3724 2019-02-06 16:28:23Z kanani ! Correct double-used log_point_s unit ! ! 3700 2019-01-26 17:03:42Z knoop ! nopointer option removed ! 1053 2012-11-13 17:11:03Z hoffmann ! initial revision ! ! Description: ! ------------ !> Calculate bulk cloud microphysics. !--------------------------------------------------------------------------------------------------! MODULE bulk_cloud_model_mod USE advec_s_bc_mod, & ONLY: advec_s_bc USE advec_s_pw_mod, & ONLY: advec_s_pw USE advec_s_up_mod, & ONLY: advec_s_up USE advec_ws, & ONLY: advec_s_ws USE arrays_3d, & ONLY: diss_s_qc, diss_s_qr, diss_s_nc, diss_s_nr, & diss_l_qc, diss_l_qr, diss_l_nc, diss_l_nr, & diss_l_qi, diss_l_ni, diss_s_qi, diss_s_ni, & ddzu, diss, dzu, dzw, hyp, hyrho, & exner, zu, tnc_m, tnr_m, tqc_m, tqr_m, tend, rdf_sc, & flux_l_qc, flux_l_qr, flux_l_nc, flux_l_nr, & flux_l_qi, flux_l_ni, flux_s_qi, flux_s_ni, & flux_s_qc, flux_s_qr, flux_s_nc, flux_s_nr, & nc, nc_1, nc_2, nc_3, nc_p, nr, nr_1, nr_2, nr_3, nr_p, & ni, ni_1, ni_2, ni_3, ni_p, tni_m, & 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, & qi, qi_1, qi_2, qi_3, qi_p, tqi_m USE averaging, & ONLY: nc_av,ni_av, nr_av, prr_av, qc_av, ql_av, qi_av, qr_av USE basic_constants_and_equations_mod, & ONLY: c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, magnus_ice, magnus_tl, & exner_function, exner_function_invers, ideal_gas_law_rho, & ideal_gas_law_rho_pt, barometric_formula, rd_d_rv, l_s, & ls_d_cp, & molecular_weight_of_solute, & molecular_weight_of_water, pi, rho_l, rho_s, r_d, r_v, vanthoff USE control_parameters, & ONLY: bc_dirichlet_l, & bc_dirichlet_n, & bc_dirichlet_r, & bc_dirichlet_s, & bc_radiation_l, & bc_radiation_n, & bc_radiation_r, & bc_radiation_s, & debug_output, & dt_3d, dt_do2d_xy, intermediate_timestep_count, & intermediate_timestep_count_max, large_scale_forcing, & loop_optimization, simulated_time, & lsf_surf, pt_surface, restart_data_format_output, rho_surface, & surface_pressure, & time_do2d_xy, message_string, initializing_actions, & ws_scheme_sca, scalar_advec, timestep_scheme, tsc USE cpulog, & ONLY: cpu_log, log_point, log_point_s USE diffusion_s_mod, & ONLY: diffusion_s USE grid_variables, & ONLY: dx, dy USE indices, & ONLY: advc_flags_s, & nbgp, nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt, & topo_top_ind, & wall_flags_total_0 USE kinds USE pegrid, & ONLY: threads_per_task USE restart_data_mpi_io_mod, & ONLY: rd_mpi_io_check_array, rrd_mpi_io, wrd_mpi_io USE statistics, & ONLY: sums_wsncs_ws_l, sums_wsnis_ws_l, sums_wsnrs_ws_l, sums_wsqcs_ws_l, sums_wsqis_ws_l,& sums_wsqrs_ws_l, weight_pres, weight_substep USE surface_mod, & ONLY : bc_h, & surf_bulk_cloud_model, & surf_microphysics_morrison, surf_microphysics_seifert, & surf_microphysics_ice_phase, & surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & surf_usm_v 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 :: call_microphysics_at_all_substeps = .FALSE. !< namelist parameter LOGICAL :: cloud_water_sedimentation = .FALSE. !< cloud water sedimentation LOGICAL :: collision_turbulence = .FALSE. !< turbulence effects LOGICAL :: curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory LOGICAL :: ice_crystal_sedimentation = .FALSE. !< flag for ice crystal sedimentation LOGICAL :: limiter_sedimentation = .TRUE. !< sedimentation limiter LOGICAL :: microphysics_ice_phase = .FALSE. !< use ice microphysics 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_morrison_no_rain = .FALSE. !< use 2-moment Morrison LOGICAL :: microphysics_sat_adjust = .FALSE. !< use saturation adjust bulk scheme? LOGICAL :: microphysics_seifert = .FALSE. !< use 2-moment Seifert and Beheng !< scheme LOGICAL :: precipitation = .FALSE. !< namelist parameter LOGICAL :: ventilation_effect = .TRUE. !< ventilation effect 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_sedimentation = 2.0_wp !< Courant number of sedimentation process 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) :: 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) :: eps_sb = 1.0E-10_wp !< threshold in two-moments scheme REAL(wp) :: dt_precipitation = 100.0_wp !< timestep precipitation (s) REAL(wp) :: e_s !< saturation water vapor pressure REAL(wp) :: e_si !< saturation water vapor pressure over ice REAL(wp) :: eps_mr = 0.0_wp !< threshold for morrison scheme REAL(wp) :: in_init = 1000.0_wp !< initial number of potential ice nucleii 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) :: na_init = 100.0E6_wp !< Total particle/aerosol concentration (cm-3) REAL(wp) :: nc_const = 70.0E6_wp !< cloud droplet concentration REAL(wp) :: pirho_l !< pi * rho_l / 6.0 REAL(wp) :: prec_time_const = 0.001_wp !< coef. in Kessler scheme (s-1) REAL(wp) :: precipitation_amount_interval = 9999999.9_wp !< namelist parameter REAL(wp) :: ql_crit = 0.0005_wp !< coef. in Kessler scheme (kg kg-1) REAL(wp) :: q_s !< saturation mixing ratio REAL(wp) :: q_si !< saturation mixing ratio over ice REAL(wp) :: sat !< supersaturation REAL(wp) :: sat_ice !< supersaturation over ice REAL(wp) :: start_ice_microphysics = 0.0_wp !< time from which on ice microhysics are executed REAL(wp) :: sed_qc_const !< const. for sedimentation of cloud water REAL(wp) :: schmidt_p_1d3=0.8921121_wp !< Schmidt number**0.33333, 0.71**0.33333 REAL(wp) :: sigma_bulk = 2.0_wp !< width of aerosol spectrum REAL(wp) :: sigma_gc = 1.3_wp !< geometric standard deviation cloud droplets REAL(wp) :: t_l !< liquid-(ice) water temperature 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) :: ximin = 4.42E-14_wp !< minimum mass of ice crystal (kg) (D~10µ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) SAVE PRIVATE PUBLIC bcm_parin, & bcm_check_parameters, & bcm_check_data_output, & bcm_check_data_output_pr, & bcm_init_arrays, & bcm_init, & bcm_header, & bcm_actions, & bcm_non_advective_processes, & bcm_exchange_horiz, & bcm_prognostic_equations, & bcm_boundary_conditions, & 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, & calc_liquid_water_content PUBLIC bulk_cloud_model, & call_microphysics_at_all_substeps, & cloud_water_sedimentation, & cloud_scheme, & collision_turbulence, & dt_precipitation, & ice_crystal_sedimentation, & in_init, & microphysics_morrison, & microphysics_morrison_no_rain, & microphysics_sat_adjust, & microphysics_seifert, & microphysics_ice_phase, & na_init, & nc_const, & precipitation, & sigma_gc, & start_ice_microphysics 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_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_header MODULE PROCEDURE bcm_header END INTERFACE bcm_header INTERFACE bcm_actions MODULE PROCEDURE bcm_actions MODULE PROCEDURE bcm_actions_ij END INTERFACE bcm_actions INTERFACE bcm_non_advective_processes MODULE PROCEDURE bcm_non_advective_processes MODULE PROCEDURE bcm_non_advective_processes_ij END INTERFACE bcm_non_advective_processes INTERFACE bcm_exchange_horiz MODULE PROCEDURE bcm_exchange_horiz END INTERFACE bcm_exchange_horiz INTERFACE bcm_prognostic_equations MODULE PROCEDURE bcm_prognostic_equations MODULE PROCEDURE bcm_prognostic_equations_ij END INTERFACE bcm_prognostic_equations INTERFACE bcm_boundary_conditions MODULE PROCEDURE bcm_boundary_conditions END INTERFACE bcm_boundary_conditions INTERFACE bcm_swap_timelevel MODULE PROCEDURE bcm_swap_timelevel END INTERFACE bcm_swap_timelevel 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_rrd_global MODULE PROCEDURE bcm_rrd_global_ftn MODULE PROCEDURE bcm_rrd_global_mpi END INTERFACE bcm_rrd_global INTERFACE bcm_rrd_local MODULE PROCEDURE bcm_rrd_local_ftn MODULE PROCEDURE bcm_rrd_local_mpi 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 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 bulk cloud 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, & bulk_cloud_model, & c_sedimentation, & call_microphysics_at_all_substeps, & cloud_scheme, & cloud_water_sedimentation, & collision_turbulence, & curvature_solution_effects_bulk, & dry_aerosol_radius, & ice_crystal_sedimentation, & in_init, & limiter_sedimentation, & microphysics_ice_phase, & na_init, & nc_const, & sigma_bulk, & start_ice_microphysics, & ventilation_effect line = ' ' ! !-- Try to find bulk cloud module namelist 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 bulk cloud module is switched on !bulk_cloud_model = .TRUE. 10 CONTINUE END SUBROUTINE bcm_parin !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Check parameters routine for bulk cloud module !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_check_parameters IMPLICIT NONE ! !-- Check cloud scheme !-- This scheme considers only saturation adjustment, i.e. water vapor surplus is converted into !-- liquid water. No other microphysical processes are considered. IF ( cloud_scheme == 'saturation_adjust' ) THEN microphysics_sat_adjust = .TRUE. microphysics_seifert = .FALSE. microphysics_kessler = .FALSE. precipitation = .FALSE. microphysics_morrison_no_rain = .FALSE. ! !-- This scheme includes all process of the seifert beheng scheme (2001,2006). Especially rain !-- processes are considered with prognostic quantities of qr and nr. !-- Cloud droplet concentration is assumed to be constant and qc is diagnostic. !-- Technical remark: The switch 'microphysics_seifert' allocates fields of qr and nr and enables !-- all rain processes. ELSEIF ( cloud_scheme == 'seifert_beheng' ) THEN microphysics_sat_adjust = .FALSE. microphysics_seifert = .TRUE. microphysics_kessler = .FALSE. microphysics_morrison = .FALSE. precipitation = .TRUE. microphysics_morrison_no_rain = .FALSE. ! !-- The kessler scheme is a simplified scheme without any prognostic quantities for microphyical !-- variables but considering autoconversion. ELSEIF ( cloud_scheme == 'kessler' ) THEN microphysics_sat_adjust = .FALSE. microphysics_seifert = .FALSE. microphysics_kessler = .TRUE. microphysics_morrison = .FALSE. precipitation = .TRUE. microphysics_morrison_no_rain = .FALSE. ! !-- The morrison scheme is an extension of the seifer beheng scheme including also relevant !-- processes for cloud droplet size particles such as activation and an diagnostic mehtod for !-- diffusional growth. !-- I.e. here all processes of Seifert and Beheng as well as of the morrision scheme are used. !-- Therefore, ztis includes prognostic quantities for qc and nc. !-- Technical remark: The switch 'microphysics_morrison' allocates fields of qc and nc and !-- enables diagnostic diffusional growth and activation. ELSEIF ( cloud_scheme == 'morrison' ) THEN microphysics_sat_adjust = .FALSE. microphysics_seifert = .TRUE. microphysics_kessler = .FALSE. microphysics_morrison = .TRUE. precipitation = .TRUE. microphysics_morrison_no_rain = .FALSE. ! !-- The 'morrision_no_rain' scheme includes only processes of morrision scheme without the rain !-- processes of seifert beheng. Therfore, the prog. quantities of qr and nr remain unallocated. !-- This might be appropiate for cloud in which the size distribution is narrow, e.g. fog. ELSEIF ( cloud_scheme == 'morrison_no_rain' ) THEN microphysics_sat_adjust = .FALSE. microphysics_seifert = .FALSE. microphysics_kessler = .FALSE. microphysics_morrison = .TRUE. microphysics_morrison_no_rain = .TRUE. precipitation = .FALSE. 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 solution for circular dependency problem surf_bulk_cloud_model = bulk_cloud_model surf_microphysics_morrison = microphysics_morrison surf_microphysics_seifert = microphysics_seifert surf_microphysics_ice_phase = microphysics_ice_phase ! !-- 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 bulk cloud 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 ( 'ni' ) IF ( .NOT. microphysics_ice_phase ) THEN message_string = 'output of "' // TRIM( var ) // '" ' // 'requires ' // & 'microphysics_ice_phase = ".TRUE."' 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 ( 'qi' ) IF ( .NOT. microphysics_ice_phase ) THEN message_string = 'output of "' // TRIM( var ) // '" ' // 'requires ' // & 'microphysics_ice_phase = ".TRUE."' CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 ) ENDIF 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 bulk cloud 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 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 = 123 dopr_index(var_count) = pr_index dopr_unit = '1/m3' unit = dopr_unit hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 ) CASE ( 'ni' ) IF ( .NOT. microphysics_ice_phase ) THEN message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // & ' is not implemented for' // ' microphysics_ice_phase = ".F."' CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 ) ENDIF pr_index = 124 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 ( 'qi' ) IF ( .NOT. microphysics_ice_phase ) THEN message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // & ' is not implemented for' // ' microphysics_ice_phase = ".F."' CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 ) ENDIF pr_index = 125 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 = 74 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 bulk cloud module arrays and define pointers !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_init_arrays USE indices, & ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt IMPLICIT NONE ! !-- Liquid water content ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ! !-- 3D-cloud water content IF ( .NOT. microphysics_morrison ) THEN ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 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 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 IF ( microphysics_seifert ) THEN ! !-- 3D-rain water content, rain drop concentration arrays 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 IF ( microphysics_ice_phase ) THEN ! !-- 3D-cloud drop water content, cloud drop concentration arrays ALLOCATE( ni_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & ni_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & ni_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qi_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qi_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & qi_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF IF ( ws_scheme_sca ) THEN IF ( microphysics_morrison ) THEN ALLOCATE( sums_wsqcs_ws_l(nzb:nzt+1,0:threads_per_task-1) ) ALLOCATE( sums_wsncs_ws_l(nzb:nzt+1,0:threads_per_task-1) ) sums_wsqcs_ws_l = 0.0_wp sums_wsncs_ws_l = 0.0_wp ENDIF IF ( microphysics_seifert ) THEN ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) ) ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) ) sums_wsqrs_ws_l = 0.0_wp sums_wsnrs_ws_l = 0.0_wp ENDIF IF ( microphysics_ice_phase ) THEN ALLOCATE( sums_wsqis_ws_l(nzb:nzt+1,0:threads_per_task-1) ) ALLOCATE( sums_wsnis_ws_l(nzb:nzt+1,0:threads_per_task-1) ) sums_wsqis_ws_l = 0.0_wp sums_wsnis_ws_l = 0.0_wp ENDIF ENDIF ! !-- Arrays needed for reasons of speed optimization for cache version. !-- For the vector version the buffer arrays are not necessary, because the the fluxes can !-- swapped directly inside the loops of the advection routines. IF ( loop_optimization /= 'vector' ) THEN IF ( ws_scheme_sca ) THEN IF ( microphysics_morrison ) THEN ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1), & diss_s_qc(nzb+1:nzt,0:threads_per_task-1), & flux_s_nc(nzb+1:nzt,0:threads_per_task-1), & diss_s_nc(nzb+1:nzt,0:threads_per_task-1) ) ALLOCATE( flux_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & flux_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) ENDIF IF ( microphysics_seifert ) THEN ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1), & diss_s_qr(nzb+1:nzt,0:threads_per_task-1), & flux_s_nr(nzb+1:nzt,0:threads_per_task-1), & diss_s_nr(nzb+1:nzt,0:threads_per_task-1) ) ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) ENDIF IF ( microphysics_ice_phase ) THEN ALLOCATE( flux_s_qi(nzb+1:nzt,0:threads_per_task-1), & diss_s_qi(nzb+1:nzt,0:threads_per_task-1), & flux_s_ni(nzb+1:nzt,0:threads_per_task-1), & diss_s_ni(nzb+1:nzt,0:threads_per_task-1) ) ALLOCATE( flux_l_qi(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_qi(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & flux_l_ni(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_ni(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) ENDIF ENDIF ENDIF ! !-- 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 IF ( microphysics_ice_phase ) THEN qi => qi_1; qi_p => qi_2; tqi_m => qi_3 ni => ni_1; ni_p => ni_2; tni_m => ni_3 ENDIF END SUBROUTINE bcm_init_arrays !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Initialization of the bulk cloud module !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_init IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< IF ( debug_output ) CALL debug_message( 'bcm_init', 'start' ) IF ( bulk_cloud_model ) THEN IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN ! !-- 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 ! !-- Initialize the remaining quantities IF ( microphysics_ice_phase ) THEN DO i = nxlg, nxrg DO j = nysg, nyng qi(:,j,i) = 0.0_wp ni(:,j,i) = 0.0_wp ENDDO ENDDO ENDIF ! !-- Liquid water content and precipitation amount are zero at beginning of the simulation. ql = 0.0_wp qc = 0.0_wp precipitation_amount = 0.0_wp prr = 0.0_wp ! !-- 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 IF ( microphysics_ice_phase ) THEN tqi_m = 0.0_wp tni_m = 0.0_wp qi_p = qi ni_p = ni ENDIF ENDIF ! Only if not read_restart_data ! !-- 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 IF ( debug_output ) CALL debug_message( 'bcm_init', 'end' ) ELSE IF ( debug_output ) CALL debug_message( 'bcm_init skipped', 'end' ) ENDIF END SUBROUTINE bcm_init !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Header output for bulk cloud module !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_header ( io ) IMPLICIT NONE INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file ! !-- Write bulk cloud 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 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: ! ------------ !> Control of microphysics for all grid points !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_actions( location ) CHARACTER (LEN=*), INTENT(IN) :: location !< call location string SELECT CASE ( location ) CASE ( 'before_timestep' ) IF ( ws_scheme_sca ) THEN IF ( microphysics_morrison ) THEN sums_wsqcs_ws_l = 0.0_wp sums_wsncs_ws_l = 0.0_wp ENDIF IF ( microphysics_seifert ) THEN sums_wsqrs_ws_l = 0.0_wp sums_wsnrs_ws_l = 0.0_wp ENDIF IF ( microphysics_ice_phase ) THEN sums_wsqis_ws_l = 0.0_wp sums_wsnis_ws_l = 0.0_wp ENDIF ENDIF CASE DEFAULT CONTINUE END SELECT END SUBROUTINE bcm_actions !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Control of microphysics for grid points i,j !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_actions_ij( i, j, location ) CHARACTER (LEN=*), INTENT(IN) :: location !< call location string INTEGER(iwp) :: dummy !< call location string INTEGER(iwp), INTENT(IN) :: i !< grid index in x-direction INTEGER(iwp), INTENT(IN) :: j !< grid index in y-direction IF ( bulk_cloud_model ) dummy = i + j SELECT CASE ( location ) CASE ( 'before_timestep' ) IF ( ws_scheme_sca ) THEN IF ( microphysics_morrison ) THEN sums_wsqcs_ws_l = 0.0_wp sums_wsncs_ws_l = 0.0_wp ENDIF IF ( microphysics_seifert ) THEN sums_wsqrs_ws_l = 0.0_wp sums_wsnrs_ws_l = 0.0_wp ENDIF IF ( microphysics_ice_phase ) THEN sums_wsqis_ws_l = 0.0_wp sums_wsnis_ws_l = 0.0_wp ENDIF ENDIF CASE DEFAULT CONTINUE END SELECT END SUBROUTINE bcm_actions_ij !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Control of microphysics for all grid points !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_non_advective_processes CALL cpu_log( log_point(51), 'microphysics', 'start' ) IF ( .NOT. microphysics_sat_adjust .AND. ( intermediate_timestep_count == 1 .OR. & call_microphysics_at_all_substeps ) ) & THEN 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. !-- Here the the simple kessler scheme is used. IF ( microphysics_kessler ) THEN CALL autoconversion_kessler IF ( cloud_water_sedimentation ) CALL sedimentation_cloud ! !-- Here the seifert beheng scheme is used. Cloud concentration is assumed to a constant value !-- an qc a diagnostic value. ELSEIF ( microphysics_seifert .AND. .NOT. microphysics_morrison ) THEN CALL adjust_cloud IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics ) THEN CALL ice_nucleation CALL ice_deposition ENDIF CALL autoconversion CALL accretion CALL selfcollection_breakup CALL evaporation_rain CALL sedimentation_rain IF ( cloud_water_sedimentation ) CALL sedimentation_cloud IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics ) THEN CALL adjust_ice IF ( ice_crystal_sedimentation ) CALL sedimentation_ice ENDIF ! !-- Here the morrison scheme is used. No rain processes are considered and qr and nr are not !-- allocated. ELSEIF ( microphysics_morrison_no_rain .AND. .NOT. microphysics_seifert ) THEN CALL activation CALL condensation IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics ) THEN CALL adjust_ice CALL ice_nucleation CALL ice_deposition ENDIF IF ( cloud_water_sedimentation ) CALL sedimentation_cloud ! !-- Here the full morrison scheme is used and all processes of Seifert and Beheng are included ELSEIF ( microphysics_morrison .AND. microphysics_seifert ) THEN CALL adjust_cloud CALL activation CALL condensation IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics ) THEN CALL adjust_ice CALL ice_nucleation CALL ice_deposition ENDIF 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 ENDIF CALL cpu_log( log_point(51), 'microphysics', 'stop' ) END SUBROUTINE bcm_non_advective_processes !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Control of microphysics for grid points i,j !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_non_advective_processes_ij( i, j ) INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< IF ( .NOT. microphysics_sat_adjust .AND. ( intermediate_timestep_count == 1 .OR. & call_microphysics_at_all_substeps ) ) & THEN 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 !-- Here the the simple kessler scheme is used. IF( microphysics_kessler ) THEN CALL autoconversion_kessler_ij( i,j ) IF ( cloud_water_sedimentation ) CALL sedimentation_cloud_ij( i,j ) ! !-- Here the seifert beheng scheme is used. Cloud concentration is assumed to a constant value !-- an qc a diagnostic value. ELSEIF ( microphysics_seifert .AND. .NOT. microphysics_morrison ) THEN CALL adjust_cloud_ij( i,j ) IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics ) THEN CALL ice_nucleation_ij( i,j ) CALL ice_deposition_ij( i,j ) ENDIF CALL autoconversion_ij( i,j ) CALL accretion_ij( i,j ) CALL selfcollection_breakup_ij( i,j ) CALL evaporation_rain_ij( i,j ) CALL sedimentation_rain_ij( i,j ) IF ( cloud_water_sedimentation ) CALL sedimentation_cloud_ij( i,j ) IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics ) THEN CALL adjust_ice_ij ( i,j ) IF ( ice_crystal_sedimentation ) CALL sedimentation_ice_ij ( i,j ) ENDIF ! !-- Here the morrison scheme is used. No rain processes are considered and qr and nr are not !-- allocated. ELSEIF ( microphysics_morrison_no_rain .AND. .NOT. microphysics_seifert ) THEN CALL activation_ij( i,j ) CALL condensation_ij( i,j ) IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics ) THEN CALL adjust_ice_ij ( i,j ) CALL ice_nucleation_ij( i,j ) CALL ice_deposition_ij( i,j ) ENDIF IF ( cloud_water_sedimentation ) CALL sedimentation_cloud_ij( i,j ) ! !-- Here the full morrison scheme is used and all processes of Seifert and Beheng are included ELSEIF ( microphysics_morrison .AND. microphysics_seifert ) THEN CALL adjust_cloud_ij( i,j ) CALL activation_ij( i,j ) CALL condensation_ij( i,j ) IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics ) THEN CALL adjust_ice_ij ( i,j ) CALL ice_nucleation_ij( i,j ) CALL ice_deposition_ij( i,j ) ENDIF CALL autoconversion_ij( i,j ) CALL accretion_ij( i,j ) CALL selfcollection_breakup_ij( i,j ) CALL evaporation_rain_ij( i,j ) CALL sedimentation_rain_ij( i,j ) IF ( cloud_water_sedimentation ) CALL sedimentation_cloud_ij( i,j ) ENDIF CALL calc_precipitation_amount_ij( i,j ) ENDIF END SUBROUTINE bcm_non_advective_processes_ij !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Control of microphysics for all grid points !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_exchange_horiz USE exchange_horiz_mod, & ONLY: exchange_horiz IF ( .NOT. microphysics_sat_adjust .AND. ( intermediate_timestep_count == 1 .OR. & call_microphysics_at_all_substeps ) ) & THEN IF ( microphysics_morrison ) THEN CALL exchange_horiz( nc, nbgp ) CALL exchange_horiz( qc, nbgp ) ENDIF IF ( microphysics_seifert ) THEN CALL exchange_horiz( qr, nbgp ) CALL exchange_horiz( nr, nbgp ) ENDIF IF ( microphysics_ice_phase ) THEN CALL exchange_horiz( qi, nbgp ) CALL exchange_horiz( ni, nbgp ) ENDIF CALL exchange_horiz( q, nbgp ) CALL exchange_horiz( pt, nbgp ) ENDIF END SUBROUTINE bcm_exchange_horiz !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Control of microphysics for all grid points !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_prognostic_equations INTEGER(iwp) :: i !< grid index in x-direction INTEGER(iwp) :: j !< grid index in y-direction INTEGER(iwp) :: k !< grid index in z-direction REAL(wp) :: sbt !< ! !-- If required, calculate prognostic equations for cloud water content and cloud drop !-- concentration. IF ( microphysics_morrison ) THEN CALL cpu_log( log_point(67), 'qc-equation', 'start' ) ! !-- Calculate prognostic equation for cloud water content sbt = tsc(2) IF ( scalar_advec == 'bc-scheme' ) THEN IF ( timestep_scheme(1:5) /= 'runge' ) THEN ! !-- Bott-Chlond scheme always uses Euler time step. Thus: sbt = 1.0_wp ENDIF tend = 0.0_wp CALL advec_s_bc( qc, 'qc' ) ENDIF ! !-- qc-tendency terms with no communication IF ( scalar_advec /= 'bc-scheme' ) THEN tend = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( advc_flags_s, qc, 'qc', & bc_dirichlet_l .OR. bc_radiation_l, & bc_dirichlet_n .OR. bc_radiation_n, & bc_dirichlet_r .OR. bc_radiation_r, & bc_dirichlet_s .OR. bc_radiation_s ) ELSE CALL advec_s_pw( qc ) ENDIF ELSE CALL advec_s_up( qc ) ENDIF ENDIF CALL diffusion_s( qc, & surf_def_h(0)%qcsws, surf_def_h(1)%qcsws, & surf_def_h(2)%qcsws, & surf_lsm_h%qcsws, surf_usm_h%qcsws, & surf_def_v(0)%qcsws, surf_def_v(1)%qcsws, & surf_def_v(2)%qcsws, surf_def_v(3)%qcsws, & surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws, & surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws, & surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws, & surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws ) ! !-- Prognostic equation for cloud water content DO i = nxl, nxr DO j = nys, nyn !following directive is required to vectorize on Intel19 !DIR$ IVDEP DO k = nzb+1, nzt qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & tsc(3) * tqc_m(k,j,i) ) & - tsc(5) * rdf_sc(k) * & qc(k,j,i) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & ) IF ( qc_p(k,j,i) < 0.0_wp ) qc_p(k,j,i) = 0.0_wp ENDDO ENDDO ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tqc_m(k,j,i) = tend(k,j,i) ENDDO ENDDO ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tqc_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tqc_m(k,j,i) ENDDO ENDDO ENDDO ENDIF ENDIF CALL cpu_log( log_point(67), 'qc-equation', 'stop' ) CALL cpu_log( log_point(68), 'nc-equation', 'start' ) ! !-- Calculate prognostic equation for cloud drop concentration sbt = tsc(2) IF ( scalar_advec == 'bc-scheme' ) THEN IF ( timestep_scheme(1:5) /= 'runge' ) THEN ! !-- Bott-Chlond scheme always uses Euler time step. Thus: sbt = 1.0_wp ENDIF tend = 0.0_wp CALL advec_s_bc( nc, 'nc' ) ENDIF ! !-- nc-tendency terms with no communication IF ( scalar_advec /= 'bc-scheme' ) THEN tend = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( advc_flags_s, nc, 'nc', & bc_dirichlet_l .OR. bc_radiation_l, & bc_dirichlet_n .OR. bc_radiation_n, & bc_dirichlet_r .OR. bc_radiation_r, & bc_dirichlet_s .OR. bc_radiation_s ) ELSE CALL advec_s_pw( nc ) ENDIF ELSE CALL advec_s_up( nc ) ENDIF ENDIF CALL diffusion_s( nc, & surf_def_h(0)%ncsws, surf_def_h(1)%ncsws, & surf_def_h(2)%ncsws, & surf_lsm_h%ncsws, surf_usm_h%ncsws, & surf_def_v(0)%ncsws, surf_def_v(1)%ncsws, & surf_def_v(2)%ncsws, surf_def_v(3)%ncsws, & surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws, & surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws, & surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws, & surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws ) ! !-- Prognostic equation for cloud drop concentration DO i = nxl, nxr DO j = nys, nyn !following directive is required to vectorize on Intel19 !DIR$ IVDEP DO k = nzb+1, nzt nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & tsc(3) * tnc_m(k,j,i) ) & - tsc(5) * rdf_sc(k) * & nc(k,j,i) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & ) IF ( nc_p(k,j,i) < 0.0_wp ) nc_p(k,j,i) = 0.0_wp ENDDO ENDDO ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tnc_m(k,j,i) = tend(k,j,i) ENDDO ENDDO ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tnc_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tnc_m(k,j,i) ENDDO ENDDO ENDDO ENDIF ENDIF CALL cpu_log( log_point(68), 'nc-equation', 'stop' ) ENDIF ! !-- If required, calculate prognostic equations for ice crystal content and ice crystal !-- concentration IF ( microphysics_ice_phase ) THEN CALL cpu_log( log_point(70), 'qi-equation', 'start' ) ! !-- Calculate prognostic equation for ice crystal content sbt = tsc(2) IF ( scalar_advec == 'bc-scheme' ) THEN IF ( timestep_scheme(1:5) /= 'runge' ) THEN ! !-- Bott-Chlond scheme always uses Euler time step. Thus: sbt = 1.0_wp ENDIF tend = 0.0_wp CALL advec_s_bc( qi, 'qi' ) ENDIF ! !-- qi-tendency terms with no communication IF ( scalar_advec /= 'bc-scheme' ) THEN tend = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( advc_flags_s, qi, 'qi', & bc_dirichlet_l .OR. bc_radiation_l, & bc_dirichlet_n .OR. bc_radiation_n, & bc_dirichlet_r .OR. bc_radiation_r, & bc_dirichlet_s .OR. bc_radiation_s ) ELSE CALL advec_s_pw( qi ) ENDIF ELSE CALL advec_s_up( qi ) ENDIF ENDIF CALL diffusion_s( qi, & surf_def_h(0)%qisws, surf_def_h(1)%qisws, & surf_def_h(2)%qisws, & surf_lsm_h%qisws, surf_usm_h%qisws, & surf_def_v(0)%qisws, surf_def_v(1)%qisws, & surf_def_v(2)%qisws, surf_def_v(3)%qisws, & surf_lsm_v(0)%qisws, surf_lsm_v(1)%qisws, & surf_lsm_v(2)%qisws, surf_lsm_v(3)%qisws, & surf_usm_v(0)%qisws, surf_usm_v(1)%qisws, & surf_usm_v(2)%qisws, surf_usm_v(3)%qisws ) ! !-- Prognostic equation for ice crystal mixing ratio DO i = nxl, nxr DO j = nys, nyn !following directive is required to vectorize on Intel19 !DIR$ IVDEP DO k = nzb+1, nzt qi_p(k,j,i) = qi(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & tsc(3) * tqi_m(k,j,i) ) & - tsc(5) * rdf_sc(k) * & qi(k,j,i) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & ) IF ( qi_p(k,j,i) < 0.0_wp ) qi_p(k,j,i) = 0.0_wp ENDDO ENDDO ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tqi_m(k,j,i) = tend(k,j,i) ENDDO ENDDO ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tqi_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tqi_m(k,j,i) ENDDO ENDDO ENDDO ENDIF ENDIF CALL cpu_log( log_point(70), 'qi-equation', 'stop' ) CALL cpu_log( log_point(69), 'ni-equation', 'start' ) ! !-- Calculate prognostic equation for ice crystal concentration sbt = tsc(2) IF ( scalar_advec == 'bc-scheme' ) THEN IF ( timestep_scheme(1:5) /= 'runge' ) THEN ! !-- Bott-Chlond scheme always uses Euler time step. Thus: sbt = 1.0_wp ENDIF tend = 0.0_wp CALL advec_s_bc( ni, 'ni' ) ENDIF ! !-- ni-tendency terms with no communication IF ( scalar_advec /= 'bc-scheme' ) THEN tend = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( advc_flags_s, ni, 'ni', & bc_dirichlet_l .OR. bc_radiation_l, & bc_dirichlet_n .OR. bc_radiation_n, & bc_dirichlet_r .OR. bc_radiation_r, & bc_dirichlet_s .OR. bc_radiation_s ) ELSE CALL advec_s_pw( ni ) ENDIF ELSE CALL advec_s_up( ni ) ENDIF ENDIF CALL diffusion_s( ni, & surf_def_h(0)%nisws, surf_def_h(1)%nisws, & surf_def_h(2)%nisws, & surf_lsm_h%nisws, surf_usm_h%nisws, & surf_def_v(0)%nisws, surf_def_v(1)%nisws, & surf_def_v(2)%nisws, surf_def_v(3)%nisws, & surf_lsm_v(0)%nisws, surf_lsm_v(1)%nisws, & surf_lsm_v(2)%nisws, surf_lsm_v(3)%nisws, & surf_usm_v(0)%nisws, surf_usm_v(1)%nisws, & surf_usm_v(2)%nisws, surf_usm_v(3)%nisws ) ! !-- Prognostic equation for ice crystal concentration DO i = nxl, nxr DO j = nys, nyn !following directive is required to vectorize on Intel19 !DIR$ IVDEP DO k = nzb+1, nzt ni_p(k,j,i) = ni(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & tsc(3) * tni_m(k,j,i) ) & - tsc(5) * rdf_sc(k) * & ni(k,j,i) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & ) IF ( ni_p(k,j,i) < 0.0_wp ) ni_p(k,j,i) = 0.0_wp ENDDO ENDDO ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tni_m(k,j,i) = tend(k,j,i) ENDDO ENDDO ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tni_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tni_m(k,j,i) ENDDO ENDDO ENDDO ENDIF ENDIF CALL cpu_log( log_point(69), 'ni-equation', 'stop' ) ENDIF ! !-- If required, calculate prognostic equations for rain water content and rain drop !-- concentration. IF ( microphysics_seifert ) THEN CALL cpu_log( log_point(52), 'qr-equation', 'start' ) ! !-- Calculate prognostic equation for rain water content sbt = tsc(2) IF ( scalar_advec == 'bc-scheme' ) THEN IF ( timestep_scheme(1:5) /= 'runge' ) THEN ! !-- Bott-Chlond scheme always uses Euler time step. Thus: sbt = 1.0_wp ENDIF tend = 0.0_wp CALL advec_s_bc( qr, 'qr' ) ENDIF ! !-- qr-tendency terms with no communication IF ( scalar_advec /= 'bc-scheme' ) THEN tend = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( advc_flags_s, qr, 'qr', & bc_dirichlet_l .OR. bc_radiation_l, & bc_dirichlet_n .OR. bc_radiation_n, & bc_dirichlet_r .OR. bc_radiation_r, & bc_dirichlet_s .OR. bc_radiation_s ) ELSE CALL advec_s_pw( qr ) ENDIF ELSE CALL advec_s_up( qr ) ENDIF ENDIF CALL diffusion_s( qr, & surf_def_h(0)%qrsws, surf_def_h(1)%qrsws, & surf_def_h(2)%qrsws, & surf_lsm_h%qrsws, surf_usm_h%qrsws, & surf_def_v(0)%qrsws, surf_def_v(1)%qrsws, & surf_def_v(2)%qrsws, surf_def_v(3)%qrsws, & surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws, & surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws, & surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws, & surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws ) ! !-- Prognostic equation for rain water content DO i = nxl, nxr DO j = nys, nyn !following directive is required to vectorize on Intel19 !DIR$ IVDEP DO k = nzb+1, nzt qr_p(k,j,i) = qr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & tsc(3) * tqr_m(k,j,i) ) & - tsc(5) * rdf_sc(k) * & qr(k,j,i) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & ) IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp ENDDO ENDDO ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tqr_m(k,j,i) = tend(k,j,i) ENDDO ENDDO ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tqr_m(k,j,i) ENDDO ENDDO ENDDO ENDIF ENDIF CALL cpu_log( log_point(52), 'qr-equation', 'stop' ) CALL cpu_log( log_point(53), 'nr-equation', 'start' ) ! !-- Calculate prognostic equation for rain drop concentration sbt = tsc(2) IF ( scalar_advec == 'bc-scheme' ) THEN IF ( timestep_scheme(1:5) /= 'runge' ) THEN ! !-- Bott-Chlond scheme always uses Euler time step. Thus: sbt = 1.0_wp ENDIF tend = 0.0_wp CALL advec_s_bc( nr, 'nr' ) ENDIF ! !-- nr-tendency terms with no communication IF ( scalar_advec /= 'bc-scheme' ) THEN tend = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( advc_flags_s, nr, 'nr', & bc_dirichlet_l .OR. bc_radiation_l, & bc_dirichlet_n .OR. bc_radiation_n, & bc_dirichlet_r .OR. bc_radiation_r, & bc_dirichlet_s .OR. bc_radiation_s ) ELSE CALL advec_s_pw( nr ) ENDIF ELSE CALL advec_s_up( nr ) ENDIF ENDIF CALL diffusion_s( nr, & surf_def_h(0)%nrsws, surf_def_h(1)%nrsws, & surf_def_h(2)%nrsws, & surf_lsm_h%nrsws, surf_usm_h%nrsws, & surf_def_v(0)%nrsws, surf_def_v(1)%nrsws, & surf_def_v(2)%nrsws, surf_def_v(3)%nrsws, & surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws, & surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws, & surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws, & surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws ) ! !-- Prognostic equation for rain drop concentration DO i = nxl, nxr DO j = nys, nyn !following directive is required to vectorize on Intel19 !DIR$ IVDEP DO k = nzb+1, nzt nr_p(k,j,i) = nr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & tsc(3) * tnr_m(k,j,i) ) & - tsc(5) * rdf_sc(k) * & nr(k,j,i) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & ) IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp ENDDO ENDDO ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tnr_m(k,j,i) = tend(k,j,i) ENDDO ENDDO ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) & + 5.3125_wp * tnr_m(k,j,i) ENDDO ENDDO ENDDO ENDIF ENDIF CALL cpu_log( log_point(53), 'nr-equation', 'stop' ) ENDIF END SUBROUTINE bcm_prognostic_equations !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Control of microphysics for grid points i,j !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_prognostic_equations_ij( i, j, i_omp_start, tn ) INTEGER(iwp), INTENT(IN) :: i !< grid index in x-direction INTEGER(iwp), INTENT(IN) :: i_omp_start !< first loop index of i-loop in !< prognostic_equations INTEGER(iwp), INTENT(IN) :: j !< grid index in y-direction INTEGER(iwp) :: k !< grid index in z-direction INTEGER(iwp), INTENT(IN) :: tn !< task number of openmp task ! !-- If required, calculate prognostic equations for cloud water content and cloud drop !-- concentration. IF ( microphysics_morrison ) THEN ! !-- Calculate prognostic equation for cloud water content tend(:,j,i) = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) & THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( advc_flags_s, i, j, qc, 'qc', flux_s_qc, & diss_s_qc, flux_l_qc, diss_l_qc, & i_omp_start, tn, & bc_dirichlet_l .OR. bc_radiation_l, & bc_dirichlet_n .OR. bc_radiation_n, & bc_dirichlet_r .OR. bc_radiation_r, & bc_dirichlet_s .OR. bc_radiation_s ) ELSE CALL advec_s_pw( i, j, qc ) ENDIF ELSE CALL advec_s_up( i, j, qc ) ENDIF CALL diffusion_s( i, j, qc, & surf_def_h(0)%qcsws, surf_def_h(1)%qcsws, & surf_def_h(2)%qcsws, & surf_lsm_h%qcsws, surf_usm_h%qcsws, & surf_def_v(0)%qcsws, surf_def_v(1)%qcsws, & surf_def_v(2)%qcsws, surf_def_v(3)%qcsws, & surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws, & surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws, & surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws, & surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws ) ! !-- Prognostic equation for cloud water content DO k = nzb+1, nzt qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * & ( tsc(2) * tend(k,j,i) + & tsc(3) * tqc_m(k,j,i) ) & - tsc(5) * rdf_sc(k) & * qc(k,j,i) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & ) IF ( qc_p(k,j,i) < 0.0_wp ) qc_p(k,j,i) = 0.0_wp ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO k = nzb+1, nzt tqc_m(k,j,i) = tend(k,j,i) ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO k = nzb+1, nzt tqc_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tqc_m(k,j,i) ENDDO ENDIF ENDIF ! !-- Calculate prognostic equation for cloud drop concentration. tend(:,j,i) = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( advc_flags_s, i, j, nc, 'nc', flux_s_nc, & diss_s_nc, flux_l_nc, diss_l_nc, & i_omp_start, tn, & bc_dirichlet_l .OR. bc_radiation_l, & bc_dirichlet_n .OR. bc_radiation_n, & bc_dirichlet_r .OR. bc_radiation_r, & bc_dirichlet_s .OR. bc_radiation_s ) ELSE CALL advec_s_pw( i, j, nc ) ENDIF ELSE CALL advec_s_up( i, j, nc ) ENDIF CALL diffusion_s( i, j, nc, & surf_def_h(0)%ncsws, surf_def_h(1)%ncsws, & surf_def_h(2)%ncsws, & surf_lsm_h%ncsws, surf_usm_h%ncsws, & surf_def_v(0)%ncsws, surf_def_v(1)%ncsws, & surf_def_v(2)%ncsws, surf_def_v(3)%ncsws, & surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws, & surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws, & surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws, & surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws ) ! !-- Prognostic equation for cloud drop concentration DO k = nzb+1, nzt nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * & ( tsc(2) * tend(k,j,i) + & tsc(3) * tnc_m(k,j,i) ) & - tsc(5) * rdf_sc(k) & * nc(k,j,i) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & ) IF ( nc_p(k,j,i) < 0.0_wp ) nc_p(k,j,i) = 0.0_wp ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO k = nzb+1, nzt tnc_m(k,j,i) = tend(k,j,i) ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO k = nzb+1, nzt tnc_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tnc_m(k,j,i) ENDDO ENDIF ENDIF ENDIF ! !-- If required, calculate prognostic equations for ice crystal mixing ratio and ice crystal !-- concentration. IF ( microphysics_ice_phase ) THEN ! !-- Calculate prognostic equation for ice crystal mixing ratio tend(:,j,i) = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) & THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( advc_flags_s, i, j, qi, 'qi', flux_s_qi, & diss_s_qi, flux_l_qi, diss_l_qi, & i_omp_start, tn, & bc_dirichlet_l .OR. bc_radiation_l, & bc_dirichlet_n .OR. bc_radiation_n, & bc_dirichlet_r .OR. bc_radiation_r, & bc_dirichlet_s .OR. bc_radiation_s ) ELSE CALL advec_s_pw( i, j, qi ) ENDIF ELSE CALL advec_s_up( i, j, qi ) ENDIF CALL diffusion_s( i, j, qi, & surf_def_h(0)%qisws, surf_def_h(1)%qisws, & surf_def_h(2)%qisws, & surf_lsm_h%qisws, surf_usm_h%qisws, & surf_def_v(0)%qisws, surf_def_v(1)%qisws, & surf_def_v(2)%qisws, surf_def_v(3)%qisws, & surf_lsm_v(0)%qisws, surf_lsm_v(1)%qisws, & surf_lsm_v(2)%qisws, surf_lsm_v(3)%qisws, & surf_usm_v(0)%qisws, surf_usm_v(1)%qisws, & surf_usm_v(2)%qisws, surf_usm_v(3)%qisws ) ! !-- Prognostic equation for ice crystal mixing ratio DO k = nzb+1, nzt qi_p(k,j,i) = qi(k,j,i) + ( dt_3d * & ( tsc(2) * tend(k,j,i) + & tsc(3) * tqi_m(k,j,i) ) & - tsc(5) * rdf_sc(k) & * qi(k,j,i) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & ) IF ( qi_p(k,j,i) < 0.0_wp ) qi_p(k,j,i) = 0.0_wp ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO k = nzb+1, nzt tqi_m(k,j,i) = tend(k,j,i) ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO k = nzb+1, nzt tqi_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tqi_m(k,j,i) ENDDO ENDIF ENDIF ! !-- Calculate prognostic equation for ice crystal concentration. tend(:,j,i) = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( advc_flags_s, i, j, ni, 'ni', flux_s_ni, & diss_s_ni, flux_l_ni, diss_l_ni, & i_omp_start, tn, & bc_dirichlet_l .OR. bc_radiation_l, & bc_dirichlet_n .OR. bc_radiation_n, & bc_dirichlet_r .OR. bc_radiation_r, & bc_dirichlet_s .OR. bc_radiation_s ) ELSE CALL advec_s_pw( i, j, ni ) ENDIF ELSE CALL advec_s_up( i, j, ni ) ENDIF CALL diffusion_s( i, j, ni, & surf_def_h(0)%nisws, surf_def_h(1)%nisws, & surf_def_h(2)%nisws, & surf_lsm_h%nisws, surf_usm_h%nisws, & surf_def_v(0)%nisws, surf_def_v(1)%nisws, & surf_def_v(2)%nisws, surf_def_v(3)%nisws, & surf_lsm_v(0)%nisws, surf_lsm_v(1)%nisws, & surf_lsm_v(2)%nisws, surf_lsm_v(3)%nisws, & surf_usm_v(0)%nisws, surf_usm_v(1)%nisws, & surf_usm_v(2)%nisws, surf_usm_v(3)%nisws ) ! !-- Prognostic equation for ice crystal concentration DO k = nzb+1, nzt ni_p(k,j,i) = ni(k,j,i) + ( dt_3d * & ( tsc(2) * tend(k,j,i) + & tsc(3) * tni_m(k,j,i) ) & - tsc(5) * rdf_sc(k) & * ni(k,j,i) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & ) IF ( ni_p(k,j,i) < 0.0_wp ) ni_p(k,j,i) = 0.0_wp ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO k = nzb+1, nzt tni_m(k,j,i) = tend(k,j,i) ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO k = nzb+1, nzt tni_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tni_m(k,j,i) ENDDO ENDIF ENDIF ENDIF ! !-- If required, calculate prognostic equations for rain water content and rain drop !-- concentration IF ( microphysics_seifert ) THEN ! !-- Calculate prognostic equation for rain water content tend(:,j,i) = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) & THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( advc_flags_s, i, j, qr, 'qr', flux_s_qr, & diss_s_qr, flux_l_qr, diss_l_qr, & i_omp_start, tn, & bc_dirichlet_l .OR. bc_radiation_l, & bc_dirichlet_n .OR. bc_radiation_n, & bc_dirichlet_r .OR. bc_radiation_r, & bc_dirichlet_s .OR. bc_radiation_s ) ELSE CALL advec_s_pw( i, j, qr ) ENDIF ELSE CALL advec_s_up( i, j, qr ) ENDIF CALL diffusion_s( i, j, qr, & surf_def_h(0)%qrsws, surf_def_h(1)%qrsws, & surf_def_h(2)%qrsws, & surf_lsm_h%qrsws, surf_usm_h%qrsws, & surf_def_v(0)%qrsws, surf_def_v(1)%qrsws, & surf_def_v(2)%qrsws, surf_def_v(3)%qrsws, & surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws, & surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws, & surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws, & surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws ) ! !-- Prognostic equation for rain water content DO k = nzb+1, nzt qr_p(k,j,i) = qr(k,j,i) + ( dt_3d * & ( tsc(2) * tend(k,j,i) + & tsc(3) * tqr_m(k,j,i) ) & - tsc(5) * rdf_sc(k) & * qr(k,j,i) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & ) IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO k = nzb+1, nzt tqr_m(k,j,i) = tend(k,j,i) ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO k = nzb+1, nzt tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tqr_m(k,j,i) ENDDO ENDIF ENDIF ! !-- Calculate prognostic equation for rain drop concentration. tend(:,j,i) = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( advc_flags_s, i, j, nr, 'nr', flux_s_nr, & diss_s_nr, flux_l_nr, diss_l_nr, & i_omp_start, tn, & bc_dirichlet_l .OR. bc_radiation_l, & bc_dirichlet_n .OR. bc_radiation_n, & bc_dirichlet_r .OR. bc_radiation_r, & bc_dirichlet_s .OR. bc_radiation_s ) ELSE CALL advec_s_pw( i, j, nr ) ENDIF ELSE CALL advec_s_up( i, j, nr ) ENDIF CALL diffusion_s( i, j, nr, & surf_def_h(0)%nrsws, surf_def_h(1)%nrsws, & surf_def_h(2)%nrsws, & surf_lsm_h%nrsws, surf_usm_h%nrsws, & surf_def_v(0)%nrsws, surf_def_v(1)%nrsws, & surf_def_v(2)%nrsws, surf_def_v(3)%nrsws, & surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws, & surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws, & surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws, & surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws ) ! !-- Prognostic equation for rain drop concentration DO k = nzb+1, nzt nr_p(k,j,i) = nr(k,j,i) + ( dt_3d * & ( tsc(2) * tend(k,j,i) + & tsc(3) * tnr_m(k,j,i) ) & - tsc(5) * rdf_sc(k) & * nr(k,j,i) & ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & ) IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO k = nzb+1, nzt tnr_m(k,j,i) = tend(k,j,i) ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO k = nzb+1, nzt tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tnr_m(k,j,i) ENDDO ENDIF ENDIF ENDIF END SUBROUTINE bcm_prognostic_equations_ij !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Swapping of timelevels !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_swap_timelevel ( mod_count ) IMPLICIT NONE INTEGER, INTENT(IN) :: mod_count IF ( bulk_cloud_model ) THEN 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 IF ( microphysics_ice_phase ) THEN qi => qi_1; qi_p => qi_2 ni => ni_1; ni_p => ni_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 IF ( microphysics_ice_phase ) THEN qi => qi_2; qi_p => qi_1 ni => ni_2; ni_p => ni_1 ENDIF END SELECT ENDIF END SUBROUTINE bcm_swap_timelevel !--------------------------------------------------------------------------------------------------! ! Description: Boundary conditions of the bulk cloud module variables !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_boundary_conditions IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< INTEGER(iwp) :: m !< INTEGER(iwp) :: l !< IF ( microphysics_morrison ) THEN ! !-- Surface conditions cloud water (Dirichlet). !-- Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate !-- belongs to the atmospheric grid point, therefore, set qr_p and nr_p at upward (k-1) and !-- downward-facing (k+1) walls. DO l = 0, 1 !$OMP PARALLEL DO PRIVATE( i, j, k ) DO m = 1, bc_h(l)%ns i = bc_h(l)%i(m) j = bc_h(l)%j(m) k = bc_h(l)%k(m) qc_p(k+bc_h(l)%koff,j,i) = 0.0_wp nc_p(k+bc_h(l)%koff,j,i) = 0.0_wp ENDDO ENDDO ! !-- Top boundary condition for cloud water (Dirichlet) qc_p(nzt+1,:,:) = 0.0_wp nc_p(nzt+1,:,:) = 0.0_wp ENDIF IF ( microphysics_ice_phase ) THEN ! !-- Surface conditions ice crysral (Dirichlet). !-- Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate !-- belongs to the atmospheric grid point, therefore, set qr_p and nr_p at upward (k-1) and !-- downward-facing (k+1) walls DO l = 0, 1 !$OMP PARALLEL DO PRIVATE( i, j, k ) DO m = 1, bc_h(l)%ns i = bc_h(l)%i(m) j = bc_h(l)%j(m) k = bc_h(l)%k(m) qi_p(k+bc_h(l)%koff,j,i) = 0.0_wp ni_p(k+bc_h(l)%koff,j,i) = 0.0_wp ENDDO ENDDO ! !-- Top boundary condition for ice crystal (Dirichlet) qi_p(nzt+1,:,:) = 0.0_wp ni_p(nzt+1,:,:) = 0.0_wp ENDIF IF ( microphysics_seifert ) THEN ! !-- Surface conditions rain water (Dirichlet). !-- Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate !-- belongs to the atmospheric grid point, therefore, set qr_p and nr_p at upward (k-1) and !-- downward-facing (k+1) walls DO l = 0, 1 !$OMP PARALLEL DO PRIVATE( i, j, k ) DO m = 1, bc_h(l)%ns i = bc_h(l)%i(m) j = bc_h(l)%j(m) k = bc_h(l)%k(m) qr_p(k+bc_h(l)%koff,j,i) = 0.0_wp nr_p(k+bc_h(l)%koff,j,i) = 0.0_wp ENDDO ENDDO ! !-- Top boundary condition for rain water (Dirichlet) qr_p(nzt+1,:,:) = 0.0_wp nr_p(nzt+1,:,:) = 0.0_wp ENDIF ! !-- Lateral boundary conditions for scalar quantities at the outflow. !-- Lateral oundary conditions for TKE and dissipation are set in tcm_boundary_conds. IF ( bc_radiation_s ) THEN IF ( microphysics_morrison ) THEN qc_p(:,nys-1,:) = qc_p(:,nys,:) nc_p(:,nys-1,:) = nc_p(:,nys,:) ENDIF IF ( microphysics_seifert ) THEN qr_p(:,nys-1,:) = qr_p(:,nys,:) nr_p(:,nys-1,:) = nr_p(:,nys,:) ENDIF IF ( microphysics_ice_phase ) THEN qi_p(:,nys-1,:) = qi_p(:,nys,:) ni_p(:,nys-1,:) = ni_p(:,nys,:) ENDIF ELSEIF ( bc_radiation_n ) THEN IF ( microphysics_morrison ) THEN qc_p(:,nyn+1,:) = qc_p(:,nyn,:) nc_p(:,nyn+1,:) = nc_p(:,nyn,:) ENDIF IF ( microphysics_seifert ) THEN qr_p(:,nyn+1,:) = qr_p(:,nyn,:) nr_p(:,nyn+1,:) = nr_p(:,nyn,:) ENDIF IF ( microphysics_ice_phase ) THEN qi_p(:,nyn+1,:) = qi_p(:,nyn,:) ni_p(:,nyn+1,:) = ni_p(:,nyn,:) ENDIF ELSEIF ( bc_radiation_l ) THEN IF ( microphysics_morrison ) THEN qc_p(:,:,nxl-1) = qc_p(:,:,nxl) nc_p(:,:,nxl-1) = nc_p(:,:,nxl) ENDIF IF ( microphysics_seifert ) THEN qr_p(:,:,nxl-1) = qr_p(:,:,nxl) nr_p(:,:,nxl-1) = nr_p(:,:,nxl) ENDIF IF ( microphysics_ice_phase ) THEN qi_p(:,:,nxl-1) = qi_p(:,:,nxl) ni_p(:,:,nxl-1) = ni_p(:,:,nxl) ENDIF ELSEIF ( bc_radiation_r ) THEN IF ( microphysics_morrison ) THEN qc_p(:,:,nxr+1) = qc_p(:,:,nxr) nc_p(:,:,nxr+1) = nc_p(:,:,nxr) ENDIF IF ( microphysics_seifert ) THEN qr_p(:,:,nxr+1) = qr_p(:,:,nxr) nr_p(:,:,nxr+1) = nr_p(:,:,nxr) ENDIF IF ( microphysics_ice_phase ) THEN qi_p(:,:,nxr+1) = qi_p(:,:,nxr) ni_p(:,:,nxr+1) = ni_p(:,:,nxr) ENDIF ENDIF END SUBROUTINE bcm_boundary_conditions !--------------------------------------------------------------------------------------------------! ! ! 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 REAL(wp), PARAMETER :: fill_value = -999.0_wp !< value for the _FillValue attribute 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 :: resorted !< flag if output is already resorted LOGICAL, INTENT(INOUT) :: found !< flag if output variable is found LOGICAL, INTENT(INOUT) :: two_d !< flag parameter that indicates 2D variables !< (horizontal cross sections) 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 ( 'ni_xy', 'ni_xz', 'ni_yz' ) IF ( av == 0 ) THEN to_be_resorted => ni ELSE IF ( .NOT. ALLOCATED( ni_av ) ) THEN ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ni_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => ni_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 ( 'qi_xy', 'qi_xz', 'qi_yz' ) IF ( av == 0 ) THEN to_be_resorted => qi ELSE IF ( .NOT. ALLOCATED( qi_av ) ) THEN ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) qi_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => qi_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_total_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 :: resorted !< flag if output is already resorted LOGICAL, INTENT(INOUT) :: found !< flag if output variable is found 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 ( 'ni' ) IF ( av == 0 ) THEN to_be_resorted => ni ELSE IF ( .NOT. ALLOCATED( ni_av ) ) THEN ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ni_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => ni_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 ( 'qi' ) IF ( av == 0 ) THEN to_be_resorted => qi ELSE IF ( .NOT. ALLOCATED( qi_av ) ) THEN ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) qi_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => qi_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_total_0(k,j,i), flag_nr ) & ) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE bcm_data_output_3d !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Read module-specific global restart data (Fortran binary format). !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_rrd_global_ftn( 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 ( 'na_init' ) READ ( 13 ) na_init CASE ( 'dry_aerosol_radius' ) READ ( 13 ) dry_aerosol_radius CASE ( 'sigma_bulk' ) READ ( 13 ) sigma_bulk CASE ( 'aerosol_bulk' ) READ ( 13 ) aerosol_bulk CASE ( 'curvature_solution_effects_bulk' ) READ ( 13 ) curvature_solution_effects_bulk CASE ( 'microphysics_ice_phase' ) READ ( 13 ) microphysics_ice_phase CASE ( 'ice_crystal_sedimentation' ) READ ( 13 ) ice_crystal_sedimentation CASE ( 'in_init' ) READ ( 13 ) in_init CASE ( 'start_ice_microphysics' ) READ ( 13 ) start_ice_microphysics CASE DEFAULT found = .FALSE. END SELECT END SUBROUTINE bcm_rrd_global_ftn !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Read module-specific global restart data (MPI-IO). !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_rrd_global_mpi CALL rrd_mpi_io( 'c_sedimentation', c_sedimentation ) CALL rrd_mpi_io( 'bulk_cloud_model', bulk_cloud_model ) CALL rrd_mpi_io( 'cloud_scheme', cloud_scheme ) CALL rrd_mpi_io( 'cloud_water_sedimentation', cloud_water_sedimentation ) CALL rrd_mpi_io( 'collision_turbulence', collision_turbulence ) CALL rrd_mpi_io( 'limiter_sedimentation', limiter_sedimentation ) CALL rrd_mpi_io( 'nc_const', nc_const ) CALL rrd_mpi_io( 'precipitation', precipitation ) CALL rrd_mpi_io( 'ventilation_effect', ventilation_effect ) CALL rrd_mpi_io( 'na_init', na_init ) CALL rrd_mpi_io( 'dry_aerosol_radius', dry_aerosol_radius ) CALL rrd_mpi_io( 'sigma_bulk', sigma_bulk ) CALL rrd_mpi_io( 'aerosol_bulk', aerosol_bulk ) CALL rrd_mpi_io( 'curvature_solution_effects_bulk', curvature_solution_effects_bulk ) CALL rrd_mpi_io( 'start_ice_microphysics', start_ice_microphysics ) CALL rrd_mpi_io( 'microphysics_ice_phase', microphysics_ice_phase ) CALL rrd_mpi_io( 'in_init', in_init ) CALL rrd_mpi_io( 'ice_crystal_sedimentation', ice_crystal_sedimentation ) END SUBROUTINE bcm_rrd_global_mpi !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Read module-specific local restart data arrays (Fortran binary format). !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_rrd_local_ftn( 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) :: 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 !< 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 ( 'ni' ) IF ( k == 1 ) READ ( 13 ) tmp_3d ni(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'ni_av' ) IF ( .NOT. ALLOCATED( ni_av ) ) THEN ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF IF ( k == 1 ) READ ( 13 ) tmp_3d ni_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 ( 'qi' ) IF ( k == 1 ) READ ( 13 ) tmp_3d qi(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) CASE ( 'qi_av' ) IF ( .NOT. ALLOCATED( qi_av ) ) THEN ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF IF ( k == 1 ) READ ( 13 ) tmp_3d qi_av(:,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_ftn !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Read module-specific local restart data arrays (MPI-IO). !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_rrd_local_mpi LOGICAL :: array_found !< CALL rd_mpi_io_check_array( 'prr' , found = array_found ) IF ( array_found ) THEN IF ( .NOT. ALLOCATED( prr ) ) ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) CALL rrd_mpi_io( 'prr', prr ) ENDIF CALL rd_mpi_io_check_array( 'prr_av' , found = array_found ) IF ( array_found ) THEN IF ( .NOT. ALLOCATED( prr_av ) ) ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) CALL rrd_mpi_io( 'prr_av', prr_av ) ENDIF CALL rd_mpi_io_check_array( 'precipitation_amount' , found = array_found ) IF ( array_found ) THEN IF ( .NOT. ALLOCATED( precipitation_amount ) ) THEN ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg) ) ENDIF CALL rrd_mpi_io( 'precipitation_amount', precipitation_amount ) ENDIF CALL rrd_mpi_io( 'ql', ql ) CALL rd_mpi_io_check_array( 'ql_av' , found = array_found ) IF ( array_found ) THEN IF ( .NOT. ALLOCATED( ql_av ) ) ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) CALL rrd_mpi_io( 'ql_av', ql_av ) ENDIF CALL rrd_mpi_io( 'qc', qc ) CALL rd_mpi_io_check_array( 'qc_av' , found = array_found ) IF ( array_found ) THEN IF ( .NOT. ALLOCATED( qc_av ) ) ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) CALL rrd_mpi_io( 'qc_av', qc_av ) ENDIF IF ( microphysics_morrison ) THEN CALL rrd_mpi_io( 'nc', nc ) CALL rd_mpi_io_check_array( 'nc_av' , found = array_found ) IF ( array_found ) THEN IF ( .NOT. ALLOCATED( nc_av ) ) ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) CALL rrd_mpi_io( 'nc_av', nc_av ) ENDIF ENDIF IF ( microphysics_seifert ) THEN CALL rrd_mpi_io( 'nr', nr ) CALL rd_mpi_io_check_array( 'nr_av' , found = array_found ) IF ( array_found ) THEN IF ( .NOT. ALLOCATED( nr_av ) ) ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) CALL rrd_mpi_io( 'nr_av', nr_av ) ENDIF CALL rrd_mpi_io( 'qr', qr ) CALL rd_mpi_io_check_array( 'qr_av' , found = array_found ) IF ( array_found ) THEN IF ( .NOT. ALLOCATED( qr_av ) ) ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) CALL rrd_mpi_io( 'qr_av', qr_av ) ENDIF ENDIF IF ( microphysics_ice_phase ) THEN CALL rrd_mpi_io( 'ni', ni ) CALL rd_mpi_io_check_array( 'ni_av' , found = array_found ) IF ( array_found ) THEN IF ( .NOT. ALLOCATED( ni_av ) ) ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) CALL rrd_mpi_io( 'ni_av', ni_av ) ENDIF CALL rrd_mpi_io( 'qi', qi ) CALL rd_mpi_io_check_array( 'qi_av' , found = array_found ) IF ( array_found ) THEN IF ( .NOT. ALLOCATED( qi_av ) ) ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) CALL rrd_mpi_io( 'qi_av', qi_av ) ENDIF ENDIF END SUBROUTINE bcm_rrd_local_mpi !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> This routine writes the respective restart data for the bulk cloud module. !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_wrd_global IMPLICIT NONE IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN 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 CALL wrd_write_string( 'na_init' ) WRITE ( 14 ) na_init CALL wrd_write_string( 'dry_aerosol_radius' ) WRITE ( 14 ) dry_aerosol_radius CALL wrd_write_string( 'sigma_bulk' ) WRITE ( 14 ) sigma_bulk CALL wrd_write_string( 'aerosol_bulk' ) WRITE ( 14 ) aerosol_bulk CALL wrd_write_string( 'curvature_solution_effects_bulk' ) WRITE ( 14 ) curvature_solution_effects_bulk CALL wrd_write_string( 'start_ice_microphysics' ) WRITE ( 14 ) start_ice_microphysics CALL wrd_write_string( 'microphysics_ice_phase' ) WRITE ( 14 ) microphysics_ice_phase CALL wrd_write_string( 'in_init' ) WRITE ( 14 ) in_init CALL wrd_write_string( 'ice_crystal_sedimentation' ) WRITE ( 14 ) ice_crystal_sedimentation ELSEIF ( TRIM( restart_data_format_output(1:3) ) == 'mpi' ) THEN CALL wrd_mpi_io( 'c_sedimentation', c_sedimentation ) CALL wrd_mpi_io( 'bulk_cloud_model', bulk_cloud_model ) CALL wrd_mpi_io( 'cloud_scheme', cloud_scheme ) CALL wrd_mpi_io( 'cloud_water_sedimentation', cloud_water_sedimentation ) CALL wrd_mpi_io( 'collision_turbulence', collision_turbulence ) CALL wrd_mpi_io( 'limiter_sedimentation', limiter_sedimentation ) CALL wrd_mpi_io( 'nc_const', nc_const ) CALL wrd_mpi_io( 'precipitation', precipitation ) CALL wrd_mpi_io( 'ventilation_effect', ventilation_effect ) CALL wrd_mpi_io( 'na_init', na_init ) CALL wrd_mpi_io( 'dry_aerosol_radius', dry_aerosol_radius ) CALL wrd_mpi_io( 'sigma_bulk', sigma_bulk ) CALL wrd_mpi_io( 'aerosol_bulk', aerosol_bulk ) CALL wrd_mpi_io( 'curvature_solution_effects_bulk', curvature_solution_effects_bulk ) CALL wrd_mpi_io( 'start_ice_microphysics', start_ice_microphysics ) CALL wrd_mpi_io( 'microphysics_ice_phase', microphysics_ice_phase ) CALL wrd_mpi_io( 'in_init', in_init ) CALL wrd_mpi_io( 'ice_crystal_sedimentation', ice_crystal_sedimentation ) ENDIF END SUBROUTINE bcm_wrd_global !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> This routine writes the respective restart data for the bulk cloud module. !--------------------------------------------------------------------------------------------------! SUBROUTINE bcm_wrd_local IMPLICIT NONE IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN 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_ice_phase ) THEN CALL wrd_write_string( 'ni' ) WRITE ( 14 ) ni IF ( ALLOCATED( ni_av ) ) THEN CALL wrd_write_string( 'ni_av' ) WRITE ( 14 ) ni_av ENDIF CALL wrd_write_string( 'qi' ) WRITE ( 14 ) qi IF ( ALLOCATED( qi_av ) ) THEN CALL wrd_write_string( 'qi_av' ) WRITE ( 14 ) qi_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 ELSEIF ( TRIM( restart_data_format_output(1:3) ) == 'mpi' ) THEN IF ( ALLOCATED( prr ) ) CALL wrd_mpi_io( 'prr', prr ) IF ( ALLOCATED( prr_av ) ) CALL wrd_mpi_io( 'prr_av', prr_av ) IF ( ALLOCATED( precipitation_amount ) ) THEN CALL wrd_mpi_io( 'precipitation_amount', precipitation_amount ) ENDIF CALL wrd_mpi_io( 'ql', ql ) IF ( ALLOCATED( ql_av ) ) CALL wrd_mpi_io( 'ql_av', ql_av ) CALL wrd_mpi_io( 'qc', qc ) IF ( ALLOCATED( qc_av ) ) CALL wrd_mpi_io( 'qc_av', qc_av ) IF ( microphysics_morrison ) THEN CALL wrd_mpi_io( 'nc', nc ) IF ( ALLOCATED( nc_av ) ) CALL wrd_mpi_io( 'nc_av', nc_av ) ENDIF IF ( microphysics_seifert ) THEN CALL wrd_mpi_io( 'nr', nr ) IF ( ALLOCATED( nr_av ) ) CALL wrd_mpi_io( 'nr_av', nr_av ) CALL wrd_mpi_io( 'qr', qr ) IF ( ALLOCATED( qr_av ) ) CALL wrd_mpi_io( 'qr_av', qr_av ) ENDIF IF ( microphysics_ice_phase ) THEN CALL wrd_mpi_io( 'ni', ni ) IF ( ALLOCATED( ni_av ) ) CALL wrd_mpi_io( 'ni_av', ni_av ) CALL wrd_mpi_io( 'qi', qi ) IF ( ALLOCATED( qi_av ) ) CALL wrd_mpi_io( 'qi_av', qi_av ) ENDIF ENDIF END SUBROUTINE bcm_wrd_local !--------------------------------------------------------------------------------------------------! ! 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(50), 'adjust_cloud', 'start' ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) IF ( .NOT. microphysics_morrison_no_rain ) THEN 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 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(50), 'adjust_cloud', 'stop' ) END SUBROUTINE adjust_cloud !--------------------------------------------------------------------------------------------------! ! 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_total_0(k,j,i), 0 ) ) IF ( .NOT. microphysics_morrison_no_rain ) THEN 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 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 END SUBROUTINE adjust_cloud_ij !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Adjust number of ice crystal to avoid nonlinear effects in sedimentation and evaporation of ice !> crytals due to too small or too big weights of ice crytals (Stevens and Seifert, 2008). !--------------------------------------------------------------------------------------------------! SUBROUTINE adjust_ice 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(97), 'adjust_ice', 'start' ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) IF ( qi(k,j,i) <= ximin ) THEN qi(k,j,i) = 0.0_wp ni(k,j,i) = 0.0_wp ELSE IF ( ni(k,j,i) * ximin > qi(k,j,i) * hyrho(k) ) THEN ni(k,j,i) = qi(k,j,i) * hyrho(k) / ximin * flag ENDIF ENDIF ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(97), 'adjust_ice', 'stop' ) END SUBROUTINE adjust_ice !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Adjust number of of ice crystal to avoid nonlinear effects in sedimentation and evaporation of !> ice crystals due to too small or too big weights of ice crytals (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_ice_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_total_0(k,j,i), 0 ) ) IF ( qi(k,j,i) <= ximin ) THEN qi(k,j,i) = 0.0_wp ni(k,j,i) = 0.0_wp ELSE IF ( ni(k,j,i) * ximin > qi(k,j,i) * hyrho(k) ) THEN ni(k,j,i) = qi(k,j,i) * hyrho(k) / ximin * flag ENDIF ENDIF ENDDO END SUBROUTINE adjust_ice_ij !--------------------------------------------------------------------------------------------------! ! 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) :: flag !< flag to indicate first grid level above 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 !< CALL cpu_log( log_point_s(65), 'activation', 'start' ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation located in subroutine 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 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) :: beta_act !< REAL(wp) :: bfactor !< 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) :: s_0 !< REAL(wp) :: sat_max !< REAL(wp) :: sigma !< REAL(wp) :: sigma_act !< DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation 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_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_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 ) / 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 ice nucleation by applying the deposition-condensation formula as given by !> Meyers et al 1992 and as described in Seifert and Beheng 2006 !--------------------------------------------------------------------------------------------------! SUBROUTINE ice_nucleation INTEGER(iwp) :: i !< loop index INTEGER(iwp) :: j !< loop index INTEGER(iwp) :: k !< loop index LOGICAL :: isdac = .TRUE. REAL(wp) :: a_m92 = -0.639_wp !< parameter for nucleation REAL(wp) :: b_m92 = 12.96_wp !< parameter for nucleation REAL(wp) :: flag !< flag to indicate first grid level REAL(wp) :: n_in !< number of ice nucleii REAL(wp) :: nucle = 0.0_wp !< nucleation rate CALL cpu_log( log_point_s(93), 'ice_nucleation', 'start' ) IF ( isdac ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation located in subroutine CALL supersaturation_ice ( i, j, k ) nucle = 0.0_wp IF ( sat_ice >= 0.05_wp .OR. ql(k,j,i) >= 0.001E-3_wp ) THEN ! !-- Calculate ice nucleation nucle = MAX( ( in_init - ni(k,j,i) ) / dt_micro, 0.0_wp ) ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), in_init) ENDIF ENDDO ENDDO ENDDO ELSE DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation located in subroutine CALL supersaturation_ice ( i, j, k ) nucle = 0.0_wp IF ( sat_ice > 0.0 ) THEN ! !-- Calculate ice nucleation n_in = in_init * EXP( a_m92 + b_m92 * sat_ice ) nucle = MAX( ( n_in - ni(k,j,i) ) / dt_micro, 0.0_wp ) ENDIF ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), 1.0E10_wp ) ENDDO ENDDO ENDDO ENDIF CALL cpu_log( log_point_s(93), 'ice_nucleation', 'stop' ) END SUBROUTINE ice_nucleation !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Calculate ice nucleation by applying the deposition-condensation formula as given by !> Meyers et al 1992 and as described in Seifert and Beheng 2006. !--------------------------------------------------------------------------------------------------! SUBROUTINE ice_nucleation_ij( i, j ) INTEGER(iwp) :: i !< loop index INTEGER(iwp) :: j !< loop index INTEGER(iwp) :: k !< loop index LOGICAL :: isdac = .TRUE. REAL(wp) :: a_m92 = -0.639_wp !< parameter for nucleation REAL(wp) :: b_m92 = 12.96_wp !< parameter for nucleation REAL(wp) :: flag !< flag to indicate first grid level REAL(wp) :: n_in !< number of ice nucleii REAL(wp) :: nucle = 0.0_wp !< nucleation rate IF ( isdac ) THEN DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation located in subroutine CALL supersaturation_ice ( i, j, k ) nucle = 0.0_wp IF ( sat_ice >= 0.05_wp .OR. ql(k,j,i) >= 0.001E-3_wp ) THEN ! !-- Calculate ice nucleation nucle = MAX( ( in_init - ni(k,j,i) ) / dt_micro, 0.0_wp ) ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), in_init) ENDIF ENDDO ELSE DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation located in subroutine CALL supersaturation_ice ( i, j, k ) nucle = 0.0_wp IF ( sat_ice > 0.0 ) THEN ! !-- Calculate ice nucleation n_in = in_init * EXP( a_m92 + b_m92 * sat_ice ) nucle = MAX( ( n_in - ni(k,j,i) ) / dt_micro, 0.0_wp ) ENDIF ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), 1.0E10_wp ) ENDDO ENDIF END SUBROUTINE ice_nucleation_ij !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Calculate condensation rate for cloud water content (after Khairoutdinov and Kogan, 2000). !--------------------------------------------------------------------------------------------------! SUBROUTINE condensation IMPLICIT NONE INTEGER(iwp) :: i !< loop index INTEGER(iwp) :: j !< loop index INTEGER(iwp) :: k !< loop index REAL(wp) :: alpha_rc = 1.0_wp !< Tuning parameter see (Seifert and Stevens, 2010) REAL(wp) :: cond !< condensation rate REAL(wp) :: cond_max !< maximum condensation rate REAL(wp) :: evap !< evaporation rate REAL(wp) :: flag !< flag to indicate first grid level above surface REAL(wp) :: g_fac !< factor 1 / Fk + Fd REAL(wp) :: nu = 1.0_wp !< Shape parameter of gernerlized gama distribution REAL(wp) :: rc !< mean cloud droplets radius assuming gamma distribution REAL(wp) :: temp !< actual temperature CALL cpu_log( log_point_s(66), 'condensation', 'start' ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation CALL supersaturation ( i, j, k ) ! !-- Actual temperature, t_l is calculated directly before in supersaturation IF ( microphysics_ice_phase ) THEN temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i) ELSE temp = t_l + lv_d_cp * ql(k,j,i) ENDIF 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 ! !-- Calculating mean radius of cloud droplets assuming gamma distribution with shape !-- parameter nu=1 (Seifert and Beheng, 2006). Tuning factor alpha_rc (introduced !-- in Seifert and Stevens, 2010 ) is switched off. Minimum radius is set to 1µm !-- following (Seifert and Beheng, 2006, Kogan and Khairoutdinov, 2000, !-- Seifert and Stevens, 2010) rc = MAX( alpha_rc * GAMMA( nu + 1.33_wp ) / GAMMA( nu + 1.0_wp ) * & ( 3.0_wp * qc(k,j,i) / & ( 4.0_wp * pi * rho_l * ( nu + 2.0_wp ) * nc(k,j,i) ) & )**0.33_wp, 1.0E-6_wp ) ! !-- Condensation needs only to be calculated in supersaturated regions IF ( sat > 0.0_wp ) THEN ! !-- Condensation rate of cloud water content after KK scheme. !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128, Morrison and Grabowski, !-- 2007, and Seifert and Stevens, 2010) cond = 4.0_wp * pi * nc(k,j,i) * g_fac * sat * rc / hyrho(k) IF ( microphysics_seifert ) THEN cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i) ELSEIF ( microphysics_morrison_no_rain ) THEN cond_max = q(k,j,i) - q_s - qc(k,j,i) ENDIF 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 = 4.0_wp * pi * nc(k,j,i) * g_fac * sat * rc / 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: ! ------------ !> Calculate condensation rate for cloud water content (after Khairoutdinov and Kogan, 2000). !--------------------------------------------------------------------------------------------------! SUBROUTINE condensation_ij( i, j ) IMPLICIT NONE INTEGER(iwp) :: i !< loop index INTEGER(iwp) :: j !< loop index INTEGER(iwp) :: k !< loop index REAL(wp) :: alpha_rc = 1.0_wp !< Tuning parameter see (Seifert and Stevens, 2010) REAL(wp) :: cond !< condensation rate REAL(wp) :: cond_max !< maximum condensation rate REAL(wp) :: evap !< evaporation rate REAL(wp) :: flag !< flag to indicate first grid level above surface REAL(wp) :: g_fac !< factor 1 / Fk + Fd REAL(wp) :: nu = 1.0_wp !< Shape parameter of gernerlized gama distribution REAL(wp) :: rc !< mean cloud droplets radius assuming gamma distribution REAL(wp) :: temp !< actual temperature DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation CALL supersaturation ( i, j, k ) ! !-- Actual temperature, t_l is calculated directly before in supersaturation IF ( microphysics_ice_phase ) THEN temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i) ELSE temp = t_l + lv_d_cp * ql(k,j,i) ENDIF 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 ! !-- Calculating mean radius of cloud droplets assuming gamma distribution with shape !-- parameter nu=1 (Seifert and Beheng, 2006). Tuning factor alpha_rc (introduced !-- in Seifert and Stevens, 2010 ) is switched off. Minimum radius is set to 1µm following !-- (Seifert and Beheng, 2006, Kogan and Khairoutdinov, 2000, Seifert and Stevens, 2010) rc = MAX( alpha_rc * GAMMA( nu + 1.33_wp ) / GAMMA( nu + 1.0_wp ) * & ( 3.0_wp * qc(k,j,i) / ( 4.0_wp * pi * rho_l * ( nu + 2.0_wp ) * nc(k,j,i) ) & )**0.33_wp, 1.0E-6_wp ) ! !-- Condensation needs only to be calculated in supersaturated regions IF ( sat > 0.0_wp ) THEN ! !-- Condensation rate of cloud water content after KK scheme. !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128, Morrison and Grabowski, 2007 !-- and Seifert and Stevens, 2010) cond = 4.0_wp * pi * nc(k,j,i) * g_fac * sat * rc / hyrho(k) IF ( microphysics_seifert ) THEN cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i) ELSEIF ( microphysics_morrison_no_rain ) THEN cond_max = q(k,j,i) - q_s - qc(k,j,i) ENDIF 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 = 4.0_wp * pi * nc(k,j,i) * g_fac * sat * rc / 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: ! ------------ !> Calculate the growth of ice particles by water vapor deposition (after Seifert and Beheng, 2006). !--------------------------------------------------------------------------------------------------! SUBROUTINE ice_deposition INTEGER(iwp) :: i !< loop index INTEGER(iwp) :: j !< loop index INTEGER(iwp) :: k !< loop index REAL(wp) :: ac = 0.09_wp !< parameter for ice capacitance REAL(wp) :: bc = 0.33_wp !< parameter for ice capacitance REAL(wp) :: deposition_rate !< depositions rate REAL(wp) :: deposition_rate_max !< maximum deposition rate REAL(wp) :: fac_gamma = 0.76_wp !< parameter to describe spectral distribution, here !< following gamma size distribution with µ =1/3 and nu=0 REAL(wp) :: flag !< flag to indicate first grid level above REAL(wp) :: gfac_dep !< factor REAL(wp) :: sublimation_rate !< sublimations rate REAL(wp) :: temp !< actual temperature REAL(wp) :: xi !< mean mass of ice crystal CALL cpu_log( log_point_s(95), 'ice deposition', 'start' ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation over a plane ice surface CALL supersaturation_ice ( i, j, k ) ! !-- Actual temperature: temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i) IF ( temp >= 273.15_wp ) CYCLE ! !-- calculating gfac_dep ( 1/ (Fk + Fd) ) see e.g. Rogers and Yau, 1989 gfac_dep = 1.0_wp / ( ( l_s / ( r_v * temp ) - 1.0_wp ) * & l_s / ( thermal_conductivity_l * temp ) & + r_v * temp / ( diff_coeff_l * e_si ) & ) ! !-- If there is nothing nucleated, than there is also no deposition (above -38°C) IF ( ni(k,j,i) <= 0.0_wp ) CYCLE ! !-- calculate mean mass of ice crystal xi = 0.0_wp xi = MAX( ( qi(k,j,i) * hyrho(k) / ni(k,j,i)), ximin ) ! !-- Condensation needs only to be calculated in supersaturated regions (regarding ice) IF ( sat_ice > 0.0_wp ) THEN ! !-- Calculate deposition rate assuming ice crystal shape as prescribed in !-- Ovchinnikov et al., 2014 and a gamma size distribution according to !-- Seifert and Beheng with µ =1/3 and nu=0 deposition_rate = 4.0_wp * pi * sat_ice * gfac_dep * fac_gamma * ac * xi**bc * & ni(k,j,i) IF ( microphysics_seifert ) THEN deposition_rate_max = q(k,j,i) - q_si - qr(k,j,i) - qi(k,j,i) ELSEIF ( microphysics_morrison_no_rain ) THEN deposition_rate_max = q(k,j,i) - q_si - qc(k,j,i) - qi(k,j,i) ENDIF deposition_rate = MIN( deposition_rate, deposition_rate_max / dt_micro ) qi(k,j,i) = qi(k,j,i) + deposition_rate * dt_micro * flag ELSEIF ( sat_ice < 0.0_wp ) THEN sublimation_rate = 4.0_wp * pi * sat_ice * gfac_dep * fac_gamma * ac * xi**bc * & ni(k,j,i) sublimation_rate = MAX( sublimation_rate, - qi(k,j,i) / dt_micro ) qi(k,j,i) = qi(k,j,i) + sublimation_rate * dt_micro * flag ENDIF ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(95), 'ice deposition', 'stop' ) END SUBROUTINE ice_deposition !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Calculate condensation rate for cloud water content (after Khairoutdinov and Kogan, 2000). !--------------------------------------------------------------------------------------------------! SUBROUTINE ice_deposition_ij( i, j ) INTEGER(iwp) :: i !< loop index INTEGER(iwp) :: j !< loop index INTEGER(iwp) :: k !< loop index REAL(wp) :: ac = 0.09_wp !< parameter for ice capacitance REAL(wp) :: bc = 0.33_wp !< parameter for ice capacitance REAL(wp) :: deposition_rate !< depositions rate REAL(wp) :: deposition_rate_max !< maximum deposition rate REAL(wp) :: fac_gamma = 0.76_wp !< parameter to describe spectral distribution, here !< following gamma size distribution with nu=1, v=0 REAL(wp) :: flag !< flag to indicate first grid level above REAL(wp) :: gfac_dep !< factor REAL(wp) :: sublimation_rate !< sublimations rate REAL(wp) :: temp !< actual temperature REAL(wp) :: xi !< mean mass of ice crystal DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation over a plane ice surface CALL supersaturation_ice (i,j,k) ! !-- Actual temperature: temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i) IF ( temp >= 273.15_wp ) CYCLE ! !-- calculating gfac_dep ( 1/ (Fk + Fd) ) see e.g. !-- Rogers and Yau, 1989 gfac_dep = 1.0_wp / ( ( l_s / ( r_v * temp ) - 1.0_wp ) * & l_s / ( thermal_conductivity_l * temp ) & + r_v * temp / ( diff_coeff_l * e_si ) & ) ! !-- If there is nothing nucleated, than there is also no deposition (above -38°C) IF ( ni(k,j,i) <= 0.0_wp ) CYCLE ! !-- calculate mean mass of ice crystal xi = 0.0_wp xi = MAX( ( qi(k,j,i) * hyrho(k) / ni(k,j,i)), ximin ) ! !-- Condensation needs only to be calculated in supersaturated regions (regarding ice) IF ( sat_ice > 0.0_wp ) THEN ! !-- Calculate deposition rate assuming ice crystal shape as prescribed in !-- Ovchinnikov et al., 2014 and a gamma size distribution according to !-- Seifert and Beheng with µ =1/3 and nu=0 deposition_rate = 4.0_wp * pi * sat_ice * gfac_dep * fac_gamma * ac * xi**bc * & ni(k,j,i) IF ( microphysics_seifert ) THEN deposition_rate_max = q(k,j,i) - q_si - qr(k,j,i) - qi(k,j,i) ELSEIF ( microphysics_morrison_no_rain ) THEN deposition_rate_max = q(k,j,i) - q_si - qc(k,j,i) - qi(k,j,i) ENDIF deposition_rate = MIN( deposition_rate, deposition_rate_max / dt_micro ) qi(k,j,i) = qi(k,j,i) + deposition_rate * dt_micro * flag ELSEIF ( sat_ice < 0.0_wp ) THEN sublimation_rate = 4.0_wp * pi * sat_ice * gfac_dep * fac_gamma * ac * xi**bc * & ni(k,j,i) sublimation_rate = MAX( sublimation_rate, - qi(k,j,i) / dt_micro ) qi(k,j,i) = qi(k,j,i) + sublimation_rate * dt_micro * flag ENDIF ENDDO END SUBROUTINE ice_deposition_ij !--------------------------------------------------------------------------------------------------! ! 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(47), 'autoconversion', 'start' ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_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(47), 'autoconversion', 'stop' ) END SUBROUTINE autoconversion !--------------------------------------------------------------------------------------------------! ! 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_total_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 = 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 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 = nxl, nxr DO j = nys, nyn ! !-- Determine vertical index of topography top k_wall = topo_top_ind(j,i,0) DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_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: ! ------------ !> 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 = topo_top_ind(j,i,0) DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_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). !--------------------------------------------------------------------------------------------------! 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 = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_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: ! ------------ !> 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_total_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 END SUBROUTINE accretion_ij !--------------------------------------------------------------------------------------------------! ! 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 = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_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: ! ------------ !> 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_total_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. !--------------------------------------------------------------------------------------------------! 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 = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) IF ( qr(k,j,i) > eps_sb ) THEN ! !-- Call calculation of supersaturation 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: ! ------------ !> 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) :: dr !< 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) :: 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_total_0(k,j,i), 0 ) ) IF ( qr(k,j,i) > eps_sb ) THEN ! !-- Call calculation of supersaturation 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 END SUBROUTINE evaporation_rain_ij !--------------------------------------------------------------------------------------------------! ! 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_nc !< REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !< 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 = nxl, nxr DO j = nys, nyn DO k = nzt, nzb+1, -1 ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_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: ! ------------ !> 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_total_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 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: ! ------------ !> Sedimentation of ice crystals !--------------------------------------------------------------------------------------------------! SUBROUTINE sedimentation_ice INTEGER(iwp) :: i !< loop index INTEGER(iwp) :: j !< loop index INTEGER(iwp) :: k !< loop index REAL(wp) :: av = 6.39_wp !< parameter for calculating fall speed REAL(wp) :: bv = 0.1666_wp !< parameter (1/6) REAL(wp) :: factor_sed_gamma_k0 = 0.76_wp !< factor for zeroth moment and µ =1/3 and nu=0 REAL(wp) :: factor_sed_gamma_k1 = 1.61_wp !< factor for first moment and µ =1/3 and nu=0 REAL(wp) :: flag !< flag to indicate first grid level REAL(wp) :: xi = 0.0_wp !< mean mass of ice crystal REAL(wp) :: vi = 0.0_wp !< mean fall speed of ice crystal REAL(wp), DIMENSION(nzb:nzt+1) :: sed_ni !< sedimentation rate zeroth moment REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qi !< sedimentation rate fist moment CALL cpu_log( log_point_s(96), 'sed_ice', 'start' ) sed_qi(nzt+1) = 0.0_wp sed_ni(nzt+1) = 0.0_wp DO i = nxl, nxr DO j = nys, nyn DO k = nzt, nzb+1, -1 IF ( ni(k,j,i) <= 0.0_wp ) THEN xi = 0.0_wp ELSE !-- Calculate mean mass of ice crystal xi = MAX( (hyrho(k) * qi(k,j,i) / ni(k,j,i)), ximin) ENDIF ! !-- Calculate fall speed of ice crystal vi = av * xi**bv ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) ! !-- Calculate sedimentation rate for each grid box, factors are calculated using !> explanation missing eventually IF ( qi(k,j,i) > eps_sb .AND. ni(k,j,i) >= 0.0_wp ) THEN sed_qi(k) = qi(k,j,i) * vi * factor_sed_gamma_k1 * flag sed_ni(k) = ni(k,j,i) * vi * factor_sed_gamma_k0 * flag ELSE sed_qi(k) = 0.0_wp sed_ni(k) = 0.0_wp ENDIF ! !-- Calculate sedimentation: divergence of sedimentation flux sed_qi(k) = MIN( sed_qi(k), hyrho(k) * dzu(k+1) * q(k,j,i) / dt_micro + sed_qi(k+1)& ) * flag q(k,j,i) = q(k,j,i) + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) & / hyrho(k) * dt_micro * flag qi(k,j,i) = qi(k,j,i) + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) & / hyrho(k) * dt_micro * flag ni(k,j,i) = ni(k,j,i) + ( sed_ni(k+1) - sed_ni(k) ) * ddzu(k+1) & / hyrho(k) * dt_micro * flag pt(k,j,i) = pt(k,j,i) - ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) & / hyrho(k) * l_s / c_p * 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_qi(k) / hyrho(k) * & weight_substep(intermediate_timestep_count) * flag ELSE prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) * flag ENDIF ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(96), 'sed_ice', 'stop' ) END SUBROUTINE sedimentation_ice !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Sedimentation of ice crystals !--------------------------------------------------------------------------------------------------! SUBROUTINE sedimentation_ice_ij( i, j ) INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< REAL(wp) :: av = 6.39_wp !< parameter for calculating fall speed REAL(wp) :: bv = 0.1666_wp !< REAL(wp) :: factor_sed_gamma_k0 = 0.76_wp REAL(wp) :: factor_sed_gamma_k1 = 1.61_wp REAL(wp) :: flag !< flag to indicate first grid level REAL(wp) :: xi = 0.0_wp !< mean mass of ice crystal REAL(wp) :: vi = 0.0_wp !< mean fall speed of ice crystal REAL(wp), DIMENSION(nzb:nzt+1) :: sed_ni !< REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qi !< sed_qi(nzt+1) = 0.0_wp sed_ni(nzt+1) = 0.0_wp DO k = nzt, nzb+1, -1 IF ( ni(k,j,i) <= 0.0_wp ) THEN xi = 0.0_wp ELSE !-- Calculate mean mass of ice crystal xi = MAX( (hyrho(k) * qi(k,j,i) / ni(k,j,i)), ximin) ENDIF ! !-- calculate fall speed of ice crystal vi = av * xi**bv ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) ! !-- Calculate sedimentation rate for each grid box, factors are !-- calculated using IF ( qi(k,j,i) > eps_sb .AND. ni(k,j,i) >= 0.0_wp ) THEN sed_qi(k) = qi(k,j,i) * vi * factor_sed_gamma_k1 * flag sed_ni(k) = ni(k,j,i) * vi * factor_sed_gamma_k0 * flag ELSE sed_qi(k) = 0.0_wp sed_ni(k) = 0.0_wp ENDIF ! !-- calculate fall speed of ice crystal vi = av * xi**bv ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) IF ( qi(k,j,i) > eps_sb .AND. ni(k,j,i) >= 0.0_wp ) THEN sed_qi(k) = qi(k,j,i) * vi * factor_sed_gamma_k1 * flag sed_ni(k) = ni(k,j,i) * vi * factor_sed_gamma_k0 * flag ELSE sed_qi(k) = 0.0_wp sed_ni(k) = 0.0_wp ENDIF sed_qi(k) = MIN( sed_qi(k), hyrho(k) * dzu(k+1) * q(k,j,i) / dt_micro + sed_qi(k+1) & ) * flag q(k,j,i) = q(k,j,i) + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) / & hyrho(k) * dt_micro * flag qi(k,j,i) = qi(k,j,i) + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) / & hyrho(k) * dt_micro * flag ni(k,j,i) = ni(k,j,i) + ( sed_ni(k+1) - sed_ni(k) ) * ddzu(k+1) / & hyrho(k) * dt_micro * flag pt(k,j,i) = pt(k,j,i) - ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) / & hyrho(k) * l_s / c_p * 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_qi(k) / hyrho(k) * & weight_substep(intermediate_timestep_count) * flag ELSE prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) * flag ENDIF ENDDO END SUBROUTINE sedimentation_ice_ij !--------------------------------------------------------------------------------------------------! ! 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 = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_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_total_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_total_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_total_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 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_total_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_total_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_total_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_total_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: ! ------------ !> 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: ! ------------ !> 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_tl( 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 ) IF ( .NOT. microphysics_ice_phase ) THEN IF ( microphysics_seifert ) THEN sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp ELSEIF ( microphysics_morrison_no_rain ) THEN sat = ( q(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp ENDIF ELSE IF ( microphysics_seifert ) THEN sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) - qi(k,j,i) ) / q_s - 1.0_wp ELSEIF ( microphysics_morrison_no_rain ) THEN sat = ( q(k,j,i) - qc(k,j,i) - qi(k,j,i) ) / q_s - 1.0_wp ENDIF ENDIF END SUBROUTINE supersaturation !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Computation of the diagnostic supersaturation sat, actual liquid water temperature t_l and !> saturation water vapor mixing ratio q_s !--------------------------------------------------------------------------------------------------! SUBROUTINE supersaturation_ice ( i, j, k ) INTEGER(iwp) :: i !< running index INTEGER(iwp) :: j !< running index INTEGER(iwp) :: k !< running index REAL(wp) :: e_a !< water vapor pressure REAL(wp) :: temp !< actual temperature ! !-- Actual liquid water temperature: t_l = exner(k) * pt(k,j,i) ! !-- Actual temperature: temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i) ! !-- Calculate water vapor saturation pressure with formular using actual temperature e_si = magnus_ice( temp ) ! !-- Computation of ice saturation mixing ratio: q_si = rd_d_rv * e_si / ( hyp(k) - e_si ) ! !-- Current vapor pressure e_a = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) - qi(k,j,i) ) * hyp(k) / & ( ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) - qi(k,j,i) ) + rd_d_rv ) ! !-- Supersaturation: !-- Not in case of microphysics_kessler or microphysics_sat_adjust since qr is unallocated sat_ice = e_a / e_si - 1.0_wp END SUBROUTINE supersaturation_ice !--------------------------------------------------------------------------------------------------! ! 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 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_total_0(k,j,i), 0 ) ) ! !-- Call calculation of supersaturation located CALL supersaturation( i, j, k ) ! !-- Compute the liquid water content IF ( .NOT. microphysics_ice_phase ) THEN IF ( microphysics_seifert .AND. .NOT. microphysics_morrison ) THEN !-- Seifert and Beheng scheme: saturation adjustment 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 ! !-- Morrison scheme: explicit condensation (see above) ELSEIF ( microphysics_morrison .AND. microphysics_seifert ) THEN ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag ! !-- Morrison without rain: explicit condensation ELSEIF ( microphysics_morrison .AND. .NOT. microphysics_seifert ) THEN ql(k,j,i) = qc(k,j,i) ! !-- Kessler and saturation adjustment scheme 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 !-- Calculations of liquid water content in case of mixed-phase cloud microphyics ELSE IF ( microphysics_seifert .AND. .NOT. microphysics_morrison ) & THEN ! !-- Seifert and Beheng scheme: saturation adjustment IF ( ( q(k,j,i) - q_s - qr(k,j,i) - qi(k,j,i) ) > 0.0_wp ) THEN qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) - qi(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) + qi(k,j,i) ) ) THEN q(k,j,i) = qr(k,j,i) + qi(k,j,i) ENDIF qc(k,j,i) = 0.0_wp ql(k,j,i) = qr(k,j,i) * flag ENDIF ! !-- Morrison scheme: explicit condensation (see above) ELSEIF ( microphysics_morrison .AND. microphysics_seifert ) THEN ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag ! !-- Morrison without rain: explicit condensation ELSEIF ( microphysics_morrison .AND. .NOT. microphysics_seifert ) THEN ql(k,j,i) = qc(k,j,i) ! !-- Kessler and saturation adjustment scheme 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 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), 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 /) !< REAL(wp) :: gamm !< REAL(wp) :: ser !< REAL(wp) :: tmp !< REAL(wp) :: x_gamm !< REAL(wp) :: xx !< REAL(wp) :: y_gamm !< 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