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

Last change on this file since 2215 was 2156, checked in by hoffmann, 7 years ago

last commit documented

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