source: palm/trunk/SOURCE/bulk_cloud_model_mod.f90 @ 3623

Last change on this file since 3623 was 3622, checked in by schwenkel, 5 years ago

Important bugfix in case of restart runs

  • Property svn:keywords set to Id
File size: 162.5 KB
Line 
1!> @file bulk_cloud_model_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
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.
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!
17! Copyright 1997-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: bulk_cloud_model_mod.f90 3622 2018-12-12 09:52:53Z forkel $
27! Important bugfix in case of restart runs.
28!
29! 3507 2018-11-08 14:28:46Z schwenkel
30! Minor bugfixes for bcm cache version and initializing prr
31!
32! 3452 2018-10-30 13:13:34Z schwenkel
33! Bugfix for profiles output
34!
35! 3445 2018-10-29 12:23:02Z schwenkel
36! Minor bugfix and use of subroutine for supersaturation calculation in case
37! of cache version
38!
39! 3383 2018-10-19 14:22:58Z knoop
40! Modularization of all bulk cloud physics code components
41!
42! unused variables removed
43!
44! 3026 2018-05-22 10:30:53Z schwenkel
45! Changed the name specific humidity to mixing ratio, since we are computing
46! mixing ratios.
47!
48! 2718 2018-01-02 08:49:38Z maronga
49! Corrected "Former revisions" section
50!
51! 2701 2017-12-15 15:40:50Z suehring
52! Changes from last commit documented
53!
54! 2698 2017-12-14 18:46:24Z suehring
55! Bugfix in get_topography_top_index
56!
57! 2696 2017-12-14 17:12:51Z kanani
58! Change in file header (GPL part)
59!
60! 2608 2017-11-13 14:04:26Z schwenkel
61! Calculation of supersaturation in external module (diagnostic_quantities_mod).
62! Change: correct calculation of saturation specific humidity to saturation
63! mixing ratio (the factor of 0.378 vanishes).
64!
65! 2522 2017-10-05 14:20:37Z schwenkel
66! Minor bugfix
67!
68! 2375 2017-08-29 14:10:28Z schwenkel
69! Improved aerosol initilization and some minor bugfixes
70! for droplet sedimenation
71!
72! 2318 2017-07-20 17:27:44Z suehring
73! Get topography top index via Function call
74!
75! 2317 2017-07-20 17:27:19Z suehring
76! s1 changed to log_sigma
77!
78! 2292 2017-06-20 09:51:42Z schwenkel
79! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
80! includes two more prognostic equations for cloud drop concentration (nc)
81! and cloud water content (qc).
82! - The process of activation is parameterized with a simple Twomey
83!   activion scheme or with considering solution and curvature
84!   effects (Khvorostyanov and Curry ,2006).
85! - The saturation adjustment scheme is replaced by the parameterization
86!   of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
87! - All other microphysical processes of Seifert and Beheng are used.
88!   Additionally, in those processes the reduction of cloud number concentration
89!   is considered.
90!
91! 2233 2017-05-30 18:08:54Z suehring
92!
93! 2232 2017-05-30 17:47:52Z suehring
94! Adjustments to new topography and surface concept
95!
96! 2155 2017-02-21 09:57:40Z hoffmann
97! Bugfix in the calculation of microphysical quantities on ghost points.
98!
99! 2031 2016-10-21 15:11:58Z knoop
100! renamed variable rho to rho_ocean
101!
102! 2000 2016-08-20 18:09:15Z knoop
103! Forced header and separation lines into 80 columns
104!
105! 1850 2016-04-08 13:29:27Z maronga
106! Module renamed
107! Adapted for modularization of microphysics.
108!
109! 1845 2016-04-08 08:29:13Z raasch
110! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
111! point (instead of one point above surface)
112!
113! 1831 2016-04-07 13:15:51Z hoffmann
114! turbulence renamed collision_turbulence,
115! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
116! avaialble for microphysics_kessler.
117!
118! 1822 2016-04-07 07:49:42Z hoffmann
119! Unused variables removed.
120! Kessler scheme integrated.
121!
122! 1691 2015-10-26 16:17:44Z maronga
123! Added new routine calc_precipitation_amount. The routine now allows to account
124! for precipitation due to sedimenation of cloud (fog) droplets
125!
126! 1682 2015-10-07 23:56:08Z knoop
127! Code annotations made doxygen readable
128!
129! 1646 2015-09-02 16:00:10Z hoffmann
130! Bugfix: Wrong computation of d_mean.
131!
132! 1361 2014-04-16 15:17:48Z hoffmann
133! Bugfix in sedimentation_rain: Index corrected.
134! Vectorized version of adjust_cloud added.
135! Little reformatting of the code.
136!
137! 1353 2014-04-08 15:21:23Z heinze
138! REAL constants provided with KIND-attribute
139!
140! 1346 2014-03-27 13:18:20Z heinze
141! Bugfix: REAL constants provided with KIND-attribute especially in call of
142! intrinsic function like MAX, MIN, SIGN
143!
144! 1334 2014-03-25 12:21:40Z heinze
145! Bugfix: REAL constants provided with KIND-attribute
146!
147! 1322 2014-03-20 16:38:49Z raasch
148! REAL constants defined as wp-kind
149!
150! 1320 2014-03-20 08:40:49Z raasch
151! ONLY-attribute added to USE-statements,
152! kind-parameters added to all INTEGER and REAL declaration statements,
153! kinds are defined in new module kinds,
154! comment fields (!:) to be used for variable explanations added to
155! all variable declaration statements
156!
157! 1241 2013-10-30 11:36:58Z heinze
158! hyp and rho_ocean have to be calculated at each time step if data from external
159! file LSF_DATA are used
160!
161! 1115 2013-03-26 18:16:16Z hoffmann
162! microphyical tendencies are calculated in bcm_actions in an optimized
163! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
164!
165! 1106 2013-03-04 05:31:38Z raasch
166! small changes in code formatting
167!
168! 1092 2013-02-02 11:24:22Z raasch
169! unused variables removed
170! file put under GPL
171!
172! 1065 2012-11-22 17:42:36Z hoffmann
173! Sedimentation process implemented according to Stevens and Seifert (2008).
174! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
175! and Stevens, 2010).
176!
177! 1053 2012-11-13 17:11:03Z hoffmann
178! initial revision
179!
180! Description:
181! ------------
182!> Calculate bulk cloud microphysics.
183!------------------------------------------------------------------------------!
184 MODULE bulk_cloud_model_mod
185
186    USE arrays_3d,                                                             &
187#if defined (__nopointer)
188        ONLY:  ddzu, diss, dzu, dzw, hyp, hyrho,                               &
189               nc,                   nc_p, nr,                   nr_p,         &
190               precipitation_amount, prr, pt, d_exner, pt_init, q, ql,         &
191               qc,                   qc_p, qr,                   qr_p,         &
192               exner, zu, tnc_m, tnr_m, tqc_m, tqr_m
193#else
194        ONLY:  ddzu, diss, dzu, dzw, hyp, hyrho,                               &
195               nc, nc_1, nc_2, nc_3, nc_p, nr, nr_1, nr_2, nr_3, nr_p,         &
196               precipitation_amount, prr, pt, d_exner, pt_init, q, ql, ql_1,   &
197               qc, qc_1, qc_2, qc_3, qc_p, qr, qr_1, qr_2, qr_3, qr_p,         &
198               exner, zu, tnc_m, tnr_m, tqc_m, tqr_m
199#endif
200
201    USE averaging,                                                             &
202        ONLY:  nc_av, nr_av, prr_av, qc_av, ql_av, qr_av
203
204    USE basic_constants_and_equations_mod,                                     &
205        ONLY:  c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, molecular_weight_of_solute,&
206               molecular_weight_of_water, pi, rho_l, rho_s, r_d, r_v, vanthoff,&
207               exner_function, exner_function_invers, ideal_gas_law_rho,       &
208               ideal_gas_law_rho_pt, barometric_formula, rd_d_rv
209
210    USE control_parameters,                                                    &
211        ONLY:  dt_3d, dt_do2d_xy, intermediate_timestep_count,                 &
212               intermediate_timestep_count_max, large_scale_forcing,           &
213               lsf_surf, pt_surface, rho_surface, surface_pressure,            &
214               time_do2d_xy, message_string, initializing_actions
215
216    USE cpulog,                                                                &
217        ONLY:  cpu_log, log_point_s
218
219    USE grid_variables,                                                        &
220        ONLY:  dx, dy
221
222    USE indices,                                                               &
223        ONLY:  nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt,           &
224               wall_flags_0
225
226    USE  kinds
227
228    USE statistics,                                                            &
229        ONLY:  weight_pres, weight_substep
230
231    USE surface_mod,                                                           &
232        ONLY :  bc_h, get_topography_top_index_ji, surf_bulk_cloud_model,      &
233                surf_microphysics_morrison, surf_microphysics_seifert
234
235    IMPLICIT NONE
236
237    CHARACTER (LEN=20)   ::  aerosol_bulk = 'nacl'                        !< namelist parameter
238    CHARACTER (LEN=20)   ::  cloud_scheme = 'saturation_adjust'           !< namelist parameter
239
240    LOGICAL ::  aerosol_nacl =.TRUE.                             !< nacl aerosol for bulk scheme
241    LOGICAL ::  aerosol_c3h4o4 =.FALSE.                          !< malonic acid aerosol for bulk scheme
242    LOGICAL ::  aerosol_nh4no3 =.FALSE.                          !< malonic acid aerosol for bulk scheme
243
244    LOGICAL ::  bulk_cloud_model = .FALSE.                       !< namelist parameter
245
246    LOGICAL ::  cloud_water_sedimentation = .FALSE.       !< cloud water sedimentation
247    LOGICAL ::  curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
248    LOGICAL ::  limiter_sedimentation = .TRUE.            !< sedimentation limiter
249    LOGICAL ::  collision_turbulence = .FALSE.            !< turbulence effects
250    LOGICAL ::  ventilation_effect = .TRUE.               !< ventilation effect
251
252    LOGICAL ::  call_microphysics_at_all_substeps = .FALSE.      !< namelist parameter
253    LOGICAL ::  microphysics_sat_adjust = .FALSE.                !< use saturation adjust bulk scheme?
254    LOGICAL ::  microphysics_kessler = .FALSE.                   !< use kessler bulk scheme?
255    LOGICAL ::  microphysics_morrison = .FALSE.                  !< use 2-moment Morrison (add. prog. eq. for nc and qc)
256    LOGICAL ::  microphysics_seifert = .FALSE.                   !< use 2-moment Seifert and Beheng scheme
257    LOGICAL ::  precipitation = .FALSE.                          !< namelist parameter
258
259    REAL(wp) ::  precipitation_amount_interval = 9999999.9_wp    !< namelist parameter
260
261    REAL(wp) ::  a_1 = 8.69E-4_wp          !< coef. in turb. parametrization (cm-2 s3)
262    REAL(wp) ::  a_2 = -7.38E-5_wp         !< coef. in turb. parametrization (cm-2 s3)
263    REAL(wp) ::  a_3 = -1.40E-2_wp         !< coef. in turb. parametrization
264    REAL(wp) ::  a_term = 9.65_wp          !< coef. for terminal velocity (m s-1)
265    REAL(wp) ::  a_vent = 0.78_wp          !< coef. for ventilation effect
266    REAL(wp) ::  b_1 = 11.45E-6_wp         !< coef. in turb. parametrization (m)
267    REAL(wp) ::  b_2 = 9.68E-6_wp          !< coef. in turb. parametrization (m)
268    REAL(wp) ::  b_3 = 0.62_wp             !< coef. in turb. parametrization
269    REAL(wp) ::  b_term = 9.8_wp           !< coef. for terminal velocity (m s-1)
270    REAL(wp) ::  b_vent = 0.308_wp         !< coef. for ventilation effect
271    REAL(wp) ::  beta_cc = 3.09E-4_wp      !< coef. in turb. parametrization (cm-2 s3)
272    REAL(wp) ::  c_1 = 4.82E-6_wp          !< coef. in turb. parametrization (m)
273    REAL(wp) ::  c_2 = 4.8E-6_wp           !< coef. in turb. parametrization (m)
274    REAL(wp) ::  c_3 = 0.76_wp             !< coef. in turb. parametrization
275    REAL(wp) ::  c_const = 0.93_wp         !< const. in Taylor-microscale Reynolds number
276    REAL(wp) ::  c_evap = 0.7_wp           !< constant in evaporation
277    REAL(wp) ::  c_term = 600.0_wp         !< coef. for terminal velocity (m-1)
278    REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
279    REAL(wp) ::  eps_sb = 1.0E-10_wp       !< threshold in two-moments scheme
280    REAL(wp) ::  eps_mr = 0.0_wp           !< threshold for morrison scheme
281    REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
282    REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
283    REAL(wp) ::  k_rr = 7.12_wp            !< const. rain-rain kernel (m3 kg-1 s-1)
284    REAL(wp) ::  k_br = 1000.0_wp          !< const. in breakup parametrization (m-1)
285    REAL(wp) ::  k_st = 1.2E8_wp           !< const. in drizzle parametrization (m-1 s-1)
286    REAL(wp) ::  kin_vis_air = 1.4086E-5_wp  !< kin. viscosity of air (m2 s-1)
287    REAL(wp) ::  prec_time_const = 0.001_wp  !< coef. in Kessler scheme (s-1)
288    REAL(wp) ::  ql_crit = 0.0005_wp       !< coef. in Kessler scheme (kg kg-1)
289    REAL(wp) ::  schmidt_p_1d3=0.8921121_wp  !< Schmidt number**0.33333, 0.71**0.33333
290    REAL(wp) ::  sigma_gc = 1.3_wp         !< geometric standard deviation cloud droplets
291    REAL(wp) ::  thermal_conductivity_l = 2.43E-2_wp  !< therm. cond. air (J m-1 s-1 K-1)
292    REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
293    REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
294    REAL(wp) ::  xamin = 5.24E-19_wp       !< average aerosol mass (kg) (~ 0.05µm)
295    REAL(wp) ::  xcmin = 4.18E-15_wp       !< minimum cloud drop size (kg) (~ 1µm)
296    REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
297    REAL(wp) ::  xrmax = 5.0E-6_wp         !< maximum rain drop site (kg)
298
299    REAL(wp) ::  c_sedimentation = 2.0_wp        !< Courant number of sedimentation process
300    REAL(wp) ::  dpirho_l                        !< 6.0 / ( pi * rho_l )
301    REAL(wp) ::  dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius
302    REAL(wp) ::  dt_micro                        !< microphysics time step
303    REAL(wp) ::  sigma_bulk = 2.0_wp             !< width of aerosol spectrum
304    REAL(wp) ::  na_init = 100.0E6_wp            !< Total particle/aerosol concentration (cm-3)
305    REAL(wp) ::  nc_const = 70.0E6_wp            !< cloud droplet concentration
306    REAL(wp) ::  dt_precipitation = 100.0_wp     !< timestep precipitation (s)
307    REAL(wp) ::  sed_qc_const                    !< const. for sedimentation of cloud water
308    REAL(wp) ::  pirho_l                         !< pi * rho_l / 6.0;
309
310    REAL(wp) ::  e_s     !< saturation water vapor pressure
311    REAL(wp) ::  q_s     !< saturation mixing ratio
312    REAL(wp) ::  sat     !< supersaturation
313    REAL(wp) ::  t_l     !< actual temperature
314
315    SAVE
316
317    PRIVATE
318
319    PUBLIC bcm_parin, &
320           bcm_check_parameters, &
321           bcm_check_data_output, &
322           bcm_check_data_output_pr, &
323           bcm_header, &
324           bcm_init_arrays, &
325           bcm_init, &
326           bcm_3d_data_averaging, &
327           bcm_data_output_2d, &
328           bcm_data_output_3d, &
329           bcm_swap_timelevel, &
330           bcm_rrd_global, &
331           bcm_rrd_local, &
332           bcm_wrd_global, &
333           bcm_wrd_local, &
334           bcm_actions, &
335           calc_liquid_water_content
336
337    PUBLIC call_microphysics_at_all_substeps, &
338           cloud_water_sedimentation, &
339           bulk_cloud_model, &
340           cloud_scheme, &
341           collision_turbulence, &
342           dt_precipitation, &
343           microphysics_morrison, &
344           microphysics_sat_adjust, &
345           microphysics_seifert, &
346           na_init, &
347           nc_const, &
348           precipitation, &
349           sigma_gc
350
351
352    INTERFACE bcm_parin
353       MODULE PROCEDURE bcm_parin
354    END INTERFACE bcm_parin
355
356    INTERFACE bcm_check_parameters
357       MODULE PROCEDURE bcm_check_parameters
358    END INTERFACE bcm_check_parameters
359
360    INTERFACE bcm_check_data_output
361       MODULE PROCEDURE bcm_check_data_output
362    END INTERFACE bcm_check_data_output
363
364    INTERFACE bcm_check_data_output_pr
365       MODULE PROCEDURE bcm_check_data_output_pr
366    END INTERFACE bcm_check_data_output_pr
367
368    INTERFACE bcm_header
369       MODULE PROCEDURE bcm_header
370    END INTERFACE bcm_header
371
372    INTERFACE bcm_init_arrays
373       MODULE PROCEDURE bcm_init_arrays
374    END INTERFACE bcm_init_arrays
375
376    INTERFACE bcm_init
377       MODULE PROCEDURE bcm_init
378    END INTERFACE bcm_init
379
380    INTERFACE bcm_3d_data_averaging
381       MODULE PROCEDURE bcm_3d_data_averaging
382    END INTERFACE bcm_3d_data_averaging
383
384    INTERFACE bcm_data_output_2d
385       MODULE PROCEDURE bcm_data_output_2d
386    END INTERFACE bcm_data_output_2d
387
388    INTERFACE bcm_data_output_3d
389       MODULE PROCEDURE bcm_data_output_3d
390    END INTERFACE bcm_data_output_3d
391
392    INTERFACE bcm_swap_timelevel
393       MODULE PROCEDURE bcm_swap_timelevel
394    END INTERFACE bcm_swap_timelevel
395
396    INTERFACE bcm_rrd_global
397       MODULE PROCEDURE bcm_rrd_global
398    END INTERFACE bcm_rrd_global
399
400    INTERFACE bcm_rrd_local
401       MODULE PROCEDURE bcm_rrd_local
402    END INTERFACE bcm_rrd_local
403
404    INTERFACE bcm_wrd_global
405       MODULE PROCEDURE bcm_wrd_global
406    END INTERFACE bcm_wrd_global
407
408    INTERFACE bcm_wrd_local
409       MODULE PROCEDURE bcm_wrd_local
410    END INTERFACE bcm_wrd_local
411
412    INTERFACE bcm_actions
413       MODULE PROCEDURE bcm_actions
414       MODULE PROCEDURE bcm_actions_ij
415    END INTERFACE bcm_actions
416
417    INTERFACE adjust_cloud
418       MODULE PROCEDURE adjust_cloud
419       MODULE PROCEDURE adjust_cloud_ij
420    END INTERFACE adjust_cloud
421
422    INTERFACE activation
423       MODULE PROCEDURE activation
424       MODULE PROCEDURE activation_ij
425    END INTERFACE activation
426
427    INTERFACE condensation
428       MODULE PROCEDURE condensation
429       MODULE PROCEDURE condensation_ij
430    END INTERFACE condensation
431
432    INTERFACE autoconversion
433       MODULE PROCEDURE autoconversion
434       MODULE PROCEDURE autoconversion_ij
435    END INTERFACE autoconversion
436
437    INTERFACE autoconversion_kessler
438       MODULE PROCEDURE autoconversion_kessler
439       MODULE PROCEDURE autoconversion_kessler_ij
440    END INTERFACE autoconversion_kessler
441
442    INTERFACE accretion
443       MODULE PROCEDURE accretion
444       MODULE PROCEDURE accretion_ij
445    END INTERFACE accretion
446
447    INTERFACE selfcollection_breakup
448       MODULE PROCEDURE selfcollection_breakup
449       MODULE PROCEDURE selfcollection_breakup_ij
450    END INTERFACE selfcollection_breakup
451
452    INTERFACE evaporation_rain
453       MODULE PROCEDURE evaporation_rain
454       MODULE PROCEDURE evaporation_rain_ij
455    END INTERFACE evaporation_rain
456
457    INTERFACE sedimentation_cloud
458       MODULE PROCEDURE sedimentation_cloud
459       MODULE PROCEDURE sedimentation_cloud_ij
460    END INTERFACE sedimentation_cloud
461
462    INTERFACE sedimentation_rain
463       MODULE PROCEDURE sedimentation_rain
464       MODULE PROCEDURE sedimentation_rain_ij
465    END INTERFACE sedimentation_rain
466
467    INTERFACE calc_precipitation_amount
468       MODULE PROCEDURE calc_precipitation_amount
469       MODULE PROCEDURE calc_precipitation_amount_ij
470    END INTERFACE calc_precipitation_amount
471
472    INTERFACE supersaturation
473       MODULE PROCEDURE supersaturation
474    END INTERFACE supersaturation
475
476    INTERFACE calc_liquid_water_content
477       MODULE PROCEDURE calc_liquid_water_content
478    END INTERFACE calc_liquid_water_content
479
480 CONTAINS
481
482
483!------------------------------------------------------------------------------!
484! Description:
485! ------------
486!> Parin for &bulk_cloud_parameters for the bulk cloud module
487!------------------------------------------------------------------------------!
488    SUBROUTINE bcm_parin
489
490
491       IMPLICIT NONE
492
493       CHARACTER (LEN=80)  ::  line  !< dummy string that contains the current line of the parameter file
494
495       NAMELIST /bulk_cloud_parameters/  &
496          aerosol_bulk, &
497          c_sedimentation, &
498          call_microphysics_at_all_substeps, &
499          bulk_cloud_model, &
500          cloud_scheme, &
501          cloud_water_sedimentation, &
502          collision_turbulence, &
503          curvature_solution_effects_bulk, &
504          dry_aerosol_radius, &
505          limiter_sedimentation, &
506          na_init, &
507          nc_const, &
508          precipitation, &
509          precipitation_amount_interval, &
510          sigma_bulk, &
511          ventilation_effect
512
513       line = ' '
514!
515!--    Try to find bulk cloud module namelist
516       REWIND ( 11 )
517       line = ' '
518       DO   WHILE ( INDEX( line, '&bulk_cloud_parameters' ) == 0 )
519          READ ( 11, '(A)', END=10 )  line
520       ENDDO
521       BACKSPACE ( 11 )
522!
523!--    Read user-defined namelist
524       READ ( 11, bulk_cloud_parameters )
525!
526!--    Set flag that indicates that the bulk cloud module is switched on
527       !bulk_cloud_model = .TRUE.
528
52910     CONTINUE
530
531
532    END SUBROUTINE bcm_parin
533
534
535!------------------------------------------------------------------------------!
536! Description:
537! ------------
538!> Check parameters routine for bulk cloud module
539!------------------------------------------------------------------------------!
540    SUBROUTINE bcm_check_parameters
541
542
543       IMPLICIT NONE
544!
545!--    Check cloud scheme
546       IF ( cloud_scheme == 'saturation_adjust' )  THEN
547          microphysics_sat_adjust = .TRUE.
548          microphysics_seifert    = .FALSE.
549          microphysics_kessler    = .FALSE.
550          precipitation           = .FALSE.
551       ELSEIF ( cloud_scheme == 'seifert_beheng' )  THEN
552          microphysics_sat_adjust = .FALSE.
553          microphysics_seifert    = .TRUE.
554          microphysics_kessler    = .FALSE.
555          microphysics_morrison  = .FALSE.
556          precipitation           = .TRUE.
557       ELSEIF ( cloud_scheme == 'kessler' )  THEN
558          microphysics_sat_adjust = .FALSE.
559          microphysics_seifert    = .FALSE.
560          microphysics_kessler    = .TRUE.
561          microphysics_morrison   = .FALSE.
562          precipitation           = .TRUE.
563       ELSEIF ( cloud_scheme == 'morrison' )  THEN
564          microphysics_sat_adjust = .FALSE.
565          microphysics_seifert    = .TRUE.
566          microphysics_kessler    = .FALSE.
567          microphysics_morrison   = .TRUE.
568          precipitation           = .TRUE.
569       ELSE
570          message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
571                           TRIM( cloud_scheme ) // '"'
572          CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
573       ENDIF
574
575
576
577!
578!--    Set the default value for the integration interval of precipitation amount
579       IF ( microphysics_seifert  .OR.  microphysics_kessler )  THEN
580          IF ( precipitation_amount_interval == 9999999.9_wp )  THEN
581             precipitation_amount_interval = dt_do2d_xy
582          ELSE
583             IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
584                WRITE( message_string, * )  'precipitation_amount_interval = ',   &
585                    precipitation_amount_interval, ' must not be larger than ',   &
586                    'dt_do2d_xy = ', dt_do2d_xy
587                CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
588             ENDIF
589          ENDIF
590       ENDIF
591
592       ! TODO: find better sollution for circular dependency problem
593       surf_bulk_cloud_model = bulk_cloud_model
594       surf_microphysics_morrison = microphysics_morrison
595       surf_microphysics_seifert = microphysics_seifert
596
597!
598!--    Check aerosol
599       IF ( aerosol_bulk == 'nacl' )  THEN
600          aerosol_nacl   = .TRUE.
601          aerosol_c3h4o4 = .FALSE.
602          aerosol_nh4no3 = .FALSE.
603       ELSEIF ( aerosol_bulk == 'c3h4o4' )  THEN
604          aerosol_nacl   = .FALSE.
605          aerosol_c3h4o4 = .TRUE.
606          aerosol_nh4no3 = .FALSE.
607       ELSEIF ( aerosol_bulk == 'nh4no3' )  THEN
608          aerosol_nacl   = .FALSE.
609          aerosol_c3h4o4 = .FALSE.
610          aerosol_nh4no3 = .TRUE.
611       ELSE
612          message_string = 'unknown aerosol = "' // TRIM( aerosol_bulk ) // '"'
613          CALL message( 'check_parameters', 'PA0469', 1, 2, 0, 6, 0 )
614       ENDIF
615
616
617    END SUBROUTINE bcm_check_parameters
618
619!------------------------------------------------------------------------------!
620! Description:
621! ------------
622!> Check data output for bulk cloud module
623!------------------------------------------------------------------------------!
624    SUBROUTINE bcm_check_data_output( var, unit )
625
626       IMPLICIT NONE
627
628       CHARACTER (LEN=*) ::  unit  !<
629       CHARACTER (LEN=*) ::  var   !<
630
631       SELECT CASE ( TRIM( var ) )
632
633          CASE ( 'nc' )
634             IF ( .NOT.  microphysics_morrison )  THEN
635                message_string = 'output of "' // TRIM( var ) // '" ' //       &
636                                 'requires ' //                                &
637                                 'cloud_scheme = "morrison"'
638                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
639             ENDIF
640             unit = '1/m3'
641
642          CASE ( 'nr' )
643             IF ( .NOT.  microphysics_seifert )  THEN
644                message_string = 'output of "' // TRIM( var ) // '" ' //       &
645                                 'requires ' //                                &
646                                 'cloud_scheme = "seifert_beheng"'
647                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
648             ENDIF
649             unit = '1/m3'
650
651          CASE ( 'prr' )
652             IF ( microphysics_sat_adjust )  THEN
653                message_string = 'output of "' // TRIM( var ) // '" ' //       &
654                                 'is not available for ' //                    &
655                                 'cloud_scheme = "saturation_adjust"'
656                CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
657             ENDIF
658             unit = 'kg/kg m/s'
659
660          CASE ( 'qc' )
661             unit = 'kg/kg'
662
663          CASE ( 'qr' )
664             IF ( .NOT.  microphysics_seifert ) THEN
665                message_string = 'output of "' // TRIM( var ) // '" ' //       &
666                                 'requires ' //                                &
667                                 'cloud_scheme = "seifert_beheng"'
668                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
669             ENDIF
670             unit = 'kg/kg'
671
672          CASE ( 'pra*' )
673             IF ( .NOT. microphysics_kessler  .AND.                            &
674                  .NOT. microphysics_seifert )  THEN
675                message_string = 'output of "' // TRIM( var ) // '" ' //       &
676                                 'requires ' //                                &
677                                 'cloud_scheme = "kessler" or "seifert_beheng"'
678                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
679             ENDIF
680! TODO: find sollution (maybe connected to flow_statistics redesign?)
681!              IF ( j == 1 )  THEN
682!                 message_string = 'temporal averaging of precipitation ' //     &
683!                           'amount "' // TRIM( var ) // '" is not possible'
684!                 CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
685!              ENDIF
686             unit = 'mm'
687
688          CASE ( 'prr*' )
689             IF ( .NOT. microphysics_kessler  .AND.                            &
690                  .NOT. microphysics_seifert )  THEN
691                message_string = 'output of "' // TRIM( var ) // '"' //        &
692                         ' requires' //                                        &
693                         ' cloud_scheme = "kessler" or "seifert_beheng"'
694                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
695             ENDIF
696             unit = 'mm/s'
697
698          CASE DEFAULT
699             unit = 'illegal'
700
701       END SELECT
702
703
704    END SUBROUTINE bcm_check_data_output
705
706
707!------------------------------------------------------------------------------!
708! Description:
709! ------------
710!> Check data output of profiles for bulk cloud module
711!------------------------------------------------------------------------------!
712    SUBROUTINE bcm_check_data_output_pr( variable, var_count, unit, dopr_unit )
713
714       USE arrays_3d,                                                          &
715           ONLY: zu
716
717       USE control_parameters,                                                 &
718           ONLY: data_output_pr
719
720       USE profil_parameter,                                                   &
721           ONLY: dopr_index
722
723       USE statistics,                                                         &
724           ONLY: hom, statistic_regions, pr_palm
725
726       IMPLICIT NONE
727
728       CHARACTER (LEN=*) ::  unit      !<
729       CHARACTER (LEN=*) ::  variable  !<
730       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
731
732       INTEGER(iwp) ::  var_count      !<
733       INTEGER(iwp) ::  pr_index       !<
734
735       SELECT CASE ( TRIM( variable ) )
736
737! TODO: make index generic: pr_index = pr_palm+1
738
739          CASE ( 'nc' )
740             IF ( .NOT.  microphysics_morrison )  THEN
741                message_string = 'data_output_pr = ' //                        &
742                                 TRIM( data_output_pr(var_count) ) //          &
743                                 ' is not implemented for' //                  &
744                                 ' cloud_scheme /= morrison'
745                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
746             ENDIF
747             pr_index = 123
748             dopr_index(var_count) = pr_index
749             dopr_unit     = '1/m3'
750             unit = dopr_unit
751             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
752
753          CASE ( 'nr' )
754             IF ( .NOT.  microphysics_seifert )  THEN
755                message_string = 'data_output_pr = ' //                        &
756                                 TRIM( data_output_pr(var_count) ) //          &
757                                 ' is not implemented for' //                  &
758                                 ' cloud_scheme /= seifert_beheng'
759                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
760             ENDIF
761             pr_index = 73
762             dopr_index(var_count) = pr_index
763             dopr_unit     = '1/m3'
764             unit = dopr_unit
765             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
766
767          CASE ( 'prr' )
768             IF ( microphysics_sat_adjust )  THEN
769                message_string = 'data_output_pr = ' //                        &
770                                 TRIM( data_output_pr(var_count) ) //          &
771                                 ' is not available for' //                    &
772                                 ' cloud_scheme = saturation_adjust'
773                CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
774             ENDIF
775             pr_index = 76
776             dopr_index(var_count) = pr_index
777             dopr_unit     = 'kg/kg m/s'
778             unit = dopr_unit
779             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
780          CASE ( 'qc' )
781             pr_index = 75
782             dopr_index(var_count) = pr_index
783             dopr_unit     = 'kg/kg'
784             unit = dopr_unit
785             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
786
787          CASE ( 'qr' )
788             IF ( .NOT.  microphysics_seifert )  THEN
789                message_string = 'data_output_pr = ' //                        &
790                                 TRIM( data_output_pr(var_count) ) //          &
791                                 ' is not implemented for' //                  &
792                                 ' cloud_scheme /= seifert_beheng'
793                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
794             ENDIF
795             pr_index = 74
796             dopr_index(var_count) = pr_index
797             dopr_unit     = 'kg/kg'
798             unit = dopr_unit
799             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
800
801          CASE DEFAULT
802             unit = 'illegal'
803
804       END SELECT
805
806    END SUBROUTINE bcm_check_data_output_pr
807
808
809!------------------------------------------------------------------------------!
810! Description:
811! ------------
812!> Allocate bulk cloud module arrays and define pointers
813!------------------------------------------------------------------------------!
814    SUBROUTINE bcm_init_arrays
815
816       USE indices,                                                            &
817           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
818
819
820       IMPLICIT NONE
821
822       INTEGER(iwp) ::  i !<
823       INTEGER(iwp) ::  j !<
824       INTEGER(iwp) ::  k !<
825!
826!--          Liquid water content
827#if defined( __nopointer )
828       ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
829#else
830       ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
831#endif
832
833!
834!--    3D-cloud water content
835       IF ( .NOT. microphysics_morrison )  THEN
836#if defined( __nopointer )
837          ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
838#else
839          ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
840#endif
841       ENDIF
842!
843!--    Precipitation amount and rate (only needed if output is switched)
844       ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg) )
845
846!
847!--    3d-precipitation rate
848       ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
849
850       IF ( microphysics_morrison )  THEN
851!
852!--       3D-cloud drop water content, cloud drop concentration arrays
853#if defined( __nopointer )
854          ALLOCATE( nc(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                      &
855                    nc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
856                    qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                      &
857                    qc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
858                    tnc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                   &
859                    tqc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
860#else
861          ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
862                    nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
863                    nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
864                    qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
865                    qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
866                    qc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
867#endif
868       ENDIF
869
870       IF ( microphysics_seifert )  THEN
871!
872!--       3D-rain water content, rain drop concentration arrays
873#if defined( __nopointer )
874          ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                      &
875                    nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
876                    qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                      &
877                    qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
878                    tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                   &
879                    tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
880#else
881          ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
882                    nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
883                    nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
884                    qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
885                    qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
886                    qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
887#endif
888       ENDIF
889
890#if ! defined( __nopointer )
891!
892!--    Initial assignment of the pointers
893       ql => ql_1
894       IF ( .NOT. microphysics_morrison )  THEN
895          qc => qc_1
896       ENDIF
897       IF ( microphysics_morrison )  THEN
898          qc => qc_1;  qc_p  => qc_2;  tqc_m  => qc_3
899          nc => nc_1;  nc_p  => nc_2;  tnc_m  => nc_3
900       ENDIF
901       IF ( microphysics_seifert )  THEN
902          qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
903          nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
904       ENDIF
905#endif
906
907
908    END SUBROUTINE bcm_init_arrays
909
910
911!------------------------------------------------------------------------------!
912! Description:
913! ------------
914!> Initialization of the bulk cloud module
915!------------------------------------------------------------------------------!
916    SUBROUTINE bcm_init !( dots_label, dots_unit, dots_num, dots_max )
917
918       IMPLICIT NONE
919
920       INTEGER(iwp) ::  i !<
921       INTEGER(iwp) ::  j !<
922       INTEGER(iwp) ::  k !<
923
924!        INTEGER(iwp) ::  dots_num
925!        INTEGER(iwp) ::  dots_max
926!        CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit
927!        CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label
928
929       CALL location_message( 'initializing bulk cloud module', .FALSE. )
930
931       IF ( bulk_cloud_model )  THEN
932          IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
933!
934!--          Initialize the remaining quantities
935             IF ( microphysics_morrison )  THEN
936                DO  i = nxlg, nxrg
937                   DO  j = nysg, nyng
938                      qc(:,j,i) = 0.0_wp
939                      nc(:,j,i) = 0.0_wp
940                   ENDDO
941                ENDDO
942             ENDIF
943
944             IF ( microphysics_seifert )  THEN
945                DO  i = nxlg, nxrg
946                   DO  j = nysg, nyng
947                      qr(:,j,i) = 0.0_wp
948                      nr(:,j,i) = 0.0_wp
949                   ENDDO
950                ENDDO
951             ENDIF
952!
953!--          Liquid water content and precipitation amount
954!--          are zero at beginning of the simulation
955             ql = 0.0_wp
956             qc = 0.0_wp
957             precipitation_amount = 0.0_wp
958             prr = 0.0_wp
959!
960!--          Initialize old and new time levels.
961             IF ( microphysics_morrison )  THEN
962                tqc_m = 0.0_wp
963                tnc_m = 0.0_wp
964                qc_p  = qc
965                nc_p  = nc
966             ENDIF
967             IF ( microphysics_seifert )  THEN
968                tqr_m = 0.0_wp
969                tnr_m = 0.0_wp
970                qr_p  = qr
971                nr_p  = nr
972             ENDIF
973          ENDIF ! Only if not read_restart_data
974!
975!--       constant for the sedimentation of cloud water (2-moment cloud physics)
976          sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l )                &
977                             )**( 2.0_wp / 3.0_wp ) *                             &
978                         EXP( 5.0_wp * LOG( sigma_gc )**2 )
979
980!
981!--       Calculate timestep according to precipitation
982          IF ( microphysics_seifert )  THEN
983             dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /      &
984                                w_precipitation
985          ENDIF
986
987!
988!--       Set constants for certain aerosol type
989          IF ( microphysics_morrison )  THEN
990             IF ( aerosol_nacl ) THEN
991                molecular_weight_of_solute = 0.05844_wp
992                rho_s                      = 2165.0_wp
993                vanthoff                   = 2.0_wp
994             ELSEIF ( aerosol_c3h4o4 ) THEN
995                molecular_weight_of_solute = 0.10406_wp
996                rho_s                      = 1600.0_wp
997                vanthoff                   = 1.37_wp
998             ELSEIF ( aerosol_nh4no3 ) THEN
999                molecular_weight_of_solute = 0.08004_wp
1000                rho_s                      = 1720.0_wp
1001                vanthoff                   = 2.31_wp
1002             ENDIF
1003          ENDIF
1004
1005!
1006!--       Pre-calculate frequently calculated fractions of pi and rho_l
1007          pirho_l  = pi * rho_l / 6.0_wp
1008          dpirho_l = 1.0_wp / pirho_l
1009
1010          CALL location_message( 'finished', .TRUE. )
1011
1012       ELSE
1013
1014          CALL location_message( 'skipped', .TRUE. )
1015
1016       ENDIF
1017
1018    END SUBROUTINE bcm_init
1019
1020
1021!------------------------------------------------------------------------------!
1022! Description:
1023! ------------
1024!> Swapping of timelevels
1025!------------------------------------------------------------------------------!
1026    SUBROUTINE bcm_swap_timelevel ( mod_count )
1027
1028       IMPLICIT NONE
1029
1030       INTEGER, INTENT(IN) :: mod_count
1031
1032       IF ( bulk_cloud_model )  THEN
1033
1034#if defined( __nopointer )
1035          IF ( microphysics_morrison )  THEN
1036             qc = qc_p
1037             nc = nc_p
1038          ENDIF
1039          IF ( microphysics_seifert )  THEN
1040             qr = qr_p
1041             nr = nr_p
1042          ENDIF
1043#else
1044          SELECT CASE ( mod_count )
1045
1046             CASE ( 0 )
1047
1048                IF ( microphysics_morrison )  THEN
1049                   qc => qc_1;    qc_p => qc_2
1050                   nc => nc_1;    nc_p => nc_2
1051                ENDIF
1052                IF ( microphysics_seifert )  THEN
1053                   qr => qr_1;    qr_p => qr_2
1054                   nr => nr_1;    nr_p => nr_2
1055                ENDIF
1056
1057             CASE ( 1 )
1058
1059                IF ( microphysics_morrison )  THEN
1060                   qc => qc_2;    qc_p => qc_1
1061                   nc => nc_2;    nc_p => nc_1
1062                ENDIF
1063                IF ( microphysics_seifert )  THEN
1064                   qr => qr_2;    qr_p => qr_1
1065                   nr => nr_2;    nr_p => nr_1
1066                ENDIF
1067
1068          END SELECT
1069#endif
1070
1071       ENDIF
1072
1073    END SUBROUTINE bcm_swap_timelevel
1074
1075
1076!------------------------------------------------------------------------------!
1077! Description:
1078! ------------
1079!> Header output for bulk cloud module
1080!------------------------------------------------------------------------------!
1081    SUBROUTINE bcm_header ( io )
1082
1083
1084       IMPLICIT NONE
1085
1086       INTEGER(iwp), INTENT(IN) ::  io  !< Unit of the output file
1087
1088!
1089!--    Write bulk cloud module header
1090       WRITE ( io, 1 )
1091
1092       WRITE ( io, 2 )
1093       WRITE ( io, 3 )
1094
1095       IF ( microphysics_kessler )  THEN
1096          WRITE ( io, 4 ) 'Kessler-Scheme'
1097       ENDIF
1098
1099       IF ( microphysics_seifert )  THEN
1100          WRITE ( io, 4 ) 'Seifert-Beheng-Scheme'
1101          IF ( cloud_water_sedimentation )  WRITE ( io, 5 )
1102          IF ( collision_turbulence )  WRITE ( io, 6 )
1103          IF ( ventilation_effect )  WRITE ( io, 7 )
1104          IF ( limiter_sedimentation )  WRITE ( io, 8 )
1105       ENDIF
1106
1107       WRITE ( io, 20 )
1108       WRITE ( io, 21 ) surface_pressure
1109       WRITE ( io, 22 ) r_d
1110       WRITE ( io, 23 ) rho_surface
1111       WRITE ( io, 24 ) c_p
1112       WRITE ( io, 25 ) l_v
1113
1114       IF ( microphysics_seifert )  THEN
1115          WRITE ( io, 26 ) 1.0E-6_wp * nc_const
1116          WRITE ( io, 27 ) c_sedimentation
1117       ENDIF
1118
1119
11201   FORMAT ( //' Bulk cloud module information:'/ &
1121               ' ------------------------------------------'/ )
11222   FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/ &
1123             '    total water content is used.' )
11243   FORMAT ( '--> Condensation is parameterized via 0% - or 100% scheme.' )
11254   FORMAT ( '--> Precipitation parameterization via ', A )
1126
11275   FORMAT ( '--> Cloud water sedimentation parameterization via Stokes law' )
11286   FORMAT ( '--> Turbulence effects on precipitation process' )
11297   FORMAT ( '--> Ventilation effects on evaporation of rain drops' )
11308   FORMAT ( '--> Slope limiter used for sedimentation process' )
1131
113220  FORMAT ( '--> Essential parameters:' )
113321  FORMAT ( '       Surface pressure             :   p_0   = ', F7.2, ' hPa')
113422  FORMAT ( '       Gas constant                 :   R     = ', F5.1, ' J/(kg K)')
113523  FORMAT ( '       Density of air               :   rho_0 = ', F6.3, ' kg/m**3')
113624  FORMAT ( '       Specific heat cap.           :   c_p   = ', F6.1, ' J/(kg K)')
113725  FORMAT ( '       Vapourization heat           :   L_v   = ', E9.2, ' J/kg')
113826  FORMAT ( '       Droplet density              :   N_c   = ', F6.1, ' 1/cm**3' )
113927  FORMAT ( '       Sedimentation Courant number :   C_s   = ', F4.1 )
1140
1141
1142    END SUBROUTINE bcm_header
1143
1144
1145!------------------------------------------------------------------------------!
1146!
1147! Description:
1148! ------------
1149!> Subroutine for averaging 3D data
1150!------------------------------------------------------------------------------!
1151    SUBROUTINE bcm_3d_data_averaging( mode, variable )
1152
1153       USE control_parameters,                                                 &
1154           ONLY:  average_count_3d
1155
1156       IMPLICIT NONE
1157
1158       CHARACTER (LEN=*) ::  mode     !<
1159       CHARACTER (LEN=*) ::  variable !<
1160
1161       INTEGER(iwp) ::  i       !< local index
1162       INTEGER(iwp) ::  j       !< local index
1163       INTEGER(iwp) ::  k       !< local index
1164
1165       IF ( mode == 'allocate' )  THEN
1166
1167          SELECT CASE ( TRIM( variable ) )
1168
1169             CASE ( 'nc' )
1170                IF ( .NOT. ALLOCATED( nc_av ) )  THEN
1171                   ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1172                ENDIF
1173                nc_av = 0.0_wp
1174
1175             CASE ( 'nr' )
1176                IF ( .NOT. ALLOCATED( nr_av ) )  THEN
1177                   ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1178                ENDIF
1179                nr_av = 0.0_wp
1180
1181             CASE ( 'prr' )
1182                IF ( .NOT. ALLOCATED( prr_av ) )  THEN
1183                   ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1184                ENDIF
1185                prr_av = 0.0_wp
1186
1187             CASE ( 'qc' )
1188                IF ( .NOT. ALLOCATED( qc_av ) )  THEN
1189                   ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1190                ENDIF
1191                qc_av = 0.0_wp
1192
1193             CASE ( 'ql' )
1194                IF ( .NOT. ALLOCATED( ql_av ) )  THEN
1195                   ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1196                ENDIF
1197                ql_av = 0.0_wp
1198
1199             CASE ( 'qr' )
1200                IF ( .NOT. ALLOCATED( qr_av ) )  THEN
1201                   ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1202                ENDIF
1203                qr_av = 0.0_wp
1204
1205             CASE DEFAULT
1206                CONTINUE
1207
1208          END SELECT
1209
1210       ELSEIF ( mode == 'sum' )  THEN
1211
1212          SELECT CASE ( TRIM( variable ) )
1213
1214             CASE ( 'nc' )
1215                IF ( ALLOCATED( nc_av ) ) THEN
1216                   DO  i = nxlg, nxrg
1217                      DO  j = nysg, nyng
1218                         DO  k = nzb, nzt+1
1219                            nc_av(k,j,i) = nc_av(k,j,i) + nc(k,j,i)
1220                         ENDDO
1221                      ENDDO
1222                   ENDDO
1223                ENDIF
1224
1225             CASE ( 'nr' )
1226                IF ( ALLOCATED( nr_av ) ) THEN
1227                   DO  i = nxlg, nxrg
1228                      DO  j = nysg, nyng
1229                         DO  k = nzb, nzt+1
1230                            nr_av(k,j,i) = nr_av(k,j,i) + nr(k,j,i)
1231                         ENDDO
1232                      ENDDO
1233                   ENDDO
1234                ENDIF
1235
1236             CASE ( 'prr' )
1237                IF ( ALLOCATED( prr_av ) ) THEN
1238                   DO  i = nxlg, nxrg
1239                      DO  j = nysg, nyng
1240                         DO  k = nzb, nzt+1
1241                            prr_av(k,j,i) = prr_av(k,j,i) + prr(k,j,i)
1242                         ENDDO
1243                      ENDDO
1244                   ENDDO
1245                ENDIF
1246
1247             CASE ( 'qc' )
1248                IF ( ALLOCATED( qc_av ) ) THEN
1249                   DO  i = nxlg, nxrg
1250                      DO  j = nysg, nyng
1251                         DO  k = nzb, nzt+1
1252                            qc_av(k,j,i) = qc_av(k,j,i) + qc(k,j,i)
1253                         ENDDO
1254                      ENDDO
1255                   ENDDO
1256                ENDIF
1257
1258             CASE ( 'ql' )
1259                IF ( ALLOCATED( ql_av ) ) THEN
1260                   DO  i = nxlg, nxrg
1261                      DO  j = nysg, nyng
1262                         DO  k = nzb, nzt+1
1263                            ql_av(k,j,i) = ql_av(k,j,i) + ql(k,j,i)
1264                         ENDDO
1265                      ENDDO
1266                   ENDDO
1267                ENDIF
1268
1269             CASE ( 'qr' )
1270                IF ( ALLOCATED( qr_av ) ) THEN
1271                   DO  i = nxlg, nxrg
1272                      DO  j = nysg, nyng
1273                         DO  k = nzb, nzt+1
1274                            qr_av(k,j,i) = qr_av(k,j,i) + qr(k,j,i)
1275                         ENDDO
1276                      ENDDO
1277                   ENDDO
1278                ENDIF
1279
1280             CASE DEFAULT
1281                CONTINUE
1282
1283          END SELECT
1284
1285       ELSEIF ( mode == 'average' )  THEN
1286
1287          SELECT CASE ( TRIM( variable ) )
1288
1289             CASE ( 'nc' )
1290                IF ( ALLOCATED( nc_av ) ) THEN
1291                   DO  i = nxlg, nxrg
1292                      DO  j = nysg, nyng
1293                         DO  k = nzb, nzt+1
1294                            nc_av(k,j,i) = nc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1295                         ENDDO
1296                      ENDDO
1297                   ENDDO
1298                ENDIF
1299
1300             CASE ( 'nr' )
1301                IF ( ALLOCATED( nr_av ) ) THEN
1302                   DO  i = nxlg, nxrg
1303                      DO  j = nysg, nyng
1304                         DO  k = nzb, nzt+1
1305                            nr_av(k,j,i) = nr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1306                         ENDDO
1307                      ENDDO
1308                   ENDDO
1309                ENDIF
1310
1311             CASE ( 'prr' )
1312                IF ( ALLOCATED( prr_av ) ) THEN
1313                   DO  i = nxlg, nxrg
1314                      DO  j = nysg, nyng
1315                         DO  k = nzb, nzt+1
1316                            prr_av(k,j,i) = prr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1317                         ENDDO
1318                      ENDDO
1319                   ENDDO
1320                ENDIF
1321
1322             CASE ( 'qc' )
1323                IF ( ALLOCATED( qc_av ) ) THEN
1324                   DO  i = nxlg, nxrg
1325                      DO  j = nysg, nyng
1326                         DO  k = nzb, nzt+1
1327                            qc_av(k,j,i) = qc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1328                         ENDDO
1329                      ENDDO
1330                   ENDDO
1331                ENDIF
1332
1333             CASE ( 'ql' )
1334                IF ( ALLOCATED( ql_av ) ) THEN
1335                   DO  i = nxlg, nxrg
1336                      DO  j = nysg, nyng
1337                         DO  k = nzb, nzt+1
1338                            ql_av(k,j,i) = ql_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1339                         ENDDO
1340                      ENDDO
1341                   ENDDO
1342                ENDIF
1343
1344             CASE ( 'qr' )
1345                IF ( ALLOCATED( qr_av ) ) THEN
1346                   DO  i = nxlg, nxrg
1347                      DO  j = nysg, nyng
1348                         DO  k = nzb, nzt+1
1349                            qr_av(k,j,i) = qr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1350                         ENDDO
1351                      ENDDO
1352                   ENDDO
1353                ENDIF
1354
1355             CASE DEFAULT
1356                CONTINUE
1357
1358          END SELECT
1359
1360       ENDIF
1361
1362    END SUBROUTINE bcm_3d_data_averaging
1363
1364
1365!------------------------------------------------------------------------------!
1366! Description:
1367! ------------
1368!> Define 2D output variables.
1369!------------------------------------------------------------------------------!
1370 SUBROUTINE bcm_data_output_2d( av, variable, found, grid, mode, local_pf,     &
1371                                two_d, nzb_do, nzt_do )
1372
1373
1374    IMPLICIT NONE
1375
1376    CHARACTER (LEN=*), INTENT(INOUT) ::  grid       !< name of vertical grid
1377    CHARACTER (LEN=*), INTENT(IN) ::  mode       !< either 'xy', 'xz' or 'yz'
1378    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< name of variable
1379
1380    INTEGER(iwp), INTENT(IN) ::  av        !< flag for (non-)average output
1381    INTEGER(iwp), INTENT(IN) ::  nzb_do    !< vertical output index (bottom)
1382    INTEGER(iwp), INTENT(IN) ::  nzt_do    !< vertical output index (top)
1383
1384    INTEGER(iwp) ::  flag_nr   !< number of masking flag
1385
1386    INTEGER(iwp) ::  i         !< loop index along x-direction
1387    INTEGER(iwp) ::  j         !< loop index along y-direction
1388    INTEGER(iwp) ::  k         !< loop index along z-direction
1389
1390    LOGICAL, INTENT(INOUT) ::  found   !< flag if output variable is found
1391    LOGICAL, INTENT(INOUT) ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
1392    LOGICAL ::  resorted  !< flag if output is already resorted
1393
1394    REAL(wp), PARAMETER ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
1395
1396    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) ::  local_pf !< local
1397       !< array to which output data is resorted to
1398
1399    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
1400
1401    found = .TRUE.
1402    resorted = .FALSE.
1403!
1404!-- Set masking flag for topography for not resorted arrays
1405    flag_nr = 0    ! 0 = scalar, 1 = u, 2 = v, 3 = w
1406
1407    SELECT CASE ( TRIM( variable ) )
1408
1409       CASE ( 'nc_xy', 'nc_xz', 'nc_yz' )
1410          IF ( av == 0 )  THEN
1411             to_be_resorted => nc
1412          ELSE
1413             IF ( .NOT. ALLOCATED( nc_av ) ) THEN
1414                ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1415                nc_av = REAL( fill_value, KIND = wp )
1416             ENDIF
1417             to_be_resorted => nc_av
1418          ENDIF
1419          IF ( mode == 'xy' ) grid = 'zu'
1420
1421       CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
1422          IF ( av == 0 )  THEN
1423             to_be_resorted => nr
1424          ELSE
1425             IF ( .NOT. ALLOCATED( nr_av ) ) THEN
1426                ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1427                nr_av = REAL( fill_value, KIND = wp )
1428             ENDIF
1429             to_be_resorted => nr_av
1430          ENDIF
1431          IF ( mode == 'xy' ) grid = 'zu'
1432
1433       CASE ( 'pra*_xy' )        ! 2d-array / integral quantity => no av
1434!                CALL exchange_horiz_2d( precipitation_amount )
1435          DO  i = nxl, nxr
1436             DO  j = nys, nyn
1437                local_pf(i,j,nzb+1) =  precipitation_amount(j,i)
1438             ENDDO
1439          ENDDO
1440          precipitation_amount = 0.0_wp   ! reset for next integ. interval
1441          resorted = .TRUE.
1442          two_d = .TRUE.
1443          IF ( mode == 'xy' ) grid = 'zu1'
1444
1445       CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
1446          IF ( av == 0 )  THEN
1447!                   CALL exchange_horiz( prr, nbgp )
1448             DO  i = nxl, nxr
1449                DO  j = nys, nyn
1450                   DO  k = nzb, nzt+1
1451                      local_pf(i,j,k) = prr(k,j,i) * hyrho(nzb+1)
1452                   ENDDO
1453                ENDDO
1454             ENDDO
1455          ELSE
1456             IF ( .NOT. ALLOCATED( prr_av ) ) THEN
1457                ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1458                prr_av = REAL( fill_value, KIND = wp )
1459             ENDIF
1460!                   CALL exchange_horiz( prr_av, nbgp )
1461             DO  i = nxl, nxr
1462                DO  j = nys, nyn
1463                   DO  k = nzb, nzt+1
1464                      local_pf(i,j,k) = prr_av(k,j,i) * hyrho(nzb+1)
1465                   ENDDO
1466                ENDDO
1467             ENDDO
1468          ENDIF
1469          resorted = .TRUE.
1470          IF ( mode == 'xy' ) grid = 'zu'
1471
1472       CASE ( 'qc_xy', 'qc_xz', 'qc_yz' )
1473          IF ( av == 0 )  THEN
1474             to_be_resorted => qc
1475          ELSE
1476             IF ( .NOT. ALLOCATED( qc_av ) ) THEN
1477                ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1478                qc_av = REAL( fill_value, KIND = wp )
1479             ENDIF
1480             to_be_resorted => qc_av
1481          ENDIF
1482          IF ( mode == 'xy' ) grid = 'zu'
1483
1484       CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
1485          IF ( av == 0 )  THEN
1486             to_be_resorted => ql
1487          ELSE
1488             IF ( .NOT. ALLOCATED( ql_av ) ) THEN
1489                ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1490                ql_av = REAL( fill_value, KIND = wp )
1491             ENDIF
1492             to_be_resorted => ql_av
1493          ENDIF
1494          IF ( mode == 'xy' ) grid = 'zu'
1495
1496       CASE ( 'qr_xy', 'qr_xz', 'qr_yz' )
1497          IF ( av == 0 )  THEN
1498             to_be_resorted => qr
1499          ELSE
1500             IF ( .NOT. ALLOCATED( qr_av ) ) THEN
1501                ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1502                qr_av = REAL( fill_value, KIND = wp )
1503             ENDIF
1504             to_be_resorted => qr_av
1505          ENDIF
1506          IF ( mode == 'xy' ) grid = 'zu'
1507
1508       CASE DEFAULT
1509          found = .FALSE.
1510          grid  = 'none'
1511
1512    END SELECT
1513
1514    IF ( found .AND. .NOT. resorted )  THEN
1515       DO  i = nxl, nxr
1516          DO  j = nys, nyn
1517             DO  k = nzb_do, nzt_do
1518                local_pf(i,j,k) = MERGE( &
1519                                     to_be_resorted(k,j,i),                    &
1520                                     REAL( fill_value, KIND = wp ),            &
1521                                     BTEST( wall_flags_0(k,j,i), flag_nr )     &
1522                                  )
1523             ENDDO
1524          ENDDO
1525       ENDDO
1526    ENDIF
1527
1528 END SUBROUTINE bcm_data_output_2d
1529
1530
1531!------------------------------------------------------------------------------!
1532! Description:
1533! ------------
1534!> Define 3D output variables.
1535!------------------------------------------------------------------------------!
1536 SUBROUTINE bcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
1537
1538
1539    IMPLICIT NONE
1540
1541    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< name of variable
1542
1543    INTEGER(iwp), INTENT(IN) ::  av        !< flag for (non-)average output
1544    INTEGER(iwp), INTENT(IN) ::  nzb_do    !< lower limit of the data output (usually 0)
1545    INTEGER(iwp), INTENT(IN) ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
1546
1547    INTEGER(iwp) ::  flag_nr   !< number of masking flag
1548
1549    INTEGER(iwp) ::  i         !< loop index along x-direction
1550    INTEGER(iwp) ::  j         !< loop index along y-direction
1551    INTEGER(iwp) ::  k         !< loop index along z-direction
1552
1553    LOGICAL, INTENT(INOUT) ::  found     !< flag if output variable is found
1554    LOGICAL ::  resorted  !< flag if output is already resorted
1555
1556    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
1557
1558    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) ::  local_pf   !< local
1559       !< array to which output data is resorted to
1560
1561    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
1562
1563    found = .TRUE.
1564    resorted = .FALSE.
1565!
1566!-- Set masking flag for topography for not resorted arrays
1567    flag_nr = 0    ! 0 = scalar, 1 = u, 2 = v, 3 = w
1568
1569    SELECT CASE ( TRIM( variable ) )
1570
1571       CASE ( 'nc' )
1572          IF ( av == 0 )  THEN
1573             to_be_resorted => nc
1574          ELSE
1575             IF ( .NOT. ALLOCATED( nc_av ) ) THEN
1576                ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1577                nc_av = REAL( fill_value, KIND = wp )
1578             ENDIF
1579             to_be_resorted => nc_av
1580          ENDIF
1581
1582       CASE ( 'nr' )
1583          IF ( av == 0 )  THEN
1584             to_be_resorted => nr
1585          ELSE
1586             IF ( .NOT. ALLOCATED( nr_av ) ) THEN
1587                ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1588                nr_av = REAL( fill_value, KIND = wp )
1589             ENDIF
1590             to_be_resorted => nr_av
1591          ENDIF
1592
1593       CASE ( 'prr' )
1594          IF ( av == 0 )  THEN
1595             DO  i = nxl, nxr
1596                DO  j = nys, nyn
1597                   DO  k = nzb_do, nzt_do
1598                      local_pf(i,j,k) = prr(k,j,i)
1599                   ENDDO
1600                ENDDO
1601             ENDDO
1602          ELSE
1603             IF ( .NOT. ALLOCATED( prr_av ) ) THEN
1604                ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1605                prr_av = REAL( fill_value, KIND = wp )
1606             ENDIF
1607             DO  i = nxl, nxr
1608                DO  j = nys, nyn
1609                   DO  k = nzb_do, nzt_do
1610                      local_pf(i,j,k) = prr_av(k,j,i)
1611                   ENDDO
1612                ENDDO
1613             ENDDO
1614          ENDIF
1615          resorted = .TRUE.
1616
1617       CASE ( 'qc' )
1618          IF ( av == 0 )  THEN
1619             to_be_resorted => qc
1620          ELSE
1621             IF ( .NOT. ALLOCATED( qc_av ) ) THEN
1622                ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1623                qc_av = REAL( fill_value, KIND = wp )
1624             ENDIF
1625             to_be_resorted => qc_av
1626          ENDIF
1627
1628       CASE ( 'ql' )
1629          IF ( av == 0 )  THEN
1630             to_be_resorted => ql
1631          ELSE
1632             IF ( .NOT. ALLOCATED( ql_av ) ) THEN
1633                ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1634                ql_av = REAL( fill_value, KIND = wp )
1635             ENDIF
1636             to_be_resorted => ql_av
1637          ENDIF
1638
1639       CASE ( 'qr' )
1640          IF ( av == 0 )  THEN
1641             to_be_resorted => qr
1642          ELSE
1643             IF ( .NOT. ALLOCATED( qr_av ) ) THEN
1644                ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1645                qr_av = REAL( fill_value, KIND = wp )
1646             ENDIF
1647             to_be_resorted => qr_av
1648          ENDIF
1649
1650       CASE DEFAULT
1651          found = .FALSE.
1652
1653    END SELECT
1654
1655
1656    IF ( found .AND. .NOT. resorted )  THEN
1657       DO  i = nxl, nxr
1658          DO  j = nys, nyn
1659             DO  k = nzb_do, nzt_do
1660                local_pf(i,j,k) = MERGE(                                       &
1661                                     to_be_resorted(k,j,i),                    &
1662                                     REAL( fill_value, KIND = wp ),            &
1663                                     BTEST( wall_flags_0(k,j,i), flag_nr )     &
1664                                  )
1665             ENDDO
1666          ENDDO
1667       ENDDO
1668    ENDIF
1669
1670 END SUBROUTINE bcm_data_output_3d
1671
1672
1673!------------------------------------------------------------------------------!
1674! Description:
1675! ------------
1676!> This routine reads the respective restart data for the bulk cloud module.
1677!------------------------------------------------------------------------------!
1678    SUBROUTINE bcm_rrd_global( found )
1679
1680
1681       USE control_parameters,                                                 &
1682           ONLY: length, restart_string
1683
1684
1685       IMPLICIT NONE
1686
1687       LOGICAL, INTENT(OUT)  ::  found
1688
1689
1690       found = .TRUE.
1691
1692       SELECT CASE ( restart_string(1:length) )
1693
1694          CASE ( 'c_sedimentation' )
1695             READ ( 13 )  c_sedimentation
1696
1697          CASE ( 'bulk_cloud_model' )
1698             READ ( 13 )  bulk_cloud_model
1699
1700          CASE ( 'cloud_scheme' )
1701             READ ( 13 )  cloud_scheme
1702
1703          CASE ( 'cloud_water_sedimentation' )
1704             READ ( 13 )  cloud_water_sedimentation
1705
1706          CASE ( 'collision_turbulence' )
1707             READ ( 13 )  collision_turbulence
1708
1709          CASE ( 'limiter_sedimentation' )
1710             READ ( 13 )  limiter_sedimentation
1711
1712          CASE ( 'nc_const' )
1713             READ ( 13 )  nc_const
1714
1715          CASE ( 'precipitation' )
1716             READ ( 13 ) precipitation
1717
1718          CASE ( 'ventilation_effect' )
1719             READ ( 13 )  ventilation_effect
1720
1721          CASE ( 'na_init' )
1722             READ ( 13 )  na_init
1723
1724          CASE ( 'dry_aerosol_radius' )
1725             READ ( 13 )  dry_aerosol_radius
1726
1727          CASE ( 'sigma_bulk' )
1728             READ ( 13 )  sigma_bulk
1729
1730          CASE ( 'aerosol_bulk' )
1731             READ ( 13 )  aerosol_bulk
1732
1733          CASE ( 'curvature_solution_effects_bulk' )
1734             READ ( 13 )  curvature_solution_effects_bulk
1735
1736
1737!          CASE ( 'global_paramter' )
1738!             READ ( 13 )  global_parameter
1739!          CASE ( 'global_array' )
1740!             IF ( .NOT. ALLOCATED( global_array ) )  ALLOCATE( global_array(1:10) )
1741!             READ ( 13 )  global_array
1742
1743          CASE DEFAULT
1744
1745             found = .FALSE.
1746
1747       END SELECT
1748
1749
1750    END SUBROUTINE bcm_rrd_global
1751
1752
1753!------------------------------------------------------------------------------!
1754! Description:
1755! ------------
1756!> This routine reads the respective restart data for the bulk cloud module.
1757!------------------------------------------------------------------------------!
1758    SUBROUTINE bcm_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,       &
1759                              nxr_on_file, nynf, nync, nyn_on_file, nysf,      &
1760                              nysc, nys_on_file, tmp_2d, tmp_3d, found )
1761
1762
1763       USE control_parameters
1764
1765       USE indices
1766
1767       USE pegrid
1768
1769
1770       IMPLICIT NONE
1771
1772       INTEGER(iwp) ::  i               !<
1773       INTEGER(iwp) ::  k               !<
1774       INTEGER(iwp) ::  nxlc            !<
1775       INTEGER(iwp) ::  nxlf            !<
1776       INTEGER(iwp) ::  nxl_on_file     !<
1777       INTEGER(iwp) ::  nxrc            !<
1778       INTEGER(iwp) ::  nxrf            !<
1779       INTEGER(iwp) ::  nxr_on_file     !<
1780       INTEGER(iwp) ::  nync            !<
1781       INTEGER(iwp) ::  nynf            !<
1782       INTEGER(iwp) ::  nyn_on_file     !<
1783       INTEGER(iwp) ::  nysc            !<
1784       INTEGER(iwp) ::  nysf            !<
1785       INTEGER(iwp) ::  nys_on_file     !<
1786
1787       LOGICAL, INTENT(OUT)  ::  found
1788
1789       REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d   !<
1790       REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
1791
1792!
1793!-- Here the reading of user-defined restart data follows:
1794!-- Sample for user-defined output
1795
1796
1797       found = .TRUE.
1798
1799       SELECT CASE ( restart_string(1:length) )
1800
1801          CASE ( 'nc' )
1802             IF ( k == 1 )  READ ( 13 )  tmp_3d
1803             nc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1804                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1805
1806          CASE ( 'nc_av' )
1807             IF ( .NOT. ALLOCATED( nc_av ) )  THEN
1808                ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1809             ENDIF
1810             IF ( k == 1 )  READ ( 13 )  tmp_3d
1811             nc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1812                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1813
1814          CASE ( 'nr' )
1815             IF ( k == 1 )  READ ( 13 )  tmp_3d
1816             nr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1817                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1818
1819          CASE ( 'nr_av' )
1820             IF ( .NOT. ALLOCATED( nr_av ) )  THEN
1821                ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1822             ENDIF
1823             IF ( k == 1 )  READ ( 13 )  tmp_3d
1824             nr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1825                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1826
1827          CASE ( 'precipitation_amount' )
1828             IF ( k == 1 )  READ ( 13 )  tmp_2d
1829             precipitation_amount(nysc-nbgp:nync+nbgp,                   &
1830                                  nxlc-nbgp:nxrc+nbgp)  =                &
1831                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1832
1833          CASE ( 'prr' )
1834             IF ( .NOT. ALLOCATED( prr ) )  THEN
1835                ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1836             ENDIF
1837             IF ( k == 1 )  READ ( 13 )  tmp_3d
1838             prr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =            &
1839                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1840
1841          CASE ( 'prr_av' )
1842             IF ( .NOT. ALLOCATED( prr_av ) )  THEN
1843                ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1844             ENDIF
1845             IF ( k == 1 )  READ ( 13 )  tmp_3d
1846             prr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
1847                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1848
1849          CASE ( 'qc' )
1850             IF ( k == 1 )  READ ( 13 )  tmp_3d
1851             qc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1852                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1853
1854          CASE ( 'qc_av' )
1855             IF ( .NOT. ALLOCATED( qc_av ) )  THEN
1856                ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1857             ENDIF
1858             IF ( k == 1 )  READ ( 13 )  tmp_3d
1859             qc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1860                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1861
1862          CASE ( 'ql' )
1863             IF ( k == 1 )  READ ( 13 )  tmp_3d
1864             ql(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1865                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1866
1867          CASE ( 'ql_av' )
1868             IF ( .NOT. ALLOCATED( ql_av ) )  THEN
1869                ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1870             ENDIF
1871             IF ( k == 1 )  READ ( 13 )  tmp_3d
1872             ql_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1873                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1874
1875          CASE ( 'qr' )
1876             IF ( k == 1 )  READ ( 13 )  tmp_3d
1877             qr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1878                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1879
1880          CASE ( 'qr_av' )
1881             IF ( .NOT. ALLOCATED( qr_av ) )  THEN
1882                ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1883             ENDIF
1884             IF ( k == 1 )  READ ( 13 )  tmp_3d
1885             qr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1886                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1887!
1888          CASE DEFAULT
1889
1890             found = .FALSE.
1891
1892          END SELECT
1893
1894
1895    END SUBROUTINE bcm_rrd_local
1896
1897
1898!------------------------------------------------------------------------------!
1899! Description:
1900! ------------
1901!> This routine writes the respective restart data for the bulk cloud module.
1902!------------------------------------------------------------------------------!
1903    SUBROUTINE bcm_wrd_global
1904
1905
1906       IMPLICIT NONE
1907
1908       CALL wrd_write_string( 'c_sedimentation' )
1909       WRITE ( 14 )  c_sedimentation
1910
1911       CALL wrd_write_string( 'bulk_cloud_model' )
1912       WRITE ( 14 )  bulk_cloud_model
1913
1914       CALL wrd_write_string( 'cloud_scheme' )
1915       WRITE ( 14 )  cloud_scheme
1916
1917       CALL wrd_write_string( 'cloud_water_sedimentation' )
1918       WRITE ( 14 )  cloud_water_sedimentation
1919
1920       CALL wrd_write_string( 'collision_turbulence' )
1921       WRITE ( 14 )  collision_turbulence
1922
1923       CALL wrd_write_string( 'limiter_sedimentation' )
1924       WRITE ( 14 )  limiter_sedimentation
1925
1926       CALL wrd_write_string( 'nc_const' )
1927       WRITE ( 14 )  nc_const
1928
1929       CALL wrd_write_string( 'precipitation' )
1930       WRITE ( 14 )  precipitation
1931
1932       CALL wrd_write_string( 'ventilation_effect' )
1933       WRITE ( 14 )  ventilation_effect
1934
1935       CALL wrd_write_string( 'na_init' )
1936       WRITE ( 14 )  na_init
1937
1938       CALL wrd_write_string( 'dry_aerosol_radius' )
1939       WRITE ( 14 )  dry_aerosol_radius
1940
1941       CALL wrd_write_string( 'sigma_bulk' )
1942       WRITE ( 14 )  sigma_bulk
1943
1944       CALL wrd_write_string( 'aerosol_bulk' )
1945       WRITE ( 14 )  aerosol_bulk
1946
1947       CALL wrd_write_string( 'curvature_solution_effects_bulk' )
1948       WRITE ( 14 )  curvature_solution_effects_bulk
1949
1950
1951! needs preceeding allocation if array
1952!       CALL wrd_write_string( 'global_parameter' )
1953!       WRITE ( 14 )  global_parameter
1954
1955!       IF ( ALLOCATED( inflow_damping_factor ) )  THEN
1956!          CALL wrd_write_string( 'inflow_damping_factor' )
1957!          WRITE ( 14 )  inflow_damping_factor
1958!       ENDIF
1959
1960
1961    END SUBROUTINE bcm_wrd_global
1962
1963
1964!------------------------------------------------------------------------------!
1965! Description:
1966! ------------
1967!> This routine writes the respective restart data for the bulk cloud module.
1968!------------------------------------------------------------------------------!
1969    SUBROUTINE bcm_wrd_local
1970
1971
1972       IMPLICIT NONE
1973
1974       IF ( ALLOCATED( prr ) )  THEN
1975          CALL wrd_write_string( 'prr' )
1976          WRITE ( 14 )  prr
1977       ENDIF
1978
1979       IF ( ALLOCATED( prr_av ) )  THEN
1980          CALL wrd_write_string( 'prr_av' )
1981          WRITE ( 14 )  prr_av
1982       ENDIF
1983
1984       IF ( ALLOCATED( precipitation_amount ) )  THEN
1985          CALL wrd_write_string( 'precipitation_amount' )
1986          WRITE ( 14 )  precipitation_amount
1987       ENDIF
1988
1989       CALL wrd_write_string( 'ql' )
1990       WRITE ( 14 )  ql
1991
1992       IF ( ALLOCATED( ql_av ) )  THEN
1993          CALL wrd_write_string( 'ql_av' )
1994          WRITE ( 14 )  ql_av
1995       ENDIF
1996
1997       CALL wrd_write_string( 'qc' )
1998       WRITE ( 14 )  qc
1999
2000       IF ( ALLOCATED( qc_av ) )  THEN
2001          CALL wrd_write_string( 'qc_av' )
2002          WRITE ( 14 )  qc_av
2003       ENDIF
2004
2005       IF ( microphysics_morrison )  THEN
2006
2007          CALL wrd_write_string( 'nc' )
2008          WRITE ( 14 )  nc
2009
2010          IF ( ALLOCATED( nc_av ) )  THEN
2011             CALL wrd_write_string( 'nc_av' )
2012             WRITE ( 14 )  nc_av
2013          ENDIF
2014
2015       ENDIF
2016
2017       IF ( microphysics_seifert )  THEN
2018
2019          CALL wrd_write_string( 'nr' )
2020          WRITE ( 14 )  nr
2021
2022          IF ( ALLOCATED( nr_av ) )  THEN
2023             CALL wrd_write_string( 'nr_av' )
2024             WRITE ( 14 )  nr_av
2025          ENDIF
2026
2027          CALL wrd_write_string( 'qr' )
2028          WRITE ( 14 )  qr
2029
2030          IF ( ALLOCATED( qr_av ) )  THEN
2031             CALL wrd_write_string( 'qr_av' )
2032             WRITE ( 14 )  qr_av
2033          ENDIF
2034
2035       ENDIF
2036
2037
2038    END SUBROUTINE bcm_wrd_local
2039
2040
2041!------------------------------------------------------------------------------!
2042! Description:
2043! ------------
2044!> Control of microphysics for all grid points
2045!------------------------------------------------------------------------------!
2046    SUBROUTINE bcm_actions
2047
2048       IMPLICIT NONE
2049
2050       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
2051!
2052!--       Calculate vertical profile of the hydrostatic pressure (hyp)
2053          hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
2054          d_exner = exner_function_invers(hyp)
2055          exner = 1.0_wp / exner_function_invers(hyp)
2056          hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
2057!
2058!--       Compute reference density
2059          rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
2060       ENDIF
2061
2062!
2063!--    Compute length of time step
2064       IF ( call_microphysics_at_all_substeps )  THEN
2065          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
2066       ELSE
2067          dt_micro = dt_3d
2068       ENDIF
2069
2070!
2071!--    Reset precipitation rate
2072       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
2073
2074!
2075!--    Compute cloud physics
2076       IF ( microphysics_kessler )  THEN
2077
2078          CALL autoconversion_kessler
2079          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
2080
2081       ELSEIF ( microphysics_seifert )  THEN
2082
2083          CALL adjust_cloud
2084          IF ( microphysics_morrison )  CALL activation
2085          IF ( microphysics_morrison )  CALL condensation
2086          CALL autoconversion
2087          CALL accretion
2088          CALL selfcollection_breakup
2089          CALL evaporation_rain
2090          CALL sedimentation_rain
2091          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
2092
2093       ENDIF
2094
2095       CALL calc_precipitation_amount
2096
2097    END SUBROUTINE bcm_actions
2098
2099!------------------------------------------------------------------------------!
2100! Description:
2101! ------------
2102!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
2103!> evaporation of rain drops due to too small or too big weights
2104!> of rain drops (Stevens and Seifert, 2008).
2105!------------------------------------------------------------------------------!
2106    SUBROUTINE adjust_cloud
2107
2108       IMPLICIT NONE
2109
2110       INTEGER(iwp) ::  i                 !<
2111       INTEGER(iwp) ::  j                 !<
2112       INTEGER(iwp) ::  k                 !<
2113
2114       REAL(wp) ::  flag                  !< flag to indicate first grid level above
2115
2116       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
2117
2118       DO  i = nxlg, nxrg
2119          DO  j = nysg, nyng
2120             DO  k = nzb+1, nzt
2121!
2122!--             Predetermine flag to mask topography
2123                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2124
2125                IF ( qr(k,j,i) <= eps_sb )  THEN
2126                   qr(k,j,i) = 0.0_wp
2127                   nr(k,j,i) = 0.0_wp
2128                ELSE
2129                   IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
2130                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
2131                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
2132                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
2133                   ENDIF
2134                ENDIF
2135
2136                IF ( microphysics_morrison ) THEN
2137                   IF ( qc(k,j,i) <= eps_sb )  THEN
2138                      qc(k,j,i) = 0.0_wp
2139                      nc(k,j,i) = 0.0_wp
2140                   ELSE
2141                      IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
2142                         nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin * flag
2143                      ENDIF
2144                   ENDIF
2145                ENDIF
2146
2147             ENDDO
2148          ENDDO
2149       ENDDO
2150
2151       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
2152
2153    END SUBROUTINE adjust_cloud
2154
2155!------------------------------------------------------------------------------!
2156! Description:
2157! ------------
2158!> Calculate number of activated condensation nucleii after simple activation
2159!> scheme of Twomey, 1959.
2160!------------------------------------------------------------------------------!
2161    SUBROUTINE activation
2162
2163       IMPLICIT NONE
2164
2165       INTEGER(iwp) ::  i                 !<
2166       INTEGER(iwp) ::  j                 !<
2167       INTEGER(iwp) ::  k                 !<
2168
2169       REAL(wp)     ::  activ             !<
2170       REAL(wp)     ::  afactor           !<
2171       REAL(wp)     ::  beta_act          !<
2172       REAL(wp)     ::  bfactor           !<
2173       REAL(wp)     ::  k_act             !<
2174       REAL(wp)     ::  n_act             !<
2175       REAL(wp)     ::  n_ccn             !<
2176       REAL(wp)     ::  s_0               !<
2177       REAL(wp)     ::  sat_max           !<
2178       REAL(wp)     ::  sigma             !<
2179       REAL(wp)     ::  sigma_act         !<
2180
2181       REAL(wp) ::  flag                  !< flag to indicate first grid level above
2182
2183       CALL cpu_log( log_point_s(65), 'activation', 'start' )
2184
2185       DO  i = nxlg, nxrg
2186          DO  j = nysg, nyng
2187             DO  k = nzb+1, nzt
2188!
2189!--             Predetermine flag to mask topography
2190                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2191
2192!
2193!--             Call calculation of supersaturation located in subroutine
2194                CALL supersaturation ( i, j, k )
2195!
2196!--             Prescribe parameters for activation
2197!--             (see: Bott + Trautmann, 2002, Atm. Res., 64)
2198                k_act  = 0.7_wp
2199                activ  = 0.0_wp
2200
2201
2202                IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
2203!
2204!--                Compute the number of activated Aerosols
2205!--                (see: Twomey, 1959, Pure and applied Geophysics, 43)
2206                   n_act     = na_init * sat**k_act
2207!
2208!--                Compute the number of cloud droplets
2209!--                (see: Morrison + Grabowski, 2007, JAS, 64)
2210!                  activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro
2211
2212!
2213!--                Compute activation rate after Khairoutdinov and Kogan
2214!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
2215                   sat_max = 1.0_wp / 100.0_wp
2216                   activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN       &
2217                      ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /    &
2218                       dt_micro
2219                ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN
2220!
2221!--                Curvature effect (afactor) with surface tension
2222!--                parameterization by Straka (2009)
2223                   sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
2224                   afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
2225!
2226!--                Solute effect (bfactor)
2227                   bfactor = vanthoff * molecular_weight_of_water *            &
2228                       rho_s / ( molecular_weight_of_solute * rho_l )
2229
2230!
2231!--                Prescribe power index that describes the soluble fraction
2232!--                of an aerosol particle (beta)
2233!--                (see: Morrison + Grabowski, 2007, JAS, 64)
2234                   beta_act  = 0.5_wp
2235                   sigma_act = sigma_bulk**( 1.0_wp + beta_act )
2236!
2237!--                Calculate mean geometric supersaturation (s_0) with
2238!--                parameterization by Khvorostyanov and Curry (2006)
2239                   s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *    &
2240                       ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
2241
2242!
2243!--                Calculate number of activated CCN as a function of
2244!--                supersaturation and taking Koehler theory into account
2245!--                (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
2246                   n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(              &
2247                      LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
2248                   activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
2249                ENDIF
2250
2251                nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)
2252
2253             ENDDO
2254          ENDDO
2255       ENDDO
2256
2257       CALL cpu_log( log_point_s(65), 'activation', 'stop' )
2258
2259    END SUBROUTINE activation
2260
2261
2262!------------------------------------------------------------------------------!
2263! Description:
2264! ------------
2265!> Calculate condensation rate for cloud water content (after Khairoutdinov and
2266!> Kogan, 2000).
2267!------------------------------------------------------------------------------!
2268    SUBROUTINE condensation
2269
2270       IMPLICIT NONE
2271
2272       INTEGER(iwp) ::  i                 !<
2273       INTEGER(iwp) ::  j                 !<
2274       INTEGER(iwp) ::  k                 !<
2275
2276       REAL(wp)     ::  cond              !<
2277       REAL(wp)     ::  cond_max          !<
2278       REAL(wp)     ::  dc                !<
2279       REAL(wp)     ::  evap              !<
2280       REAL(wp)     ::  g_fac             !<
2281       REAL(wp)     ::  nc_0              !<
2282       REAL(wp)     ::  temp              !<
2283       REAL(wp)     ::  xc                !<
2284
2285       REAL(wp) ::  flag                  !< flag to indicate first grid level above
2286
2287       CALL cpu_log( log_point_s(66), 'condensation', 'start' )
2288
2289       DO  i = nxlg, nxrg
2290          DO  j = nysg, nyng
2291             DO  k = nzb+1, nzt
2292!
2293!--             Predetermine flag to mask topography
2294                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2295!
2296!--             Call calculation of supersaturation
2297                CALL supersaturation ( i, j, k )
2298!
2299!--             Actual temperature:
2300                temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
2301
2302                g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
2303                                    l_v / ( thermal_conductivity_l * temp )    &
2304                                    + r_v * temp / ( diff_coeff_l * e_s )      &
2305                                    )
2306!
2307!--             Mean weight of cloud drops
2308                IF ( nc(k,j,i) <= 0.0_wp) CYCLE
2309                xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
2310!
2311!--             Weight averaged diameter of cloud drops:
2312                dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2313!
2314!--             Integral diameter of cloud drops
2315                nc_0 = nc(k,j,i) * dc
2316!
2317!--             Condensation needs only to be calculated in supersaturated regions
2318                IF ( sat > 0.0_wp )  THEN
2319!
2320!--                Condensation rate of cloud water content
2321!--                after KK scheme.
2322!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
2323                   cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2324                   cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
2325                   cond      = MIN( cond, cond_max / dt_micro )
2326
2327                   qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
2328                ELSEIF ( sat < 0.0_wp ) THEN
2329                   evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2330                   evap      = MAX( evap, -qc(k,j,i) / dt_micro )
2331
2332                   qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag
2333                ENDIF
2334                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
2335                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin
2336                ENDIF
2337             ENDDO
2338          ENDDO
2339       ENDDO
2340
2341       CALL cpu_log( log_point_s(66), 'condensation', 'stop' )
2342
2343    END SUBROUTINE condensation
2344
2345
2346!------------------------------------------------------------------------------!
2347! Description:
2348! ------------
2349!> Autoconversion rate (Seifert and Beheng, 2006).
2350!------------------------------------------------------------------------------!
2351    SUBROUTINE autoconversion
2352
2353       IMPLICIT NONE
2354
2355       INTEGER(iwp) ::  i                 !<
2356       INTEGER(iwp) ::  j                 !<
2357       INTEGER(iwp) ::  k                 !<
2358
2359       REAL(wp)     ::  alpha_cc          !<
2360       REAL(wp)     ::  autocon           !<
2361       REAL(wp)     ::  dissipation       !<
2362       REAL(wp)     ::  flag              !< flag to mask topography grid points
2363       REAL(wp)     ::  k_au              !<
2364       REAL(wp)     ::  l_mix             !<
2365       REAL(wp)     ::  nc_auto           !<
2366       REAL(wp)     ::  nu_c              !<
2367       REAL(wp)     ::  phi_au            !<
2368       REAL(wp)     ::  r_cc              !<
2369       REAL(wp)     ::  rc                !<
2370       REAL(wp)     ::  re_lambda         !<
2371       REAL(wp)     ::  sigma_cc          !<
2372       REAL(wp)     ::  tau_cloud         !<
2373       REAL(wp)     ::  xc                !<
2374
2375       CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
2376
2377       DO  i = nxlg, nxrg
2378          DO  j = nysg, nyng
2379             DO  k = nzb+1, nzt
2380!
2381!--             Predetermine flag to mask topography
2382                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2383
2384                IF ( microphysics_morrison ) THEN
2385                   nc_auto = nc(k,j,i)
2386                ELSE
2387                   nc_auto = nc_const
2388                ENDIF
2389
2390                IF ( qc(k,j,i) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
2391
2392                   k_au = k_cc / ( 20.0_wp * x0 )
2393!
2394!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
2395!--                (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
2396                   tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) +         &
2397                                    qc(k,j,i) ), 0.0_wp )
2398!
2399!--                Universal function for autoconversion process
2400!--                (Seifert and Beheng, 2006):
2401                   phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
2402                            ( 1.0_wp - tau_cloud**0.68_wp )**3
2403!
2404!--                Shape parameter of gamma distribution (Geoffroy et al., 2010):
2405!--                (Use constant nu_c = 1.0_wp instead?)
2406                   nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
2407!
2408!--                Mean weight of cloud droplets:
2409                   xc = MAX( hyrho(k) * qc(k,j,i) / nc_auto, xcmin) 
2410!
2411!--                Parameterized turbulence effects on autoconversion (Seifert,
2412!--                Nuijens and Stevens, 2010)
2413                   IF ( collision_turbulence )  THEN
2414!
2415!--                   Weight averaged radius of cloud droplets:
2416                      rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2417
2418                      alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
2419                      r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
2420                      sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
2421!
2422!--                   Mixing length (neglecting distance to ground and
2423!--                   stratification)
2424                      l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
2425!
2426!--                   Limit dissipation rate according to Seifert, Nuijens and
2427!--                   Stevens (2010)
2428                      dissipation = MIN( 0.06_wp, diss(k,j,i) )
2429!
2430!--                   Compute Taylor-microscale Reynolds number:
2431                      re_lambda = 6.0_wp / 11.0_wp *                           &
2432                                  ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *   &
2433                                  SQRT( 15.0_wp / kin_vis_air ) *              &
2434                                  dissipation**( 1.0_wp / 6.0_wp )
2435!
2436!--                   The factor of 1.0E4 is needed to convert the dissipation
2437!--                   rate from m2 s-3 to cm2 s-3.
2438                      k_au = k_au * ( 1.0_wp +                                 &
2439                                      dissipation * 1.0E4_wp *                 &
2440                                      ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
2441                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &
2442                                                                      r_cc ) / &
2443                                                        sigma_cc )**2          &
2444                                                      ) + beta_cc              &
2445                                      )                                        &
2446                                    )
2447                   ENDIF
2448!
2449!--                Autoconversion rate (Seifert and Beheng, 2006):
2450                   autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /    &
2451                             ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *     &
2452                             ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * &
2453                             rho_surface
2454                   autocon = MIN( autocon, qc(k,j,i) / dt_micro )
2455
2456                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
2457                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
2458                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro  &
2459                                                              * flag
2460                   IF ( microphysics_morrison ) THEN
2461                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *         &
2462                                    autocon / x0 * hyrho(k) * dt_micro * flag )
2463                   ENDIF
2464
2465                ENDIF
2466
2467             ENDDO
2468          ENDDO
2469       ENDDO
2470
2471       CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
2472
2473    END SUBROUTINE autoconversion
2474
2475
2476!------------------------------------------------------------------------------!
2477! Description:
2478! ------------
2479!> Autoconversion process (Kessler, 1969).
2480!------------------------------------------------------------------------------!
2481    SUBROUTINE autoconversion_kessler
2482
2483
2484       IMPLICIT NONE
2485
2486       INTEGER(iwp) ::  i      !<
2487       INTEGER(iwp) ::  j      !<
2488       INTEGER(iwp) ::  k      !<
2489       INTEGER(iwp) ::  k_wall !< topgraphy top index
2490
2491       REAL(wp)    ::  dqdt_precip !<
2492       REAL(wp)    ::  flag        !< flag to mask topography grid points
2493
2494       DO  i = nxlg, nxrg
2495          DO  j = nysg, nyng
2496!
2497!--          Determine vertical index of topography top
2498             k_wall = get_topography_top_index_ji( j, i, 's' )
2499             DO  k = nzb+1, nzt
2500!
2501!--             Predetermine flag to mask topography
2502                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2503
2504                IF ( qc(k,j,i) > ql_crit )  THEN
2505                   dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
2506                ELSE
2507                   dqdt_precip = 0.0_wp
2508                ENDIF
2509
2510                qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
2511                q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
2512                pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp *     &
2513                                        d_exner(k)              * flag
2514
2515!
2516!--             Compute the rain rate (stored on surface grid point)
2517                prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
2518
2519             ENDDO
2520          ENDDO
2521       ENDDO
2522
2523   END SUBROUTINE autoconversion_kessler
2524
2525
2526!------------------------------------------------------------------------------!
2527! Description:
2528! ------------
2529!> Accretion rate (Seifert and Beheng, 2006).
2530!------------------------------------------------------------------------------!
2531    SUBROUTINE accretion
2532
2533       IMPLICIT NONE
2534
2535       INTEGER(iwp) ::  i                 !<
2536       INTEGER(iwp) ::  j                 !<
2537       INTEGER(iwp) ::  k                 !<
2538
2539       REAL(wp)     ::  accr              !<
2540       REAL(wp)     ::  flag              !< flag to mask topography grid points
2541       REAL(wp)     ::  k_cr              !<
2542       REAL(wp)     ::  nc_accr           !<
2543       REAL(wp)     ::  phi_ac            !<
2544       REAL(wp)     ::  tau_cloud         !<
2545       REAL(wp)     ::  xc                !<
2546
2547
2548       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
2549
2550       DO  i = nxlg, nxrg
2551          DO  j = nysg, nyng
2552             DO  k = nzb+1, nzt
2553!
2554!--             Predetermine flag to mask topography
2555                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2556
2557                IF ( microphysics_morrison ) THEN
2558                   nc_accr = nc(k,j,i)
2559                ELSE
2560                   nc_accr = nc_const
2561                ENDIF
2562
2563                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )    &
2564                                             .AND.  ( nc_accr > eps_mr ) ) THEN
2565!
2566!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
2567                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
2568!
2569!--                Universal function for accretion process (Seifert and
2570!--                Beheng, 2001):
2571                   phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
2572
2573!
2574!--                Mean weight of cloud drops
2575                   xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
2576!
2577!--                Parameterized turbulence effects on autoconversion (Seifert,
2578!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
2579!--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
2580                   IF ( collision_turbulence )  THEN
2581                      k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                      &
2582                                       MIN( 600.0_wp,                          &
2583                                            diss(k,j,i) * 1.0E4_wp )**0.25_wp  &
2584                                     )
2585                   ELSE
2586                      k_cr = k_cr0
2587                   ENDIF
2588!
2589!--                Accretion rate (Seifert and Beheng, 2006):
2590                   accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *              &
2591                          SQRT( rho_surface * hyrho(k) )
2592                   accr = MIN( accr, qc(k,j,i) / dt_micro )
2593
2594                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
2595                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
2596                   IF ( microphysics_morrison )  THEN
2597                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i),                  &
2598                                         accr / xc * hyrho(k) * dt_micro * flag)
2599                   ENDIF
2600
2601                ENDIF
2602
2603             ENDDO
2604          ENDDO
2605       ENDDO
2606
2607       CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
2608
2609    END SUBROUTINE accretion
2610
2611
2612!------------------------------------------------------------------------------!
2613! Description:
2614! ------------
2615!> Collisional breakup rate (Seifert, 2008).
2616!------------------------------------------------------------------------------!
2617    SUBROUTINE selfcollection_breakup
2618
2619       IMPLICIT NONE
2620
2621       INTEGER(iwp) ::  i                 !<
2622       INTEGER(iwp) ::  j                 !<
2623       INTEGER(iwp) ::  k                 !<
2624
2625       REAL(wp)     ::  breakup           !<
2626       REAL(wp)     ::  dr                !<
2627       REAL(wp)     ::  flag              !< flag to mask topography grid points
2628       REAL(wp)     ::  phi_br            !<
2629       REAL(wp)     ::  selfcoll          !<
2630
2631       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
2632
2633       DO  i = nxlg, nxrg
2634          DO  j = nysg, nyng
2635             DO  k = nzb+1, nzt
2636!
2637!--             Predetermine flag to mask topography
2638                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2639
2640                IF ( qr(k,j,i) > eps_sb )  THEN
2641!
2642!--                Selfcollection rate (Seifert and Beheng, 2001):
2643                   selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) *                   &
2644                              SQRT( hyrho(k) * rho_surface )
2645!
2646!--                Weight averaged diameter of rain drops:
2647                   dr = ( hyrho(k) * qr(k,j,i) /                               &
2648                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
2649!
2650!--                Collisional breakup rate (Seifert, 2008):
2651                   IF ( dr >= 0.3E-3_wp )  THEN
2652                      phi_br  = k_br * ( dr - 1.1E-3_wp )
2653                      breakup = selfcoll * ( phi_br + 1.0_wp )
2654                   ELSE
2655                      breakup = 0.0_wp
2656                   ENDIF
2657
2658                   selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
2659                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
2660
2661                ENDIF
2662             ENDDO
2663          ENDDO
2664       ENDDO
2665
2666       CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
2667
2668    END SUBROUTINE selfcollection_breakup
2669
2670
2671!------------------------------------------------------------------------------!
2672! Description:
2673! ------------
2674!> Evaporation of precipitable water. Condensation is neglected for
2675!> precipitable water.
2676!------------------------------------------------------------------------------!
2677    SUBROUTINE evaporation_rain
2678
2679       IMPLICIT NONE
2680
2681       INTEGER(iwp) ::  i                 !<
2682       INTEGER(iwp) ::  j                 !<
2683       INTEGER(iwp) ::  k                 !<
2684
2685       REAL(wp)     ::  dr                !<
2686       REAL(wp)     ::  evap              !<
2687       REAL(wp)     ::  evap_nr           !<
2688       REAL(wp)     ::  f_vent            !<
2689       REAL(wp)     ::  flag              !< flag to mask topography grid points
2690       REAL(wp)     ::  g_evap            !<
2691       REAL(wp)     ::  lambda_r          !<
2692       REAL(wp)     ::  mu_r              !<
2693       REAL(wp)     ::  mu_r_2            !<
2694       REAL(wp)     ::  mu_r_5d2          !<
2695       REAL(wp)     ::  nr_0              !<
2696       REAL(wp)     ::  temp              !<
2697       REAL(wp)     ::  xr                !<
2698
2699       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
2700
2701       DO  i = nxlg, nxrg
2702          DO  j = nysg, nyng
2703             DO  k = nzb+1, nzt
2704!
2705!--             Predetermine flag to mask topography
2706                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2707
2708                IF ( qr(k,j,i) > eps_sb )  THEN
2709
2710!
2711!--                Call calculation of supersaturation 
2712                   CALL supersaturation ( i, j, k )
2713!
2714!--                Evaporation needs only to be calculated in subsaturated regions
2715                   IF ( sat < 0.0_wp )  THEN
2716!
2717!--                   Actual temperature:
2718                      temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
2719
2720                      g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *     &
2721                                          l_v / ( thermal_conductivity_l * temp ) &
2722                                          + r_v * temp / ( diff_coeff_l * e_s )   &
2723                                        )
2724!
2725!--                   Mean weight of rain drops
2726                      xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
2727!
2728!--                   Weight averaged diameter of rain drops:
2729                      dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
2730!
2731!--                   Compute ventilation factor and intercept parameter
2732!--                   (Seifert and Beheng, 2006; Seifert, 2008):
2733                      IF ( ventilation_effect )  THEN
2734!
2735!--                      Shape parameter of gamma distribution (Milbrandt and Yau,
2736!--                      2005; Stevens and Seifert, 2008):
2737                         mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *          &
2738                                                          ( dr - 1.4E-3_wp ) ) )
2739!
2740!--                      Slope parameter of gamma distribution (Seifert, 2008):
2741                         lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *  &
2742                                      ( mu_r + 1.0_wp )                        &
2743                                    )**( 1.0_wp / 3.0_wp ) / dr
2744
2745                         mu_r_2   = mu_r + 2.0_wp
2746                         mu_r_5d2 = mu_r + 2.5_wp
2747
2748                         f_vent = a_vent * gamm( mu_r_2 ) *                     &
2749                                  lambda_r**( -mu_r_2 ) + b_vent *              &
2750                                  schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
2751                                  gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) *  &
2752                                  ( 1.0_wp -                                    &
2753                                    0.5_wp * ( b_term / a_term ) *              &
2754                                    ( lambda_r / ( c_term + lambda_r )          &
2755                                    )**mu_r_5d2 -                               &
2756                                    0.125_wp * ( b_term / a_term )**2 *         &
2757                                    ( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
2758                                    )**mu_r_5d2 -                               &
2759                                    0.0625_wp * ( b_term / a_term )**3 *        &
2760                                    ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
2761                                    )**mu_r_5d2 -                               &
2762                                    0.0390625_wp * ( b_term / a_term )**4 *     &
2763                                    ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
2764                                    )**mu_r_5d2                                 &
2765                                  )
2766
2767                         nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
2768                                  gamm( mu_r + 1.0_wp )
2769                      ELSE
2770                         f_vent = 1.0_wp
2771                         nr_0   = nr(k,j,i) * dr
2772                      ENDIF
2773   !
2774   !--                Evaporation rate of rain water content (Seifert and
2775   !--                Beheng, 2006):
2776                      evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
2777                                hyrho(k)
2778                      evap    = MAX( evap, -qr(k,j,i) / dt_micro )
2779                      evap_nr = MAX( c_evap * evap / xr * hyrho(k),            &
2780                                     -nr(k,j,i) / dt_micro )
2781
2782                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro * flag
2783                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
2784
2785                   ENDIF
2786                ENDIF
2787
2788             ENDDO
2789          ENDDO
2790       ENDDO
2791
2792       CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
2793
2794    END SUBROUTINE evaporation_rain
2795
2796
2797!------------------------------------------------------------------------------!
2798! Description:
2799! ------------
2800!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
2801!------------------------------------------------------------------------------!
2802    SUBROUTINE sedimentation_cloud
2803
2804
2805       IMPLICIT NONE
2806
2807       INTEGER(iwp) ::  i             !<
2808       INTEGER(iwp) ::  j             !<
2809       INTEGER(iwp) ::  k             !<
2810
2811       REAL(wp) ::  flag              !< flag to mask topography grid points
2812       REAL(wp) ::  nc_sedi           !<
2813
2814       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
2815       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc !<
2816
2817
2818       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
2819
2820       sed_qc(nzt+1) = 0.0_wp
2821       sed_nc(nzt+1) = 0.0_wp
2822
2823       DO  i = nxlg, nxrg
2824          DO  j = nysg, nyng
2825             DO  k = nzt, nzb+1, -1
2826!
2827!--             Predetermine flag to mask topography
2828                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2829
2830                IF ( microphysics_morrison ) THEN
2831                   nc_sedi = nc(k,j,i)
2832                ELSE
2833                   nc_sedi = nc_const
2834                ENDIF
2835
2836!
2837!--             Sedimentation fluxes for number concentration are only calculated
2838!--             for cloud_scheme = 'morrison'
2839                IF ( microphysics_morrison ) THEN
2840                   IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
2841                      sed_nc(k) = sed_qc_const *                               &
2842                              ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
2843                              ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
2844                   ELSE
2845                      sed_nc(k) = 0.0_wp
2846                   ENDIF
2847
2848                   sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *           &
2849                                    nc(k,j,i) / dt_micro + sed_nc(k+1)         &
2850                                  ) * flag
2851
2852                   nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *       &
2853                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
2854                ENDIF
2855
2856                IF ( qc(k,j,i) > eps_sb .AND.  nc_sedi > eps_mr )  THEN
2857                   sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *  &
2858                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
2859                                                                           flag
2860                ELSE
2861                   sed_qc(k) = 0.0_wp
2862                ENDIF
2863
2864                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
2865                                            dt_micro + sed_qc(k+1)             &
2866                               ) * flag
2867
2868                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
2869                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
2870                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
2871                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
2872                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
2873                                        ddzu(k+1) / hyrho(k) * lv_d_cp *       &
2874                                        d_exner(k) * dt_micro            * flag
2875
2876!
2877!--             Compute the precipitation rate due to cloud (fog) droplets
2878                IF ( call_microphysics_at_all_substeps )  THEN
2879                   prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
2880                                * weight_substep(intermediate_timestep_count)  &
2881                                * flag
2882                ELSE
2883                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
2884                ENDIF
2885
2886             ENDDO
2887          ENDDO
2888       ENDDO
2889
2890       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
2891
2892    END SUBROUTINE sedimentation_cloud
2893
2894
2895!------------------------------------------------------------------------------!
2896! Description:
2897! ------------
2898!> Computation of sedimentation flux. Implementation according to Stevens
2899!> and Seifert (2008). Code is based on UCLA-LES.
2900!------------------------------------------------------------------------------!
2901    SUBROUTINE sedimentation_rain
2902
2903       IMPLICIT NONE
2904
2905       INTEGER(iwp) ::  i             !< running index x direction
2906       INTEGER(iwp) ::  j             !< running index y direction
2907       INTEGER(iwp) ::  k             !< running index z direction
2908       INTEGER(iwp) ::  k_run         !<
2909       INTEGER(iwp) ::  m             !< running index surface elements
2910       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
2911       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
2912
2913       REAL(wp)     ::  c_run                      !<
2914       REAL(wp)     ::  d_max                      !<
2915       REAL(wp)     ::  d_mean                     !<
2916       REAL(wp)     ::  d_min                      !<
2917       REAL(wp)     ::  dr                         !<
2918       REAL(wp)     ::  flux                       !<
2919       REAL(wp)     ::  flag                       !< flag to mask topography grid points
2920       REAL(wp)     ::  lambda_r                   !<
2921       REAL(wp)     ::  mu_r                       !<
2922       REAL(wp)     ::  z_run                      !<
2923
2924       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
2925       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
2926       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
2927       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
2928       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
2929       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
2930       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
2931       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
2932
2933       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
2934
2935!
2936!--    Compute velocities
2937       DO  i = nxlg, nxrg
2938          DO  j = nysg, nyng
2939             DO  k = nzb+1, nzt
2940!
2941!--             Predetermine flag to mask topography
2942                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2943
2944                IF ( qr(k,j,i) > eps_sb )  THEN
2945!
2946!--                Weight averaged diameter of rain drops:
2947                   dr = ( hyrho(k) * qr(k,j,i) /                               &
2948                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
2949!
2950!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
2951!--                Stevens and Seifert, 2008):
2952                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
2953                                                     ( dr - 1.4E-3_wp ) ) )
2954!
2955!--                Slope parameter of gamma distribution (Seifert, 2008):
2956                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
2957                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
2958
2959                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
2960                                               a_term - b_term * ( 1.0_wp +    &
2961                                                  c_term /                     &
2962                                                  lambda_r )**( -1.0_wp *      &
2963                                                  ( mu_r + 1.0_wp ) )          &
2964                                              )                                &
2965                                ) * flag
2966
2967                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
2968                                               a_term - b_term * ( 1.0_wp +    &
2969                                                  c_term /                     &
2970                                                  lambda_r )**( -1.0_wp *      &
2971                                                  ( mu_r + 4.0_wp ) )          &
2972                                             )                                 &
2973                                ) * flag
2974                ELSE
2975                   w_nr(k) = 0.0_wp
2976                   w_qr(k) = 0.0_wp
2977                ENDIF
2978             ENDDO
2979!
2980!--          Adjust boundary values using surface data type.
2981!--          Upward-facing
2982             surf_s = bc_h(0)%start_index(j,i)
2983             surf_e = bc_h(0)%end_index(j,i)
2984             DO  m = surf_s, surf_e
2985                k         = bc_h(0)%k(m)
2986                w_nr(k-1) = w_nr(k)
2987                w_qr(k-1) = w_qr(k)
2988             ENDDO
2989!
2990!--          Downward-facing
2991             surf_s = bc_h(1)%start_index(j,i)
2992             surf_e = bc_h(1)%end_index(j,i)
2993             DO  m = surf_s, surf_e
2994                k         = bc_h(1)%k(m)
2995                w_nr(k+1) = w_nr(k)
2996                w_qr(k+1) = w_qr(k)
2997             ENDDO
2998!
2999!--          Model top boundary value
3000             w_nr(nzt+1) = 0.0_wp
3001             w_qr(nzt+1) = 0.0_wp
3002!
3003!--          Compute Courant number
3004             DO  k = nzb+1, nzt
3005!
3006!--             Predetermine flag to mask topography
3007                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3008
3009                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
3010                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
3011                          dt_micro * ddzu(k) * flag
3012                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
3013                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
3014                          dt_micro * ddzu(k) * flag
3015             ENDDO
3016!
3017!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
3018             IF ( limiter_sedimentation )  THEN
3019
3020                DO k = nzb+1, nzt
3021!
3022!--                Predetermine flag to mask topography
3023                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3024
3025                   d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
3026                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
3027                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
3028
3029                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
3030                                                              2.0_wp * d_max,  &
3031                                                              ABS( d_mean ) )  &
3032                                                      * flag
3033
3034                   d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
3035                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
3036                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
3037
3038                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
3039                                                              2.0_wp * d_max,  &
3040                                                              ABS( d_mean ) )
3041                ENDDO
3042
3043             ELSE
3044
3045                nr_slope = 0.0_wp
3046                qr_slope = 0.0_wp
3047
3048             ENDIF
3049
3050             sed_nr(nzt+1) = 0.0_wp
3051             sed_qr(nzt+1) = 0.0_wp
3052!
3053!--          Compute sedimentation flux
3054             DO  k = nzt, nzb+1, -1
3055!
3056!--             Predetermine flag to mask topography
3057                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3058!
3059!--             Sum up all rain drop number densities which contribute to the flux
3060!--             through k-1/2
3061                flux  = 0.0_wp
3062                z_run = 0.0_wp ! height above z(k)
3063                k_run = k
3064                c_run = MIN( 1.0_wp, c_nr(k) )
3065                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
3066                   flux  = flux + hyrho(k_run) *                               &
3067                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
3068                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)  &
3069                                              * flag
3070                   z_run = z_run + dzu(k_run) * flag
3071                   k_run = k_run + 1          * flag
3072                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )    &
3073                                              * flag
3074                ENDDO
3075!
3076!--             It is not allowed to sediment more rain drop number density than
3077!--             available
3078                flux = MIN( flux,                                              &
3079                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
3080                            dt_micro                                           &
3081                          )
3082
3083                sed_nr(k) = flux / dt_micro * flag
3084                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
3085                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
3086!
3087!--             Sum up all rain water content which contributes to the flux
3088!--             through k-1/2
3089                flux  = 0.0_wp
3090                z_run = 0.0_wp ! height above z(k)
3091                k_run = k
3092                c_run = MIN( 1.0_wp, c_qr(k) )
3093
3094                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
3095
3096                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
3097                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
3098                                  0.5_wp ) * c_run * dzu(k_run) * flag
3099                   z_run = z_run + dzu(k_run)                   * flag
3100                   k_run = k_run + 1                            * flag
3101                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )    &
3102                                                                * flag
3103
3104                ENDDO
3105!
3106!--             It is not allowed to sediment more rain water content than
3107!--             available
3108                flux = MIN( flux,                                              &
3109                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
3110                            dt_micro                                           &
3111                          )
3112
3113                sed_qr(k) = flux / dt_micro * flag
3114
3115                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
3116                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
3117                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
3118                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
3119                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
3120                                        ddzu(k+1) / hyrho(k) * lv_d_cp *       &
3121                                        d_exner(k) * dt_micro            * flag
3122!
3123!--             Compute the rain rate
3124                IF ( call_microphysics_at_all_substeps )  THEN
3125                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
3126                                * weight_substep(intermediate_timestep_count) &
3127                                * flag
3128                ELSE
3129                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
3130                ENDIF
3131
3132             ENDDO
3133          ENDDO
3134       ENDDO
3135
3136       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
3137
3138    END SUBROUTINE sedimentation_rain
3139
3140
3141!------------------------------------------------------------------------------!
3142! Description:
3143! ------------
3144!> Computation of the precipitation amount due to gravitational settling of
3145!> rain and cloud (fog) droplets
3146!------------------------------------------------------------------------------!
3147    SUBROUTINE calc_precipitation_amount
3148
3149       IMPLICIT NONE
3150
3151       INTEGER(iwp) ::  i             !< running index x direction
3152       INTEGER(iwp) ::  j             !< running index y direction
3153       INTEGER(iwp) ::  k             !< running index y direction
3154       INTEGER(iwp) ::  m             !< running index surface elements
3155
3156       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
3157            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
3158            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
3159       THEN
3160!
3161!--       Run over all upward-facing surface elements, i.e. non-natural,
3162!--       natural and urban
3163          DO  m = 1, bc_h(0)%ns
3164             i = bc_h(0)%i(m)
3165             j = bc_h(0)%j(m)
3166             k = bc_h(0)%k(m)
3167             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
3168                                               prr(k,j,i) * hyrho(k) * dt_3d
3169          ENDDO
3170
3171       ENDIF
3172
3173    END SUBROUTINE calc_precipitation_amount
3174
3175
3176!------------------------------------------------------------------------------!
3177! Description:
3178! ------------
3179!> Control of microphysics for grid points i,j
3180!------------------------------------------------------------------------------!
3181
3182    SUBROUTINE bcm_actions_ij( i, j )
3183
3184       IMPLICIT NONE
3185
3186       INTEGER(iwp) ::  i                 !<
3187       INTEGER(iwp) ::  j                 !<
3188
3189       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
3190!
3191!--       Calculate vertical profile of the hydrostatic pressure (hyp)
3192          hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
3193          d_exner = exner_function_invers(hyp)
3194          exner = 1.0_wp / exner_function_invers(hyp)
3195          hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
3196!
3197!--       Compute reference density
3198          rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
3199       ENDIF
3200
3201!
3202!--    Compute length of time step
3203       IF ( call_microphysics_at_all_substeps )  THEN
3204          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
3205       ELSE
3206          dt_micro = dt_3d
3207       ENDIF
3208!
3209!--    Reset precipitation rate
3210       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
3211
3212!
3213!--    Compute cloud physics
3214       IF( microphysics_kessler )  THEN
3215
3216          CALL autoconversion_kessler( i,j )
3217          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
3218
3219       ELSEIF ( microphysics_seifert )  THEN
3220
3221          CALL adjust_cloud( i,j )
3222          IF ( microphysics_morrison ) CALL activation( i,j )
3223          IF ( microphysics_morrison ) CALL condensation( i,j )
3224          CALL autoconversion( i,j )
3225          CALL accretion( i,j )
3226          CALL selfcollection_breakup( i,j )
3227          CALL evaporation_rain( i,j )
3228          CALL sedimentation_rain( i,j )
3229          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
3230
3231       ENDIF
3232
3233       CALL calc_precipitation_amount( i,j )
3234
3235    END SUBROUTINE bcm_actions_ij
3236
3237!------------------------------------------------------------------------------!
3238! Description:
3239! ------------
3240!> Adjust number of raindrops to avoid nonlinear effects in
3241!> sedimentation and evaporation of rain drops due to too small or
3242!> too big weights of rain drops (Stevens and Seifert, 2008).
3243!> The same procedure is applied to cloud droplets if they are determined
3244!> prognostically. Call for grid point i,j
3245!------------------------------------------------------------------------------!
3246    SUBROUTINE adjust_cloud_ij( i, j )
3247
3248       IMPLICIT NONE
3249
3250       INTEGER(iwp) ::  i                 !<
3251       INTEGER(iwp) ::  j                 !<
3252       INTEGER(iwp) ::  k                 !<
3253
3254       REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
3255
3256       DO  k = nzb+1, nzt
3257!
3258!--       Predetermine flag to mask topography
3259          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3260
3261          IF ( qr(k,j,i) <= eps_sb )  THEN
3262             qr(k,j,i) = 0.0_wp
3263             nr(k,j,i) = 0.0_wp
3264          ELSE
3265!
3266!--          Adjust number of raindrops to avoid nonlinear effects in
3267!--          sedimentation and evaporation of rain drops due to too small or
3268!--          too big weights of rain drops (Stevens and Seifert, 2008).
3269             IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
3270                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
3271             ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
3272                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
3273             ENDIF
3274
3275          ENDIF
3276
3277          IF ( microphysics_morrison ) THEN
3278             IF ( qc(k,j,i) <= eps_sb )  THEN
3279                qc(k,j,i) = 0.0_wp
3280                nc(k,j,i) = 0.0_wp
3281             ELSE
3282                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
3283                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xamin * flag
3284                ENDIF
3285             ENDIF
3286          ENDIF
3287
3288       ENDDO
3289
3290    END SUBROUTINE adjust_cloud_ij
3291
3292!------------------------------------------------------------------------------!
3293! Description:
3294! ------------
3295!> Calculate number of activated condensation nucleii after simple activation
3296!> scheme of Twomey, 1959.
3297!------------------------------------------------------------------------------!
3298    SUBROUTINE activation_ij( i, j )
3299
3300       IMPLICIT NONE
3301
3302       INTEGER(iwp) ::  i                 !<
3303       INTEGER(iwp) ::  j                 !<
3304       INTEGER(iwp) ::  k                 !<
3305
3306       REAL(wp)     ::  activ             !<
3307       REAL(wp)     ::  afactor           !<
3308       REAL(wp)     ::  beta_act          !<
3309       REAL(wp)     ::  bfactor           !<
3310       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3311       REAL(wp)     ::  k_act             !<
3312       REAL(wp)     ::  n_act             !<
3313       REAL(wp)     ::  n_ccn             !<
3314       REAL(wp)     ::  s_0               !<
3315       REAL(wp)     ::  sat_max           !<
3316       REAL(wp)     ::  sigma             !<
3317       REAL(wp)     ::  sigma_act         !<
3318
3319       DO  k = nzb+1, nzt
3320!
3321!--       Predetermine flag to mask topography
3322          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3323!
3324!--       Call calculation of supersaturation 
3325          CALL supersaturation ( i, j, k )
3326!
3327!--       Prescribe parameters for activation
3328!--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
3329          k_act  = 0.7_wp
3330          activ  = 0.0_wp
3331
3332          IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
3333!
3334!--          Compute the number of activated Aerosols
3335!--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
3336             n_act     = na_init * sat**k_act
3337!
3338!--          Compute the number of cloud droplets
3339!--          (see: Morrison + Grabowski, 2007, JAS, 64)
3340!            activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
3341
3342!
3343!--          Compute activation rate after Khairoutdinov and Kogan
3344!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
3345             sat_max = 0.8_wp / 100.0_wp
3346             activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN             &
3347                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /         &
3348                  dt_micro
3349
3350             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
3351          ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk )  THEN
3352!
3353!--          Curvature effect (afactor) with surface tension
3354!--          parameterization by Straka (2009)
3355             sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
3356             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
3357!
3358!--          Solute effect (bfactor)
3359             bfactor = vanthoff * molecular_weight_of_water *                  &
3360                  rho_s / ( molecular_weight_of_solute * rho_l )
3361
3362!
3363!--          Prescribe power index that describes the soluble fraction
3364!--          of an aerosol particle (beta).
3365!--          (see: Morrison + Grabowski, 2007, JAS, 64)
3366             beta_act  = 0.5_wp
3367             sigma_act = sigma_bulk**( 1.0_wp + beta_act )
3368!
3369!--          Calculate mean geometric supersaturation (s_0) with
3370!--          parameterization by Khvorostyanov and Curry (2006)
3371             s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *         &
3372               ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
3373
3374!
3375!--          Calculate number of activated CCN as a function of
3376!--          supersaturation and taking Koehler theory into account
3377!--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
3378             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                    &
3379                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
3380             activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
3381
3382             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)
3383          ENDIF
3384
3385       ENDDO
3386
3387    END SUBROUTINE activation_ij
3388
3389!------------------------------------------------------------------------------!
3390! Description:
3391! ------------
3392!> Calculate condensation rate for cloud water content (after Khairoutdinov and
3393!> Kogan, 2000).
3394!------------------------------------------------------------------------------!
3395    SUBROUTINE condensation_ij( i, j )
3396
3397       IMPLICIT NONE
3398
3399       INTEGER(iwp) ::  i                 !<
3400       INTEGER(iwp) ::  j                 !<
3401       INTEGER(iwp) ::  k                 !<
3402
3403       REAL(wp)     ::  cond              !<
3404       REAL(wp)     ::  cond_max          !<
3405       REAL(wp)     ::  dc                !<
3406       REAL(wp)     ::  evap              !<
3407       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3408       REAL(wp)     ::  g_fac             !<
3409       REAL(wp)     ::  nc_0              !<
3410       REAL(wp)     ::  temp              !<
3411       REAL(wp)     ::  xc                !<
3412
3413
3414       DO  k = nzb+1, nzt
3415!
3416!--       Predetermine flag to mask topography
3417          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3418!
3419!--       Call calculation of supersaturation 
3420          CALL supersaturation ( i, j, k )
3421!
3422!--       Actual temperature:
3423          temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
3424
3425          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
3426                              l_v / ( thermal_conductivity_l * temp )    &
3427                              + r_v * temp / ( diff_coeff_l * e_s )      &
3428                            )
3429!
3430!--       Mean weight of cloud drops
3431          IF ( nc(k,j,i) <= 0.0_wp) CYCLE
3432          xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
3433!
3434!--       Weight averaged diameter of cloud drops:
3435          dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
3436!
3437!--       Integral diameter of cloud drops
3438          nc_0 = nc(k,j,i) * dc
3439!
3440!--       Condensation needs only to be calculated in supersaturated regions
3441          IF ( sat > 0.0_wp )  THEN
3442!
3443!--          Condensation rate of cloud water content
3444!--          after KK scheme.
3445!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
3446             cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
3447             cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
3448             cond      = MIN( cond, cond_max / dt_micro )
3449
3450             qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
3451          ELSEIF ( sat < 0.0_wp ) THEN
3452             evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
3453             evap      = MAX( evap, -qc(k,j,i) / dt_micro )
3454
3455             qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag
3456          ENDIF
3457       ENDDO
3458
3459    END SUBROUTINE condensation_ij
3460
3461
3462!------------------------------------------------------------------------------!
3463! Description:
3464! ------------
3465!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
3466!------------------------------------------------------------------------------!
3467    SUBROUTINE autoconversion_ij( i, j )
3468
3469       IMPLICIT NONE
3470
3471       INTEGER(iwp) ::  i                 !<
3472       INTEGER(iwp) ::  j                 !<
3473       INTEGER(iwp) ::  k                 !<
3474
3475       REAL(wp)     ::  alpha_cc          !<
3476       REAL(wp)     ::  autocon           !<
3477       REAL(wp)     ::  dissipation       !<
3478       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3479       REAL(wp)     ::  k_au              !<
3480       REAL(wp)     ::  l_mix             !<
3481       REAL(wp)     ::  nc_auto           !<
3482       REAL(wp)     ::  nu_c              !<
3483       REAL(wp)     ::  phi_au            !<
3484       REAL(wp)     ::  r_cc              !<
3485       REAL(wp)     ::  rc                !<
3486       REAL(wp)     ::  re_lambda         !<
3487       REAL(wp)     ::  sigma_cc          !<
3488       REAL(wp)     ::  tau_cloud         !<
3489       REAL(wp)     ::  xc                !<
3490
3491       DO  k = nzb+1, nzt
3492!
3493!--       Predetermine flag to mask topography
3494          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3495          IF ( microphysics_morrison ) THEN
3496             nc_auto = nc(k,j,i)
3497          ELSE
3498             nc_auto = nc_const
3499          ENDIF
3500
3501          IF ( qc(k,j,i) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
3502
3503             k_au = k_cc / ( 20.0_wp * x0 )
3504!
3505!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
3506!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
3507             tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ),     &
3508                              0.0_wp )
3509!
3510!--          Universal function for autoconversion process
3511!--          (Seifert and Beheng, 2006):
3512             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
3513!
3514!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
3515!--          (Use constant nu_c = 1.0_wp instead?)
3516             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
3517!
3518!--          Mean weight of cloud droplets:
3519             xc = hyrho(k) * qc(k,j,i) / nc_auto
3520!
3521!--          Parameterized turbulence effects on autoconversion (Seifert,
3522!--          Nuijens and Stevens, 2010)
3523             IF ( collision_turbulence )  THEN
3524!
3525!--             Weight averaged radius of cloud droplets:
3526                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
3527
3528                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
3529                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
3530                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
3531!
3532!--             Mixing length (neglecting distance to ground and stratification)
3533                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
3534!
3535!--             Limit dissipation rate according to Seifert, Nuijens and
3536!--             Stevens (2010)
3537                dissipation = MIN( 0.06_wp, diss(k,j,i) )
3538!
3539!--             Compute Taylor-microscale Reynolds number:
3540                re_lambda = 6.0_wp / 11.0_wp *                                 &
3541                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
3542                            SQRT( 15.0_wp / kin_vis_air ) *                    &
3543                            dissipation**( 1.0_wp / 6.0_wp )
3544!
3545!--             The factor of 1.0E4 is needed to convert the dissipation rate
3546!--             from m2 s-3 to cm2 s-3.
3547                k_au = k_au * ( 1.0_wp +                                       &
3548                                dissipation * 1.0E4_wp *                       &
3549                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
3550                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
3551                                                  sigma_cc )**2                &
3552                                                ) + beta_cc                    &
3553                                )                                              &
3554                              )
3555             ENDIF
3556!
3557!--          Autoconversion rate (Seifert and Beheng, 2006):
3558             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
3559                       ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *            &
3560                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
3561                       rho_surface
3562             autocon = MIN( autocon, qc(k,j,i) / dt_micro )
3563
3564             qr(k,j,i) = qr(k,j,i) + autocon * dt_micro                 * flag
3565             qc(k,j,i) = qc(k,j,i) - autocon * dt_micro                 * flag
3566             nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag
3567             IF ( microphysics_morrison ) THEN
3568                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *                &
3569                                    autocon / x0 * hyrho(k) * dt_micro * flag )
3570             ENDIF
3571
3572          ENDIF
3573
3574       ENDDO
3575
3576    END SUBROUTINE autoconversion_ij
3577
3578!------------------------------------------------------------------------------!
3579! Description:
3580! ------------
3581!> Autoconversion process (Kessler, 1969).
3582!------------------------------------------------------------------------------!
3583    SUBROUTINE autoconversion_kessler_ij( i, j )
3584
3585
3586       IMPLICIT NONE
3587
3588       INTEGER(iwp) ::  i      !<
3589       INTEGER(iwp) ::  j      !<
3590       INTEGER(iwp) ::  k      !<
3591       INTEGER(iwp) ::  k_wall !< topography top index
3592
3593       REAL(wp)    ::  dqdt_precip !<
3594       REAL(wp)    ::  flag              !< flag to indicate first grid level above surface
3595
3596!
3597!--    Determine vertical index of topography top
3598       k_wall = get_topography_top_index_ji( j, i, 's' )
3599       DO  k = nzb+1, nzt
3600!
3601!--       Predetermine flag to mask topography
3602          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3603
3604          IF ( qc(k,j,i) > ql_crit )  THEN
3605             dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
3606          ELSE
3607             dqdt_precip = 0.0_wp
3608          ENDIF
3609
3610          qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
3611          q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
3612          pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k)  &
3613                                                                 * flag
3614
3615!
3616!--       Compute the rain rate (stored on surface grid point)
3617          prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
3618
3619       ENDDO
3620
3621    END SUBROUTINE autoconversion_kessler_ij
3622
3623!------------------------------------------------------------------------------!
3624! Description:
3625! ------------
3626!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
3627!------------------------------------------------------------------------------!
3628    SUBROUTINE accretion_ij( i, j )
3629
3630       IMPLICIT NONE
3631
3632       INTEGER(iwp) ::  i                 !<
3633       INTEGER(iwp) ::  j                 !<
3634       INTEGER(iwp) ::  k                 !<
3635
3636       REAL(wp)     ::  accr              !<
3637       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3638       REAL(wp)     ::  k_cr              !<
3639       REAL(wp)     ::  nc_accr           !<
3640       REAL(wp)     ::  phi_ac            !<
3641       REAL(wp)     ::  tau_cloud         !<
3642       REAL(wp)     ::  xc                !<
3643
3644
3645       DO  k = nzb+1, nzt
3646!
3647!--       Predetermine flag to mask topography
3648          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3649          IF ( microphysics_morrison ) THEN
3650             nc_accr = nc(k,j,i)
3651          ELSE
3652             nc_accr = nc_const
3653          ENDIF
3654
3655          IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )  .AND.  &
3656               ( nc_accr > eps_mr ) )  THEN
3657!
3658!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
3659             tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
3660!
3661!--          Universal function for accretion process
3662!--          (Seifert and Beheng, 2001):
3663             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
3664
3665!
3666!--          Mean weight of cloud drops
3667             xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
3668!
3669!--          Parameterized turbulence effects on autoconversion (Seifert,
3670!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
3671!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
3672             IF ( collision_turbulence )  THEN
3673                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
3674                                 MIN( 600.0_wp,                                &
3675                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
3676                               )
3677             ELSE
3678                k_cr = k_cr0
3679             ENDIF
3680!
3681!--          Accretion rate (Seifert and Beheng, 2006):
3682             accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *                      &
3683                    SQRT( rho_surface * hyrho(k) )
3684             accr = MIN( accr, qc(k,j,i) / dt_micro )
3685
3686             qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
3687             qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
3688             IF ( microphysics_morrison )  THEN
3689                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc *               &
3690                                             hyrho(k) * dt_micro * flag        &
3691                                         )
3692             ENDIF
3693
3694
3695          ENDIF
3696
3697       ENDDO
3698
3699    END SUBROUTINE accretion_ij
3700
3701
3702!------------------------------------------------------------------------------!
3703! Description:
3704! ------------
3705!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
3706!------------------------------------------------------------------------------!
3707    SUBROUTINE selfcollection_breakup_ij( i, j )
3708
3709       IMPLICIT NONE
3710
3711       INTEGER(iwp) ::  i                 !<
3712       INTEGER(iwp) ::  j                 !<
3713       INTEGER(iwp) ::  k                 !<
3714
3715       REAL(wp)     ::  breakup           !<
3716       REAL(wp)     ::  dr                !<
3717       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3718       REAL(wp)     ::  phi_br            !<
3719       REAL(wp)     ::  selfcoll          !<
3720
3721       DO  k = nzb+1, nzt
3722!
3723!--       Predetermine flag to mask topography
3724          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3725
3726          IF ( qr(k,j,i) > eps_sb )  THEN
3727!
3728!--          Selfcollection rate (Seifert and Beheng, 2001):
3729             selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * SQRT( hyrho(k) * rho_surface )
3730!
3731!--          Weight averaged diameter of rain drops:
3732             dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
3733!
3734!--          Collisional breakup rate (Seifert, 2008):
3735             IF ( dr >= 0.3E-3_wp )  THEN
3736                phi_br  = k_br * ( dr - 1.1E-3_wp )
3737                breakup = selfcoll * ( phi_br + 1.0_wp )
3738             ELSE
3739                breakup = 0.0_wp
3740             ENDIF
3741
3742             selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
3743             nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
3744
3745          ENDIF
3746       ENDDO
3747
3748    END SUBROUTINE selfcollection_breakup_ij
3749
3750
3751!------------------------------------------------------------------------------!
3752! Description:
3753! ------------
3754!> Evaporation of precipitable water. Condensation is neglected for
3755!> precipitable water. Call for grid point i,j
3756!------------------------------------------------------------------------------!
3757    SUBROUTINE evaporation_rain_ij( i, j )
3758
3759       IMPLICIT NONE
3760
3761       INTEGER(iwp) ::  i                 !<
3762       INTEGER(iwp) ::  j                 !<
3763       INTEGER(iwp) ::  k                 !<
3764
3765       REAL(wp)     ::  dr                !<
3766       REAL(wp)     ::  evap              !<
3767       REAL(wp)     ::  evap_nr           !<
3768       REAL(wp)     ::  f_vent            !<
3769       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3770       REAL(wp)     ::  g_evap            !<
3771       REAL(wp)     ::  lambda_r          !<
3772       REAL(wp)     ::  mu_r              !<
3773       REAL(wp)     ::  mu_r_2            !<
3774       REAL(wp)     ::  mu_r_5d2          !<
3775       REAL(wp)     ::  nr_0              !<
3776       REAL(wp)     ::  temp              !<
3777       REAL(wp)     ::  xr                !<
3778
3779       DO  k = nzb+1, nzt
3780!
3781!--       Predetermine flag to mask topography
3782          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3783
3784          IF ( qr(k,j,i) > eps_sb )  THEN
3785!
3786!--          Call calculation of supersaturation 
3787             CALL supersaturation ( i, j, k )
3788!
3789!--          Evaporation needs only to be calculated in subsaturated regions
3790             IF ( sat < 0.0_wp )  THEN
3791!
3792!--             Actual temperature:
3793                temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
3794
3795                g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /  &
3796                                    ( thermal_conductivity_l * temp ) +        &
3797                                    r_v * temp / ( diff_coeff_l * e_s )        &
3798                                  )
3799!
3800!--             Mean weight of rain drops
3801                xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
3802!
3803!--             Weight averaged diameter of rain drops:
3804                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
3805!
3806!--             Compute ventilation factor and intercept parameter
3807!--             (Seifert and Beheng, 2006; Seifert, 2008):
3808                IF ( ventilation_effect )  THEN
3809!
3810!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
3811!--                Stevens and Seifert, 2008):
3812                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
3813!
3814!--                Slope parameter of gamma distribution (Seifert, 2008):
3815                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
3816                                ( mu_r + 1.0_wp )                              &
3817                              )**( 1.0_wp / 3.0_wp ) / dr
3818
3819                   mu_r_2   = mu_r + 2.0_wp
3820                   mu_r_5d2 = mu_r + 2.5_wp
3821
3822                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  &
3823                            b_vent * schmidt_p_1d3 *                           &
3824                            SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
3825                            lambda_r**( -mu_r_5d2 ) *                          &
3826                            ( 1.0_wp -                                         &
3827                              0.5_wp * ( b_term / a_term ) *                   &
3828                              ( lambda_r / ( c_term + lambda_r )               &
3829                              )**mu_r_5d2 -                                    &
3830                              0.125_wp * ( b_term / a_term )**2 *              &
3831                              ( lambda_r / ( 2.0_wp * c_term + lambda_r )      &
3832                              )**mu_r_5d2 -                                    &
3833                              0.0625_wp * ( b_term / a_term )**3 *             &
3834                              ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
3835                              )**mu_r_5d2 -                                    &
3836                              0.0390625_wp * ( b_term / a_term )**4 *          &
3837                              ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
3838                              )**mu_r_5d2                                      &
3839                            )
3840
3841                   nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /           &
3842                            gamm( mu_r + 1.0_wp )
3843                ELSE
3844                   f_vent = 1.0_wp
3845                   nr_0   = nr(k,j,i) * dr
3846                ENDIF
3847!
3848!--             Evaporation rate of rain water content (Seifert and Beheng, 2006):
3849                evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
3850                evap    = MAX( evap, -qr(k,j,i) / dt_micro )
3851                evap_nr = MAX( c_evap * evap / xr * hyrho(k),                  &
3852                               -nr(k,j,i) / dt_micro )
3853
3854                qr(k,j,i) = qr(k,j,i) + evap    * dt_micro * flag
3855                nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
3856
3857             ENDIF
3858          ENDIF
3859
3860       ENDDO
3861
3862    END SUBROUTINE evaporation_rain_ij
3863
3864
3865!------------------------------------------------------------------------------!
3866! Description:
3867! ------------
3868!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
3869!> Call for grid point i,j
3870!------------------------------------------------------------------------------!
3871    SUBROUTINE sedimentation_cloud_ij( i, j )
3872
3873       IMPLICIT NONE
3874
3875       INTEGER(iwp) ::  i             !<
3876       INTEGER(iwp) ::  j             !<
3877       INTEGER(iwp) ::  k             !<
3878
3879       REAL(wp)     ::  flag    !< flag to indicate first grid level above surface
3880       REAL(wp)     ::  nc_sedi !<
3881
3882       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc  !<
3883       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc  !<
3884
3885       sed_qc(nzt+1) = 0.0_wp
3886       sed_nc(nzt+1) = 0.0_wp
3887
3888
3889       DO  k = nzt, nzb+1, -1
3890!
3891!--       Predetermine flag to mask topography
3892          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3893          IF ( microphysics_morrison ) THEN
3894             nc_sedi = nc(k,j,i)
3895          ELSE
3896             nc_sedi = nc_const
3897          ENDIF
3898!
3899!--       Sedimentation fluxes for number concentration are only calculated
3900!--       for cloud_scheme = 'morrison'
3901          IF ( microphysics_morrison ) THEN
3902             IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
3903                sed_nc(k) = sed_qc_const *                                     &
3904                            ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *     &
3905                            ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
3906             ELSE
3907                sed_nc(k) = 0.0_wp
3908             ENDIF
3909
3910             sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *                 &
3911                              nc(k,j,i) / dt_micro + sed_nc(k+1)                &
3912                            ) * flag
3913
3914             nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *               &
3915                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
3916          ENDIF
3917
3918          IF ( qc(k,j,i) > eps_sb  .AND.  nc_sedi > eps_mr )  THEN
3919             sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *        &
3920                         ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )  * flag
3921          ELSE
3922             sed_qc(k) = 0.0_wp
3923          ENDIF
3924
3925          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /          &
3926                                      dt_micro + sed_qc(k+1)                   &
3927                         ) * flag
3928
3929          q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
3930                                hyrho(k) * dt_micro * flag
3931          qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
3932                                hyrho(k) * dt_micro * flag
3933          pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
3934                                hyrho(k) * lv_d_cp * d_exner(k) * dt_micro     &
3935                                                                * flag
3936
3937!
3938!--       Compute the precipitation rate of cloud (fog) droplets
3939          IF ( call_microphysics_at_all_substeps )  THEN
3940             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *                  &
3941                              weight_substep(intermediate_timestep_count) * flag
3942          ELSE
3943             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
3944          ENDIF
3945
3946       ENDDO
3947
3948    END SUBROUTINE sedimentation_cloud_ij
3949
3950
3951!------------------------------------------------------------------------------!
3952! Description:
3953! ------------
3954!> Computation of sedimentation flux. Implementation according to Stevens
3955!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
3956!------------------------------------------------------------------------------!
3957    SUBROUTINE sedimentation_rain_ij( i, j )
3958
3959       IMPLICIT NONE
3960
3961       INTEGER(iwp) ::  i             !< running index x direction
3962       INTEGER(iwp) ::  j             !< running index y direction
3963       INTEGER(iwp) ::  k             !< running index z direction
3964       INTEGER(iwp) ::  k_run         !<
3965       INTEGER(iwp) ::  m             !< running index surface elements
3966       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
3967       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
3968
3969       REAL(wp)     ::  c_run                      !<
3970       REAL(wp)     ::  d_max                      !<
3971       REAL(wp)     ::  d_mean                     !<
3972       REAL(wp)     ::  d_min                      !<
3973       REAL(wp)     ::  dr                         !<
3974       REAL(wp)     ::  flux                       !<
3975       REAL(wp)     ::  flag                       !< flag to indicate first grid level above surface
3976       REAL(wp)     ::  lambda_r                   !<
3977       REAL(wp)     ::  mu_r                       !<
3978       REAL(wp)     ::  z_run                      !<
3979
3980       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
3981       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
3982       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
3983       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
3984       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
3985       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
3986       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
3987       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
3988
3989!
3990!--    Compute velocities
3991       DO  k = nzb+1, nzt
3992!
3993!--       Predetermine flag to mask topography
3994          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3995
3996          IF ( qr(k,j,i) > eps_sb )  THEN
3997!
3998!--          Weight averaged diameter of rain drops:
3999             dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
4000!
4001!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
4002!--          Stevens and Seifert, 2008):
4003             mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
4004!
4005!--          Slope parameter of gamma distribution (Seifert, 2008):
4006             lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *              &
4007                          ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
4008
4009             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
4010                                         a_term - b_term * ( 1.0_wp +          &
4011                                            c_term / lambda_r )**( -1.0_wp *   &
4012                                            ( mu_r + 1.0_wp ) )                &
4013                                        )                                      &
4014                          ) * flag
4015             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
4016                                         a_term - b_term * ( 1.0_wp +          &
4017                                            c_term / lambda_r )**( -1.0_wp *   &
4018                                            ( mu_r + 4.0_wp ) )                &
4019                                       )                                       &
4020                          ) * flag
4021          ELSE
4022             w_nr(k) = 0.0_wp
4023             w_qr(k) = 0.0_wp
4024          ENDIF
4025       ENDDO
4026!
4027!--    Adjust boundary values using surface data type.
4028!--    Upward facing non-natural
4029       surf_s = bc_h(0)%start_index(j,i)
4030       surf_e = bc_h(0)%end_index(j,i)
4031       DO  m = surf_s, surf_e
4032          k         = bc_h(0)%k(m)
4033          w_nr(k-1) = w_nr(k)
4034          w_qr(k-1) = w_qr(k)
4035       ENDDO
4036!
4037!--    Downward facing non-natural
4038       surf_s = bc_h(1)%start_index(j,i)
4039       surf_e = bc_h(1)%end_index(j,i)
4040       DO  m = surf_s, surf_e
4041          k         = bc_h(1)%k(m)
4042          w_nr(k+1) = w_nr(k)
4043          w_qr(k+1) = w_qr(k)
4044       ENDDO
4045!
4046!--    Neumann boundary condition at model top
4047       w_nr(nzt+1) = 0.0_wp
4048       w_qr(nzt+1) = 0.0_wp
4049!
4050!--    Compute Courant number
4051       DO  k = nzb+1, nzt
4052!
4053!--       Predetermine flag to mask topography
4054          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4055
4056          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
4057                    dt_micro * ddzu(k) * flag
4058          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
4059                    dt_micro * ddzu(k) * flag
4060       ENDDO
4061!
4062!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
4063       IF ( limiter_sedimentation )  THEN
4064
4065          DO k = nzb+1, nzt
4066!
4067!--          Predetermine flag to mask topography
4068             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4069
4070             d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
4071             d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
4072             d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
4073
4074             qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
4075                                                        2.0_wp * d_max,        &
4076                                                        ABS( d_mean ) ) * flag
4077
4078             d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
4079             d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
4080             d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
4081
4082             nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
4083                                                        2.0_wp * d_max,        &
4084                                                        ABS( d_mean ) ) * flag
4085          ENDDO
4086
4087       ELSE
4088
4089          nr_slope = 0.0_wp
4090          qr_slope = 0.0_wp
4091
4092       ENDIF
4093
4094       sed_nr(nzt+1) = 0.0_wp
4095       sed_qr(nzt+1) = 0.0_wp
4096!
4097!--    Compute sedimentation flux
4098       DO  k = nzt, nzb+1, -1
4099!
4100!--       Predetermine flag to mask topography
4101          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4102!
4103!--       Sum up all rain drop number densities which contribute to the flux
4104!--       through k-1/2
4105          flux  = 0.0_wp
4106          z_run = 0.0_wp ! height above z(k)
4107          k_run = k
4108          c_run = MIN( 1.0_wp, c_nr(k) )
4109          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
4110             flux  = flux + hyrho(k_run) *                                     &
4111                     ( nr(k_run,j,i) + nr_slope(k_run) * ( 1.0_wp - c_run ) *  &
4112                     0.5_wp ) * c_run * dzu(k_run) * flag
4113             z_run = z_run + dzu(k_run)            * flag
4114             k_run = k_run + 1                     * flag
4115             c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag
4116          ENDDO
4117!
4118!--       It is not allowed to sediment more rain drop number density than
4119!--       available
4120          flux = MIN( flux,                                                    &
4121                      hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * dt_micro )
4122
4123          sed_nr(k) = flux / dt_micro * flag
4124          nr(k,j,i)  = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
4125                                    hyrho(k) * dt_micro * flag
4126!
4127!--       Sum up all rain water content which contributes to the flux
4128!--       through k-1/2
4129          flux  = 0.0_wp
4130          z_run = 0.0_wp ! height above z(k)
4131          k_run = k
4132          c_run = MIN( 1.0_wp, c_qr(k) )
4133
4134          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
4135
4136             flux  = flux + hyrho(k_run) *                                     &
4137                     ( qr(k_run,j,i) + qr_slope(k_run) * ( 1.0_wp - c_run ) *  &
4138                     0.5_wp ) * c_run * dzu(k_run) * flag
4139             z_run = z_run + dzu(k_run)            * flag
4140             k_run = k_run + 1                     * flag
4141             c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag
4142
4143          ENDDO
4144!
4145!--       It is not allowed to sediment more rain water content than available
4146          flux = MIN( flux,                                                    &
4147                      hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * dt_micro )
4148
4149          sed_qr(k) = flux / dt_micro * flag
4150
4151          qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
4152                                hyrho(k) * dt_micro * flag
4153          q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
4154                                hyrho(k) * dt_micro * flag
4155          pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
4156                                hyrho(k) * lv_d_cp * d_exner(k) * dt_micro     &
4157                                                                * flag
4158!
4159!--       Compute the rain rate
4160          IF ( call_microphysics_at_all_substeps )  THEN
4161             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)                    &
4162                          * weight_substep(intermediate_timestep_count) * flag
4163          ELSE
4164             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
4165          ENDIF
4166
4167       ENDDO
4168
4169    END SUBROUTINE sedimentation_rain_ij
4170
4171
4172!------------------------------------------------------------------------------!
4173! Description:
4174! ------------
4175!> This subroutine computes the precipitation amount due to gravitational
4176!> settling of rain and cloud (fog) droplets
4177!------------------------------------------------------------------------------!
4178    SUBROUTINE calc_precipitation_amount_ij( i, j )
4179
4180       IMPLICIT NONE
4181
4182       INTEGER(iwp) ::  i             !< running index x direction
4183       INTEGER(iwp) ::  j             !< running index y direction
4184       INTEGER(iwp) ::  k             !< running index z direction
4185       INTEGER(iwp) ::  m             !< running index surface elements
4186       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
4187       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
4188
4189       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
4190            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
4191            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
4192       THEN
4193
4194          surf_s = bc_h(0)%start_index(j,i)
4195          surf_e = bc_h(0)%end_index(j,i)
4196          DO  m = surf_s, surf_e
4197             k                         = bc_h(0)%k(m)
4198             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
4199                                               prr(k,j,i) * hyrho(k) * dt_3d
4200          ENDDO
4201
4202       ENDIF
4203
4204    END SUBROUTINE calc_precipitation_amount_ij
4205
4206
4207!------------------------------------------------------------------------------!
4208! Description:
4209! ------------
4210!> Computation of the diagnostic supersaturation sat, actual liquid water
4211!< temperature t_l and saturation water vapor mixing ratio q_s
4212!------------------------------------------------------------------------------!
4213    SUBROUTINE supersaturation ( i,j,k )
4214
4215       IMPLICIT NONE
4216
4217       INTEGER(iwp) ::  i   !< running index
4218       INTEGER(iwp) ::  j   !< running index
4219       INTEGER(iwp) ::  k   !< running index
4220
4221       REAL(wp) ::  alpha   !< correction factor
4222!
4223!--    Actual liquid water temperature:
4224       t_l = exner(k) * pt(k,j,i)
4225!
4226!--    Calculate water vapor saturation pressure
4227       e_s = magnus( t_l )
4228!
4229!--    Computation of saturation mixing ratio:
4230       q_s   = rd_d_rv * e_s / ( hyp(k) - e_s )
4231!
4232!--    Correction factor
4233       alpha = rd_d_rv * lv_d_rd * lv_d_cp / ( t_l * t_l )
4234!
4235!--    Correction of the approximated value
4236!--    (see: Cuijpers + Duynkerke, 1993, JAS, 23)
4237       q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / ( 1.0_wp + alpha * q_s )
4238!
4239!--    Supersaturation:
4240!--    Not in case of microphysics_kessler or microphysics_sat_adjust
4241!--    since qr is unallocated
4242       IF ( .NOT. microphysics_kessler  .AND.                                  &
4243            .NOT.  microphysics_sat_adjust )  THEN
4244          sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
4245       ENDIF
4246
4247    END SUBROUTINE supersaturation
4248
4249
4250!------------------------------------------------------------------------------!
4251! Description:
4252! ------------
4253!> Calculation of the liquid water content (0%-or-100%-scheme). This scheme is
4254!> used by the one and the two moment cloud physics scheme. Using the two moment
4255!> scheme, this calculation results in the cloud water content.
4256!------------------------------------------------------------------------------!
4257    SUBROUTINE calc_liquid_water_content
4258
4259
4260
4261       IMPLICIT NONE
4262
4263       INTEGER(iwp) ::  i !<
4264       INTEGER(iwp) ::  j !<
4265       INTEGER(iwp) ::  k !<
4266
4267       REAL(wp)     ::  flag !< flag to indicate first grid level above surface
4268
4269
4270       DO  i = nxlg, nxrg
4271          DO  j = nysg, nyng
4272             DO  k = nzb+1, nzt
4273!
4274!--             Predetermine flag to mask topography
4275                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4276
4277   !
4278   !--          Call calculation of supersaturation located
4279                CALL supersaturation( i, j, k )
4280
4281   !
4282   !--          Compute the liquid water content
4283                IF ( microphysics_seifert  .AND. .NOT. microphysics_morrison ) &
4284                THEN
4285                   IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp )  THEN
4286                      qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) ) * flag
4287                      ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag
4288                   ELSE
4289                      IF ( q(k,j,i) < qr(k,j,i) )  q(k,j,i) = qr(k,j,i)
4290                      qc(k,j,i) = 0.0_wp
4291                      ql(k,j,i) = qr(k,j,i) * flag
4292                   ENDIF
4293                ELSEIF ( microphysics_morrison )  THEN
4294                   ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag
4295                ELSE
4296                   IF ( ( q(k,j,i) - q_s ) > 0.0_wp )  THEN
4297                      qc(k,j,i) = ( q(k,j,i) - q_s ) * flag
4298                      ql(k,j,i) = qc(k,j,i) * flag
4299                   ELSE
4300                      qc(k,j,i) = 0.0_wp
4301                      ql(k,j,i) = 0.0_wp
4302                   ENDIF
4303                ENDIF
4304             ENDDO
4305          ENDDO
4306       ENDDO
4307
4308    END SUBROUTINE calc_liquid_water_content
4309
4310!------------------------------------------------------------------------------!
4311! Description:
4312! ------------
4313!> This function computes the gamma function (Press et al., 1992).
4314!> The gamma function is needed for the calculation of the evaporation
4315!> of rain drops.
4316!------------------------------------------------------------------------------!
4317    FUNCTION gamm( xx )
4318
4319       IMPLICIT NONE
4320
4321       INTEGER(iwp) ::  j            !<
4322
4323       REAL(wp)     ::  gamm         !<
4324       REAL(wp)     ::  ser          !<
4325       REAL(wp)     ::  tmp          !<
4326       REAL(wp)     ::  x_gamm       !<
4327       REAL(wp)     ::  xx           !<
4328       REAL(wp)     ::  y_gamm       !<
4329
4330
4331       REAL(wp), PARAMETER  ::  stp = 2.5066282746310005_wp               !<
4332       REAL(wp), PARAMETER  ::  cof(6) = (/ 76.18009172947146_wp,      &
4333                                           -86.50532032941677_wp,      &
4334                                            24.01409824083091_wp,      &
4335                                            -1.231739572450155_wp,     &
4336                                             0.1208650973866179E-2_wp, &
4337                                            -0.5395239384953E-5_wp /)     !<
4338
4339       x_gamm = xx
4340       y_gamm = x_gamm
4341       tmp = x_gamm + 5.5_wp
4342       tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
4343       ser = 1.000000000190015_wp
4344
4345       DO  j = 1, 6
4346          y_gamm = y_gamm + 1.0_wp
4347          ser    = ser + cof( j ) / y_gamm
4348       ENDDO
4349
4350!
4351!--    Until this point the algorithm computes the logarithm of the gamma
4352!--    function. Hence, the exponential function is used.
4353!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
4354       gamm = EXP( tmp ) * stp * ser / x_gamm
4355
4356       RETURN
4357
4358    END FUNCTION gamm
4359
4360 END MODULE bulk_cloud_model_mod
Note: See TracBrowser for help on using the repository browser.