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

Last change on this file since 2749 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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