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

Last change on this file since 3348 was 3274, checked in by knoop, 6 years ago

Modularization of all bulk cloud physics code components

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