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

Last change on this file since 2031 was 2031, checked in by knoop, 7 years ago

Renamed variable rho to rho_ocean, rho_init to rho_ocean_init and rho_av to rho_ocean_av

  • Property svn:keywords set to Id
File size: 86.8 KB
RevLine 
[1850]1!> @file microphysics_mod.f90
[2000]2!------------------------------------------------------------------------------!
[1093]3! This file is part of PALM.
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!
[1818]17! Copyright 1997-2016 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1093]19!
[1000]20! Current revisions:
[1092]21! ------------------
[2031]22! renamed variable rho to rho_ocean
[1851]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: microphysics_mod.f90 2031 2016-10-21 15:11:58Z knoop $
27!
[2001]28! 2000 2016-08-20 18:09:15Z knoop
29! Forced header and separation lines into 80 columns
30!
[1851]31! 1850 2016-04-08 13:29:27Z maronga
32! Module renamed
33! Adapted for modularization of microphysics.
34!
[1846]35! 1845 2016-04-08 08:29:13Z raasch
36! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
37! point (instead of one point above surface)
38!
[1832]39! 1831 2016-04-07 13:15:51Z hoffmann
40! turbulence renamed collision_turbulence,
41! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
42! avaialble for microphysics_kessler.
43!
[1823]44! 1822 2016-04-07 07:49:42Z hoffmann
45! Unused variables removed.
46! Kessler scheme integrated.
47!
[1692]48! 1691 2015-10-26 16:17:44Z maronga
49! Added new routine calc_precipitation_amount. The routine now allows to account
50! for precipitation due to sedimenation of cloud (fog) droplets
51!
[1683]52! 1682 2015-10-07 23:56:08Z knoop
53! Code annotations made doxygen readable
54!
[1647]55! 1646 2015-09-02 16:00:10Z hoffmann
56! Bugfix: Wrong computation of d_mean.
57!
[1362]58! 1361 2014-04-16 15:17:48Z hoffmann
59! Bugfix in sedimentation_rain: Index corrected.
60! Vectorized version of adjust_cloud added.
61! Little reformatting of the code.
62!
[1354]63! 1353 2014-04-08 15:21:23Z heinze
64! REAL constants provided with KIND-attribute
65!
[1347]66! 1346 2014-03-27 13:18:20Z heinze
67! Bugfix: REAL constants provided with KIND-attribute especially in call of
68! intrinsic function like MAX, MIN, SIGN
69!
[1335]70! 1334 2014-03-25 12:21:40Z heinze
71! Bugfix: REAL constants provided with KIND-attribute
72!
[1323]73! 1322 2014-03-20 16:38:49Z raasch
74! REAL constants defined as wp-kind
75!
[1321]76! 1320 2014-03-20 08:40:49Z raasch
[1320]77! ONLY-attribute added to USE-statements,
78! kind-parameters added to all INTEGER and REAL declaration statements,
79! kinds are defined in new module kinds,
80! comment fields (!:) to be used for variable explanations added to
81! all variable declaration statements
[1000]82!
[1242]83! 1241 2013-10-30 11:36:58Z heinze
[2031]84! hyp and rho_ocean have to be calculated at each time step if data from external
[1242]85! file LSF_DATA are used
86!
[1116]87! 1115 2013-03-26 18:16:16Z hoffmann
88! microphyical tendencies are calculated in microphysics_control in an optimized
89! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
90!
[1107]91! 1106 2013-03-04 05:31:38Z raasch
92! small changes in code formatting
93!
[1093]94! 1092 2013-02-02 11:24:22Z raasch
95! unused variables removed
96! file put under GPL
97!
[1066]98! 1065 2012-11-22 17:42:36Z hoffmann
99! Sedimentation process implemented according to Stevens and Seifert (2008).
[1115]100! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
[1066]101! and Stevens, 2010).
102!
[1054]103! 1053 2012-11-13 17:11:03Z hoffmann
104! initial revision
[1000]105!
106! Description:
107! ------------
[1849]108!> Calculate bilk cloud microphysics.
[1000]109!------------------------------------------------------------------------------!
[1682]110 MODULE microphysics_mod
[1000]111
[1849]112    USE  kinds
113
114    IMPLICIT NONE
115
116    LOGICAL ::  cloud_water_sedimentation = .FALSE.  !< cloud water sedimentation
117    LOGICAL ::  limiter_sedimentation = .TRUE.       !< sedimentation limiter
118    LOGICAL ::  collision_turbulence = .FALSE.       !< turbulence effects
119    LOGICAL ::  ventilation_effect = .TRUE.          !< ventilation effect
120
121    REAL(wp) ::  a_1 = 8.69E-4_wp          !< coef. in turb. parametrization (cm-2 s3)
122    REAL(wp) ::  a_2 = -7.38E-5_wp         !< coef. in turb. parametrization (cm-2 s3)
123    REAL(wp) ::  a_3 = -1.40E-2_wp         !< coef. in turb. parametrization
124    REAL(wp) ::  a_term = 9.65_wp          !< coef. for terminal velocity (m s-1)
125    REAL(wp) ::  a_vent = 0.78_wp          !< coef. for ventilation effect
126    REAL(wp) ::  b_1 = 11.45E-6_wp         !< coef. in turb. parametrization (m)
127    REAL(wp) ::  b_2 = 9.68E-6_wp          !< coef. in turb. parametrization (m)
128    REAL(wp) ::  b_3 = 0.62_wp             !< coef. in turb. parametrization
129    REAL(wp) ::  b_term = 9.8_wp           !< coef. for terminal velocity (m s-1)
130    REAL(wp) ::  b_vent = 0.308_wp         !< coef. for ventilation effect
131    REAL(wp) ::  beta_cc = 3.09E-4_wp      !< coef. in turb. parametrization (cm-2 s3)
132    REAL(wp) ::  c_1 = 4.82E-6_wp          !< coef. in turb. parametrization (m)
133    REAL(wp) ::  c_2 = 4.8E-6_wp           !< coef. in turb. parametrization (m)
134    REAL(wp) ::  c_3 = 0.76_wp             !< coef. in turb. parametrization
135    REAL(wp) ::  c_const = 0.93_wp         !< const. in Taylor-microscale Reynolds number
136    REAL(wp) ::  c_evap = 0.7_wp           !< constant in evaporation
137    REAL(wp) ::  c_term = 600.0_wp         !< coef. for terminal velocity (m-1)
138    REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
139    REAL(wp) ::  eps_sb = 1.0E-20_wp       !< threshold in two-moments scheme
140    REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
141    REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
142    REAL(wp) ::  k_rr = 7.12_wp            !< const. rain-rain kernel (m3 kg-1 s-1)
143    REAL(wp) ::  k_br = 1000.0_wp          !< const. in breakup parametrization (m-1)
144    REAL(wp) ::  k_st = 1.2E8_wp           !< const. in drizzle parametrization (m-1 s-1)
145    REAL(wp) ::  kappa_rr = 60.7_wp        !< const. in collision kernel (kg-1/3)
146    REAL(wp) ::  kin_vis_air = 1.4086E-5_wp  !< kin. viscosity of air (m2 s-1)
147    REAL(wp) ::  prec_time_const = 0.001_wp  !< coef. in Kessler scheme (s-1)
148    REAL(wp) ::  ql_crit = 0.0005_wp       !< coef. in Kessler scheme (kg kg-1)
149    REAL(wp) ::  schmidt_p_1d3=0.8921121_wp  !< Schmidt number**0.33333, 0.71**0.33333
150    REAL(wp) ::  sigma_gc = 1.3_wp         !< geometric standard deviation cloud droplets
151    REAL(wp) ::  thermal_conductivity_l = 2.43E-2_wp  !< therm. cond. air (J m-1 s-1 K-1)
152    REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
153    REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
154    REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
155    REAL(wp) ::  xrmax = 5.0E-6_wp         !< maximum rain drop site (kg)
156
157    REAL(wp) ::  c_sedimentation = 2.0_wp  !< Courant number of sedimentation process
158    REAL(wp) ::  dpirho_l                  !< 6.0 / ( pi * rho_l )
159    REAL(wp) ::  dt_micro                  !< microphysics time step
160    REAL(wp) ::  nc_const = 70.0E6_wp      !< cloud droplet concentration
161    REAL(wp) ::  dt_precipitation = 100.0_wp !< timestep precipitation (s)
162    REAL(wp) ::  sed_qc_const              !< const. for sedimentation of cloud water
163    REAL(wp) ::  pirho_l                   !< pi * rho_l / 6.0;
164
165    REAL(wp), DIMENSION(:), ALLOCATABLE ::  nc_1d  !<
166    REAL(wp), DIMENSION(:), ALLOCATABLE ::  nr_1d  !<
167    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_1d  !<
168    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qc_1d  !<
169    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qr_1d  !<
170    REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_1d   !<
171
172    SAVE
173
[1000]174    PRIVATE
[1849]175    PUBLIC microphysics_control, microphysics_init
[1000]176
[1849]177    PUBLIC cloud_water_sedimentation, collision_turbulence, c_sedimentation,   &
178           dt_precipitation, limiter_sedimentation, nc_const, sigma_gc,        &
179           ventilation_effect
180
[1115]181    INTERFACE microphysics_control
182       MODULE PROCEDURE microphysics_control
183       MODULE PROCEDURE microphysics_control_ij
184    END INTERFACE microphysics_control
[1022]185
[1115]186    INTERFACE adjust_cloud
187       MODULE PROCEDURE adjust_cloud
188       MODULE PROCEDURE adjust_cloud_ij
189    END INTERFACE adjust_cloud
190
[1000]191    INTERFACE autoconversion
192       MODULE PROCEDURE autoconversion
193       MODULE PROCEDURE autoconversion_ij
194    END INTERFACE autoconversion
195
[1822]196    INTERFACE autoconversion_kessler
197       MODULE PROCEDURE autoconversion_kessler
198       MODULE PROCEDURE autoconversion_kessler_ij
199    END INTERFACE autoconversion_kessler
200
[1000]201    INTERFACE accretion
202       MODULE PROCEDURE accretion
203       MODULE PROCEDURE accretion_ij
204    END INTERFACE accretion
[1005]205
206    INTERFACE selfcollection_breakup
207       MODULE PROCEDURE selfcollection_breakup
208       MODULE PROCEDURE selfcollection_breakup_ij
209    END INTERFACE selfcollection_breakup
[1012]210
211    INTERFACE evaporation_rain
212       MODULE PROCEDURE evaporation_rain
213       MODULE PROCEDURE evaporation_rain_ij
214    END INTERFACE evaporation_rain
215
216    INTERFACE sedimentation_cloud
217       MODULE PROCEDURE sedimentation_cloud
218       MODULE PROCEDURE sedimentation_cloud_ij
219    END INTERFACE sedimentation_cloud
[1000]220 
[1012]221    INTERFACE sedimentation_rain
222       MODULE PROCEDURE sedimentation_rain
223       MODULE PROCEDURE sedimentation_rain_ij
224    END INTERFACE sedimentation_rain
225
[1691]226    INTERFACE calc_precipitation_amount
227       MODULE PROCEDURE calc_precipitation_amount
228       MODULE PROCEDURE calc_precipitation_amount_ij
229    END INTERFACE calc_precipitation_amount
230
[1000]231 CONTAINS
[1849]232!------------------------------------------------------------------------------!
233! Description:
234! ------------
235!> Initialization of bulk microphysics
236!------------------------------------------------------------------------------!
237    SUBROUTINE microphysics_init
[1000]238
[1849]239       USE arrays_3d,                                                          &
240           ONLY:  dzu
[1000]241
[1849]242       USE constants,                                                          &
243           ONLY:  pi
244
245       USE cloud_parameters,                                                   &
246           ONLY:  rho_l
247
248       USE control_parameters,                                                 &
249           ONLY:  microphysics_seifert
250
251       USE indices,                                                            &
252           ONLY:  nzb, nzt
253
254       IMPLICIT NONE
255
256!
257!--    constant for the sedimentation of cloud water (2-moment cloud physics)
258       sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l )                &
259                          )**( 2.0_wp / 3.0_wp ) *                             &
260                      EXP( 5.0_wp * LOG( sigma_gc )**2 )
261
262!
263!--    Calculate timestep according to precipitation
264       IF ( microphysics_seifert )  THEN
265          dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /      &
266                             w_precipitation
267       ENDIF
268
269!
270!--    Pre-calculate frequently calculated fractions of pi and rho_l
271       pirho_l  = pi * rho_l / 6.0_wp
272       dpirho_l = 1.0_wp / pirho_l
273
274!
275!--    Allocate 1D microphysics arrays
276       ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1),         &
277                  qc_1d(nzb:nzt+1) )
278
279       IF ( microphysics_seifert )  THEN
280          ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
281       ENDIF
282
283!
284!--    Initialze nc_1d with nc_const
285       nc_1d = nc_const
286
287    END SUBROUTINE microphysics_init
288
289
[1000]290!------------------------------------------------------------------------------!
[1682]291! Description:
292! ------------
[1849]293!> Control of microphysics for all grid points
[1000]294!------------------------------------------------------------------------------!
[1115]295    SUBROUTINE microphysics_control
[1022]296
[1361]297       USE arrays_3d,                                                          &
[1849]298           ONLY:  hyp, pt_init, prr, zu
[1361]299
300       USE cloud_parameters,                                                   &
[1849]301           ONLY:  cp, hyrho, pt_d_t, r_d, t_d_pt
[1361]302
303       USE control_parameters,                                                 &
[1849]304           ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
305                  intermediate_timestep_count, large_scale_forcing,            &
[1822]306                  lsf_surf, microphysics_kessler, microphysics_seifert,        &
307                  pt_surface, rho_surface,surface_pressure
[1361]308
309       USE indices,                                                            &
310           ONLY:  nzb, nzt
311
[1320]312       USE kinds
[1115]313
[1361]314       USE statistics,                                                         &
315           ONLY:  weight_pres
316
[1115]317       IMPLICIT NONE
318
[1682]319       INTEGER(iwp) ::  k                 !<
[1115]320
[1682]321       REAL(wp)     ::  t_surface         !<
[1361]322
323       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
324!
325!--       Calculate:
326!--       pt / t : ratio of potential and actual temperature (pt_d_t)
327!--       t / pt : ratio of actual and potential temperature (t_d_pt)
328!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
329          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
330          DO  k = nzb, nzt+1
331             hyp(k)    = surface_pressure * 100.0_wp * &
332                         ( ( t_surface - g / cp * zu(k) ) /                    &
333                         t_surface )**(1.0_wp / 0.286_wp)
334             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
335             t_d_pt(k) = 1.0_wp / pt_d_t(k)
336             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
[1115]337          ENDDO
[1822]338
[1361]339!
340!--       Compute reference density
341          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
342       ENDIF
[1115]343
[1361]344!
345!--    Compute length of time step
346       IF ( call_microphysics_at_all_substeps )  THEN
347          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
348       ELSE
349          dt_micro = dt_3d
350       ENDIF
351
352!
[1822]353!--    Reset precipitation rate
354       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
355
356!
[1361]357!--    Compute cloud physics
[1822]358       IF ( microphysics_kessler )  THEN
359
360          CALL autoconversion_kessler
[1831]361          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
[1822]362
363       ELSEIF ( microphysics_seifert )  THEN
364
[1361]365          CALL adjust_cloud
366          CALL autoconversion
367          CALL accretion
368          CALL selfcollection_breakup
369          CALL evaporation_rain
370          CALL sedimentation_rain
[1831]371          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
[1361]372
[1691]373       ENDIF
374
[1822]375       CALL calc_precipitation_amount
376
[1115]377    END SUBROUTINE microphysics_control
378
[1682]379!------------------------------------------------------------------------------!
380! Description:
381! ------------
382!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
383!> evaporation of rain drops due to too small or too big weights
384!> of rain drops (Stevens and Seifert, 2008).
385!------------------------------------------------------------------------------!
[1115]386    SUBROUTINE adjust_cloud
387
[1361]388       USE arrays_3d,                                                          &
389           ONLY:  qr, nr
390
391       USE cloud_parameters,                                                   &
[1849]392           ONLY:  hyrho
[1361]393
394       USE cpulog,                                                             &
395           ONLY:  cpu_log, log_point_s
396
397       USE indices,                                                            &
[1822]398           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]399
[1320]400       USE kinds
[1022]401
402       IMPLICIT NONE
403
[1682]404       INTEGER(iwp) ::  i                 !<
405       INTEGER(iwp) ::  j                 !<
406       INTEGER(iwp) ::  k                 !<
[1022]407
[1361]408       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
409
[1022]410       DO  i = nxl, nxr
411          DO  j = nys, nyn
[1115]412             DO  k = nzb_s_inner(j,i)+1, nzt
[1361]413                IF ( qr(k,j,i) <= eps_sb )  THEN
414                   qr(k,j,i) = 0.0_wp
415                   nr(k,j,i) = 0.0_wp
416                ELSE
417                   IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
418                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
419                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
420                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
421                   ENDIF
422                ENDIF
[1022]423             ENDDO
424          ENDDO
425       ENDDO
426
[1361]427       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
428
[1115]429    END SUBROUTINE adjust_cloud
[1022]430
[1106]431
[1682]432!------------------------------------------------------------------------------!
433! Description:
434! ------------
435!> Autoconversion rate (Seifert and Beheng, 2006).
436!------------------------------------------------------------------------------!
[1000]437    SUBROUTINE autoconversion
438
[1361]439       USE arrays_3d,                                                          &
440           ONLY:  diss, dzu, nr, qc, qr
441
442       USE cloud_parameters,                                                   &
[1849]443           ONLY:  hyrho
[1361]444
445       USE control_parameters,                                                 &
[1849]446           ONLY:  rho_surface
[1361]447
448       USE cpulog,                                                             &
449           ONLY:  cpu_log, log_point_s
450
451       USE grid_variables,                                                     &
452           ONLY:  dx, dy
453
454       USE indices,                                                            &
[1822]455           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]456
[1320]457       USE kinds
[1000]458
459       IMPLICIT NONE
460
[1682]461       INTEGER(iwp) ::  i                 !<
462       INTEGER(iwp) ::  j                 !<
463       INTEGER(iwp) ::  k                 !<
[1000]464
[1682]465       REAL(wp)     ::  alpha_cc          !<                   
466       REAL(wp)     ::  autocon           !<
467       REAL(wp)     ::  dissipation       !<
468       REAL(wp)     ::  k_au              !<
469       REAL(wp)     ::  l_mix             !<
470       REAL(wp)     ::  nu_c              !<
471       REAL(wp)     ::  phi_au            !<
472       REAL(wp)     ::  r_cc              !<
473       REAL(wp)     ::  rc                !<
474       REAL(wp)     ::  re_lambda         !<
475       REAL(wp)     ::  sigma_cc          !<
476       REAL(wp)     ::  tau_cloud         !<
477       REAL(wp)     ::  xc                !<
[1361]478
479       CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
480
[1000]481       DO  i = nxl, nxr
482          DO  j = nys, nyn
[1115]483             DO  k = nzb_s_inner(j,i)+1, nzt
[1000]484
[1361]485                IF ( qc(k,j,i) > eps_sb )  THEN
486
487                   k_au = k_cc / ( 20.0_wp * x0 )
488!
489!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
490!--                (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
491                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
492!
493!--                Universal function for autoconversion process
494!--                (Seifert and Beheng, 2006):
495                   phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
496                            ( 1.0_wp - tau_cloud**0.68_wp )**3
497!
498!--                Shape parameter of gamma distribution (Geoffroy et al., 2010):
499!--                (Use constant nu_c = 1.0_wp instead?)
500                   nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
501!
502!--                Mean weight of cloud droplets:
503                   xc = hyrho(k) * qc(k,j,i) / nc_const
504!
505!--                Parameterized turbulence effects on autoconversion (Seifert,
506!--                Nuijens and Stevens, 2010)
[1831]507                   IF ( collision_turbulence )  THEN
[1361]508!
509!--                   Weight averaged radius of cloud droplets:
510                      rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
511
512                      alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
513                      r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
514                      sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
515!
516!--                   Mixing length (neglecting distance to ground and
517!--                   stratification)
518                      l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
519!
520!--                   Limit dissipation rate according to Seifert, Nuijens and
521!--                   Stevens (2010)
522                      dissipation = MIN( 0.06_wp, diss(k,j,i) )
523!
524!--                   Compute Taylor-microscale Reynolds number:
525                      re_lambda = 6.0_wp / 11.0_wp *                           &
526                                  ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *   &
527                                  SQRT( 15.0_wp / kin_vis_air ) *              &
528                                  dissipation**( 1.0_wp / 6.0_wp )
529!
530!--                   The factor of 1.0E4 is needed to convert the dissipation
531!--                   rate from m2 s-3 to cm2 s-3.
532                      k_au = k_au * ( 1.0_wp +                                 &
533                                      dissipation * 1.0E4_wp *                 &
534                                      ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
535                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &     
536                                                                      r_cc ) / &
537                                                        sigma_cc )**2          &
538                                                      ) + beta_cc              &
539                                      )                                        &
540                                    )
541                   ENDIF
542!
543!--                Autoconversion rate (Seifert and Beheng, 2006):
544                   autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /    &
545                             ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *     &
546                             ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * &
547                             rho_surface
548                   autocon = MIN( autocon, qc(k,j,i) / dt_micro )
549
550                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro
551                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro 
552                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro
553
554                ENDIF
555
[1000]556             ENDDO
557          ENDDO
558       ENDDO
559
[1361]560       CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
561
[1000]562    END SUBROUTINE autoconversion
563
[1106]564
[1682]565!------------------------------------------------------------------------------!
566! Description:
567! ------------
[1822]568!> Autoconversion process (Kessler, 1969).
569!------------------------------------------------------------------------------!
570    SUBROUTINE autoconversion_kessler
571
572       USE arrays_3d,                                                          &
[1849]573           ONLY:  dzw, pt, prr, q, qc
[1822]574
575       USE cloud_parameters,                                                   &
[1849]576           ONLY:  l_d_cp, pt_d_t
[1822]577
578       USE indices,                                                            &
[1845]579           ONLY:  nxl, nxr, nyn, nys, nzb_s_inner, nzt
[1822]580
581       USE kinds
582
583
584       IMPLICIT NONE
585
586       INTEGER(iwp) ::  i !<
587       INTEGER(iwp) ::  j !<
588       INTEGER(iwp) ::  k !<
589
590       REAL(wp)    ::  dqdt_precip !<
591
592       DO  i = nxl, nxr
593          DO  j = nys, nyn
[1845]594             DO  k = nzb_s_inner(j,i)+1, nzt
[1822]595
596                IF ( qc(k,j,i) > ql_crit )  THEN
597                   dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
598                ELSE
599                   dqdt_precip = 0.0_wp
600                ENDIF
601
602                qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro
603                q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro
[1845]604                pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp *      &
605                                        pt_d_t(k)
[1822]606
607!
[1845]608!--             Compute the rain rate (stored on surface grid point)
609                prr(nzb_s_inner(j,i),j,i) = prr(nzb_s_inner(j,i),j,i) +        &
610                                            dqdt_precip * dzw(k)
[1822]611
612             ENDDO
613          ENDDO
614       ENDDO
615
616   END SUBROUTINE autoconversion_kessler
617
618
619!------------------------------------------------------------------------------!
620! Description:
621! ------------
[1682]622!> Accretion rate (Seifert and Beheng, 2006).
623!------------------------------------------------------------------------------!
[1005]624    SUBROUTINE accretion
[1000]625
[1361]626       USE arrays_3d,                                                          &
627           ONLY:  diss, qc, qr
628
629       USE cloud_parameters,                                                   &
[1849]630           ONLY:  hyrho
[1361]631
632       USE control_parameters,                                                 &
[1849]633           ONLY:  rho_surface
[1361]634
635       USE cpulog,                                                             &
636           ONLY:  cpu_log, log_point_s
637
638       USE indices,                                                            &
[1822]639           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]640
[1320]641       USE kinds
[1005]642
[1000]643       IMPLICIT NONE
644
[1682]645       INTEGER(iwp) ::  i                 !<
646       INTEGER(iwp) ::  j                 !<
647       INTEGER(iwp) ::  k                 !<
[1000]648
[1682]649       REAL(wp)     ::  accr              !<
650       REAL(wp)     ::  k_cr              !<
651       REAL(wp)     ::  phi_ac            !<
652       REAL(wp)     ::  tau_cloud         !<
[1361]653
654       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
655
[1005]656       DO  i = nxl, nxr
657          DO  j = nys, nyn
[1115]658             DO  k = nzb_s_inner(j,i)+1, nzt
[1000]659
[1361]660                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb ) )  THEN
661!
662!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
663                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ) 
664!
665!--                Universal function for accretion process (Seifert and
666!--                Beheng, 2001):
667                   phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
668!
669!--                Parameterized turbulence effects on autoconversion (Seifert,
670!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
671!--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1831]672                   IF ( collision_turbulence )  THEN
[1361]673                      k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                      &
674                                       MIN( 600.0_wp,                          &
675                                            diss(k,j,i) * 1.0E4_wp )**0.25_wp  &
676                                     )
677                   ELSE
678                      k_cr = k_cr0                       
679                   ENDIF
680!
681!--                Accretion rate (Seifert and Beheng, 2006):
682                   accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *              &
683                          SQRT( rho_surface * hyrho(k) )
684                   accr = MIN( accr, qc(k,j,i) / dt_micro )
685
686                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro 
687                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro
688
689                ENDIF
690
[1005]691             ENDDO
692          ENDDO
[1000]693       ENDDO
694
[1361]695       CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
696
[1005]697    END SUBROUTINE accretion
[1000]698
[1106]699
[1682]700!------------------------------------------------------------------------------!
701! Description:
702! ------------
703!> Collisional breakup rate (Seifert, 2008).
704!------------------------------------------------------------------------------!
[1005]705    SUBROUTINE selfcollection_breakup
[1000]706
[1361]707       USE arrays_3d,                                                          &
708           ONLY:  nr, qr
709
710       USE cloud_parameters,                                                   &
[1849]711           ONLY:  hyrho
[1361]712
713       USE control_parameters,                                                 &
[1849]714           ONLY:  rho_surface
[1361]715
716       USE cpulog,                                                             &
717           ONLY:  cpu_log, log_point_s
718
719       USE indices,                                                            &
[1822]720           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]721
[1320]722       USE kinds
[1361]723   
[1000]724       IMPLICIT NONE
725
[1682]726       INTEGER(iwp) ::  i                 !<
727       INTEGER(iwp) ::  j                 !<
728       INTEGER(iwp) ::  k                 !<
[1000]729
[1682]730       REAL(wp)     ::  breakup           !<
731       REAL(wp)     ::  dr                !<
732       REAL(wp)     ::  phi_br            !<
733       REAL(wp)     ::  selfcoll          !<
[1361]734
735       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
736
[1000]737       DO  i = nxl, nxr
738          DO  j = nys, nyn
[1115]739             DO  k = nzb_s_inner(j,i)+1, nzt
[1361]740                IF ( qr(k,j,i) > eps_sb )  THEN
741!
742!--                Selfcollection rate (Seifert and Beheng, 2001):
743                   selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) *                   &
744                              SQRT( hyrho(k) * rho_surface )
745!
746!--                Weight averaged diameter of rain drops:
747                   dr = ( hyrho(k) * qr(k,j,i) /                               &
748                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
749!
750!--                Collisional breakup rate (Seifert, 2008):
751                   IF ( dr >= 0.3E-3_wp )  THEN
752                      phi_br  = k_br * ( dr - 1.1E-3_wp )
753                      breakup = selfcoll * ( phi_br + 1.0_wp )
754                   ELSE
755                      breakup = 0.0_wp
756                   ENDIF
[1000]757
[1361]758                   selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
759                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro
760
761                ENDIF         
[1000]762             ENDDO
763          ENDDO
764       ENDDO
765
[1361]766       CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
767
[1005]768    END SUBROUTINE selfcollection_breakup
[1000]769
[1106]770
[1682]771!------------------------------------------------------------------------------!
772! Description:
773! ------------
774!> Evaporation of precipitable water. Condensation is neglected for
775!> precipitable water.
776!------------------------------------------------------------------------------!
[1012]777    SUBROUTINE evaporation_rain
[1000]778
[1361]779       USE arrays_3d,                                                          &
780           ONLY:  hyp, nr, pt, q,  qc, qr
781
782       USE cloud_parameters,                                                   &
[1849]783           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
[1361]784
785       USE constants,                                                          &
786           ONLY:  pi
787
788       USE cpulog,                                                             &
789           ONLY:  cpu_log, log_point_s
790
791       USE indices,                                                            &
[1822]792           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]793
[1320]794       USE kinds
[1012]795
796       IMPLICIT NONE
797
[1682]798       INTEGER(iwp) ::  i                 !<
799       INTEGER(iwp) ::  j                 !<
800       INTEGER(iwp) ::  k                 !<
[1361]801
[1682]802       REAL(wp)     ::  alpha             !<
803       REAL(wp)     ::  dr                !<
804       REAL(wp)     ::  e_s               !<
805       REAL(wp)     ::  evap              !<
806       REAL(wp)     ::  evap_nr           !<
807       REAL(wp)     ::  f_vent            !<
808       REAL(wp)     ::  g_evap            !<
809       REAL(wp)     ::  lambda_r          !<
810       REAL(wp)     ::  mu_r              !<
811       REAL(wp)     ::  mu_r_2            !<
812       REAL(wp)     ::  mu_r_5d2          !<
813       REAL(wp)     ::  nr_0              !<
814       REAL(wp)     ::  q_s               !<
815       REAL(wp)     ::  sat               !<
816       REAL(wp)     ::  t_l               !<
817       REAL(wp)     ::  temp              !<
818       REAL(wp)     ::  xr                !<
[1361]819
820       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
821
[1012]822       DO  i = nxl, nxr
823          DO  j = nys, nyn
[1115]824             DO  k = nzb_s_inner(j,i)+1, nzt
[1361]825                IF ( qr(k,j,i) > eps_sb )  THEN
826!
827!--                Actual liquid water temperature:
828                   t_l = t_d_pt(k) * pt(k,j,i)
829!
830!--                Saturation vapor pressure at t_l:
831                   e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
832                                          ( t_l - 35.86_wp )                   &
833                                        )
834!
835!--                Computation of saturation humidity:
836                   q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
837                   alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
838                   q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /               &
839                           ( 1.0_wp + alpha * q_s )
840!
841!--                Supersaturation:
842                   sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
843!
844!--                Evaporation needs only to be calculated in subsaturated regions
845                   IF ( sat < 0.0_wp )  THEN
846!
847!--                   Actual temperature:
848                      temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
849             
850                      g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *     &
851                                          l_v / ( thermal_conductivity_l * temp ) &
852                                          + r_v * temp / ( diff_coeff_l * e_s )   &
853                                        )
854!
855!--                   Mean weight of rain drops
856                      xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
857!
858!--                   Weight averaged diameter of rain drops:
859                      dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
860!
861!--                   Compute ventilation factor and intercept parameter
862!--                   (Seifert and Beheng, 2006; Seifert, 2008):
863                      IF ( ventilation_effect )  THEN
864!
865!--                      Shape parameter of gamma distribution (Milbrandt and Yau,
866!--                      2005; Stevens and Seifert, 2008):
867                         mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *          &
868                                                          ( dr - 1.4E-3_wp ) ) )
869!
870!--                      Slope parameter of gamma distribution (Seifert, 2008):
871                         lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *  &
872                                      ( mu_r + 1.0_wp )                        &
873                                    )**( 1.0_wp / 3.0_wp ) / dr
[1012]874
[1361]875                         mu_r_2   = mu_r + 2.0_wp
876                         mu_r_5d2 = mu_r + 2.5_wp 
877
878                         f_vent = a_vent * gamm( mu_r_2 ) *                     &
879                                  lambda_r**( -mu_r_2 ) + b_vent *              &
880                                  schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
881                                  gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) *  &
882                                  ( 1.0_wp -                                    &
883                                    0.5_wp * ( b_term / a_term ) *              &
884                                    ( lambda_r / ( c_term + lambda_r )          &
885                                    )**mu_r_5d2 -                               &
886                                    0.125_wp * ( b_term / a_term )**2 *         &
887                                    ( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
888                                    )**mu_r_5d2 -                               &
889                                    0.0625_wp * ( b_term / a_term )**3 *        &
890                                    ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
891                                    )**mu_r_5d2 -                               &
892                                    0.0390625_wp * ( b_term / a_term )**4 *     & 
893                                    ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
894                                    )**mu_r_5d2                                 &
895                                  )
896
897                         nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
898                                  gamm( mu_r + 1.0_wp ) 
899                      ELSE
900                         f_vent = 1.0_wp
901                         nr_0   = nr(k,j,i) * dr
902                      ENDIF
903   !
904   !--                Evaporation rate of rain water content (Seifert and
905   !--                Beheng, 2006):
906                      evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
907                                hyrho(k)
908                      evap    = MAX( evap, -qr(k,j,i) / dt_micro )
909                      evap_nr = MAX( c_evap * evap / xr * hyrho(k),            &
910                                     -nr(k,j,i) / dt_micro )
911
912                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro
913                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro
914
915                   ENDIF
916                ENDIF         
917
[1012]918             ENDDO
919          ENDDO
920       ENDDO
921
[1361]922       CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
923
[1012]924    END SUBROUTINE evaporation_rain
925
[1106]926
[1682]927!------------------------------------------------------------------------------!
928! Description:
929! ------------
930!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
931!------------------------------------------------------------------------------!
[1012]932    SUBROUTINE sedimentation_cloud
933
[1361]934       USE arrays_3d,                                                          &
[1849]935           ONLY:  ddzu, dzu, pt, prr, q, qc
[1361]936
937       USE cloud_parameters,                                                   &
[1849]938           ONLY:  hyrho, l_d_cp, pt_d_t
[1361]939
940       USE control_parameters,                                                 &
[1849]941           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
[1361]942
943       USE cpulog,                                                             &
944           ONLY:  cpu_log, log_point_s
945
946       USE indices,                                                            &
947           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
948
[1320]949       USE kinds
[1691]950
951       USE statistics,                                                         &
952           ONLY:  weight_substep
953   
954
[1012]955       IMPLICIT NONE
956
[1849]957       INTEGER(iwp) ::  i             !<
958       INTEGER(iwp) ::  j             !<
959       INTEGER(iwp) ::  k             !<
[1361]960
[1682]961       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
[1361]962
963       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
964
965       sed_qc(nzt+1) = 0.0_wp
966
[1012]967       DO  i = nxl, nxr
968          DO  j = nys, nyn
[1361]969             DO  k = nzt, nzb_s_inner(j,i)+1, -1
[1012]970
[1361]971                IF ( qc(k,j,i) > eps_sb )  THEN
972                   sed_qc(k) = sed_qc_const * nc_const**( -2.0_wp / 3.0_wp ) * &
973                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )
974                ELSE
975                   sed_qc(k) = 0.0_wp
976                ENDIF
977
978                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
979                                            dt_micro + sed_qc(k+1)             &
980                               )
981
982                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
983                                        ddzu(k+1) / hyrho(k) * dt_micro
984                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
985                                        ddzu(k+1) / hyrho(k) * dt_micro
986                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
987                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
988                                        pt_d_t(k) * dt_micro
989
[1691]990!
991!--             Compute the precipitation rate due to cloud (fog) droplets
[1822]992                IF ( call_microphysics_at_all_substeps )  THEN
993                   prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
994                                * weight_substep(intermediate_timestep_count)
995                ELSE
996                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
[1691]997                ENDIF
998
[1012]999             ENDDO
1000          ENDDO
1001       ENDDO
1002
[1361]1003       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
1004
[1012]1005    END SUBROUTINE sedimentation_cloud
1006
[1106]1007
[1682]1008!------------------------------------------------------------------------------!
1009! Description:
1010! ------------
1011!> Computation of sedimentation flux. Implementation according to Stevens
1012!> and Seifert (2008). Code is based on UCLA-LES.
1013!------------------------------------------------------------------------------!
[1012]1014    SUBROUTINE sedimentation_rain
1015
[1361]1016       USE arrays_3d,                                                          &
[1849]1017           ONLY:  ddzu, dzu, nr, pt, prr, q, qr
[1361]1018
1019       USE cloud_parameters,                                                   &
[1849]1020           ONLY:  hyrho, l_d_cp, pt_d_t
[1361]1021
1022       USE control_parameters,                                                 &
[1849]1023           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
[1361]1024       USE cpulog,                                                             &
1025           ONLY:  cpu_log, log_point_s
1026
1027       USE indices,                                                            &
1028           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
1029
[1320]1030       USE kinds
[1012]1031
[1361]1032       USE statistics,                                                         &
1033           ONLY:  weight_substep
1034       
[1012]1035       IMPLICIT NONE
1036
[1682]1037       INTEGER(iwp) ::  i                          !<
1038       INTEGER(iwp) ::  j                          !<
1039       INTEGER(iwp) ::  k                          !<
1040       INTEGER(iwp) ::  k_run                      !<
[1361]1041
[1682]1042       REAL(wp)     ::  c_run                      !<
1043       REAL(wp)     ::  d_max                      !<
1044       REAL(wp)     ::  d_mean                     !<
1045       REAL(wp)     ::  d_min                      !<
1046       REAL(wp)     ::  dr                         !<
1047       REAL(wp)     ::  flux                       !<
1048       REAL(wp)     ::  lambda_r                   !<
1049       REAL(wp)     ::  mu_r                       !<
1050       REAL(wp)     ::  z_run                      !<
[1361]1051
[1682]1052       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
1053       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
1054       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
1055       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
1056       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
1057       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
1058       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
1059       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
[1361]1060
1061       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
[1682]1062
[1361]1063!
1064!--    Compute velocities
[1012]1065       DO  i = nxl, nxr
1066          DO  j = nys, nyn
[1115]1067             DO  k = nzb_s_inner(j,i)+1, nzt
[1361]1068                IF ( qr(k,j,i) > eps_sb )  THEN
1069!
1070!--                Weight averaged diameter of rain drops:
1071                   dr = ( hyrho(k) * qr(k,j,i) /                               &
1072                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
1073!
1074!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1075!--                Stevens and Seifert, 2008):
1076                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
1077                                                     ( dr - 1.4E-3_wp ) ) )
1078!
1079!--                Slope parameter of gamma distribution (Seifert, 2008):
1080                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
1081                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1012]1082
[1361]1083                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
1084                                               a_term - b_term * ( 1.0_wp +    &
1085                                                  c_term /                     &
1086                                                  lambda_r )**( -1.0_wp *      &
1087                                                  ( mu_r + 1.0_wp ) )          &
1088                                              )                                &
1089                                )
1090
1091                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
1092                                               a_term - b_term * ( 1.0_wp +    &
1093                                                  c_term /                     &
1094                                                  lambda_r )**( -1.0_wp *      &
1095                                                  ( mu_r + 4.0_wp ) )          &
1096                                             )                                 &
1097                                )
1098                ELSE
1099                   w_nr(k) = 0.0_wp
1100                   w_qr(k) = 0.0_wp
1101                ENDIF
[1012]1102             ENDDO
[1361]1103!
1104!--          Adjust boundary values
1105             w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
1106             w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
1107             w_nr(nzt+1) = 0.0_wp
1108             w_qr(nzt+1) = 0.0_wp
1109!
1110!--          Compute Courant number
1111             DO  k = nzb_s_inner(j,i)+1, nzt
1112                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
1113                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
1114                          dt_micro * ddzu(k)
1115                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
1116                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
1117                          dt_micro * ddzu(k)
1118             ENDDO     
1119!
1120!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
1121             IF ( limiter_sedimentation )  THEN
1122
1123                DO k = nzb_s_inner(j,i)+1, nzt
[1646]1124                   d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
[1361]1125                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
1126                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
1127
1128                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
1129                                                              2.0_wp * d_max,  &
1130                                                              ABS( d_mean ) )
1131
[1646]1132                   d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
[1361]1133                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
1134                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
1135
1136                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
1137                                                              2.0_wp * d_max,  &
1138                                                              ABS( d_mean ) )
1139                ENDDO
1140
1141             ELSE
1142
1143                nr_slope = 0.0_wp
1144                qr_slope = 0.0_wp
1145
1146             ENDIF
1147
1148             sed_nr(nzt+1) = 0.0_wp
1149             sed_qr(nzt+1) = 0.0_wp
1150!
1151!--          Compute sedimentation flux
1152             DO  k = nzt, nzb_s_inner(j,i)+1, -1
1153!
1154!--             Sum up all rain drop number densities which contribute to the flux
1155!--             through k-1/2
1156                flux  = 0.0_wp
1157                z_run = 0.0_wp ! height above z(k)
1158                k_run = k
1159                c_run = MIN( 1.0_wp, c_nr(k) )
1160                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
1161                   flux  = flux + hyrho(k_run) *                               &
1162                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
1163                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)
1164                   z_run = z_run + dzu(k_run)
1165                   k_run = k_run + 1
1166                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
1167                ENDDO
1168!
1169!--             It is not allowed to sediment more rain drop number density than
1170!--             available
1171                flux = MIN( flux,                                              &
1172                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
1173                            dt_micro                                           &
1174                          )
1175
1176                sed_nr(k) = flux / dt_micro
1177                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
1178                                        ddzu(k+1) / hyrho(k) * dt_micro
1179!
1180!--             Sum up all rain water content which contributes to the flux
1181!--             through k-1/2
1182                flux  = 0.0_wp
1183                z_run = 0.0_wp ! height above z(k)
1184                k_run = k
1185                c_run = MIN( 1.0_wp, c_qr(k) )
1186
1187                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
1188
1189                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
1190                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
1191                                  0.5_wp ) * c_run * dzu(k_run)
1192                   z_run = z_run + dzu(k_run)
1193                   k_run = k_run + 1
1194                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
1195
1196                ENDDO
1197!
1198!--             It is not allowed to sediment more rain water content than
1199!--             available
1200                flux = MIN( flux,                                              &
1201                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
1202                            dt_micro                                           &
1203                          )
1204
1205                sed_qr(k) = flux / dt_micro
1206
1207                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
1208                                        ddzu(k+1) / hyrho(k) * dt_micro
1209                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
1210                                        ddzu(k+1) / hyrho(k) * dt_micro 
1211                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
1212                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
1213                                        pt_d_t(k) * dt_micro
1214!
1215!--             Compute the rain rate
1216                IF ( call_microphysics_at_all_substeps )  THEN
[1691]1217                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
1218                                * weight_substep(intermediate_timestep_count)
[1361]1219                ELSE
[1691]1220                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
[1361]1221                ENDIF
1222
1223             ENDDO
[1012]1224          ENDDO
1225       ENDDO
1226
[1691]1227       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
1228
1229    END SUBROUTINE sedimentation_rain
1230
1231
1232!------------------------------------------------------------------------------!
1233! Description:
1234! ------------
1235!> Computation of the precipitation amount due to gravitational settling of
1236!> rain and cloud (fog) droplets
1237!------------------------------------------------------------------------------!
1238    SUBROUTINE calc_precipitation_amount
1239
[1849]1240       USE arrays_3d,                                                          &
1241           ONLY:  precipitation_amount, prr
1242
[1691]1243       USE cloud_parameters,                                                   &
[1849]1244           ONLY:  hyrho
[1691]1245
1246       USE control_parameters,                                                 &
1247           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
1248                  intermediate_timestep_count, intermediate_timestep_count_max,&
1249                  precipitation_amount_interval, time_do2d_xy
1250
1251       USE indices,                                                            &
1252           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner
1253
1254       USE kinds
1255
1256       IMPLICIT NONE
1257
1258       INTEGER(iwp) ::  i                          !:
1259       INTEGER(iwp) ::  j                          !:
1260
1261
1262       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
1263            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
1264            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
1265       THEN
1266
[1361]1267          DO  i = nxl, nxr
1268             DO  j = nys, nyn
[1691]1269
[1361]1270                precipitation_amount(j,i) = precipitation_amount(j,i) +        &
1271                                            prr(nzb_s_inner(j,i)+1,j,i) *      &
1272                                            hyrho(nzb_s_inner(j,i)+1) * dt_3d
[1691]1273
[1361]1274             ENDDO
1275          ENDDO
1276       ENDIF
1277
[1691]1278    END SUBROUTINE calc_precipitation_amount
[1361]1279
[1012]1280
[1000]1281!------------------------------------------------------------------------------!
[1682]1282! Description:
1283! ------------
[1849]1284!> Control of microphysics for grid points i,j
[1000]1285!------------------------------------------------------------------------------!
[1022]1286
[1115]1287    SUBROUTINE microphysics_control_ij( i, j )
1288
[1320]1289       USE arrays_3d,                                                          &
[1849]1290           ONLY:  hyp, nr, pt, pt_init, prr, q, qc, qr, zu
[1115]1291
[1320]1292       USE cloud_parameters,                                                   &
[1849]1293           ONLY:  cp, hyrho, pt_d_t, r_d, t_d_pt
[1320]1294
1295       USE control_parameters,                                                 &
[1849]1296           ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
1297                  intermediate_timestep_count, large_scale_forcing,            &
[1822]1298                  lsf_surf, microphysics_seifert, microphysics_kessler,        &
1299                  pt_surface, rho_surface, surface_pressure
[1320]1300
1301       USE indices,                                                            &
1302           ONLY:  nzb, nzt
1303
1304       USE kinds
1305
1306       USE statistics,                                                         &
1307           ONLY:  weight_pres
1308
[1022]1309       IMPLICIT NONE
1310
[1682]1311       INTEGER(iwp) ::  i                 !<
1312       INTEGER(iwp) ::  j                 !<
1313       INTEGER(iwp) ::  k                 !<
[1115]1314
[1682]1315       REAL(wp)     ::  t_surface         !<
[1320]1316
[1361]1317       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
[1241]1318!
1319!--       Calculate:
1320!--       pt / t : ratio of potential and actual temperature (pt_d_t)
1321!--       t / pt : ratio of actual and potential temperature (t_d_pt)
1322!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
[1353]1323          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
[1241]1324          DO  k = nzb, nzt+1
[1353]1325             hyp(k)    = surface_pressure * 100.0_wp * &
[1361]1326                         ( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
[1353]1327             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
1328             t_d_pt(k) = 1.0_wp / pt_d_t(k)
[1241]1329             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
1330          ENDDO
1331!
1332!--       Compute reference density
[1353]1333          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
[1241]1334       ENDIF
1335
[1361]1336!
1337!--    Compute length of time step
1338       IF ( call_microphysics_at_all_substeps )  THEN
1339          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
1340       ELSE
1341          dt_micro = dt_3d
1342       ENDIF
[1241]1343
[1115]1344!
[1361]1345!--    Use 1d arrays
[1115]1346       q_1d(:)  = q(:,j,i)
1347       pt_1d(:) = pt(:,j,i)
1348       qc_1d(:) = qc(:,j,i)
1349       nc_1d(:) = nc_const
[1822]1350       IF ( microphysics_seifert )  THEN
[1115]1351          qr_1d(:) = qr(:,j,i)
1352          nr_1d(:) = nr(:,j,i)
1353       ENDIF
[1361]1354
[1115]1355!
[1822]1356!--    Reset precipitation rate
1357       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
1358
1359!
[1115]1360!--    Compute cloud physics
[1822]1361       IF( microphysics_kessler )  THEN
1362
1363          CALL autoconversion_kessler( i,j )
[1831]1364          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
[1822]1365
1366       ELSEIF ( microphysics_seifert )  THEN
1367
1368          CALL adjust_cloud( i,j )
[1115]1369          CALL autoconversion( i,j )
1370          CALL accretion( i,j )
1371          CALL selfcollection_breakup( i,j )
1372          CALL evaporation_rain( i,j )
1373          CALL sedimentation_rain( i,j )
[1831]1374          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
[1115]1375
[1691]1376       ENDIF
1377
[1822]1378       CALL calc_precipitation_amount( i,j )
1379
[1115]1380!
[1361]1381!--    Store results on the 3d arrays
1382       q(:,j,i)  = q_1d(:)
1383       pt(:,j,i) = pt_1d(:)
[1822]1384       IF ( microphysics_seifert )  THEN
[1361]1385          qr(:,j,i) = qr_1d(:)
1386          nr(:,j,i) = nr_1d(:)
[1115]1387       ENDIF
1388
1389    END SUBROUTINE microphysics_control_ij
1390
[1682]1391!------------------------------------------------------------------------------!
1392! Description:
1393! ------------
1394!> Adjust number of raindrops to avoid nonlinear effects in
1395!> sedimentation and evaporation of rain drops due to too small or
1396!> too big weights of rain drops (Stevens and Seifert, 2008).
1397!> The same procedure is applied to cloud droplets if they are determined
1398!> prognostically. Call for grid point i,j
1399!------------------------------------------------------------------------------!
[1115]1400    SUBROUTINE adjust_cloud_ij( i, j )
1401
[1320]1402       USE cloud_parameters,                                                   &
[1849]1403           ONLY:  hyrho
[1320]1404
1405       USE indices,                                                            &
[1822]1406           ONLY:  nzb_s_inner, nzt
[1320]1407
1408       USE kinds
1409
[1115]1410       IMPLICIT NONE
1411
[1682]1412       INTEGER(iwp) ::  i                 !<
1413       INTEGER(iwp) ::  j                 !<
1414       INTEGER(iwp) ::  k                 !<
1415
[1115]1416       DO  k = nzb_s_inner(j,i)+1, nzt
[1022]1417
[1361]1418          IF ( qr_1d(k) <= eps_sb )  THEN
1419             qr_1d(k) = 0.0_wp
1420             nr_1d(k) = 0.0_wp
[1065]1421          ELSE
[1022]1422!
[1048]1423!--          Adjust number of raindrops to avoid nonlinear effects in
1424!--          sedimentation and evaporation of rain drops due to too small or
[1065]1425!--          too big weights of rain drops (Stevens and Seifert, 2008).
[1361]1426             IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) )  THEN
1427                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin
1428             ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) )  THEN
1429                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax
[1048]1430             ENDIF
[1115]1431
[1022]1432          ENDIF
[1115]1433
[1022]1434       ENDDO
1435
[1115]1436    END SUBROUTINE adjust_cloud_ij
[1022]1437
[1106]1438
[1682]1439!------------------------------------------------------------------------------!
1440! Description:
1441! ------------
1442!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
1443!------------------------------------------------------------------------------!
[1005]1444    SUBROUTINE autoconversion_ij( i, j )
[1000]1445
[1320]1446       USE arrays_3d,                                                          &
[1849]1447           ONLY:  diss, dzu
[1115]1448
[1320]1449       USE cloud_parameters,                                                   &
[1849]1450           ONLY:  hyrho
[1320]1451
1452       USE control_parameters,                                                 &
[1849]1453           ONLY:  rho_surface
[1320]1454
1455       USE grid_variables,                                                     &
1456           ONLY:  dx, dy
1457
1458       USE indices,                                                            &
[1822]1459           ONLY:  nzb_s_inner, nzt
[1320]1460
1461       USE kinds
1462
[1000]1463       IMPLICIT NONE
1464
[1682]1465       INTEGER(iwp) ::  i                 !<
1466       INTEGER(iwp) ::  j                 !<
1467       INTEGER(iwp) ::  k                 !<
[1000]1468
[1682]1469       REAL(wp)     ::  alpha_cc          !<                   
1470       REAL(wp)     ::  autocon           !<
1471       REAL(wp)     ::  dissipation       !<
1472       REAL(wp)     ::  k_au              !<
1473       REAL(wp)     ::  l_mix             !<
1474       REAL(wp)     ::  nu_c              !<
1475       REAL(wp)     ::  phi_au            !<
1476       REAL(wp)     ::  r_cc              !<
1477       REAL(wp)     ::  rc                !<
1478       REAL(wp)     ::  re_lambda         !<
1479       REAL(wp)     ::  sigma_cc          !<
1480       REAL(wp)     ::  tau_cloud         !<
1481       REAL(wp)     ::  xc                !<
[1106]1482
[1115]1483       DO  k = nzb_s_inner(j,i)+1, nzt
[1000]1484
[1115]1485          IF ( qc_1d(k) > eps_sb )  THEN
[1361]1486
1487             k_au = k_cc / ( 20.0_wp * x0 )
[1012]1488!
[1048]1489!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]1490!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
1491             tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
[1012]1492!
1493!--          Universal function for autoconversion process
1494!--          (Seifert and Beheng, 2006):
[1361]1495             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
[1012]1496!
1497!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
[1353]1498!--          (Use constant nu_c = 1.0_wp instead?)
[1361]1499             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
[1012]1500!
1501!--          Mean weight of cloud droplets:
[1115]1502             xc = hyrho(k) * qc_1d(k) / nc_1d(k)
[1012]1503!
[1065]1504!--          Parameterized turbulence effects on autoconversion (Seifert,
1505!--          Nuijens and Stevens, 2010)
[1831]1506             IF ( collision_turbulence )  THEN
[1065]1507!
1508!--             Weight averaged radius of cloud droplets:
[1353]1509                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
[1065]1510
[1353]1511                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
1512                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
1513                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
[1065]1514!
1515!--             Mixing length (neglecting distance to ground and stratification)
[1334]1516                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
[1065]1517!
1518!--             Limit dissipation rate according to Seifert, Nuijens and
1519!--             Stevens (2010)
[1361]1520                dissipation = MIN( 0.06_wp, diss(k,j,i) )
[1065]1521!
1522!--             Compute Taylor-microscale Reynolds number:
[1361]1523                re_lambda = 6.0_wp / 11.0_wp *                                 &
1524                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
1525                            SQRT( 15.0_wp / kin_vis_air ) *                    &
1526                            dissipation**( 1.0_wp / 6.0_wp )
[1065]1527!
1528!--             The factor of 1.0E4 is needed to convert the dissipation rate
1529!--             from m2 s-3 to cm2 s-3.
[1361]1530                k_au = k_au * ( 1.0_wp +                                       &
1531                                dissipation * 1.0E4_wp *                       &
1532                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
1533                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
1534                                                  sigma_cc )**2                &
1535                                                ) + beta_cc                    &
1536                                )                                              &
1537                              )
[1065]1538             ENDIF
1539!
[1012]1540!--          Autoconversion rate (Seifert and Beheng, 2006):
[1361]1541             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
1542                       ( nu_c + 1.0_wp )**2 * qc_1d(k)**2 * xc**2 *            &
1543                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
[1115]1544                       rho_surface
1545             autocon = MIN( autocon, qc_1d(k) / dt_micro )
[1106]1546
[1115]1547             qr_1d(k) = qr_1d(k) + autocon * dt_micro
1548             qc_1d(k) = qc_1d(k) - autocon * dt_micro 
1549             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
1550
[1005]1551          ENDIF
[1000]1552
1553       ENDDO
1554
[1005]1555    END SUBROUTINE autoconversion_ij
1556
[1822]1557!------------------------------------------------------------------------------!
1558! Description:
1559! ------------
1560!> Autoconversion process (Kessler, 1969).
1561!------------------------------------------------------------------------------!
1562    SUBROUTINE autoconversion_kessler_ij( i, j )
[1106]1563
[1822]1564       USE arrays_3d,                                                          &
[1849]1565           ONLY:  dzw, prr
[1822]1566
1567       USE cloud_parameters,                                                   &
[1849]1568           ONLY:  l_d_cp, pt_d_t
[1822]1569
1570       USE indices,                                                            &
[1845]1571           ONLY:  nzb_s_inner, nzt
[1822]1572
1573       USE kinds
1574
1575
1576       IMPLICIT NONE
1577
1578       INTEGER(iwp) ::  i !<
1579       INTEGER(iwp) ::  j !<
1580       INTEGER(iwp) ::  k !<
1581
1582       REAL(wp)    ::  dqdt_precip !<
1583
[1845]1584       DO  k = nzb_s_inner(j,i)+1, nzt
[1822]1585
1586          IF ( qc_1d(k) > ql_crit )  THEN
1587             dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit )
1588          ELSE
1589             dqdt_precip = 0.0_wp
1590          ENDIF
1591
1592          qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro
1593          q_1d(k)  = q_1d(k)  - dqdt_precip * dt_micro
1594          pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k)
1595
1596!
[1845]1597!--       Compute the rain rate (stored on surface grid point)
1598          prr(nzb_s_inner(j,i),j,i) = prr(nzb_s_inner(j,i),j,i) +              &
1599                                      dqdt_precip * dzw(k)
[1822]1600
1601       ENDDO
1602
1603    END SUBROUTINE autoconversion_kessler_ij
1604
[1682]1605!------------------------------------------------------------------------------!
1606! Description:
1607! ------------
1608!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
1609!------------------------------------------------------------------------------!
[1005]1610    SUBROUTINE accretion_ij( i, j )
1611
[1320]1612       USE arrays_3d,                                                          &
[1849]1613           ONLY:  diss
[1115]1614
[1320]1615       USE cloud_parameters,                                                   &
[1849]1616           ONLY:  hyrho
[1320]1617
1618       USE control_parameters,                                                 &
[1849]1619           ONLY:  rho_surface
[1320]1620
1621       USE indices,                                                            &
[1822]1622           ONLY:  nzb_s_inner, nzt
[1320]1623
1624       USE kinds
1625
[1005]1626       IMPLICIT NONE
1627
[1682]1628       INTEGER(iwp) ::  i                 !<
1629       INTEGER(iwp) ::  j                 !<
1630       INTEGER(iwp) ::  k                 !<
[1005]1631
[1682]1632       REAL(wp)     ::  accr              !<
1633       REAL(wp)     ::  k_cr              !<
1634       REAL(wp)     ::  phi_ac            !<
1635       REAL(wp)     ::  tau_cloud         !<
[1320]1636
[1115]1637       DO  k = nzb_s_inner(j,i)+1, nzt
1638          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
[1012]1639!
[1048]1640!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]1641             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 
[1012]1642!
1643!--          Universal function for accretion process
[1048]1644!--          (Seifert and Beheng, 2001):
[1361]1645             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
[1012]1646!
[1065]1647!--          Parameterized turbulence effects on autoconversion (Seifert,
1648!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
[1361]1649!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1831]1650             IF ( collision_turbulence )  THEN
[1361]1651                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
1652                                 MIN( 600.0_wp,                                &
1653                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
1654                               )
[1065]1655             ELSE
1656                k_cr = k_cr0                       
1657             ENDIF
1658!
[1012]1659!--          Accretion rate (Seifert and Beheng, 2006):
[1361]1660             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
[1115]1661             accr = MIN( accr, qc_1d(k) / dt_micro )
[1106]1662
[1115]1663             qr_1d(k) = qr_1d(k) + accr * dt_micro 
1664             qc_1d(k) = qc_1d(k) - accr * dt_micro
1665
[1005]1666          ENDIF
[1106]1667
[1005]1668       ENDDO
1669
[1000]1670    END SUBROUTINE accretion_ij
1671
[1005]1672
[1682]1673!------------------------------------------------------------------------------!
1674! Description:
1675! ------------
1676!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
1677!------------------------------------------------------------------------------!
[1005]1678    SUBROUTINE selfcollection_breakup_ij( i, j )
1679
[1320]1680       USE cloud_parameters,                                                   &
[1849]1681           ONLY:  hyrho
[1320]1682
1683       USE control_parameters,                                                 &
[1849]1684           ONLY:  rho_surface
[1320]1685
1686       USE indices,                                                            &
[1822]1687           ONLY:  nzb_s_inner, nzt
[1320]1688
1689       USE kinds
[1005]1690   
1691       IMPLICIT NONE
1692
[1682]1693       INTEGER(iwp) ::  i                 !<
1694       INTEGER(iwp) ::  j                 !<
1695       INTEGER(iwp) ::  k                 !<
[1005]1696
[1682]1697       REAL(wp)     ::  breakup           !<
1698       REAL(wp)     ::  dr                !<
1699       REAL(wp)     ::  phi_br            !<
1700       REAL(wp)     ::  selfcoll          !<
[1320]1701
[1115]1702       DO  k = nzb_s_inner(j,i)+1, nzt
1703          IF ( qr_1d(k) > eps_sb )  THEN
[1012]1704!
[1115]1705!--          Selfcollection rate (Seifert and Beheng, 2001):
[1361]1706             selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
[1012]1707!
[1115]1708!--          Weight averaged diameter of rain drops:
[1334]1709             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]1710!
[1048]1711!--          Collisional breakup rate (Seifert, 2008):
[1353]1712             IF ( dr >= 0.3E-3_wp )  THEN
1713                phi_br  = k_br * ( dr - 1.1E-3_wp )
1714                breakup = selfcoll * ( phi_br + 1.0_wp )
[1005]1715             ELSE
[1353]1716                breakup = 0.0_wp
[1005]1717             ENDIF
[1048]1718
[1115]1719             selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
1720             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
[1106]1721
[1005]1722          ENDIF         
1723       ENDDO
1724
1725    END SUBROUTINE selfcollection_breakup_ij
1726
[1106]1727
[1682]1728!------------------------------------------------------------------------------!
1729! Description:
1730! ------------
1731!> Evaporation of precipitable water. Condensation is neglected for
1732!> precipitable water. Call for grid point i,j
1733!------------------------------------------------------------------------------!
[1012]1734    SUBROUTINE evaporation_rain_ij( i, j )
1735
[1320]1736       USE arrays_3d,                                                          &
[1849]1737           ONLY:  hyp
[1048]1738
[1320]1739       USE cloud_parameters,                                                   &
[1849]1740           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
[1320]1741
1742       USE constants,                                                          &
1743           ONLY:  pi
1744
1745       USE indices,                                                            &
[1822]1746           ONLY:  nzb_s_inner, nzt
[1320]1747
1748       USE kinds
1749
[1012]1750       IMPLICIT NONE
1751
[1682]1752       INTEGER(iwp) ::  i                 !<
1753       INTEGER(iwp) ::  j                 !<
1754       INTEGER(iwp) ::  k                 !<
[1012]1755
[1682]1756       REAL(wp)     ::  alpha             !<
1757       REAL(wp)     ::  dr                !<
1758       REAL(wp)     ::  e_s               !<
1759       REAL(wp)     ::  evap              !<
1760       REAL(wp)     ::  evap_nr           !<
1761       REAL(wp)     ::  f_vent            !<
1762       REAL(wp)     ::  g_evap            !<
1763       REAL(wp)     ::  lambda_r          !<
1764       REAL(wp)     ::  mu_r              !<
1765       REAL(wp)     ::  mu_r_2            !<
1766       REAL(wp)     ::  mu_r_5d2          !<
1767       REAL(wp)     ::  nr_0              !<
1768       REAL(wp)     ::  q_s               !<
1769       REAL(wp)     ::  sat               !<
1770       REAL(wp)     ::  t_l               !<
1771       REAL(wp)     ::  temp              !<
1772       REAL(wp)     ::  xr                !<
[1320]1773
[1115]1774       DO  k = nzb_s_inner(j,i)+1, nzt
1775          IF ( qr_1d(k) > eps_sb )  THEN
[1012]1776!
1777!--          Actual liquid water temperature:
[1115]1778             t_l = t_d_pt(k) * pt_1d(k)
[1012]1779!
1780!--          Saturation vapor pressure at t_l:
[1361]1781             e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /          &
1782                                    ( t_l - 35.86_wp )                         &
1783                                  )
[1012]1784!
1785!--          Computation of saturation humidity:
[1361]1786             q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
[1353]1787             alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
[1361]1788             q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
[1012]1789!
[1106]1790!--          Supersaturation:
[1361]1791             sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
[1012]1792!
[1361]1793!--          Evaporation needs only to be calculated in subsaturated regions
1794             IF ( sat < 0.0_wp )  THEN
[1012]1795!
[1361]1796!--             Actual temperature:
1797                temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
1798       
1799                g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /  &
1800                                    ( thermal_conductivity_l * temp ) +        &
1801                                    r_v * temp / ( diff_coeff_l * e_s )        &
1802                                  )
[1012]1803!
[1361]1804!--             Mean weight of rain drops
1805                xr = hyrho(k) * qr_1d(k) / nr_1d(k)
[1115]1806!
[1361]1807!--             Weight averaged diameter of rain drops:
1808                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]1809!
[1361]1810!--             Compute ventilation factor and intercept parameter
1811!--             (Seifert and Beheng, 2006; Seifert, 2008):
1812                IF ( ventilation_effect )  THEN
[1115]1813!
[1361]1814!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1815!--                Stevens and Seifert, 2008):
1816                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
1817!
1818!--                Slope parameter of gamma distribution (Seifert, 2008):
1819                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
1820                                ( mu_r + 1.0_wp )                              &
1821                              )**( 1.0_wp / 3.0_wp ) / dr
[1115]1822
[1361]1823                   mu_r_2   = mu_r + 2.0_wp
1824                   mu_r_5d2 = mu_r + 2.5_wp 
1825
1826                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  & 
1827                            b_vent * schmidt_p_1d3 *                           &
1828                            SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
1829                            lambda_r**( -mu_r_5d2 ) *                          &
1830                            ( 1.0_wp -                                         &
1831                              0.5_wp * ( b_term / a_term ) *                   &
1832                              ( lambda_r / ( c_term + lambda_r )               &
1833                              )**mu_r_5d2 -                                    &
1834                              0.125_wp * ( b_term / a_term )**2 *              &
1835                              ( lambda_r / ( 2.0_wp * c_term + lambda_r )      &
1836                              )**mu_r_5d2 -                                    &
1837                              0.0625_wp * ( b_term / a_term )**3 *             &
1838                              ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
1839                              )**mu_r_5d2 -                                    &
1840                              0.0390625_wp * ( b_term / a_term )**4 *          & 
1841                              ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
1842                              )**mu_r_5d2                                      &
1843                            )
1844
1845                   nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) /           &
1846                            gamm( mu_r + 1.0_wp ) 
1847                ELSE
1848                   f_vent = 1.0_wp
1849                   nr_0   = nr_1d(k) * dr
1850                ENDIF
[1012]1851!
[1361]1852!--             Evaporation rate of rain water content (Seifert and Beheng, 2006):
1853                evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
1854                evap    = MAX( evap, -qr_1d(k) / dt_micro )
1855                evap_nr = MAX( c_evap * evap / xr * hyrho(k),                  &
1856                               -nr_1d(k) / dt_micro )
[1106]1857
[1361]1858                qr_1d(k) = qr_1d(k) + evap    * dt_micro
1859                nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
[1115]1860
[1361]1861             ENDIF
[1012]1862          ENDIF         
[1106]1863
[1012]1864       ENDDO
1865
1866    END SUBROUTINE evaporation_rain_ij
1867
[1106]1868
[1682]1869!------------------------------------------------------------------------------!
1870! Description:
1871! ------------
1872!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
1873!> Call for grid point i,j
1874!------------------------------------------------------------------------------!
[1012]1875    SUBROUTINE sedimentation_cloud_ij( i, j )
1876
[1320]1877       USE arrays_3d,                                                          &
[1849]1878           ONLY:  ddzu, dzu, prr
[1320]1879
1880       USE cloud_parameters,                                                   &
[1849]1881           ONLY:  hyrho, l_d_cp, pt_d_t
[1320]1882
1883       USE control_parameters,                                                 &
[1849]1884           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
[1320]1885
1886       USE indices,                                                            &
1887           ONLY:  nzb, nzb_s_inner, nzt
1888
1889       USE kinds
[1012]1890       
[1691]1891       USE statistics,                                                         &
1892           ONLY:  weight_substep
1893
[1012]1894       IMPLICIT NONE
1895
[1849]1896       INTEGER(iwp) ::  i             !<
1897       INTEGER(iwp) ::  j             !<
1898       INTEGER(iwp) ::  k             !<
[1106]1899
[1682]1900       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc  !<
[1115]1901
[1353]1902       sed_qc(nzt+1) = 0.0_wp
[1012]1903
[1115]1904       DO  k = nzt, nzb_s_inner(j,i)+1, -1
1905          IF ( qc_1d(k) > eps_sb )  THEN
[1361]1906             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) *       &
1907                         ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
[1115]1908          ELSE
[1353]1909             sed_qc(k) = 0.0_wp
[1012]1910          ENDIF
[1115]1911
[1361]1912          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /          &
1913                                      dt_micro + sed_qc(k+1)                   &
1914                         )
[1115]1915
[1361]1916          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
[1115]1917                                hyrho(k) * dt_micro
[1361]1918          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      & 
[1115]1919                                hyrho(k) * dt_micro
[1361]1920          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
[1115]1921                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
1922
[1691]1923!
1924!--       Compute the precipitation rate of cloud (fog) droplets
[1822]1925          IF ( call_microphysics_at_all_substeps )  THEN
1926             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *               &
[1691]1927                             weight_substep(intermediate_timestep_count)
[1822]1928          ELSE
1929             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
[1691]1930          ENDIF
1931
[1012]1932       ENDDO
1933
1934    END SUBROUTINE sedimentation_cloud_ij
1935
[1106]1936
[1682]1937!------------------------------------------------------------------------------!
1938! Description:
1939! ------------
1940!> Computation of sedimentation flux. Implementation according to Stevens
1941!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
1942!------------------------------------------------------------------------------!
[1012]1943    SUBROUTINE sedimentation_rain_ij( i, j )
1944
[1320]1945       USE arrays_3d,                                                          &
[1849]1946           ONLY:  ddzu, dzu, prr
[1320]1947
1948       USE cloud_parameters,                                                   &
[1849]1949           ONLY:  hyrho, l_d_cp, pt_d_t
[1320]1950
1951       USE control_parameters,                                                 &
[1849]1952           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
[1320]1953
1954       USE indices,                                                            &
1955           ONLY:  nzb, nzb_s_inner, nzt
1956
1957       USE kinds
1958
1959       USE statistics,                                                         &
1960           ONLY:  weight_substep
[1012]1961       
1962       IMPLICIT NONE
1963
[1682]1964       INTEGER(iwp) ::  i                          !<
1965       INTEGER(iwp) ::  j                          !<
1966       INTEGER(iwp) ::  k                          !<
1967       INTEGER(iwp) ::  k_run                      !<
[1012]1968
[1682]1969       REAL(wp)     ::  c_run                      !<
1970       REAL(wp)     ::  d_max                      !<
1971       REAL(wp)     ::  d_mean                     !<
1972       REAL(wp)     ::  d_min                      !<
1973       REAL(wp)     ::  dr                         !<
1974       REAL(wp)     ::  flux                       !<
1975       REAL(wp)     ::  lambda_r                   !<
1976       REAL(wp)     ::  mu_r                       !<
1977       REAL(wp)     ::  z_run                      !<
[1320]1978
[1682]1979       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
1980       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
1981       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
1982       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
1983       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
1984       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
1985       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
1986       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
[1320]1987
[1012]1988!
[1065]1989!--    Compute velocities
1990       DO  k = nzb_s_inner(j,i)+1, nzt
[1115]1991          IF ( qr_1d(k) > eps_sb )  THEN
1992!
1993!--          Weight averaged diameter of rain drops:
[1334]1994             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]1995!
1996!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1997!--          Stevens and Seifert, 2008):
[1353]1998             mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
[1115]1999!
2000!--          Slope parameter of gamma distribution (Seifert, 2008):
[1361]2001             lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *              &
2002                          ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1115]2003
[1361]2004             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
2005                                         a_term - b_term * ( 1.0_wp +          &
2006                                            c_term / lambda_r )**( -1.0_wp *   &
2007                                            ( mu_r + 1.0_wp ) )                &
2008                                        )                                      &
2009                          )
2010             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
2011                                         a_term - b_term * ( 1.0_wp +          &
2012                                            c_term / lambda_r )**( -1.0_wp *   &
2013                                            ( mu_r + 4.0_wp ) )                &
2014                                       )                                       &
2015                          )
[1065]2016          ELSE
[1353]2017             w_nr(k) = 0.0_wp
2018             w_qr(k) = 0.0_wp
[1065]2019          ENDIF
2020       ENDDO
[1048]2021!
[1065]2022!--    Adjust boundary values
[1115]2023       w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
2024       w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
[1353]2025       w_nr(nzt+1) = 0.0_wp
2026       w_qr(nzt+1) = 0.0_wp
[1065]2027!
2028!--    Compute Courant number
[1115]2029       DO  k = nzb_s_inner(j,i)+1, nzt
[1361]2030          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
[1115]2031                    dt_micro * ddzu(k)
[1361]2032          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
[1115]2033                    dt_micro * ddzu(k)
2034       ENDDO     
[1065]2035!
2036!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
2037       IF ( limiter_sedimentation )  THEN
2038
[1115]2039          DO k = nzb_s_inner(j,i)+1, nzt
[1646]2040             d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
[1115]2041             d_min  = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
2042             d_max  = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
[1065]2043
[1361]2044             qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
2045                                                        2.0_wp * d_max,        &
2046                                                        ABS( d_mean ) )
[1065]2047
[1646]2048             d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
[1115]2049             d_min  = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
2050             d_max  = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
[1065]2051
[1361]2052             nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
2053                                                        2.0_wp * d_max,        &
2054                                                        ABS( d_mean ) )
[1022]2055          ENDDO
[1048]2056
[1065]2057       ELSE
[1106]2058
[1353]2059          nr_slope = 0.0_wp
2060          qr_slope = 0.0_wp
[1106]2061
[1065]2062       ENDIF
[1115]2063
[1353]2064       sed_nr(nzt+1) = 0.0_wp
2065       sed_qr(nzt+1) = 0.0_wp
[1065]2066!
2067!--    Compute sedimentation flux
[1115]2068       DO  k = nzt, nzb_s_inner(j,i)+1, -1
[1065]2069!
2070!--       Sum up all rain drop number densities which contribute to the flux
2071!--       through k-1/2
[1353]2072          flux  = 0.0_wp
2073          z_run = 0.0_wp ! height above z(k)
[1065]2074          k_run = k
[1346]2075          c_run = MIN( 1.0_wp, c_nr(k) )
[1353]2076          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
[1361]2077             flux  = flux + hyrho(k_run) *                                     &
2078                     ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) *   &
[1353]2079                     0.5_wp ) * c_run * dzu(k_run)
[1065]2080             z_run = z_run + dzu(k_run)
2081             k_run = k_run + 1
[1346]2082             c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
[1022]2083          ENDDO
2084!
[1065]2085!--       It is not allowed to sediment more rain drop number density than
2086!--       available
[1361]2087          flux = MIN( flux,                                                    &
[1115]2088                      hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
[1065]2089
[1115]2090          sed_nr(k) = flux / dt_micro
[1361]2091          nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
2092                                    hyrho(k) * dt_micro
[1065]2093!
2094!--       Sum up all rain water content which contributes to the flux
2095!--       through k-1/2
[1353]2096          flux  = 0.0_wp
2097          z_run = 0.0_wp ! height above z(k)
[1065]2098          k_run = k
[1346]2099          c_run = MIN( 1.0_wp, c_qr(k) )
[1106]2100
[1361]2101          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
[1106]2102
[1361]2103             flux  = flux + hyrho(k_run) *                                     &
2104                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) *   &
[1353]2105                     0.5_wp ) * c_run * dzu(k_run)
[1065]2106             z_run = z_run + dzu(k_run)
2107             k_run = k_run + 1
[1346]2108             c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
[1106]2109
[1065]2110          ENDDO
2111!
2112!--       It is not allowed to sediment more rain water content than available
[1361]2113          flux = MIN( flux,                                                    &
[1115]2114                      hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
[1065]2115
[1115]2116          sed_qr(k) = flux / dt_micro
2117
[1361]2118          qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
[1115]2119                                hyrho(k) * dt_micro
[1361]2120          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
[1115]2121                                hyrho(k) * dt_micro 
[1361]2122          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
[1115]2123                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
[1065]2124!
2125!--       Compute the rain rate
[1361]2126          IF ( call_microphysics_at_all_substeps )  THEN
[1691]2127             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)                    &
2128                          * weight_substep(intermediate_timestep_count)
[1361]2129          ELSE
[1691]2130             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
[1361]2131          ENDIF
2132
[1065]2133       ENDDO
[1115]2134
[1691]2135    END SUBROUTINE sedimentation_rain_ij
[1012]2136
[1691]2137
2138!------------------------------------------------------------------------------!
2139! Description:
2140! ------------
2141!> This subroutine computes the precipitation amount due to gravitational
2142!> settling of rain and cloud (fog) droplets
2143!------------------------------------------------------------------------------!
2144    SUBROUTINE calc_precipitation_amount_ij( i, j )
2145
[1849]2146       USE arrays_3d,                                                          &
2147           ONLY:  precipitation_amount, prr
2148
[1691]2149       USE cloud_parameters,                                                   &
[1849]2150           ONLY:  hyrho
[1691]2151
2152       USE control_parameters,                                                 &
2153           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
2154                  intermediate_timestep_count, intermediate_timestep_count_max,&
[1822]2155                  precipitation_amount_interval, time_do2d_xy
[1691]2156
2157       USE indices,                                                            &
2158           ONLY:  nzb_s_inner
2159
2160       USE kinds
2161
2162       IMPLICIT NONE
2163
2164       INTEGER(iwp) ::  i                          !:
2165       INTEGER(iwp) ::  j                          !:
2166
2167
2168       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
2169            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
2170            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
2171       THEN
2172
[1361]2173          precipitation_amount(j,i) = precipitation_amount(j,i) +              &
2174                                      prr(nzb_s_inner(j,i)+1,j,i) *            &
[1115]2175                                      hyrho(nzb_s_inner(j,i)+1) * dt_3d
[1048]2176       ENDIF
2177
[1691]2178    END SUBROUTINE calc_precipitation_amount_ij
[1012]2179
[1361]2180!------------------------------------------------------------------------------!
[1682]2181! Description:
2182! ------------
2183!> This function computes the gamma function (Press et al., 1992).
2184!> The gamma function is needed for the calculation of the evaporation
2185!> of rain drops.
[1361]2186!------------------------------------------------------------------------------!
[1012]2187    FUNCTION gamm( xx ) 
[1320]2188
2189       USE kinds
2190
[1012]2191       IMPLICIT NONE
[1106]2192
[1682]2193       INTEGER(iwp) ::  j            !<
[1320]2194
[1682]2195       REAL(wp)     ::  gamm         !<
2196       REAL(wp)     ::  ser          !<
2197       REAL(wp)     ::  tmp          !<
2198       REAL(wp)     ::  x_gamm       !<
2199       REAL(wp)     ::  xx           !<
2200       REAL(wp)     ::  y_gamm       !<
[1320]2201
[1849]2202
2203       REAL(wp), PARAMETER  ::  stp = 2.5066282746310005_wp               !<
2204       REAL(wp), PARAMETER  ::  cof(6) = (/ 76.18009172947146_wp,      &
2205                                           -86.50532032941677_wp,      &
2206                                            24.01409824083091_wp,      &
2207                                            -1.231739572450155_wp,     &
2208                                             0.1208650973866179E-2_wp, &
2209                                            -0.5395239384953E-5_wp /)     !<
2210
2211       x_gamm = xx
2212       y_gamm = x_gamm
[1353]2213       tmp = x_gamm + 5.5_wp
[1849]2214       tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
[1334]2215       ser = 1.000000000190015_wp
[1106]2216
2217       DO  j = 1, 6 
[1353]2218          y_gamm = y_gamm + 1.0_wp 
[1012]2219          ser    = ser + cof( j ) / y_gamm 
[1106]2220       ENDDO
2221
[1012]2222!
2223!--    Until this point the algorithm computes the logarithm of the gamma
2224!--    function. Hence, the exponential function is used. 
2225!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
2226       gamm = EXP( tmp ) * stp * ser / x_gamm 
[1106]2227
[1012]2228       RETURN
2229 
2230    END FUNCTION gamm 
2231
2232 END MODULE microphysics_mod
Note: See TracBrowser for help on using the repository browser.