source: palm/trunk/SOURCE/microphysics_mod.f90 @ 3049

Last change on this file since 3049 was 3026, checked in by schwenkel, 6 years ago

Changed the name specific humidity to mixing ratio

  • Property svn:keywords set to Id
File size: 123.1 KB
RevLine 
[1850]1!> @file microphysics_mod.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1093]4!
[2000]5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
[1093]9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
[2718]17! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1093]19!
[1000]20! Current revisions:
[1092]21! ------------------
[2318]22!
[2701]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: microphysics_mod.f90 3026 2018-05-22 10:30:53Z Giersch $
[3026]27! Changed the name specific humidity to mixing ratio, since we are computing
28! mixing ratios.
29!
30! 2718 2018-01-02 08:49:38Z maronga
[2716]31! Corrected "Former revisions" section
[2701]32!
[2716]33! 2701 2017-12-15 15:40:50Z suehring
34! Changes from last commit documented
35!
[2701]36! 2698 2017-12-14 18:46:24Z suehring
[2716]37! Bugfix in get_topography_top_index
38!
39! 2696 2017-12-14 17:12:51Z kanani
40! Change in file header (GPL part)
41!
42! 2608 2017-11-13 14:04:26Z schwenkel
[2608]43! Calculation of supersaturation in external module (diagnostic_quantities_mod).
44! Change: correct calculation of saturation specific humidity to saturation
45! mixing ratio (the factor of 0.378 vanishes).
46!
47! 2522 2017-10-05 14:20:37Z schwenkel
[2522]48! Minor bugfix
49!
50! 2375 2017-08-29 14:10:28Z schwenkel
[2375]51! Improved aerosol initilization and some minor bugfixes
52! for droplet sedimenation
53!
54! 2318 2017-07-20 17:27:44Z suehring
[2318]55! Get topography top index via Function call
56!
57! 2317 2017-07-20 17:27:19Z suehring
[2312]58! s1 changed to log_sigma
59!
60! 2292 2017-06-20 09:51:42Z schwenkel
61! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
62! includes two more prognostic equations for cloud drop concentration (nc)
63! and cloud water content (qc).
64! - The process of activation is parameterized with a simple Twomey
65!   activion scheme or with considering solution and curvature
[2292]66!   effects (Khvorostyanov and Curry ,2006).
67! - The saturation adjustment scheme is replaced by the parameterization
68!   of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
69! - All other microphysical processes of Seifert and Beheng are used.
[2312]70!   Additionally, in those processes the reduction of cloud number concentration
71!   is considered.
72!
[2292]73! 2233 2017-05-30 18:08:54Z suehring
[1321]74!
[2233]75! 2232 2017-05-30 17:47:52Z suehring
76! Adjustments to new topography and surface concept
[2312]77!
[2156]78! 2155 2017-02-21 09:57:40Z hoffmann
79! Bugfix in the calculation of microphysical quantities on ghost points.
80!
[2032]81! 2031 2016-10-21 15:11:58Z knoop
82! renamed variable rho to rho_ocean
[2155]83!
[2001]84! 2000 2016-08-20 18:09:15Z knoop
85! Forced header and separation lines into 80 columns
[2155]86!
[1851]87! 1850 2016-04-08 13:29:27Z maronga
88! Module renamed
89! Adapted for modularization of microphysics.
90!
[1846]91! 1845 2016-04-08 08:29:13Z raasch
92! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
93! point (instead of one point above surface)
94!
[1832]95! 1831 2016-04-07 13:15:51Z hoffmann
96! turbulence renamed collision_turbulence,
97! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
98! avaialble for microphysics_kessler.
99!
[1823]100! 1822 2016-04-07 07:49:42Z hoffmann
101! Unused variables removed.
102! Kessler scheme integrated.
103!
[1692]104! 1691 2015-10-26 16:17:44Z maronga
105! Added new routine calc_precipitation_amount. The routine now allows to account
106! for precipitation due to sedimenation of cloud (fog) droplets
[2155]107!
[1683]108! 1682 2015-10-07 23:56:08Z knoop
[2155]109! Code annotations made doxygen readable
[1683]110!
[1647]111! 1646 2015-09-02 16:00:10Z hoffmann
112! Bugfix: Wrong computation of d_mean.
113!
[1362]114! 1361 2014-04-16 15:17:48Z hoffmann
115! Bugfix in sedimentation_rain: Index corrected.
116! Vectorized version of adjust_cloud added.
117! Little reformatting of the code.
[2155]118!
[1354]119! 1353 2014-04-08 15:21:23Z heinze
120! REAL constants provided with KIND-attribute
[2155]121!
[1347]122! 1346 2014-03-27 13:18:20Z heinze
[2155]123! Bugfix: REAL constants provided with KIND-attribute especially in call of
[1347]124! intrinsic function like MAX, MIN, SIGN
[2155]125!
[1335]126! 1334 2014-03-25 12:21:40Z heinze
127! Bugfix: REAL constants provided with KIND-attribute
128!
[1323]129! 1322 2014-03-20 16:38:49Z raasch
130! REAL constants defined as wp-kind
131!
[1321]132! 1320 2014-03-20 08:40:49Z raasch
[1320]133! ONLY-attribute added to USE-statements,
134! kind-parameters added to all INTEGER and REAL declaration statements,
135! kinds are defined in new module kinds,
136! comment fields (!:) to be used for variable explanations added to
137! all variable declaration statements
[1000]138!
[1242]139! 1241 2013-10-30 11:36:58Z heinze
[2031]140! hyp and rho_ocean have to be calculated at each time step if data from external
[1242]141! file LSF_DATA are used
142!
[1116]143! 1115 2013-03-26 18:16:16Z hoffmann
144! microphyical tendencies are calculated in microphysics_control in an optimized
145! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
146!
[1107]147! 1106 2013-03-04 05:31:38Z raasch
148! small changes in code formatting
149!
[1093]150! 1092 2013-02-02 11:24:22Z raasch
151! unused variables removed
152! file put under GPL
153!
[1066]154! 1065 2012-11-22 17:42:36Z hoffmann
155! Sedimentation process implemented according to Stevens and Seifert (2008).
[1115]156! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
[1066]157! and Stevens, 2010).
158!
[1054]159! 1053 2012-11-13 17:11:03Z hoffmann
160! initial revision
[2155]161!
[1000]162! Description:
163! ------------
[1849]164!> Calculate bilk cloud microphysics.
[1000]165!------------------------------------------------------------------------------!
[1682]166 MODULE microphysics_mod
[1000]167
[1849]168    USE  kinds
169
170    IMPLICIT NONE
171
[2292]172    LOGICAL ::  cloud_water_sedimentation = .FALSE.       !< cloud water sedimentation
173    LOGICAL ::  curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
174    LOGICAL ::  limiter_sedimentation = .TRUE.            !< sedimentation limiter
175    LOGICAL ::  collision_turbulence = .FALSE.            !< turbulence effects
176    LOGICAL ::  ventilation_effect = .TRUE.               !< ventilation effect
[1849]177
178    REAL(wp) ::  a_1 = 8.69E-4_wp          !< coef. in turb. parametrization (cm-2 s3)
179    REAL(wp) ::  a_2 = -7.38E-5_wp         !< coef. in turb. parametrization (cm-2 s3)
180    REAL(wp) ::  a_3 = -1.40E-2_wp         !< coef. in turb. parametrization
181    REAL(wp) ::  a_term = 9.65_wp          !< coef. for terminal velocity (m s-1)
182    REAL(wp) ::  a_vent = 0.78_wp          !< coef. for ventilation effect
183    REAL(wp) ::  b_1 = 11.45E-6_wp         !< coef. in turb. parametrization (m)
184    REAL(wp) ::  b_2 = 9.68E-6_wp          !< coef. in turb. parametrization (m)
185    REAL(wp) ::  b_3 = 0.62_wp             !< coef. in turb. parametrization
186    REAL(wp) ::  b_term = 9.8_wp           !< coef. for terminal velocity (m s-1)
187    REAL(wp) ::  b_vent = 0.308_wp         !< coef. for ventilation effect
188    REAL(wp) ::  beta_cc = 3.09E-4_wp      !< coef. in turb. parametrization (cm-2 s3)
189    REAL(wp) ::  c_1 = 4.82E-6_wp          !< coef. in turb. parametrization (m)
190    REAL(wp) ::  c_2 = 4.8E-6_wp           !< coef. in turb. parametrization (m)
191    REAL(wp) ::  c_3 = 0.76_wp             !< coef. in turb. parametrization
192    REAL(wp) ::  c_const = 0.93_wp         !< const. in Taylor-microscale Reynolds number
193    REAL(wp) ::  c_evap = 0.7_wp           !< constant in evaporation
194    REAL(wp) ::  c_term = 600.0_wp         !< coef. for terminal velocity (m-1)
195    REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
[2155]196    REAL(wp) ::  eps_sb = 1.0E-10_wp       !< threshold in two-moments scheme
[2292]197    REAL(wp) ::  eps_mr = 0.0_wp           !< threshold for morrison scheme
[1849]198    REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
199    REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
200    REAL(wp) ::  k_rr = 7.12_wp            !< const. rain-rain kernel (m3 kg-1 s-1)
201    REAL(wp) ::  k_br = 1000.0_wp          !< const. in breakup parametrization (m-1)
202    REAL(wp) ::  k_st = 1.2E8_wp           !< const. in drizzle parametrization (m-1 s-1)
203    REAL(wp) ::  kappa_rr = 60.7_wp        !< const. in collision kernel (kg-1/3)
204    REAL(wp) ::  kin_vis_air = 1.4086E-5_wp  !< kin. viscosity of air (m2 s-1)
205    REAL(wp) ::  prec_time_const = 0.001_wp  !< coef. in Kessler scheme (s-1)
206    REAL(wp) ::  ql_crit = 0.0005_wp       !< coef. in Kessler scheme (kg kg-1)
207    REAL(wp) ::  schmidt_p_1d3=0.8921121_wp  !< Schmidt number**0.33333, 0.71**0.33333
208    REAL(wp) ::  sigma_gc = 1.3_wp         !< geometric standard deviation cloud droplets
209    REAL(wp) ::  thermal_conductivity_l = 2.43E-2_wp  !< therm. cond. air (J m-1 s-1 K-1)
210    REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
211    REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
[2292]212    REAL(wp) ::  xamin = 5.24E-19_wp       !< average aerosol mass (kg) (~ 0.05µm)
213    REAL(wp) ::  xcmin = 4.18E-15_wp       !< minimum cloud drop size (kg) (~ 1µm)
214    REAL(wp) ::  xcmax = 2.6E-10_wp        !< maximum cloud drop size (kg) (~ 40µm)
[1849]215    REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
216    REAL(wp) ::  xrmax = 5.0E-6_wp         !< maximum rain drop site (kg)
217
[2375]218    REAL(wp) ::  c_sedimentation = 2.0_wp        !< Courant number of sedimentation process
219    REAL(wp) ::  dpirho_l                        !< 6.0 / ( pi * rho_l )
220    REAL(wp) ::  dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius
221    REAL(wp) ::  dt_micro                        !< microphysics time step
222    REAL(wp) ::  sigma_bulk = 2.0_wp             !< width of aerosol spectrum
223    REAL(wp) ::  na_init = 100.0E6_wp            !< Total particle/aerosol concentration (cm-3)
224    REAL(wp) ::  nc_const = 70.0E6_wp            !< cloud droplet concentration
225    REAL(wp) ::  dt_precipitation = 100.0_wp     !< timestep precipitation (s)
226    REAL(wp) ::  sed_qc_const                    !< const. for sedimentation of cloud water
227    REAL(wp) ::  pirho_l                         !< pi * rho_l / 6.0;
[1849]228
229    REAL(wp), DIMENSION(:), ALLOCATABLE ::  nc_1d  !<
230    REAL(wp), DIMENSION(:), ALLOCATABLE ::  nr_1d  !<
231    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_1d  !<
232    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qc_1d  !<
233    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qr_1d  !<
234    REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_1d   !<
235
236    SAVE
237
[1000]238    PRIVATE
[1849]239    PUBLIC microphysics_control, microphysics_init
[1000]240
[2292]241    PUBLIC cloud_water_sedimentation, collision_turbulence,                    &
[2375]242           curvature_solution_effects_bulk, c_sedimentation,                   &
243           dry_aerosol_radius, dt_precipitation,                               &
244           limiter_sedimentation, na_init, nc_const, sigma_gc, sigma_bulk,     &
[1849]245           ventilation_effect
246
[2312]247
[1115]248    INTERFACE microphysics_control
249       MODULE PROCEDURE microphysics_control
250       MODULE PROCEDURE microphysics_control_ij
251    END INTERFACE microphysics_control
[1022]252
[1115]253    INTERFACE adjust_cloud
254       MODULE PROCEDURE adjust_cloud
255       MODULE PROCEDURE adjust_cloud_ij
256    END INTERFACE adjust_cloud
257
[2292]258    INTERFACE activation
259       MODULE PROCEDURE activation
260       MODULE PROCEDURE activation_ij
261    END INTERFACE activation
262
263    INTERFACE condensation
264       MODULE PROCEDURE condensation
265       MODULE PROCEDURE condensation_ij
266    END INTERFACE condensation
267
[1000]268    INTERFACE autoconversion
269       MODULE PROCEDURE autoconversion
270       MODULE PROCEDURE autoconversion_ij
271    END INTERFACE autoconversion
272
[1822]273    INTERFACE autoconversion_kessler
274       MODULE PROCEDURE autoconversion_kessler
275       MODULE PROCEDURE autoconversion_kessler_ij
276    END INTERFACE autoconversion_kessler
277
[1000]278    INTERFACE accretion
279       MODULE PROCEDURE accretion
280       MODULE PROCEDURE accretion_ij
281    END INTERFACE accretion
[1005]282
283    INTERFACE selfcollection_breakup
284       MODULE PROCEDURE selfcollection_breakup
285       MODULE PROCEDURE selfcollection_breakup_ij
286    END INTERFACE selfcollection_breakup
[1012]287
288    INTERFACE evaporation_rain
289       MODULE PROCEDURE evaporation_rain
290       MODULE PROCEDURE evaporation_rain_ij
291    END INTERFACE evaporation_rain
292
293    INTERFACE sedimentation_cloud
294       MODULE PROCEDURE sedimentation_cloud
295       MODULE PROCEDURE sedimentation_cloud_ij
296    END INTERFACE sedimentation_cloud
[2155]297
[1012]298    INTERFACE sedimentation_rain
299       MODULE PROCEDURE sedimentation_rain
300       MODULE PROCEDURE sedimentation_rain_ij
301    END INTERFACE sedimentation_rain
302
[1691]303    INTERFACE calc_precipitation_amount
304       MODULE PROCEDURE calc_precipitation_amount
305       MODULE PROCEDURE calc_precipitation_amount_ij
306    END INTERFACE calc_precipitation_amount
307
[1000]308 CONTAINS
[1849]309!------------------------------------------------------------------------------!
310! Description:
311! ------------
312!> Initialization of bulk microphysics
313!------------------------------------------------------------------------------!
314    SUBROUTINE microphysics_init
[1000]315
[1849]316       USE arrays_3d,                                                          &
317           ONLY:  dzu
[1000]318
[1849]319       USE constants,                                                          &
320           ONLY:  pi
321
322       USE cloud_parameters,                                                   &
[2375]323           ONLY:  molecular_weight_of_solute, rho_l, rho_s, vanthoff
[1849]324
325       USE control_parameters,                                                 &
[2375]326           ONLY:  aerosol_nacl, aerosol_c3h4o4 , aerosol_nh4no3,               &
327                  microphysics_morrison, microphysics_seifert
[1849]328
329       USE indices,                                                            &
330           ONLY:  nzb, nzt
331
332       IMPLICIT NONE
333
334!
335!--    constant for the sedimentation of cloud water (2-moment cloud physics)
336       sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l )                &
337                          )**( 2.0_wp / 3.0_wp ) *                             &
338                      EXP( 5.0_wp * LOG( sigma_gc )**2 )
339
340!
341!--    Calculate timestep according to precipitation
342       IF ( microphysics_seifert )  THEN
343          dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /      &
344                             w_precipitation
345       ENDIF
346
347!
[2375]348!--    Set constants for certain aerosol type
349       IF ( microphysics_morrison )  THEN
350          IF ( aerosol_nacl ) THEN
351             molecular_weight_of_solute = 0.05844_wp 
352             rho_s                      = 2165.0_wp
353             vanthoff                   = 2.0_wp
354          ELSEIF ( aerosol_c3h4o4 ) THEN
355             molecular_weight_of_solute = 0.10406_wp 
356             rho_s                      = 1600.0_wp
357             vanthoff                   = 1.37_wp
358          ELSEIF ( aerosol_nh4no3 ) THEN
359             molecular_weight_of_solute = 0.08004_wp 
360             rho_s                      = 1720.0_wp
361             vanthoff                   = 2.31_wp
362          ENDIF
363       ENDIF
364 
365!
[1849]366!--    Pre-calculate frequently calculated fractions of pi and rho_l
367       pirho_l  = pi * rho_l / 6.0_wp
368       dpirho_l = 1.0_wp / pirho_l
369
370!
371!--    Allocate 1D microphysics arrays
[2292]372       ALLOCATE ( pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1),         &
[1849]373                  qc_1d(nzb:nzt+1) )
374
[2292]375       IF ( microphysics_morrison .OR. microphysics_seifert ) THEN
376          ALLOCATE ( nc_1d(nzb:nzt+1) )
377       ENDIF
378
[1849]379       IF ( microphysics_seifert )  THEN
380          ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
381       ENDIF
382
383    END SUBROUTINE microphysics_init
384
385
[1000]386!------------------------------------------------------------------------------!
[1682]387! Description:
388! ------------
[1849]389!> Control of microphysics for all grid points
[1000]390!------------------------------------------------------------------------------!
[1115]391    SUBROUTINE microphysics_control
[1022]392
[1361]393       USE arrays_3d,                                                          &
[1849]394           ONLY:  hyp, pt_init, prr, zu
[1361]395
396       USE cloud_parameters,                                                   &
[1849]397           ONLY:  cp, hyrho, pt_d_t, r_d, t_d_pt
[1361]398
399       USE control_parameters,                                                 &
[1849]400           ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
401                  intermediate_timestep_count, large_scale_forcing,            &
[2292]402                  lsf_surf, microphysics_kessler, microphysics_morrison,       &
403                  microphysics_seifert, pt_surface,                            &
404                  rho_surface,surface_pressure
[1361]405
406       USE indices,                                                            &
407           ONLY:  nzb, nzt
408
[1320]409       USE kinds
[1115]410
[1361]411       USE statistics,                                                         &
412           ONLY:  weight_pres
413
[1115]414       IMPLICIT NONE
415
[1682]416       INTEGER(iwp) ::  k                 !<
[1115]417
[1682]418       REAL(wp)     ::  t_surface         !<
[1361]419
420       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
421!
422!--       Calculate:
423!--       pt / t : ratio of potential and actual temperature (pt_d_t)
424!--       t / pt : ratio of actual and potential temperature (t_d_pt)
425!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
426          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
427          DO  k = nzb, nzt+1
428             hyp(k)    = surface_pressure * 100.0_wp * &
429                         ( ( t_surface - g / cp * zu(k) ) /                    &
430                         t_surface )**(1.0_wp / 0.286_wp)
431             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
432             t_d_pt(k) = 1.0_wp / pt_d_t(k)
[2155]433             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
[1115]434          ENDDO
[1822]435
[1361]436!
437!--       Compute reference density
438          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
439       ENDIF
[1115]440
[1361]441!
[2155]442!--    Compute length of time step
[1361]443       IF ( call_microphysics_at_all_substeps )  THEN
444          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
445       ELSE
446          dt_micro = dt_3d
447       ENDIF
448
449!
[1822]450!--    Reset precipitation rate
451       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
452
453!
[1361]454!--    Compute cloud physics
[1822]455       IF ( microphysics_kessler )  THEN
456
457          CALL autoconversion_kessler
[1831]458          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
[1822]459
460       ELSEIF ( microphysics_seifert )  THEN
461
[1361]462          CALL adjust_cloud
[2292]463          IF ( microphysics_morrison )  CALL activation
464          IF ( microphysics_morrison )  CALL condensation
[1361]465          CALL autoconversion
466          CALL accretion
467          CALL selfcollection_breakup
468          CALL evaporation_rain
469          CALL sedimentation_rain
[1831]470          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
[1361]471
[1691]472       ENDIF
473
[1822]474       CALL calc_precipitation_amount
475
[1115]476    END SUBROUTINE microphysics_control
477
[1682]478!------------------------------------------------------------------------------!
479! Description:
480! ------------
[2155]481!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
482!> evaporation of rain drops due to too small or too big weights
[1682]483!> of rain drops (Stevens and Seifert, 2008).
484!------------------------------------------------------------------------------!
[1115]485    SUBROUTINE adjust_cloud
486
[1361]487       USE arrays_3d,                                                          &
[2292]488           ONLY:  qc, nc, qr, nr
[1361]489
490       USE cloud_parameters,                                                   &
[1849]491           ONLY:  hyrho
[1361]492
[2292]493       USE control_parameters,                                                 &
494           ONLY:  microphysics_morrison
495
[1361]496       USE cpulog,                                                             &
497           ONLY:  cpu_log, log_point_s
498
499       USE indices,                                                            &
[2317]500           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]501
[1320]502       USE kinds
[1022]503
504       IMPLICIT NONE
505
[1682]506       INTEGER(iwp) ::  i                 !<
507       INTEGER(iwp) ::  j                 !<
508       INTEGER(iwp) ::  k                 !<
[1022]509
[1361]510       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
511
[2155]512       DO  i = nxlg, nxrg
513          DO  j = nysg, nyng
[2232]514             DO  k = nzb+1, nzt
[1361]515                IF ( qr(k,j,i) <= eps_sb )  THEN
516                   qr(k,j,i) = 0.0_wp
517                   nr(k,j,i) = 0.0_wp
518                ELSE
519                   IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
[2232]520                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin *               &
[2312]521                                       MERGE( 1.0_wp, 0.0_wp,                  &
[2232]522                                              BTEST( wall_flags_0(k,j,i), 0 ) )
[1361]523                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
[2232]524                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax *               &
[2312]525                                       MERGE( 1.0_wp, 0.0_wp,                  &
[2232]526                                              BTEST( wall_flags_0(k,j,i), 0 ) )
[1361]527                   ENDIF
528                ENDIF
[2292]529                IF ( microphysics_morrison ) THEN
530                   IF ( qc(k,j,i) <= eps_sb )  THEN
531                      qc(k,j,i) = 0.0_wp
532                      nc(k,j,i) = 0.0_wp
533                   ELSE
534                      IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
535                         nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin *            &
[2312]536                                          MERGE( 1.0_wp, 0.0_wp,               &
[2292]537                                              BTEST( wall_flags_0(k,j,i), 0 ) )
538                      ENDIF
539                   ENDIF
540                ENDIF
[1022]541             ENDDO
542          ENDDO
543       ENDDO
544
[1361]545       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
546
[1115]547    END SUBROUTINE adjust_cloud
[1022]548
[2292]549!------------------------------------------------------------------------------!
550! Description:
551! ------------
552!> Calculate number of activated condensation nucleii after simple activation
553!> scheme of Twomey, 1959.
554!------------------------------------------------------------------------------!
555    SUBROUTINE activation
[1106]556
[2292]557       USE arrays_3d,                                                          &
558           ONLY:  hyp, nc, nr, pt, q,  qc, qr
559
560       USE cloud_parameters,                                                   &
[2375]561           ONLY:  hyrho, l_d_cp, l_d_r, l_v, molecular_weight_of_solute,       &
562                  molecular_weight_of_water, rho_l, rho_s, r_v, t_d_pt,        &
563                  vanthoff
[2292]564
565       USE constants,                                                          &
566           ONLY:  pi
567
568       USE cpulog,                                                             &
569           ONLY:  cpu_log, log_point_s
570
[2608]571       USE diagnostic_quantities_mod,                                          &
572           ONLY: e_s, magnus, q_s, sat, supersaturation, t_l
573
[2292]574       USE indices,                                                            &
575           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
576
577       USE kinds
578
579       USE control_parameters,                                                 &
580          ONLY: simulated_time
581
582       IMPLICIT NONE
583
584       INTEGER(iwp) ::  i                 !<
585       INTEGER(iwp) ::  j                 !<
586       INTEGER(iwp) ::  k                 !<
587
588       REAL(wp)     ::  activ             !<
589       REAL(wp)     ::  afactor           !<
590       REAL(wp)     ::  beta_act          !<
591       REAL(wp)     ::  bfactor           !<
592       REAL(wp)     ::  k_act             !<
593       REAL(wp)     ::  n_act             !<
594       REAL(wp)     ::  n_ccn             !<
595       REAL(wp)     ::  s_0               !<
596       REAL(wp)     ::  sat_max           !<
597       REAL(wp)     ::  sigma             !<
[2375]598       REAL(wp)     ::  sigma_act         !<
[2292]599
600       CALL cpu_log( log_point_s(65), 'activation', 'start' )
601
602       DO  i = nxlg, nxrg
603          DO  j = nysg, nyng
604             DO  k = nzb+1, nzt
605
606!
[2608]607!--             Call calculation of supersaturation located
608!--             in diagnostic_quantities_mod
609                CALL supersaturation ( i, j, k )
[2292]610!
611!--             Prescribe parameters for activation
612!--             (see: Bott + Trautmann, 2002, Atm. Res., 64)
[2312]613                k_act  = 0.7_wp
[2522]614                activ  = 0.0_wp
[2292]615
[2522]616
617                IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
[2292]618!
[2312]619!--                Compute the number of activated Aerosols
[2292]620!--                (see: Twomey, 1959, Pure and applied Geophysics, 43)
621                   n_act     = na_init * sat**k_act
622!
623!--                Compute the number of cloud droplets
624!--                (see: Morrison + Grabowski, 2007, JAS, 64)
625!                  activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro
626
627!
628!--                Compute activation rate after Khairoutdinov and Kogan
[2312]629!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
630                   sat_max = 1.0_wp / 100.0_wp
631                   activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN       &
632                      ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /    &
[2292]633                       dt_micro
[2522]634                ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN
[2292]635!
[2312]636!--                Curvature effect (afactor) with surface tension
[2292]637!--                parameterization by Straka (2009)
[2608]638                   sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
639                   afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
[2292]640!
641!--                Solute effect (bfactor)
642                   bfactor = vanthoff * molecular_weight_of_water *            &
643                       rho_s / ( molecular_weight_of_solute * rho_l )
644
645!
[2312]646!--                Prescribe power index that describes the soluble fraction
[2375]647!--                of an aerosol particle (beta)
[2292]648!--                (see: Morrison + Grabowski, 2007, JAS, 64)
[2375]649                   beta_act  = 0.5_wp
650                   sigma_act = sigma_bulk**( 1.0_wp + beta_act )
[2292]651!
[2312]652!--                Calculate mean geometric supersaturation (s_0) with
[2292]653!--                parameterization by Khvorostyanov and Curry (2006)
[2375]654                   s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *    &
[2292]655                       ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
656
657!
[2312]658!--                Calculate number of activated CCN as a function of
[2292]659!--                supersaturation and taking Koehler theory into account
[2312]660!--                (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
661                   n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(              &
[2375]662                      LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
[2292]663                   activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
664                ENDIF
665
666                nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
667
668             ENDDO
669          ENDDO
670       ENDDO
671
672       CALL cpu_log( log_point_s(65), 'activation', 'stop' )
673
674    END SUBROUTINE activation
675
676
[1682]677!------------------------------------------------------------------------------!
678! Description:
679! ------------
[2312]680!> Calculate condensation rate for cloud water content (after Khairoutdinov and
[2292]681!> Kogan, 2000).
682!------------------------------------------------------------------------------!
683    SUBROUTINE condensation
684
685       USE arrays_3d,                                                          &
686           ONLY:  hyp, nr, pt, q,  qc, qr, nc
687
688       USE cloud_parameters,                                                   &
689           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
690
691       USE constants,                                                          &
692           ONLY:  pi
693
694       USE cpulog,                                                             &
695           ONLY:  cpu_log, log_point_s
696
[2608]697       USE diagnostic_quantities_mod,                                          &
698           ONLY: e_s, magnus, q_s, sat, supersaturation, t_l
699
[2292]700       USE indices,                                                            &
701           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
702
703       USE kinds
704
705       USE control_parameters,                                                 &
706          ONLY: simulated_time
707
708       IMPLICIT NONE
709
710       INTEGER(iwp) ::  i                 !<
711       INTEGER(iwp) ::  j                 !<
712       INTEGER(iwp) ::  k                 !<
713
714       REAL(wp)     ::  cond              !<
715       REAL(wp)     ::  cond_max          !<
716       REAL(wp)     ::  dc                !<
717       REAL(wp)     ::  evap              !<
718       REAL(wp)     ::  evap_nc           !<
719       REAL(wp)     ::  g_fac             !<
720       REAL(wp)     ::  nc_0              !<
721       REAL(wp)     ::  temp              !<
722       REAL(wp)     ::  xc                !<
723
724       CALL cpu_log( log_point_s(66), 'condensation', 'start' )
725
726       DO  i = nxlg, nxrg
727          DO  j = nysg, nyng
728             DO  k = nzb+1, nzt
729!
[2608]730!--             Call calculation of supersaturation located
731!--             in diagnostic_quantities_mod
732                CALL supersaturation ( i, j, k )
[2292]733!
734!--             Actual temperature:
735                temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
736
737                g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
738                                    l_v / ( thermal_conductivity_l * temp )    &
739                                    + r_v * temp / ( diff_coeff_l * e_s )      &
740                                    )
741!
742!--             Mean weight of cloud drops
743                IF ( nc(k,j,i) <= 0.0_wp) CYCLE
[2312]744                xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
[2292]745!
746!--             Weight averaged diameter of cloud drops:
747                dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
748!
749!--             Integral diameter of cloud drops
[2312]750                nc_0 = nc(k,j,i) * dc
[2292]751!
752!--             Condensation needs only to be calculated in supersaturated regions
753                IF ( sat > 0.0_wp )  THEN
754!
755!--                Condensation rate of cloud water content
756!--                after KK scheme.
757!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
758                   cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
759                   cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
760                   cond      = MIN( cond, cond_max / dt_micro )
[2312]761
[2292]762                   qc(k,j,i) = qc(k,j,i) + cond * dt_micro
763                ELSEIF ( sat < 0.0_wp ) THEN
764                   evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
765                   evap      = MAX( evap, -qc(k,j,i) / dt_micro )
766
767                   qc(k,j,i) = qc(k,j,i) + evap * dt_micro
768                ENDIF
769                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
770                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin
771                ENDIF
772             ENDDO
773          ENDDO
774       ENDDO
775
776       CALL cpu_log( log_point_s(66), 'condensation', 'stop' )
777
778    END SUBROUTINE condensation
779
780
781!------------------------------------------------------------------------------!
782! Description:
783! ------------
[1682]784!> Autoconversion rate (Seifert and Beheng, 2006).
785!------------------------------------------------------------------------------!
[1000]786    SUBROUTINE autoconversion
787
[1361]788       USE arrays_3d,                                                          &
[2292]789           ONLY:  diss, dzu, nc, nr, qc, qr
[1361]790
791       USE cloud_parameters,                                                   &
[1849]792           ONLY:  hyrho
[1361]793
794       USE control_parameters,                                                 &
[2292]795           ONLY:  microphysics_morrison, rho_surface
[1361]796
797       USE cpulog,                                                             &
798           ONLY:  cpu_log, log_point_s
799
800       USE grid_variables,                                                     &
801           ONLY:  dx, dy
802
803       USE indices,                                                            &
[2317]804           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]805
[1320]806       USE kinds
[1000]807
808       IMPLICIT NONE
809
[1682]810       INTEGER(iwp) ::  i                 !<
811       INTEGER(iwp) ::  j                 !<
812       INTEGER(iwp) ::  k                 !<
[1000]813
[2155]814       REAL(wp)     ::  alpha_cc          !<
[1682]815       REAL(wp)     ::  autocon           !<
816       REAL(wp)     ::  dissipation       !<
[2232]817       REAL(wp)     ::  flag              !< flag to mask topography grid points
[1682]818       REAL(wp)     ::  k_au              !<
819       REAL(wp)     ::  l_mix             !<
[2292]820       REAL(wp)     ::  nc_auto           !<
[1682]821       REAL(wp)     ::  nu_c              !<
822       REAL(wp)     ::  phi_au            !<
823       REAL(wp)     ::  r_cc              !<
824       REAL(wp)     ::  rc                !<
825       REAL(wp)     ::  re_lambda         !<
826       REAL(wp)     ::  sigma_cc          !<
827       REAL(wp)     ::  tau_cloud         !<
828       REAL(wp)     ::  xc                !<
[1361]829
830       CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
831
[2155]832       DO  i = nxlg, nxrg
833          DO  j = nysg, nyng
[2232]834             DO  k = nzb+1, nzt
835!
836!--             Predetermine flag to mask topography
837                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1000]838
[2292]839                nc_auto = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
840
[2522]841                IF ( qc(k,j,i) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
[1361]842
843                   k_au = k_cc / ( 20.0_wp * x0 )
844!
845!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
846!--                (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
[2522]847                   tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) +         &
848                                    qc(k,j,i) ), 0.0_wp )
[1361]849!
[2155]850!--                Universal function for autoconversion process
[1361]851!--                (Seifert and Beheng, 2006):
852                   phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
853                            ( 1.0_wp - tau_cloud**0.68_wp )**3
854!
855!--                Shape parameter of gamma distribution (Geoffroy et al., 2010):
856!--                (Use constant nu_c = 1.0_wp instead?)
857                   nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
858!
859!--                Mean weight of cloud droplets:
[2522]860                   xc = MAX( hyrho(k) * qc(k,j,i) / nc_auto, xcmin) 
[1361]861!
[2155]862!--                Parameterized turbulence effects on autoconversion (Seifert,
[1361]863!--                Nuijens and Stevens, 2010)
[1831]864                   IF ( collision_turbulence )  THEN
[1361]865!
866!--                   Weight averaged radius of cloud droplets:
867                      rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
868
869                      alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
870                      r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
871                      sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
872!
[2155]873!--                   Mixing length (neglecting distance to ground and
[1361]874!--                   stratification)
875                      l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
876!
[2155]877!--                   Limit dissipation rate according to Seifert, Nuijens and
[1361]878!--                   Stevens (2010)
879                      dissipation = MIN( 0.06_wp, diss(k,j,i) )
880!
881!--                   Compute Taylor-microscale Reynolds number:
882                      re_lambda = 6.0_wp / 11.0_wp *                           &
883                                  ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *   &
884                                  SQRT( 15.0_wp / kin_vis_air ) *              &
885                                  dissipation**( 1.0_wp / 6.0_wp )
886!
[2155]887!--                   The factor of 1.0E4 is needed to convert the dissipation
[1361]888!--                   rate from m2 s-3 to cm2 s-3.
889                      k_au = k_au * ( 1.0_wp +                                 &
890                                      dissipation * 1.0E4_wp *                 &
891                                      ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
[2155]892                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &
[1361]893                                                                      r_cc ) / &
894                                                        sigma_cc )**2          &
895                                                      ) + beta_cc              &
896                                      )                                        &
897                                    )
898                   ENDIF
899!
900!--                Autoconversion rate (Seifert and Beheng, 2006):
901                   autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /    &
902                             ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *     &
903                             ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * &
904                             rho_surface
905                   autocon = MIN( autocon, qc(k,j,i) / dt_micro )
906
[2232]907                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
908                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
909                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro  &
910                                                              * flag
[2292]911                   IF ( microphysics_morrison ) THEN
[2312]912                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *         &
[2292]913                                    autocon / x0 * hyrho(k) * dt_micro * flag )
914                   ENDIF
[1361]915
916                ENDIF
917
[1000]918             ENDDO
919          ENDDO
920       ENDDO
921
[1361]922       CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
923
[1000]924    END SUBROUTINE autoconversion
925
[1106]926
[1682]927!------------------------------------------------------------------------------!
928! Description:
929! ------------
[1822]930!> Autoconversion process (Kessler, 1969).
931!------------------------------------------------------------------------------!
932    SUBROUTINE autoconversion_kessler
933
934       USE arrays_3d,                                                          &
[1849]935           ONLY:  dzw, pt, prr, q, qc
[1822]936
937       USE cloud_parameters,                                                   &
[1849]938           ONLY:  l_d_cp, pt_d_t
[1822]939
940       USE indices,                                                            &
[2317]941           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1822]942
943       USE kinds
944
[2317]945       USE surface_mod,                                                        &
[2698]946           ONLY:  get_topography_top_index_ji
[1822]947
[2317]948
[1822]949       IMPLICIT NONE
950
[2232]951       INTEGER(iwp) ::  i      !<
952       INTEGER(iwp) ::  j      !<
953       INTEGER(iwp) ::  k      !<
954       INTEGER(iwp) ::  k_wall !< topgraphy top index
[1822]955
956       REAL(wp)    ::  dqdt_precip !<
[2232]957       REAL(wp)    ::  flag        !< flag to mask topography grid points
[1822]958
[2155]959       DO  i = nxlg, nxrg
960          DO  j = nysg, nyng
[2232]961!
962!--          Determine vertical index of topography top
[2698]963             k_wall = get_topography_top_index_ji( j, i, 's' )
[2232]964             DO  k = nzb+1, nzt
965!
966!--             Predetermine flag to mask topography
967                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1822]968
969                IF ( qc(k,j,i) > ql_crit )  THEN
970                   dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
971                ELSE
972                   dqdt_precip = 0.0_wp
973                ENDIF
974
[2232]975                qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
976                q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
[1845]977                pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp *      &
[2232]978                                        pt_d_t(k)              * flag
[1822]979
980!
[1845]981!--             Compute the rain rate (stored on surface grid point)
[2232]982                prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
[1822]983
984             ENDDO
985          ENDDO
986       ENDDO
987
988   END SUBROUTINE autoconversion_kessler
989
990
991!------------------------------------------------------------------------------!
992! Description:
993! ------------
[1682]994!> Accretion rate (Seifert and Beheng, 2006).
995!------------------------------------------------------------------------------!
[1005]996    SUBROUTINE accretion
[1000]997
[1361]998       USE arrays_3d,                                                          &
[2292]999           ONLY:  diss, qc, qr, nc
[1361]1000
1001       USE cloud_parameters,                                                   &
[1849]1002           ONLY:  hyrho
[1361]1003
1004       USE control_parameters,                                                 &
[2292]1005           ONLY:  microphysics_morrison, rho_surface
[1361]1006
1007       USE cpulog,                                                             &
1008           ONLY:  cpu_log, log_point_s
1009
1010       USE indices,                                                            &
[2317]1011           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]1012
[1320]1013       USE kinds
[1005]1014
[1000]1015       IMPLICIT NONE
1016
[1682]1017       INTEGER(iwp) ::  i                 !<
1018       INTEGER(iwp) ::  j                 !<
1019       INTEGER(iwp) ::  k                 !<
[1000]1020
[1682]1021       REAL(wp)     ::  accr              !<
[2232]1022       REAL(wp)     ::  flag              !< flag to mask topography grid points
[1682]1023       REAL(wp)     ::  k_cr              !<
[2292]1024       REAL(wp)     ::  nc_accr           !<
[1682]1025       REAL(wp)     ::  phi_ac            !<
1026       REAL(wp)     ::  tau_cloud         !<
[2292]1027       REAL(wp)     ::  xc                !<
[1361]1028
[2292]1029
[1361]1030       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
1031
[2155]1032       DO  i = nxlg, nxrg
1033          DO  j = nysg, nyng
[2232]1034             DO  k = nzb+1, nzt
1035!
1036!--             Predetermine flag to mask topography
1037                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1000]1038
[2292]1039                nc_accr = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
1040
1041                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )    &
1042                                             .AND.  ( nc_accr > eps_mr ) ) THEN
[1361]1043!
1044!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
[2155]1045                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
[1361]1046!
[2155]1047!--                Universal function for accretion process (Seifert and
[1361]1048!--                Beheng, 2001):
1049                   phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
[2292]1050
[1361]1051!
[2292]1052!--                Mean weight of cloud drops
[2312]1053                   xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
[2292]1054!
[2155]1055!--                Parameterized turbulence effects on autoconversion (Seifert,
1056!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
[1361]1057!--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1831]1058                   IF ( collision_turbulence )  THEN
[1361]1059                      k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                      &
1060                                       MIN( 600.0_wp,                          &
1061                                            diss(k,j,i) * 1.0E4_wp )**0.25_wp  &
1062                                     )
1063                   ELSE
[2155]1064                      k_cr = k_cr0
[1361]1065                   ENDIF
1066!
1067!--                Accretion rate (Seifert and Beheng, 2006):
1068                   accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *              &
1069                          SQRT( rho_surface * hyrho(k) )
1070                   accr = MIN( accr, qc(k,j,i) / dt_micro )
1071
[2232]1072                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
1073                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
[2292]1074                   IF ( microphysics_morrison )  THEN
1075                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i),                  &
1076                                         accr / xc * hyrho(k) * dt_micro * flag)
1077                   ENDIF
[1361]1078
1079                ENDIF
1080
[1005]1081             ENDDO
1082          ENDDO
[1000]1083       ENDDO
1084
[1361]1085       CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
1086
[1005]1087    END SUBROUTINE accretion
[1000]1088
[1106]1089
[1682]1090!------------------------------------------------------------------------------!
1091! Description:
1092! ------------
1093!> Collisional breakup rate (Seifert, 2008).
1094!------------------------------------------------------------------------------!
[1005]1095    SUBROUTINE selfcollection_breakup
[1000]1096
[1361]1097       USE arrays_3d,                                                          &
1098           ONLY:  nr, qr
1099
1100       USE cloud_parameters,                                                   &
[1849]1101           ONLY:  hyrho
[1361]1102
1103       USE control_parameters,                                                 &
[1849]1104           ONLY:  rho_surface
[1361]1105
1106       USE cpulog,                                                             &
1107           ONLY:  cpu_log, log_point_s
1108
1109       USE indices,                                                            &
[2317]1110           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]1111
[1320]1112       USE kinds
[2155]1113
[1000]1114       IMPLICIT NONE
1115
[1682]1116       INTEGER(iwp) ::  i                 !<
1117       INTEGER(iwp) ::  j                 !<
1118       INTEGER(iwp) ::  k                 !<
[1000]1119
[1682]1120       REAL(wp)     ::  breakup           !<
1121       REAL(wp)     ::  dr                !<
[2232]1122       REAL(wp)     ::  flag              !< flag to mask topography grid points
[1682]1123       REAL(wp)     ::  phi_br            !<
1124       REAL(wp)     ::  selfcoll          !<
[1361]1125
1126       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
1127
[2155]1128       DO  i = nxlg, nxrg
1129          DO  j = nysg, nyng
[2232]1130             DO  k = nzb+1, nzt
1131!
1132!--             Predetermine flag to mask topography
1133                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1134
[1361]1135                IF ( qr(k,j,i) > eps_sb )  THEN
1136!
1137!--                Selfcollection rate (Seifert and Beheng, 2001):
1138                   selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) *                   &
1139                              SQRT( hyrho(k) * rho_surface )
1140!
1141!--                Weight averaged diameter of rain drops:
1142                   dr = ( hyrho(k) * qr(k,j,i) /                               &
1143                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
1144!
1145!--                Collisional breakup rate (Seifert, 2008):
1146                   IF ( dr >= 0.3E-3_wp )  THEN
1147                      phi_br  = k_br * ( dr - 1.1E-3_wp )
1148                      breakup = selfcoll * ( phi_br + 1.0_wp )
1149                   ELSE
1150                      breakup = 0.0_wp
1151                   ENDIF
[1000]1152
[1361]1153                   selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
[2232]1154                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
[1361]1155
[2155]1156                ENDIF
[1000]1157             ENDDO
1158          ENDDO
1159       ENDDO
1160
[1361]1161       CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
1162
[1005]1163    END SUBROUTINE selfcollection_breakup
[1000]1164
[1106]1165
[1682]1166!------------------------------------------------------------------------------!
1167! Description:
1168! ------------
[2155]1169!> Evaporation of precipitable water. Condensation is neglected for
[1682]1170!> precipitable water.
1171!------------------------------------------------------------------------------!
[1012]1172    SUBROUTINE evaporation_rain
[1000]1173
[1361]1174       USE arrays_3d,                                                          &
1175           ONLY:  hyp, nr, pt, q,  qc, qr
1176
1177       USE cloud_parameters,                                                   &
[1849]1178           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
[1361]1179
1180       USE constants,                                                          &
1181           ONLY:  pi
1182
1183       USE cpulog,                                                             &
1184           ONLY:  cpu_log, log_point_s
1185
[2608]1186       USE diagnostic_quantities_mod,                                          &
1187           ONLY: e_s, magnus, q_s, sat, supersaturation, t_l
1188
[1361]1189       USE indices,                                                            &
[2317]1190           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]1191
[1320]1192       USE kinds
[1012]1193
1194       IMPLICIT NONE
1195
[1682]1196       INTEGER(iwp) ::  i                 !<
1197       INTEGER(iwp) ::  j                 !<
1198       INTEGER(iwp) ::  k                 !<
[1361]1199
[1682]1200       REAL(wp)     ::  dr                !<
1201       REAL(wp)     ::  evap              !<
1202       REAL(wp)     ::  evap_nr           !<
1203       REAL(wp)     ::  f_vent            !<
[2232]1204       REAL(wp)     ::  flag              !< flag to mask topography grid points
[1682]1205       REAL(wp)     ::  g_evap            !<
1206       REAL(wp)     ::  lambda_r          !<
1207       REAL(wp)     ::  mu_r              !<
1208       REAL(wp)     ::  mu_r_2            !<
1209       REAL(wp)     ::  mu_r_5d2          !<
1210       REAL(wp)     ::  nr_0              !<
1211       REAL(wp)     ::  temp              !<
1212       REAL(wp)     ::  xr                !<
[1361]1213
1214       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
1215
[2155]1216       DO  i = nxlg, nxrg
1217          DO  j = nysg, nyng
[2232]1218             DO  k = nzb+1, nzt
1219!
1220!--             Predetermine flag to mask topography
1221                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1222
[1361]1223                IF ( qr(k,j,i) > eps_sb )  THEN
[2608]1224
[1361]1225!
[2608]1226!--                Call calculation of supersaturation located
1227!--                in diagnostic_quantities_mod
1228                   CALL supersaturation ( i, j, k )
[1361]1229!
1230!--                Evaporation needs only to be calculated in subsaturated regions
1231                   IF ( sat < 0.0_wp )  THEN
1232!
1233!--                   Actual temperature:
1234                      temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
[2155]1235
[1361]1236                      g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *     &
1237                                          l_v / ( thermal_conductivity_l * temp ) &
1238                                          + r_v * temp / ( diff_coeff_l * e_s )   &
1239                                        )
1240!
1241!--                   Mean weight of rain drops
1242                      xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
1243!
1244!--                   Weight averaged diameter of rain drops:
1245                      dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
1246!
[2155]1247!--                   Compute ventilation factor and intercept parameter
[1361]1248!--                   (Seifert and Beheng, 2006; Seifert, 2008):
1249                      IF ( ventilation_effect )  THEN
1250!
[2155]1251!--                      Shape parameter of gamma distribution (Milbrandt and Yau,
[1361]1252!--                      2005; Stevens and Seifert, 2008):
1253                         mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *          &
1254                                                          ( dr - 1.4E-3_wp ) ) )
1255!
1256!--                      Slope parameter of gamma distribution (Seifert, 2008):
1257                         lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *  &
1258                                      ( mu_r + 1.0_wp )                        &
1259                                    )**( 1.0_wp / 3.0_wp ) / dr
[1012]1260
[1361]1261                         mu_r_2   = mu_r + 2.0_wp
[2155]1262                         mu_r_5d2 = mu_r + 2.5_wp
[1361]1263
1264                         f_vent = a_vent * gamm( mu_r_2 ) *                     &
1265                                  lambda_r**( -mu_r_2 ) + b_vent *              &
1266                                  schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
1267                                  gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) *  &
1268                                  ( 1.0_wp -                                    &
1269                                    0.5_wp * ( b_term / a_term ) *              &
1270                                    ( lambda_r / ( c_term + lambda_r )          &
1271                                    )**mu_r_5d2 -                               &
1272                                    0.125_wp * ( b_term / a_term )**2 *         &
1273                                    ( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
1274                                    )**mu_r_5d2 -                               &
1275                                    0.0625_wp * ( b_term / a_term )**3 *        &
1276                                    ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
1277                                    )**mu_r_5d2 -                               &
[2155]1278                                    0.0390625_wp * ( b_term / a_term )**4 *     &
[1361]1279                                    ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
1280                                    )**mu_r_5d2                                 &
1281                                  )
1282
1283                         nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
[2155]1284                                  gamm( mu_r + 1.0_wp )
[1361]1285                      ELSE
1286                         f_vent = 1.0_wp
1287                         nr_0   = nr(k,j,i) * dr
1288                      ENDIF
1289   !
[2155]1290   !--                Evaporation rate of rain water content (Seifert and
[1361]1291   !--                Beheng, 2006):
1292                      evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
1293                                hyrho(k)
1294                      evap    = MAX( evap, -qr(k,j,i) / dt_micro )
1295                      evap_nr = MAX( c_evap * evap / xr * hyrho(k),            &
1296                                     -nr(k,j,i) / dt_micro )
1297
[2232]1298                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro * flag
1299                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
[1361]1300
1301                   ENDIF
[2155]1302                ENDIF
[1361]1303
[1012]1304             ENDDO
1305          ENDDO
1306       ENDDO
1307
[1361]1308       CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
1309
[1012]1310    END SUBROUTINE evaporation_rain
1311
[1106]1312
[1682]1313!------------------------------------------------------------------------------!
1314! Description:
1315! ------------
1316!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
1317!------------------------------------------------------------------------------!
[1012]1318    SUBROUTINE sedimentation_cloud
1319
[1361]1320       USE arrays_3d,                                                          &
[2292]1321           ONLY:  ddzu, dzu, nc, pt, prr, q, qc
[1361]1322
1323       USE cloud_parameters,                                                   &
[1849]1324           ONLY:  hyrho, l_d_cp, pt_d_t
[1361]1325
1326       USE control_parameters,                                                 &
[2312]1327           ONLY:  call_microphysics_at_all_substeps,                           &
[2292]1328                  intermediate_timestep_count, microphysics_morrison
[1361]1329
1330       USE cpulog,                                                             &
1331           ONLY:  cpu_log, log_point_s
1332
1333       USE indices,                                                            &
[2317]1334           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]1335
[1320]1336       USE kinds
[1691]1337
1338       USE statistics,                                                         &
1339           ONLY:  weight_substep
1340
[2155]1341
[1012]1342       IMPLICIT NONE
1343
[1849]1344       INTEGER(iwp) ::  i             !<
1345       INTEGER(iwp) ::  j             !<
1346       INTEGER(iwp) ::  k             !<
[1361]1347
[2292]1348       REAL(wp) ::  flag              !< flag to mask topography grid points
1349       REAL(wp) ::  nc_sedi           !<
1350
[2232]1351       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
[2292]1352       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc !<
[1361]1353
[2292]1354
[1361]1355       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
1356
1357       sed_qc(nzt+1) = 0.0_wp
[2375]1358       sed_nc(nzt+1) = 0.0_wp
[1361]1359
[2155]1360       DO  i = nxlg, nxrg
1361          DO  j = nysg, nyng
[2232]1362             DO  k = nzt, nzb+1, -1
1363!
1364!--             Predetermine flag to mask topography
1365                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2292]1366                nc_sedi = MERGE ( nc(k,j,i), nc_const, microphysics_morrison )
[1012]1367
[2292]1368!
1369!--             Sedimentation fluxes for number concentration are only calculated
1370!--             for cloud_scheme = 'morrison'
1371                IF ( microphysics_morrison ) THEN
1372                   IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
[2312]1373                      sed_nc(k) = sed_qc_const *                               &
[2292]1374                              ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
1375                              ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
1376                   ELSE
1377                      sed_nc(k) = 0.0_wp
1378                   ENDIF
1379
1380                   sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *           &
1381                                    nc(k,j,i) / dt_micro + sed_nc(k+1)         &
1382                                  ) * flag
1383
1384                   nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *       &
1385                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
1386                ENDIF
1387
[2375]1388                IF ( qc(k,j,i) > eps_sb .AND.  nc_sedi > eps_mr )  THEN
[2292]1389                   sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *  &
[2232]1390                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
1391                                                                           flag
[1361]1392                ELSE
1393                   sed_qc(k) = 0.0_wp
1394                ENDIF
1395
1396                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
1397                                            dt_micro + sed_qc(k+1)             &
[2232]1398                               ) * flag
[1361]1399
1400                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
[2232]1401                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]1402                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
[2232]1403                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]1404                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
1405                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
[2232]1406                                        pt_d_t(k) * dt_micro            * flag
[1361]1407
[1691]1408!
1409!--             Compute the precipitation rate due to cloud (fog) droplets
[1822]1410                IF ( call_microphysics_at_all_substeps )  THEN
1411                   prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
[2232]1412                                * weight_substep(intermediate_timestep_count)  &
1413                                * flag
[1822]1414                ELSE
[2232]1415                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
[1691]1416                ENDIF
1417
[1012]1418             ENDDO
1419          ENDDO
1420       ENDDO
1421
[1361]1422       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
1423
[1012]1424    END SUBROUTINE sedimentation_cloud
1425
[1106]1426
[1682]1427!------------------------------------------------------------------------------!
1428! Description:
1429! ------------
1430!> Computation of sedimentation flux. Implementation according to Stevens
1431!> and Seifert (2008). Code is based on UCLA-LES.
1432!------------------------------------------------------------------------------!
[1012]1433    SUBROUTINE sedimentation_rain
1434
[1361]1435       USE arrays_3d,                                                          &
[1849]1436           ONLY:  ddzu, dzu, nr, pt, prr, q, qr
[1361]1437
1438       USE cloud_parameters,                                                   &
[1849]1439           ONLY:  hyrho, l_d_cp, pt_d_t
[1361]1440
1441       USE control_parameters,                                                 &
[1849]1442           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
[1361]1443       USE cpulog,                                                             &
1444           ONLY:  cpu_log, log_point_s
1445
1446       USE indices,                                                            &
[2317]1447           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]1448
[1320]1449       USE kinds
[1012]1450
[1361]1451       USE statistics,                                                         &
1452           ONLY:  weight_substep
[2155]1453
[2232]1454       USE surface_mod,                                                        &
1455           ONLY :  bc_h
[2312]1456
[1012]1457       IMPLICIT NONE
1458
[2232]1459       INTEGER(iwp) ::  i             !< running index x direction
1460       INTEGER(iwp) ::  j             !< running index y direction
1461       INTEGER(iwp) ::  k             !< running index z direction
[2312]1462       INTEGER(iwp) ::  k_run         !<
[2232]1463       INTEGER(iwp) ::  l             !< running index of surface type
1464       INTEGER(iwp) ::  m             !< running index surface elements
1465       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
1466       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
[1361]1467
[1682]1468       REAL(wp)     ::  c_run                      !<
1469       REAL(wp)     ::  d_max                      !<
1470       REAL(wp)     ::  d_mean                     !<
1471       REAL(wp)     ::  d_min                      !<
1472       REAL(wp)     ::  dr                         !<
1473       REAL(wp)     ::  flux                       !<
[2232]1474       REAL(wp)     ::  flag                       !< flag to mask topography grid points
[1682]1475       REAL(wp)     ::  lambda_r                   !<
1476       REAL(wp)     ::  mu_r                       !<
1477       REAL(wp)     ::  z_run                      !<
[1361]1478
[1682]1479       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
1480       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
1481       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
1482       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
1483       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
1484       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
1485       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
1486       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
[1361]1487
1488       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
[1682]1489
[1361]1490!
[2155]1491!--    Compute velocities
1492       DO  i = nxlg, nxrg
1493          DO  j = nysg, nyng
[2232]1494             DO  k = nzb+1, nzt
1495!
1496!--             Predetermine flag to mask topography
1497                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1498
[1361]1499                IF ( qr(k,j,i) > eps_sb )  THEN
1500!
1501!--                Weight averaged diameter of rain drops:
1502                   dr = ( hyrho(k) * qr(k,j,i) /                               &
1503                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
1504!
1505!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1506!--                Stevens and Seifert, 2008):
1507                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
1508                                                     ( dr - 1.4E-3_wp ) ) )
1509!
1510!--                Slope parameter of gamma distribution (Seifert, 2008):
1511                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
1512                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1012]1513
[1361]1514                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
1515                                               a_term - b_term * ( 1.0_wp +    &
1516                                                  c_term /                     &
1517                                                  lambda_r )**( -1.0_wp *      &
1518                                                  ( mu_r + 1.0_wp ) )          &
1519                                              )                                &
[2232]1520                                ) * flag
[1361]1521
1522                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
1523                                               a_term - b_term * ( 1.0_wp +    &
1524                                                  c_term /                     &
1525                                                  lambda_r )**( -1.0_wp *      &
1526                                                  ( mu_r + 4.0_wp ) )          &
1527                                             )                                 &
[2232]1528                                ) * flag
[1361]1529                ELSE
1530                   w_nr(k) = 0.0_wp
1531                   w_qr(k) = 0.0_wp
1532                ENDIF
[1012]1533             ENDDO
[1361]1534!
[2312]1535!--          Adjust boundary values using surface data type.
[2232]1536!--          Upward-facing
[2312]1537             surf_s = bc_h(0)%start_index(j,i)
[2232]1538             surf_e = bc_h(0)%end_index(j,i)
1539             DO  m = surf_s, surf_e
1540                k         = bc_h(0)%k(m)
1541                w_nr(k-1) = w_nr(k)
1542                w_qr(k-1) = w_qr(k)
1543             ENDDO
1544!
1545!--          Downward-facing
[2312]1546             surf_s = bc_h(1)%start_index(j,i)
[2232]1547             surf_e = bc_h(1)%end_index(j,i)
1548             DO  m = surf_s, surf_e
1549                k         = bc_h(1)%k(m)
1550                w_nr(k+1) = w_nr(k)
1551                w_qr(k+1) = w_qr(k)
1552             ENDDO
1553!
1554!--          Model top boundary value
[1361]1555             w_nr(nzt+1) = 0.0_wp
1556             w_qr(nzt+1) = 0.0_wp
1557!
1558!--          Compute Courant number
[2232]1559             DO  k = nzb+1, nzt
1560!
1561!--             Predetermine flag to mask topography
1562                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1563
[1361]1564                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
1565                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
[2232]1566                          dt_micro * ddzu(k) * flag
[1361]1567                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
1568                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
[2232]1569                          dt_micro * ddzu(k) * flag
[2155]1570             ENDDO
[1361]1571!
1572!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
1573             IF ( limiter_sedimentation )  THEN
1574
[2232]1575                DO k = nzb+1, nzt
1576!
1577!--                Predetermine flag to mask topography
1578                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1579
[1646]1580                   d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
[1361]1581                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
1582                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
1583
1584                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
1585                                                              2.0_wp * d_max,  &
[2232]1586                                                              ABS( d_mean ) )  &
1587                                                      * flag
[1361]1588
[1646]1589                   d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
[1361]1590                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
1591                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
1592
1593                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
1594                                                              2.0_wp * d_max,  &
1595                                                              ABS( d_mean ) )
1596                ENDDO
1597
1598             ELSE
1599
1600                nr_slope = 0.0_wp
1601                qr_slope = 0.0_wp
1602
1603             ENDIF
1604
1605             sed_nr(nzt+1) = 0.0_wp
1606             sed_qr(nzt+1) = 0.0_wp
1607!
1608!--          Compute sedimentation flux
[2232]1609             DO  k = nzt, nzb+1, -1
[1361]1610!
[2232]1611!--             Predetermine flag to mask topography
1612                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1613!
[2312]1614!--             Sum up all rain drop number densities which contribute to the flux
[1361]1615!--             through k-1/2
1616                flux  = 0.0_wp
1617                z_run = 0.0_wp ! height above z(k)
1618                k_run = k
1619                c_run = MIN( 1.0_wp, c_nr(k) )
1620                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
1621                   flux  = flux + hyrho(k_run) *                               &
1622                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
[2232]1623                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)  &
1624                                              * flag
1625                   z_run = z_run + dzu(k_run) * flag
1626                   k_run = k_run + 1          * flag
1627                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )    &
1628                                              * flag
[1361]1629                ENDDO
1630!
[2155]1631!--             It is not allowed to sediment more rain drop number density than
[1361]1632!--             available
1633                flux = MIN( flux,                                              &
1634                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
1635                            dt_micro                                           &
1636                          )
1637
[2232]1638                sed_nr(k) = flux / dt_micro * flag
[1361]1639                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
[2232]1640                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]1641!
[2155]1642!--             Sum up all rain water content which contributes to the flux
[1361]1643!--             through k-1/2
1644                flux  = 0.0_wp
1645                z_run = 0.0_wp ! height above z(k)
1646                k_run = k
1647                c_run = MIN( 1.0_wp, c_qr(k) )
1648
1649                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
1650
1651                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
1652                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
[2232]1653                                  0.5_wp ) * c_run * dzu(k_run) * flag
1654                   z_run = z_run + dzu(k_run)                   * flag
1655                   k_run = k_run + 1                            * flag
1656                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )    &
1657                                                                * flag
[1361]1658
1659                ENDDO
1660!
[2155]1661!--             It is not allowed to sediment more rain water content than
[1361]1662!--             available
1663                flux = MIN( flux,                                              &
1664                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
1665                            dt_micro                                           &
1666                          )
1667
[2232]1668                sed_qr(k) = flux / dt_micro * flag
[1361]1669
1670                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
[2232]1671                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]1672                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
[2232]1673                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]1674                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
1675                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
[2232]1676                                        pt_d_t(k) * dt_micro            * flag
[1361]1677!
1678!--             Compute the rain rate
1679                IF ( call_microphysics_at_all_substeps )  THEN
[1691]1680                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
[2232]1681                                * weight_substep(intermediate_timestep_count) &
1682                                * flag
[1361]1683                ELSE
[2232]1684                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
[1361]1685                ENDIF
1686
1687             ENDDO
[1012]1688          ENDDO
1689       ENDDO
1690
[1691]1691       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
1692
1693    END SUBROUTINE sedimentation_rain
1694
1695
1696!------------------------------------------------------------------------------!
1697! Description:
1698! ------------
1699!> Computation of the precipitation amount due to gravitational settling of
1700!> rain and cloud (fog) droplets
1701!------------------------------------------------------------------------------!
1702    SUBROUTINE calc_precipitation_amount
1703
[1849]1704       USE arrays_3d,                                                          &
1705           ONLY:  precipitation_amount, prr
1706
[1691]1707       USE cloud_parameters,                                                   &
[1849]1708           ONLY:  hyrho
[1691]1709
1710       USE control_parameters,                                                 &
1711           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
1712                  intermediate_timestep_count, intermediate_timestep_count_max,&
1713                  precipitation_amount_interval, time_do2d_xy
1714
1715       USE indices,                                                            &
[2232]1716           ONLY:  nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0
[1691]1717
1718       USE kinds
1719
[2232]1720       USE surface_mod,                                                        &
1721           ONLY :  bc_h
1722
[1691]1723       IMPLICIT NONE
1724
[2232]1725       INTEGER(iwp) ::  i             !< running index x direction
1726       INTEGER(iwp) ::  j             !< running index y direction
1727       INTEGER(iwp) ::  k             !< running index y direction
1728       INTEGER(iwp) ::  m             !< running index surface elements
1729       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
1730       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
[1691]1731
1732       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
1733            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
1734            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
1735       THEN
[2232]1736!
[2312]1737!--       Run over all upward-facing surface elements, i.e. non-natural,
[2232]1738!--       natural and urban
1739          DO  m = 1, bc_h(0)%ns
[2312]1740             i = bc_h(0)%i(m)
[2232]1741             j = bc_h(0)%j(m)
1742             k = bc_h(0)%k(m)
1743             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
1744                                               prr(k,j,i) * hyrho(k) * dt_3d
1745          ENDDO
[1691]1746
[1361]1747       ENDIF
1748
[1691]1749    END SUBROUTINE calc_precipitation_amount
[1361]1750
[1012]1751
[1000]1752!------------------------------------------------------------------------------!
[1682]1753! Description:
1754! ------------
[1849]1755!> Control of microphysics for grid points i,j
[1000]1756!------------------------------------------------------------------------------!
[1022]1757
[1115]1758    SUBROUTINE microphysics_control_ij( i, j )
1759
[1320]1760       USE arrays_3d,                                                          &
[2292]1761           ONLY:  hyp, nc, nr, pt, pt_init, prr, q, qc, qr, zu
[1115]1762
[1320]1763       USE cloud_parameters,                                                   &
[1849]1764           ONLY:  cp, hyrho, pt_d_t, r_d, t_d_pt
[1320]1765
1766       USE control_parameters,                                                 &
[1849]1767           ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
1768                  intermediate_timestep_count, large_scale_forcing,            &
[2292]1769                  lsf_surf, microphysics_morrison, microphysics_seifert,       &
1770                  microphysics_kessler, pt_surface, rho_surface,               &
1771                  surface_pressure
[1320]1772
1773       USE indices,                                                            &
1774           ONLY:  nzb, nzt
1775
1776       USE kinds
1777
1778       USE statistics,                                                         &
1779           ONLY:  weight_pres
1780
[1022]1781       IMPLICIT NONE
1782
[1682]1783       INTEGER(iwp) ::  i                 !<
1784       INTEGER(iwp) ::  j                 !<
1785       INTEGER(iwp) ::  k                 !<
[1115]1786
[1682]1787       REAL(wp)     ::  t_surface         !<
[1320]1788
[1361]1789       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
[1241]1790!
1791!--       Calculate:
1792!--       pt / t : ratio of potential and actual temperature (pt_d_t)
1793!--       t / pt : ratio of actual and potential temperature (t_d_pt)
1794!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
[1353]1795          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
[1241]1796          DO  k = nzb, nzt+1
[1353]1797             hyp(k)    = surface_pressure * 100.0_wp * &
[1361]1798                         ( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
[1353]1799             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
1800             t_d_pt(k) = 1.0_wp / pt_d_t(k)
[2155]1801             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
[1241]1802          ENDDO
1803!
1804!--       Compute reference density
[1353]1805          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
[1241]1806       ENDIF
1807
[1361]1808!
[2155]1809!--    Compute length of time step
[1361]1810       IF ( call_microphysics_at_all_substeps )  THEN
1811          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
1812       ELSE
1813          dt_micro = dt_3d
1814       ENDIF
[1241]1815
[1115]1816!
[1361]1817!--    Use 1d arrays
[1115]1818       q_1d(:)  = q(:,j,i)
1819       pt_1d(:) = pt(:,j,i)
1820       qc_1d(:) = qc(:,j,i)
[1822]1821       IF ( microphysics_seifert )  THEN
[1115]1822          qr_1d(:) = qr(:,j,i)
1823          nr_1d(:) = nr(:,j,i)
1824       ENDIF
[2292]1825       IF ( microphysics_morrison )  THEN
1826          nc_1d(:) = nc(:,j,i)
1827       ENDIF
[1361]1828
[2292]1829
[1115]1830!
[1822]1831!--    Reset precipitation rate
1832       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
1833
1834!
[1115]1835!--    Compute cloud physics
[1822]1836       IF( microphysics_kessler )  THEN
1837
1838          CALL autoconversion_kessler( i,j )
[1831]1839          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
[1822]1840
1841       ELSEIF ( microphysics_seifert )  THEN
1842
1843          CALL adjust_cloud( i,j )
[2292]1844          IF ( microphysics_morrison ) CALL activation( i,j )
1845          IF ( microphysics_morrison ) CALL condensation( i,j )
[1115]1846          CALL autoconversion( i,j )
1847          CALL accretion( i,j )
1848          CALL selfcollection_breakup( i,j )
1849          CALL evaporation_rain( i,j )
1850          CALL sedimentation_rain( i,j )
[1831]1851          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
[1115]1852
[1691]1853       ENDIF
1854
[1822]1855       CALL calc_precipitation_amount( i,j )
1856
[1115]1857!
[1361]1858!--    Store results on the 3d arrays
1859       q(:,j,i)  = q_1d(:)
1860       pt(:,j,i) = pt_1d(:)
[2292]1861       IF ( microphysics_morrison ) THEN
1862          qc(:,j,i) = qc_1d(:)
1863          nc(:,j,i) = nc_1d(:)
1864       ENDIF
[1822]1865       IF ( microphysics_seifert )  THEN
[1361]1866          qr(:,j,i) = qr_1d(:)
1867          nr(:,j,i) = nr_1d(:)
[1115]1868       ENDIF
1869
1870    END SUBROUTINE microphysics_control_ij
1871
[1682]1872!------------------------------------------------------------------------------!
1873! Description:
1874! ------------
[2155]1875!> Adjust number of raindrops to avoid nonlinear effects in
[1682]1876!> sedimentation and evaporation of rain drops due to too small or
1877!> too big weights of rain drops (Stevens and Seifert, 2008).
1878!> The same procedure is applied to cloud droplets if they are determined
1879!> prognostically. Call for grid point i,j
1880!------------------------------------------------------------------------------!
[1115]1881    SUBROUTINE adjust_cloud_ij( i, j )
1882
[1320]1883       USE cloud_parameters,                                                   &
[1849]1884           ONLY:  hyrho
[1320]1885
[2292]1886       USE control_parameters,                                                 &
1887           ONLY: microphysics_morrison
1888
[1320]1889       USE indices,                                                            &
[2232]1890           ONLY:  nzb, nzt, wall_flags_0
[1320]1891
1892       USE kinds
1893
[1115]1894       IMPLICIT NONE
1895
[1682]1896       INTEGER(iwp) ::  i                 !<
1897       INTEGER(iwp) ::  j                 !<
1898       INTEGER(iwp) ::  k                 !<
1899
[2232]1900       REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
[1022]1901
[2232]1902       DO  k = nzb+1, nzt
1903!
1904!--       Predetermine flag to mask topography
1905          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1906
[1361]1907          IF ( qr_1d(k) <= eps_sb )  THEN
1908             qr_1d(k) = 0.0_wp
1909             nr_1d(k) = 0.0_wp
[1065]1910          ELSE
[1022]1911!
[2155]1912!--          Adjust number of raindrops to avoid nonlinear effects in
[1048]1913!--          sedimentation and evaporation of rain drops due to too small or
[1065]1914!--          too big weights of rain drops (Stevens and Seifert, 2008).
[1361]1915             IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) )  THEN
[2232]1916                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin * flag
[1361]1917             ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) )  THEN
[2232]1918                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax * flag
[1048]1919             ENDIF
[1115]1920
[1022]1921          ENDIF
[1115]1922
[2292]1923          IF ( microphysics_morrison ) THEN
1924             IF ( qc_1d(k) <= eps_sb )  THEN
1925                qc_1d(k) = 0.0_wp
1926                nc_1d(k) = 0.0_wp
1927             ELSE
1928                IF ( nc_1d(k) * xcmin > qc_1d(k) * hyrho(k) )  THEN
1929                   nc_1d(k) = qc_1d(k) * hyrho(k) / xamin * flag
1930                ENDIF
1931             ENDIF
1932          ENDIF
1933
[1022]1934       ENDDO
1935
[1115]1936    END SUBROUTINE adjust_cloud_ij
[1022]1937
[2292]1938!------------------------------------------------------------------------------!
1939! Description:
1940! ------------
1941!> Calculate number of activated condensation nucleii after simple activation
1942!> scheme of Twomey, 1959.
1943!------------------------------------------------------------------------------!
1944    SUBROUTINE activation_ij( i, j )
[1106]1945
[2292]1946       USE arrays_3d,                                                          &
1947           ONLY:  hyp, nr, pt, q,  qc, qr, nc
1948
1949       USE cloud_parameters,                                                   &
[2375]1950           ONLY:  hyrho, l_d_cp, l_d_r, l_v, molecular_weight_of_solute,       &
1951                  molecular_weight_of_water, rho_l, rho_s, r_v, t_d_pt,        &
1952                  vanthoff
[2292]1953
1954       USE constants,                                                          &
1955           ONLY:  pi
1956
1957       USE cpulog,                                                             &
1958           ONLY:  cpu_log, log_point_s
1959
1960       USE indices,                                                            &
1961           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
1962
1963       USE kinds
1964
1965       USE control_parameters,                                                 &
1966          ONLY: simulated_time
1967
1968       IMPLICIT NONE
1969
1970       INTEGER(iwp) ::  i                 !<
1971       INTEGER(iwp) ::  j                 !<
1972       INTEGER(iwp) ::  k                 !<
1973
1974       REAL(wp)     ::  activ             !<
1975       REAL(wp)     ::  afactor           !<
1976       REAL(wp)     ::  alpha             !<
1977       REAL(wp)     ::  beta_act          !<
1978       REAL(wp)     ::  bfactor           !<
1979       REAL(wp)     ::  e_s               !<
1980       REAL(wp)     ::  k_act             !<
1981       REAL(wp)     ::  n_act             !<
1982       REAL(wp)     ::  n_ccn             !<
1983       REAL(wp)     ::  q_s               !<
1984       REAL(wp)     ::  s_0               !<
1985       REAL(wp)     ::  sat               !<
1986       REAL(wp)     ::  sat_max           !<
1987       REAL(wp)     ::  sigma             !<
[2375]1988       REAL(wp)     ::  sigma_act         !<
[2292]1989       REAL(wp)     ::  t_int             !<
1990       REAL(wp)     ::  t_l               !<
1991
1992       DO  k = nzb+1, nzt
1993!
1994!--       Actual liquid water temperature:
1995          t_l = t_d_pt(k) * pt_1d(k)
1996
1997!
1998!--       Calculate actual temperature
1999          t_int = pt_1d(k) * ( hyp(k) / 100000.0_wp )**0.286_wp
2000!
2001!--             Saturation vapor pressure at t_l:
2002          e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /             &
2003                                 ( t_l - 35.86_wp )                            &
2004                                 )
2005!
[3026]2006!--       Computation of saturation mixing ratio:
2007          q_s   = 0.622_wp * e_s / ( hyp(k) -  e_s )
[2292]2008          alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
2009          q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) /                         &
2010                        ( 1.0_wp + alpha * q_s )
2011
2012!--       Supersaturation:
2013          sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
2014
2015!
2016!--       Prescribe parameters for activation
2017!--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
[2312]2018          k_act  = 0.7_wp
[2522]2019          activ  = 0.0_wp
[2292]2020
2021          IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
2022!
[2312]2023!--          Compute the number of activated Aerosols
[2292]2024!--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
2025             n_act     = na_init * sat**k_act
2026!
2027!--          Compute the number of cloud droplets
2028!--          (see: Morrison + Grabowski, 2007, JAS, 64)
2029!            activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
2030
2031!
2032!--          Compute activation rate after Khairoutdinov and Kogan
[2312]2033!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
2034             sat_max = 0.8_wp / 100.0_wp
2035             activ   = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN              &
2036                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) /          &
[2292]2037                  dt_micro
2038
2039             nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
2040          ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk )  THEN
2041!
[2312]2042!--          Curvature effect (afactor) with surface tension
[2292]2043!--          parameterization by Straka (2009)
2044             sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
2045             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
2046!
2047!--          Solute effect (bfactor)
2048             bfactor = vanthoff * molecular_weight_of_water *                  &
2049                  rho_s / ( molecular_weight_of_solute * rho_l )
2050
2051!
[2312]2052!--          Prescribe power index that describes the soluble fraction
[2375]2053!--          of an aerosol particle (beta).
[2292]2054!--          (see: Morrison + Grabowski, 2007, JAS, 64)
[2375]2055             beta_act  = 0.5_wp
2056             sigma_act = sigma_bulk**( 1.0_wp + beta_act )
[2292]2057!
[2312]2058!--          Calculate mean geometric supersaturation (s_0) with
[2292]2059!--          parameterization by Khvorostyanov and Curry (2006)
[2375]2060             s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *         &
[2292]2061               ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
2062
2063!
[2312]2064!--          Calculate number of activated CCN as a function of
[2292]2065!--          supersaturation and taking Koehler theory into account
2066!--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
[2375]2067             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                    &
2068                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
[2292]2069             activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
[2312]2070
2071             nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
[2292]2072          ENDIF
2073
2074       ENDDO
2075
2076    END SUBROUTINE activation_ij
2077
[1682]2078!------------------------------------------------------------------------------!
2079! Description:
2080! ------------
[2312]2081!> Calculate condensation rate for cloud water content (after Khairoutdinov and
[2292]2082!> Kogan, 2000).
2083!------------------------------------------------------------------------------!
2084    SUBROUTINE condensation_ij( i, j )
2085
2086       USE arrays_3d,                                                          &
2087           ONLY:  hyp, nr, pt, q,  qc, qr, nc
2088
2089       USE cloud_parameters,                                                   &
2090           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
2091
2092       USE constants,                                                          &
2093           ONLY:  pi
2094
2095       USE cpulog,                                                             &
2096           ONLY:  cpu_log, log_point_s
2097
2098       USE indices,                                                            &
2099           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
2100
2101       USE kinds
2102
2103       USE control_parameters,                                                 &
2104          ONLY: simulated_time
2105
2106       IMPLICIT NONE
2107
2108       INTEGER(iwp) ::  i                 !<
2109       INTEGER(iwp) ::  j                 !<
2110       INTEGER(iwp) ::  k                 !<
2111
2112       REAL(wp)     ::  alpha             !<
2113       REAL(wp)     ::  cond              !<
2114       REAL(wp)     ::  cond_max          !<
2115       REAL(wp)     ::  dc                !<
2116       REAL(wp)     ::  e_s               !<
2117       REAL(wp)     ::  evap              !<
2118       REAL(wp)     ::  evap_nc           !<
2119       REAL(wp)     ::  g_fac             !<
2120       REAL(wp)     ::  nc_0              !<
2121       REAL(wp)     ::  q_s               !<
2122       REAL(wp)     ::  sat               !<
2123       REAL(wp)     ::  t_l               !<
2124       REAL(wp)     ::  temp              !<
2125       REAL(wp)     ::  xc                !<
2126
2127
2128       DO  k = nzb+1, nzt
2129!
2130!--       Actual liquid water temperature:
2131          t_l = t_d_pt(k) * pt_1d(k)
2132!
2133!--             Saturation vapor pressure at t_l:
2134          e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
2135                                 ( t_l - 35.86_wp )                   &
2136                                 )
2137!
[3026]2138!--       Computation of saturation mixing ratio:
2139          q_s   = 0.622_wp * e_s / ( hyp(k) - e_s )
[2292]2140          alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
2141          q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) /               &
2142                        ( 1.0_wp + alpha * q_s )
2143
2144!--       Supersaturation:
2145          sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
2146
2147
2148!
2149!--       Actual temperature:
2150          temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
2151
2152          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
2153                              l_v / ( thermal_conductivity_l * temp )    &
2154                              + r_v * temp / ( diff_coeff_l * e_s )      &
2155                            )
2156!
2157!--       Mean weight of cloud drops
2158          IF ( nc_1d(k) <= 0.0_wp) CYCLE
[2312]2159          xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin)
[2292]2160!
2161!--       Weight averaged diameter of cloud drops:
2162          dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2163!
2164!--       Integral diameter of cloud drops
[2312]2165          nc_0 = nc_1d(k) * dc
[2292]2166!
2167!--       Condensation needs only to be calculated in supersaturated regions
2168          IF ( sat > 0.0_wp )  THEN
2169!
2170!--          Condensation rate of cloud water content
2171!--          after KK scheme.
2172!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
2173             cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2174             cond_max  = q_1d(k) - q_s - qc_1d(k) - qr_1d(k)
2175             cond      = MIN( cond, cond_max / dt_micro )
[2312]2176
[2292]2177             qc_1d(k) = qc_1d(k) + cond * dt_micro
2178          ELSEIF ( sat < 0.0_wp ) THEN
2179             evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2180             evap      = MAX( evap, -qc_1d(k) / dt_micro )
2181
2182             qc_1d(k) = qc_1d(k) + evap * dt_micro
2183          ENDIF
2184       ENDDO
2185
2186    END SUBROUTINE condensation_ij
2187
2188
2189!------------------------------------------------------------------------------!
2190! Description:
2191! ------------
[1682]2192!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
2193!------------------------------------------------------------------------------!
[1005]2194    SUBROUTINE autoconversion_ij( i, j )
[1000]2195
[1320]2196       USE arrays_3d,                                                          &
[1849]2197           ONLY:  diss, dzu
[1115]2198
[1320]2199       USE cloud_parameters,                                                   &
[1849]2200           ONLY:  hyrho
[1320]2201
2202       USE control_parameters,                                                 &
[2292]2203           ONLY:  microphysics_morrison, rho_surface
[1320]2204
2205       USE grid_variables,                                                     &
2206           ONLY:  dx, dy
2207
2208       USE indices,                                                            &
[2232]2209           ONLY:  nzb, nzt, wall_flags_0
[1320]2210
2211       USE kinds
2212
[1000]2213       IMPLICIT NONE
2214
[1682]2215       INTEGER(iwp) ::  i                 !<
2216       INTEGER(iwp) ::  j                 !<
2217       INTEGER(iwp) ::  k                 !<
[1000]2218
[2155]2219       REAL(wp)     ::  alpha_cc          !<
[1682]2220       REAL(wp)     ::  autocon           !<
2221       REAL(wp)     ::  dissipation       !<
[2232]2222       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
[1682]2223       REAL(wp)     ::  k_au              !<
2224       REAL(wp)     ::  l_mix             !<
[2292]2225       REAL(wp)     ::  nc_auto           !<
[1682]2226       REAL(wp)     ::  nu_c              !<
2227       REAL(wp)     ::  phi_au            !<
2228       REAL(wp)     ::  r_cc              !<
2229       REAL(wp)     ::  rc                !<
2230       REAL(wp)     ::  re_lambda         !<
2231       REAL(wp)     ::  sigma_cc          !<
2232       REAL(wp)     ::  tau_cloud         !<
2233       REAL(wp)     ::  xc                !<
[1106]2234
[2232]2235       DO  k = nzb+1, nzt
2236!
2237!--       Predetermine flag to mask topography
2238          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2292]2239          nc_auto = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
[1000]2240
[2522]2241          IF ( qc_1d(k) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
[1361]2242
2243             k_au = k_cc / ( 20.0_wp * x0 )
[1012]2244!
[1048]2245!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]2246!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
[2522]2247             tau_cloud = MAX( 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) ),     & 
2248                              0.0_wp )
[1012]2249!
[2155]2250!--          Universal function for autoconversion process
[1012]2251!--          (Seifert and Beheng, 2006):
[1361]2252             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
[1012]2253!
2254!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
[1353]2255!--          (Use constant nu_c = 1.0_wp instead?)
[1361]2256             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
[1012]2257!
2258!--          Mean weight of cloud droplets:
[2292]2259             xc = hyrho(k) * qc_1d(k) / nc_auto
[1012]2260!
[2155]2261!--          Parameterized turbulence effects on autoconversion (Seifert,
[1065]2262!--          Nuijens and Stevens, 2010)
[1831]2263             IF ( collision_turbulence )  THEN
[1065]2264!
2265!--             Weight averaged radius of cloud droplets:
[1353]2266                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
[1065]2267
[1353]2268                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
2269                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
2270                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
[1065]2271!
2272!--             Mixing length (neglecting distance to ground and stratification)
[1334]2273                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
[1065]2274!
[2155]2275!--             Limit dissipation rate according to Seifert, Nuijens and
[1065]2276!--             Stevens (2010)
[1361]2277                dissipation = MIN( 0.06_wp, diss(k,j,i) )
[1065]2278!
2279!--             Compute Taylor-microscale Reynolds number:
[1361]2280                re_lambda = 6.0_wp / 11.0_wp *                                 &
2281                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
2282                            SQRT( 15.0_wp / kin_vis_air ) *                    &
2283                            dissipation**( 1.0_wp / 6.0_wp )
[1065]2284!
2285!--             The factor of 1.0E4 is needed to convert the dissipation rate
2286!--             from m2 s-3 to cm2 s-3.
[1361]2287                k_au = k_au * ( 1.0_wp +                                       &
2288                                dissipation * 1.0E4_wp *                       &
2289                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
2290                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
2291                                                  sigma_cc )**2                &
2292                                                ) + beta_cc                    &
2293                                )                                              &
2294                              )
[1065]2295             ENDIF
2296!
[1012]2297!--          Autoconversion rate (Seifert and Beheng, 2006):
[1361]2298             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
2299                       ( nu_c + 1.0_wp )**2 * qc_1d(k)**2 * xc**2 *            &
2300                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
[1115]2301                       rho_surface
2302             autocon = MIN( autocon, qc_1d(k) / dt_micro )
[1106]2303
[2232]2304             qr_1d(k) = qr_1d(k) + autocon * dt_micro                 * flag
2305             qc_1d(k) = qc_1d(k) - autocon * dt_micro                 * flag
2306             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag
[2292]2307             IF ( microphysics_morrison ) THEN
[2312]2308                nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp *                &
[2292]2309                                    autocon / x0 * hyrho(k) * dt_micro * flag )
2310             ENDIF
[1115]2311
[1005]2312          ENDIF
[1000]2313
2314       ENDDO
2315
[1005]2316    END SUBROUTINE autoconversion_ij
2317
[1822]2318!------------------------------------------------------------------------------!
2319! Description:
2320! ------------
2321!> Autoconversion process (Kessler, 1969).
2322!------------------------------------------------------------------------------!
2323    SUBROUTINE autoconversion_kessler_ij( i, j )
[1106]2324
[1822]2325       USE arrays_3d,                                                          &
[1849]2326           ONLY:  dzw, prr
[1822]2327
2328       USE cloud_parameters,                                                   &
[1849]2329           ONLY:  l_d_cp, pt_d_t
[1822]2330
2331       USE indices,                                                            &
[2317]2332           ONLY:  nzb, nzt, wall_flags_0
[1822]2333
2334       USE kinds
2335
[2317]2336       USE surface_mod,                                                        &
[2698]2337           ONLY:  get_topography_top_index_ji
[1822]2338
[2317]2339
[1822]2340       IMPLICIT NONE
2341
[2232]2342       INTEGER(iwp) ::  i      !<
2343       INTEGER(iwp) ::  j      !<
2344       INTEGER(iwp) ::  k      !<
2345       INTEGER(iwp) ::  k_wall !< topography top index
[1822]2346
2347       REAL(wp)    ::  dqdt_precip !<
[2232]2348       REAL(wp)    ::  flag              !< flag to indicate first grid level above surface
[1822]2349
[2232]2350!
2351!--    Determine vertical index of topography top
[2698]2352       k_wall = get_topography_top_index_ji( j, i, 's' )
[2232]2353       DO  k = nzb+1, nzt
2354!
2355!--       Predetermine flag to mask topography
2356          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1822]2357
2358          IF ( qc_1d(k) > ql_crit )  THEN
2359             dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit )
2360          ELSE
2361             dqdt_precip = 0.0_wp
2362          ENDIF
2363
[2232]2364          qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro * flag
2365          q_1d(k)  = q_1d(k)  - dqdt_precip * dt_micro * flag
2366          pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k) * flag
[1822]2367
2368!
[1845]2369!--       Compute the rain rate (stored on surface grid point)
[2232]2370          prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
[1822]2371
2372       ENDDO
2373
2374    END SUBROUTINE autoconversion_kessler_ij
2375
[1682]2376!------------------------------------------------------------------------------!
2377! Description:
2378! ------------
2379!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
2380!------------------------------------------------------------------------------!
[1005]2381    SUBROUTINE accretion_ij( i, j )
2382
[1320]2383       USE arrays_3d,                                                          &
[1849]2384           ONLY:  diss
[1115]2385
[1320]2386       USE cloud_parameters,                                                   &
[1849]2387           ONLY:  hyrho
[1320]2388
2389       USE control_parameters,                                                 &
[2292]2390           ONLY:  microphysics_morrison, rho_surface
[1320]2391
2392       USE indices,                                                            &
[2232]2393           ONLY:  nzb, nzt, wall_flags_0
[1320]2394
2395       USE kinds
2396
[1005]2397       IMPLICIT NONE
2398
[1682]2399       INTEGER(iwp) ::  i                 !<
2400       INTEGER(iwp) ::  j                 !<
2401       INTEGER(iwp) ::  k                 !<
[1005]2402
[1682]2403       REAL(wp)     ::  accr              !<
[2232]2404       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
[1682]2405       REAL(wp)     ::  k_cr              !<
[2292]2406       REAL(wp)     ::  nc_accr           !<
[1682]2407       REAL(wp)     ::  phi_ac            !<
2408       REAL(wp)     ::  tau_cloud         !<
[2292]2409       REAL(wp)     ::  xc                !<
[1320]2410
[2292]2411
[2232]2412       DO  k = nzb+1, nzt
2413!
2414!--       Predetermine flag to mask topography
2415          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2292]2416          nc_accr = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
[2232]2417
[2292]2418          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb )  .AND.  &
2419               ( nc_accr > eps_mr ) )  THEN
[1012]2420!
[1048]2421!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
[2155]2422             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
[1012]2423!
[2155]2424!--          Universal function for accretion process
[1048]2425!--          (Seifert and Beheng, 2001):
[1361]2426             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
[2292]2427
[1012]2428!
[2292]2429!--          Mean weight of cloud drops
[2312]2430             xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin)
[2292]2431!
[2155]2432!--          Parameterized turbulence effects on autoconversion (Seifert,
2433!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
[1361]2434!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1831]2435             IF ( collision_turbulence )  THEN
[1361]2436                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
2437                                 MIN( 600.0_wp,                                &
2438                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
2439                               )
[1065]2440             ELSE
[2155]2441                k_cr = k_cr0
[1065]2442             ENDIF
2443!
[1012]2444!--          Accretion rate (Seifert and Beheng, 2006):
[2292]2445             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac *                      &
2446                    SQRT( rho_surface * hyrho(k) )
[1115]2447             accr = MIN( accr, qc_1d(k) / dt_micro )
[1106]2448
[2232]2449             qr_1d(k) = qr_1d(k) + accr * dt_micro * flag
2450             qc_1d(k) = qc_1d(k) - accr * dt_micro * flag
[2292]2451             IF ( microphysics_morrison )  THEN
2452                nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), accr / xc *               &
2453                                             hyrho(k) * dt_micro * flag        &
2454                                         )
2455             ENDIF
[1115]2456
[2292]2457
[1005]2458          ENDIF
[1106]2459
[1005]2460       ENDDO
2461
[1000]2462    END SUBROUTINE accretion_ij
2463
[1005]2464
[1682]2465!------------------------------------------------------------------------------!
2466! Description:
2467! ------------
2468!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
2469!------------------------------------------------------------------------------!
[1005]2470    SUBROUTINE selfcollection_breakup_ij( i, j )
2471
[1320]2472       USE cloud_parameters,                                                   &
[1849]2473           ONLY:  hyrho
[1320]2474
2475       USE control_parameters,                                                 &
[1849]2476           ONLY:  rho_surface
[1320]2477
2478       USE indices,                                                            &
[2232]2479           ONLY:  nzb, nzt, wall_flags_0
[1320]2480
2481       USE kinds
[2155]2482
[1005]2483       IMPLICIT NONE
2484
[1682]2485       INTEGER(iwp) ::  i                 !<
2486       INTEGER(iwp) ::  j                 !<
2487       INTEGER(iwp) ::  k                 !<
[1005]2488
[1682]2489       REAL(wp)     ::  breakup           !<
2490       REAL(wp)     ::  dr                !<
[2232]2491       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
[1682]2492       REAL(wp)     ::  phi_br            !<
2493       REAL(wp)     ::  selfcoll          !<
[1320]2494
[2232]2495       DO  k = nzb+1, nzt
2496!
2497!--       Predetermine flag to mask topography
2498          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2499
[1115]2500          IF ( qr_1d(k) > eps_sb )  THEN
[1012]2501!
[1115]2502!--          Selfcollection rate (Seifert and Beheng, 2001):
[1361]2503             selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
[1012]2504!
[1115]2505!--          Weight averaged diameter of rain drops:
[1334]2506             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]2507!
[1048]2508!--          Collisional breakup rate (Seifert, 2008):
[1353]2509             IF ( dr >= 0.3E-3_wp )  THEN
2510                phi_br  = k_br * ( dr - 1.1E-3_wp )
2511                breakup = selfcoll * ( phi_br + 1.0_wp )
[1005]2512             ELSE
[1353]2513                breakup = 0.0_wp
[1005]2514             ENDIF
[1048]2515
[1115]2516             selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
[2232]2517             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro * flag
[1106]2518
[2155]2519          ENDIF
[1005]2520       ENDDO
2521
2522    END SUBROUTINE selfcollection_breakup_ij
2523
[1106]2524
[1682]2525!------------------------------------------------------------------------------!
2526! Description:
2527! ------------
[2155]2528!> Evaporation of precipitable water. Condensation is neglected for
[1682]2529!> precipitable water. Call for grid point i,j
2530!------------------------------------------------------------------------------!
[1012]2531    SUBROUTINE evaporation_rain_ij( i, j )
2532
[1320]2533       USE arrays_3d,                                                          &
[1849]2534           ONLY:  hyp
[1048]2535
[1320]2536       USE cloud_parameters,                                                   &
[1849]2537           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
[1320]2538
2539       USE constants,                                                          &
2540           ONLY:  pi
2541
2542       USE indices,                                                            &
[2232]2543           ONLY:  nzb, nzt, wall_flags_0
[1320]2544
2545       USE kinds
2546
[1012]2547       IMPLICIT NONE
2548
[1682]2549       INTEGER(iwp) ::  i                 !<
2550       INTEGER(iwp) ::  j                 !<
2551       INTEGER(iwp) ::  k                 !<
[1012]2552
[1682]2553       REAL(wp)     ::  alpha             !<
2554       REAL(wp)     ::  dr                !<
2555       REAL(wp)     ::  e_s               !<
2556       REAL(wp)     ::  evap              !<
2557       REAL(wp)     ::  evap_nr           !<
2558       REAL(wp)     ::  f_vent            !<
[2232]2559       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
[1682]2560       REAL(wp)     ::  g_evap            !<
2561       REAL(wp)     ::  lambda_r          !<
2562       REAL(wp)     ::  mu_r              !<
2563       REAL(wp)     ::  mu_r_2            !<
2564       REAL(wp)     ::  mu_r_5d2          !<
2565       REAL(wp)     ::  nr_0              !<
2566       REAL(wp)     ::  q_s               !<
2567       REAL(wp)     ::  sat               !<
2568       REAL(wp)     ::  t_l               !<
2569       REAL(wp)     ::  temp              !<
2570       REAL(wp)     ::  xr                !<
[1320]2571
[2232]2572       DO  k = nzb+1, nzt
2573!
2574!--       Predetermine flag to mask topography
2575          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2576
[1115]2577          IF ( qr_1d(k) > eps_sb )  THEN
[1012]2578!
2579!--          Actual liquid water temperature:
[1115]2580             t_l = t_d_pt(k) * pt_1d(k)
[1012]2581!
2582!--          Saturation vapor pressure at t_l:
[1361]2583             e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /          &
2584                                    ( t_l - 35.86_wp )                         &
2585                                  )
[1012]2586!
[3026]2587!--          Computation of saturation mixing ratio:
2588             q_s   = 0.622_wp * e_s / ( hyp(k) - e_s )
[1353]2589             alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
[1361]2590             q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
[1012]2591!
[1106]2592!--          Supersaturation:
[1361]2593             sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
[1012]2594!
[1361]2595!--          Evaporation needs only to be calculated in subsaturated regions
2596             IF ( sat < 0.0_wp )  THEN
[1012]2597!
[1361]2598!--             Actual temperature:
2599                temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
[2155]2600
[1361]2601                g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /  &
2602                                    ( thermal_conductivity_l * temp ) +        &
2603                                    r_v * temp / ( diff_coeff_l * e_s )        &
2604                                  )
[1012]2605!
[1361]2606!--             Mean weight of rain drops
2607                xr = hyrho(k) * qr_1d(k) / nr_1d(k)
[1115]2608!
[1361]2609!--             Weight averaged diameter of rain drops:
2610                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]2611!
[2155]2612!--             Compute ventilation factor and intercept parameter
[1361]2613!--             (Seifert and Beheng, 2006; Seifert, 2008):
2614                IF ( ventilation_effect )  THEN
[1115]2615!
[1361]2616!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
2617!--                Stevens and Seifert, 2008):
2618                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
2619!
2620!--                Slope parameter of gamma distribution (Seifert, 2008):
2621                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
2622                                ( mu_r + 1.0_wp )                              &
2623                              )**( 1.0_wp / 3.0_wp ) / dr
[1115]2624
[1361]2625                   mu_r_2   = mu_r + 2.0_wp
[2155]2626                   mu_r_5d2 = mu_r + 2.5_wp
[1361]2627
[2155]2628                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  &
[1361]2629                            b_vent * schmidt_p_1d3 *                           &
2630                            SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
2631                            lambda_r**( -mu_r_5d2 ) *                          &
2632                            ( 1.0_wp -                                         &
2633                              0.5_wp * ( b_term / a_term ) *                   &
2634                              ( lambda_r / ( c_term + lambda_r )               &
2635                              )**mu_r_5d2 -                                    &
2636                              0.125_wp * ( b_term / a_term )**2 *              &
2637                              ( lambda_r / ( 2.0_wp * c_term + lambda_r )      &
2638                              )**mu_r_5d2 -                                    &
2639                              0.0625_wp * ( b_term / a_term )**3 *             &
2640                              ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
2641                              )**mu_r_5d2 -                                    &
[2155]2642                              0.0390625_wp * ( b_term / a_term )**4 *          &
[1361]2643                              ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
2644                              )**mu_r_5d2                                      &
2645                            )
2646
2647                   nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) /           &
[2155]2648                            gamm( mu_r + 1.0_wp )
[1361]2649                ELSE
2650                   f_vent = 1.0_wp
2651                   nr_0   = nr_1d(k) * dr
2652                ENDIF
[1012]2653!
[1361]2654!--             Evaporation rate of rain water content (Seifert and Beheng, 2006):
2655                evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
2656                evap    = MAX( evap, -qr_1d(k) / dt_micro )
2657                evap_nr = MAX( c_evap * evap / xr * hyrho(k),                  &
2658                               -nr_1d(k) / dt_micro )
[1106]2659
[2232]2660                qr_1d(k) = qr_1d(k) + evap    * dt_micro * flag
2661                nr_1d(k) = nr_1d(k) + evap_nr * dt_micro * flag
[1115]2662
[1361]2663             ENDIF
[2155]2664          ENDIF
[1106]2665
[1012]2666       ENDDO
2667
2668    END SUBROUTINE evaporation_rain_ij
2669
[1106]2670
[1682]2671!------------------------------------------------------------------------------!
2672! Description:
2673! ------------
[2155]2674!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
[1682]2675!> Call for grid point i,j
2676!------------------------------------------------------------------------------!
[1012]2677    SUBROUTINE sedimentation_cloud_ij( i, j )
2678
[1320]2679       USE arrays_3d,                                                          &
[1849]2680           ONLY:  ddzu, dzu, prr
[1320]2681
2682       USE cloud_parameters,                                                   &
[1849]2683           ONLY:  hyrho, l_d_cp, pt_d_t
[1320]2684
2685       USE control_parameters,                                                 &
[2292]2686           ONLY:  call_microphysics_at_all_substeps,                           &
2687                  intermediate_timestep_count, microphysics_morrison
[1320]2688
2689       USE indices,                                                            &
[2232]2690           ONLY:  nzb, nzb, nzt, wall_flags_0
[1320]2691
2692       USE kinds
[2155]2693
[1691]2694       USE statistics,                                                         &
2695           ONLY:  weight_substep
2696
[1012]2697       IMPLICIT NONE
2698
[1849]2699       INTEGER(iwp) ::  i             !<
2700       INTEGER(iwp) ::  j             !<
2701       INTEGER(iwp) ::  k             !<
[1106]2702
[2292]2703       REAL(wp)     ::  flag    !< flag to indicate first grid level above surface
2704       REAL(wp)     ::  nc_sedi !<
2705
2706       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc  !<
[2232]2707       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc  !<
[1115]2708
[1353]2709       sed_qc(nzt+1) = 0.0_wp
[2375]2710       sed_nc(nzt+1) = 0.0_wp
[1012]2711
[2375]2712
[2232]2713       DO  k = nzt, nzb+1, -1
2714!
2715!--       Predetermine flag to mask topography
2716          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2292]2717          nc_sedi = MERGE( nc_1d(k), nc_const, microphysics_morrison )
2718!
2719!--       Sedimentation fluxes for number concentration are only calculated
2720!--       for cloud_scheme = 'morrison'
2721          IF ( microphysics_morrison ) THEN
2722             IF ( qc_1d(k) > eps_sb  .AND.  nc_1d(k) > eps_mr )  THEN
[2312]2723                sed_nc(k) = sed_qc_const *                                     &
[2292]2724                            ( qc_1d(k) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *     &
2725                            ( nc_1d(k) )**( 1.0_wp / 3.0_wp )
2726             ELSE
2727                sed_nc(k) = 0.0_wp
2728             ENDIF
[2232]2729
[2292]2730             sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *                 &
2731                              nc_1d(k) / dt_micro + sed_nc(k+1)                &
2732                            ) * flag
2733
2734             nc_1d(k) = nc_1d(k) + ( sed_nc(k+1) - sed_nc(k) ) *               &
2735                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
2736          ENDIF
2737
[2375]2738          IF ( qc_1d(k) > eps_sb  .AND.  nc_sedi > eps_mr )  THEN
[2292]2739             sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *        &
[2232]2740                         ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )  * flag
[1115]2741          ELSE
[1353]2742             sed_qc(k) = 0.0_wp
[1012]2743          ENDIF
[1115]2744
[1361]2745          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /          &
2746                                      dt_micro + sed_qc(k+1)                   &
[2232]2747                         ) * flag
[1115]2748
[1361]2749          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
[2232]2750                                hyrho(k) * dt_micro * flag
[2155]2751          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
[2232]2752                                hyrho(k) * dt_micro * flag
[1361]2753          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
[2232]2754                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
[1115]2755
[1691]2756!
2757!--       Compute the precipitation rate of cloud (fog) droplets
[1822]2758          IF ( call_microphysics_at_all_substeps )  THEN
[2232]2759             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *                  &
2760                              weight_substep(intermediate_timestep_count) * flag
[1822]2761          ELSE
[2232]2762             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
[1691]2763          ENDIF
2764
[1012]2765       ENDDO
2766
2767    END SUBROUTINE sedimentation_cloud_ij
2768
[1106]2769
[1682]2770!------------------------------------------------------------------------------!
2771! Description:
2772! ------------
2773!> Computation of sedimentation flux. Implementation according to Stevens
2774!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
2775!------------------------------------------------------------------------------!
[1012]2776    SUBROUTINE sedimentation_rain_ij( i, j )
2777
[1320]2778       USE arrays_3d,                                                          &
[1849]2779           ONLY:  ddzu, dzu, prr
[1320]2780
2781       USE cloud_parameters,                                                   &
[1849]2782           ONLY:  hyrho, l_d_cp, pt_d_t
[1320]2783
2784       USE control_parameters,                                                 &
[1849]2785           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
[1320]2786
2787       USE indices,                                                            &
[2232]2788           ONLY:  nzb, nzb, nzt, wall_flags_0
[1320]2789
2790       USE kinds
2791
2792       USE statistics,                                                         &
2793           ONLY:  weight_substep
[2155]2794
[2232]2795       USE surface_mod,                                                        &
2796           ONLY :  bc_h
[2312]2797
[1012]2798       IMPLICIT NONE
2799
[2232]2800       INTEGER(iwp) ::  i             !< running index x direction
2801       INTEGER(iwp) ::  j             !< running index y direction
2802       INTEGER(iwp) ::  k             !< running index z direction
2803       INTEGER(iwp) ::  k_run         !<
2804       INTEGER(iwp) ::  m             !< running index surface elements
2805       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
2806       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
[1012]2807
[1682]2808       REAL(wp)     ::  c_run                      !<
2809       REAL(wp)     ::  d_max                      !<
2810       REAL(wp)     ::  d_mean                     !<
2811       REAL(wp)     ::  d_min                      !<
2812       REAL(wp)     ::  dr                         !<
2813       REAL(wp)     ::  flux                       !<
[2232]2814       REAL(wp)     ::  flag                       !< flag to indicate first grid level above surface
[1682]2815       REAL(wp)     ::  lambda_r                   !<
2816       REAL(wp)     ::  mu_r                       !<
2817       REAL(wp)     ::  z_run                      !<
[1320]2818
[1682]2819       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
2820       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
2821       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
2822       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
2823       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
2824       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
2825       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
2826       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
[1320]2827
[1012]2828!
[2155]2829!--    Compute velocities
[2232]2830       DO  k = nzb+1, nzt
2831!
2832!--       Predetermine flag to mask topography
2833          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2834
[1115]2835          IF ( qr_1d(k) > eps_sb )  THEN
2836!
2837!--          Weight averaged diameter of rain drops:
[1334]2838             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]2839!
2840!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
2841!--          Stevens and Seifert, 2008):
[1353]2842             mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
[1115]2843!
2844!--          Slope parameter of gamma distribution (Seifert, 2008):
[1361]2845             lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *              &
2846                          ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1115]2847
[1361]2848             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
2849                                         a_term - b_term * ( 1.0_wp +          &
2850                                            c_term / lambda_r )**( -1.0_wp *   &
2851                                            ( mu_r + 1.0_wp ) )                &
2852                                        )                                      &
[2232]2853                          ) * flag
[1361]2854             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
2855                                         a_term - b_term * ( 1.0_wp +          &
2856                                            c_term / lambda_r )**( -1.0_wp *   &
2857                                            ( mu_r + 4.0_wp ) )                &
2858                                       )                                       &
[2232]2859                          ) * flag
[1065]2860          ELSE
[1353]2861             w_nr(k) = 0.0_wp
2862             w_qr(k) = 0.0_wp
[1065]2863          ENDIF
2864       ENDDO
[1048]2865!
[2312]2866!--    Adjust boundary values using surface data type.
[2232]2867!--    Upward facing non-natural
[2312]2868       surf_s = bc_h(0)%start_index(j,i)
[2232]2869       surf_e = bc_h(0)%end_index(j,i)
2870       DO  m = surf_s, surf_e
2871          k         = bc_h(0)%k(m)
2872          w_nr(k-1) = w_nr(k)
2873          w_qr(k-1) = w_qr(k)
2874       ENDDO
2875!
2876!--    Downward facing non-natural
[2312]2877       surf_s = bc_h(1)%start_index(j,i)
[2232]2878       surf_e = bc_h(1)%end_index(j,i)
2879       DO  m = surf_s, surf_e
2880          k         = bc_h(1)%k(m)
2881          w_nr(k+1) = w_nr(k)
2882          w_qr(k+1) = w_qr(k)
2883       ENDDO
2884!
2885!--    Neumann boundary condition at model top
[1353]2886       w_nr(nzt+1) = 0.0_wp
2887       w_qr(nzt+1) = 0.0_wp
[1065]2888!
2889!--    Compute Courant number
[2232]2890       DO  k = nzb+1, nzt
2891!
2892!--       Predetermine flag to mask topography
2893          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2894
[1361]2895          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
[2232]2896                    dt_micro * ddzu(k) * flag
[1361]2897          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
[2232]2898                    dt_micro * ddzu(k) * flag
[2155]2899       ENDDO
[1065]2900!
2901!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
2902       IF ( limiter_sedimentation )  THEN
2903
[2232]2904          DO k = nzb+1, nzt
2905!
2906!--          Predetermine flag to mask topography
2907             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2908
[1646]2909             d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
[1115]2910             d_min  = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
2911             d_max  = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
[1065]2912
[1361]2913             qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
2914                                                        2.0_wp * d_max,        &
[2232]2915                                                        ABS( d_mean ) ) * flag
[1065]2916
[1646]2917             d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
[1115]2918             d_min  = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
2919             d_max  = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
[1065]2920
[1361]2921             nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
2922                                                        2.0_wp * d_max,        &
[2232]2923                                                        ABS( d_mean ) ) * flag
[1022]2924          ENDDO
[1048]2925
[1065]2926       ELSE
[1106]2927
[1353]2928          nr_slope = 0.0_wp
2929          qr_slope = 0.0_wp
[1106]2930
[1065]2931       ENDIF
[1115]2932
[1353]2933       sed_nr(nzt+1) = 0.0_wp
2934       sed_qr(nzt+1) = 0.0_wp
[1065]2935!
2936!--    Compute sedimentation flux
[2232]2937       DO  k = nzt, nzb+1, -1
[1065]2938!
[2232]2939!--       Predetermine flag to mask topography
2940          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2941!
[2312]2942!--       Sum up all rain drop number densities which contribute to the flux
[1065]2943!--       through k-1/2
[1353]2944          flux  = 0.0_wp
2945          z_run = 0.0_wp ! height above z(k)
[1065]2946          k_run = k
[1346]2947          c_run = MIN( 1.0_wp, c_nr(k) )
[1353]2948          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
[1361]2949             flux  = flux + hyrho(k_run) *                                     &
2950                     ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) *   &
[2232]2951                     0.5_wp ) * c_run * dzu(k_run) * flag
2952             z_run = z_run + dzu(k_run)            * flag
2953             k_run = k_run + 1                     * flag
2954             c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag
[1022]2955          ENDDO
2956!
[2155]2957!--       It is not allowed to sediment more rain drop number density than
[1065]2958!--       available
[1361]2959          flux = MIN( flux,                                                    &
[1115]2960                      hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
[1065]2961
[2232]2962          sed_nr(k) = flux / dt_micro * flag
[1361]2963          nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
[2232]2964                                    hyrho(k) * dt_micro * flag
[1065]2965!
[2155]2966!--       Sum up all rain water content which contributes to the flux
[1065]2967!--       through k-1/2
[1353]2968          flux  = 0.0_wp
2969          z_run = 0.0_wp ! height above z(k)
[1065]2970          k_run = k
[1346]2971          c_run = MIN( 1.0_wp, c_qr(k) )
[1106]2972
[1361]2973          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
[1106]2974
[1361]2975             flux  = flux + hyrho(k_run) *                                     &
2976                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) *   &
[2232]2977                     0.5_wp ) * c_run * dzu(k_run) * flag
2978             z_run = z_run + dzu(k_run)            * flag
2979             k_run = k_run + 1                     * flag
2980             c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag
[1106]2981
[1065]2982          ENDDO
2983!
2984!--       It is not allowed to sediment more rain water content than available
[1361]2985          flux = MIN( flux,                                                    &
[1115]2986                      hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
[1065]2987
[2232]2988          sed_qr(k) = flux / dt_micro * flag
[1115]2989
[1361]2990          qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
[2232]2991                                hyrho(k) * dt_micro * flag
[1361]2992          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
[2232]2993                                hyrho(k) * dt_micro * flag
[1361]2994          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
[2232]2995                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
[1065]2996!
2997!--       Compute the rain rate
[1361]2998          IF ( call_microphysics_at_all_substeps )  THEN
[1691]2999             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)                    &
[2232]3000                          * weight_substep(intermediate_timestep_count) * flag
[1361]3001          ELSE
[2232]3002             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
[1361]3003          ENDIF
3004
[1065]3005       ENDDO
[1115]3006
[1691]3007    END SUBROUTINE sedimentation_rain_ij
[1012]3008
[1691]3009
3010!------------------------------------------------------------------------------!
3011! Description:
3012! ------------
3013!> This subroutine computes the precipitation amount due to gravitational
3014!> settling of rain and cloud (fog) droplets
3015!------------------------------------------------------------------------------!
3016    SUBROUTINE calc_precipitation_amount_ij( i, j )
3017
[1849]3018       USE arrays_3d,                                                          &
3019           ONLY:  precipitation_amount, prr
3020
[1691]3021       USE cloud_parameters,                                                   &
[1849]3022           ONLY:  hyrho
[1691]3023
3024       USE control_parameters,                                                 &
3025           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
3026                  intermediate_timestep_count, intermediate_timestep_count_max,&
[1822]3027                  precipitation_amount_interval, time_do2d_xy
[1691]3028
3029       USE indices,                                                            &
[2232]3030           ONLY:  nzb, nzt, wall_flags_0
[1691]3031
3032       USE kinds
3033
[2232]3034       USE surface_mod,                                                        &
3035           ONLY :  bc_h
3036
[1691]3037       IMPLICIT NONE
3038
[2232]3039       INTEGER(iwp) ::  i             !< running index x direction
3040       INTEGER(iwp) ::  j             !< running index y direction
3041       INTEGER(iwp) ::  k             !< running index z direction
3042       INTEGER(iwp) ::  m             !< running index surface elements
3043       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
3044       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
[1691]3045
3046       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
3047            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
3048            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
3049       THEN
3050
[2312]3051          surf_s = bc_h(0)%start_index(j,i)
[2232]3052          surf_e = bc_h(0)%end_index(j,i)
3053          DO  m = surf_s, surf_e
3054             k                         = bc_h(0)%k(m)
3055             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
3056                                               prr(k,j,i) * hyrho(k) * dt_3d
3057          ENDDO
3058
[1048]3059       ENDIF
3060
[1691]3061    END SUBROUTINE calc_precipitation_amount_ij
[1012]3062
[1361]3063!------------------------------------------------------------------------------!
[1682]3064! Description:
3065! ------------
3066!> This function computes the gamma function (Press et al., 1992).
[2155]3067!> The gamma function is needed for the calculation of the evaporation
[1682]3068!> of rain drops.
[1361]3069!------------------------------------------------------------------------------!
[2155]3070    FUNCTION gamm( xx )
[1320]3071
3072       USE kinds
3073
[2155]3074       IMPLICIT NONE
[1106]3075
[1682]3076       INTEGER(iwp) ::  j            !<
[1320]3077
[1682]3078       REAL(wp)     ::  gamm         !<
3079       REAL(wp)     ::  ser          !<
3080       REAL(wp)     ::  tmp          !<
3081       REAL(wp)     ::  x_gamm       !<
3082       REAL(wp)     ::  xx           !<
3083       REAL(wp)     ::  y_gamm       !<
[1320]3084
[1849]3085
3086       REAL(wp), PARAMETER  ::  stp = 2.5066282746310005_wp               !<
3087       REAL(wp), PARAMETER  ::  cof(6) = (/ 76.18009172947146_wp,      &
3088                                           -86.50532032941677_wp,      &
3089                                            24.01409824083091_wp,      &
3090                                            -1.231739572450155_wp,     &
3091                                             0.1208650973866179E-2_wp, &
3092                                            -0.5395239384953E-5_wp /)     !<
3093
3094       x_gamm = xx
3095       y_gamm = x_gamm
[1353]3096       tmp = x_gamm + 5.5_wp
[1849]3097       tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
[1334]3098       ser = 1.000000000190015_wp
[1106]3099
[2155]3100       DO  j = 1, 6
3101          y_gamm = y_gamm + 1.0_wp
3102          ser    = ser + cof( j ) / y_gamm
[1106]3103       ENDDO
3104
[2155]3105!
3106!--    Until this point the algorithm computes the logarithm of the gamma
3107!--    function. Hence, the exponential function is used.
3108!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
3109       gamm = EXP( tmp ) * stp * ser / x_gamm
[1106]3110
[2155]3111       RETURN
[1012]3112
[2155]3113    END FUNCTION gamm
3114
[1012]3115 END MODULE microphysics_mod
Note: See TracBrowser for help on using the repository browser.