source: palm/trunk/SOURCE/microphysics.f90 @ 1362

Last change on this file since 1362 was 1362, checked in by hoffmann, 10 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 72.2 KB
Line 
1 MODULE microphysics_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: microphysics.f90 1362 2014-04-16 15:19:12Z hoffmann $
27!
28! 1361 2014-04-16 15:17:48Z hoffmann
29! Bugfix in sedimentation_rain: Index corrected.
30! Vectorized version of adjust_cloud added.
31! Little reformatting of the code.
32!
33! 1353 2014-04-08 15:21:23Z heinze
34! REAL constants provided with KIND-attribute
35!
36! 1346 2014-03-27 13:18:20Z heinze
37! Bugfix: REAL constants provided with KIND-attribute especially in call of
38! intrinsic function like MAX, MIN, SIGN
39!
40! 1334 2014-03-25 12:21:40Z heinze
41! Bugfix: REAL constants provided with KIND-attribute
42!
43! 1322 2014-03-20 16:38:49Z raasch
44! REAL constants defined as wp-kind
45!
46! 1320 2014-03-20 08:40:49Z raasch
47! ONLY-attribute added to USE-statements,
48! kind-parameters added to all INTEGER and REAL declaration statements,
49! kinds are defined in new module kinds,
50! comment fields (!:) to be used for variable explanations added to
51! all variable declaration statements
52!
53! 1241 2013-10-30 11:36:58Z heinze
54! hyp and rho have to be calculated at each time step if data from external
55! file LSF_DATA are used
56!
57! 1115 2013-03-26 18:16:16Z hoffmann
58! microphyical tendencies are calculated in microphysics_control in an optimized
59! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
60!
61! 1106 2013-03-04 05:31:38Z raasch
62! small changes in code formatting
63!
64! 1092 2013-02-02 11:24:22Z raasch
65! unused variables removed
66! file put under GPL
67!
68! 1065 2012-11-22 17:42:36Z hoffmann
69! Sedimentation process implemented according to Stevens and Seifert (2008).
70! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
71! and Stevens, 2010).
72!
73! 1053 2012-11-13 17:11:03Z hoffmann
74! initial revision
75!
76! Description:
77! ------------
78! Calculate cloud microphysics according to the two moment bulk
79! scheme by Seifert and Beheng (2006).
80!------------------------------------------------------------------------------!
81
82    PRIVATE
83    PUBLIC microphysics_control
84
85    INTERFACE microphysics_control
86       MODULE PROCEDURE microphysics_control
87       MODULE PROCEDURE microphysics_control_ij
88    END INTERFACE microphysics_control
89
90    INTERFACE adjust_cloud
91       MODULE PROCEDURE adjust_cloud
92       MODULE PROCEDURE adjust_cloud_ij
93    END INTERFACE adjust_cloud
94
95    INTERFACE autoconversion
96       MODULE PROCEDURE autoconversion
97       MODULE PROCEDURE autoconversion_ij
98    END INTERFACE autoconversion
99
100    INTERFACE accretion
101       MODULE PROCEDURE accretion
102       MODULE PROCEDURE accretion_ij
103    END INTERFACE accretion
104
105    INTERFACE selfcollection_breakup
106       MODULE PROCEDURE selfcollection_breakup
107       MODULE PROCEDURE selfcollection_breakup_ij
108    END INTERFACE selfcollection_breakup
109
110    INTERFACE evaporation_rain
111       MODULE PROCEDURE evaporation_rain
112       MODULE PROCEDURE evaporation_rain_ij
113    END INTERFACE evaporation_rain
114
115    INTERFACE sedimentation_cloud
116       MODULE PROCEDURE sedimentation_cloud
117       MODULE PROCEDURE sedimentation_cloud_ij
118    END INTERFACE sedimentation_cloud
119 
120    INTERFACE sedimentation_rain
121       MODULE PROCEDURE sedimentation_rain
122       MODULE PROCEDURE sedimentation_rain_ij
123    END INTERFACE sedimentation_rain
124
125 CONTAINS
126
127
128!------------------------------------------------------------------------------!
129! Call for all grid points
130!------------------------------------------------------------------------------!
131    SUBROUTINE microphysics_control
132
133       USE arrays_3d,                                                          &
134           ONLY:  hyp, nr, pt, pt_init, q, qc, qr, zu
135
136       USE cloud_parameters,                                                   &
137           ONLY:  cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
138
139       USE control_parameters,                                                 &
140           ONLY:  call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
141                  g, intermediate_timestep_count,                              &
142                  large_scale_forcing, lsf_surf, precipitation, pt_surface,    &
143                  rho_surface,surface_pressure
144
145       USE indices,                                                            &
146           ONLY:  nzb, nzt
147
148       USE kinds
149
150       USE statistics,                                                         &
151           ONLY:  weight_pres
152
153       IMPLICIT NONE
154
155       INTEGER(iwp) ::  i                 !:
156       INTEGER(iwp) ::  j                 !:
157       INTEGER(iwp) ::  k                 !:
158
159       REAL(wp)     ::  t_surface         !:
160
161       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
162!
163!--       Calculate:
164!--       pt / t : ratio of potential and actual temperature (pt_d_t)
165!--       t / pt : ratio of actual and potential temperature (t_d_pt)
166!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
167          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
168          DO  k = nzb, nzt+1
169             hyp(k)    = surface_pressure * 100.0_wp * &
170                         ( ( t_surface - g / cp * zu(k) ) /                    &
171                         t_surface )**(1.0_wp / 0.286_wp)
172             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
173             t_d_pt(k) = 1.0_wp / pt_d_t(k)
174             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
175          ENDDO
176!
177!--       Compute reference density
178          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
179       ENDIF
180
181!
182!--    Compute length of time step
183       IF ( call_microphysics_at_all_substeps )  THEN
184          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
185       ELSE
186          dt_micro = dt_3d
187       ENDIF
188
189!
190!--    Compute cloud physics
191       IF ( precipitation )  THEN
192          CALL adjust_cloud
193          CALL autoconversion
194          CALL accretion
195          CALL selfcollection_breakup
196          CALL evaporation_rain
197          CALL sedimentation_rain
198       ENDIF
199
200       IF ( drizzle )  CALL sedimentation_cloud
201
202    END SUBROUTINE microphysics_control
203
204    SUBROUTINE adjust_cloud
205
206       USE arrays_3d,                                                          &
207           ONLY:  qr, nr
208
209       USE cloud_parameters,                                                   &
210           ONLY:  eps_sb, xrmin, xrmax, hyrho, k_cc, x0
211
212       USE cpulog,                                                             &
213           ONLY:  cpu_log, log_point_s
214
215       USE indices,                                                            &
216           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
217
218       USE kinds
219
220       IMPLICIT NONE
221
222       INTEGER(iwp) ::  i                 !:
223       INTEGER(iwp) ::  j                 !:
224       INTEGER(iwp) ::  k                 !:
225
226       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
227
228       DO  i = nxl, nxr
229          DO  j = nys, nyn
230             DO  k = nzb_s_inner(j,i)+1, nzt
231                IF ( qr(k,j,i) <= eps_sb )  THEN
232                   qr(k,j,i) = 0.0_wp
233                   nr(k,j,i) = 0.0_wp
234                ELSE
235!
236!--                Adjust number of raindrops to avoid nonlinear effects in
237!--                sedimentation and evaporation of rain drops due to too small
238!--                or too big weights of rain drops (Stevens and Seifert, 2008).
239                   IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
240                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
241                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
242                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
243                   ENDIF
244
245                ENDIF
246             ENDDO
247          ENDDO
248       ENDDO
249
250       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
251
252    END SUBROUTINE adjust_cloud
253
254
255    SUBROUTINE autoconversion
256
257       USE arrays_3d,                                                          &
258           ONLY:  diss, dzu, nr, qc, qr
259
260       USE cloud_parameters,                                                   &
261           ONLY:  a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3,        &
262                  c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air,         &
263                  nc_const, x0
264
265       USE control_parameters,                                                 &
266           ONLY:  dt_micro, rho_surface, turbulence
267
268       USE cpulog,                                                             &
269           ONLY:  cpu_log, log_point_s
270
271       USE grid_variables,                                                     &
272           ONLY:  dx, dy
273
274       USE indices,                                                            &
275           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
276
277       USE kinds
278
279       IMPLICIT NONE
280
281       INTEGER(iwp) ::  i                 !:
282       INTEGER(iwp) ::  j                 !:
283       INTEGER(iwp) ::  k                 !:
284
285       REAL(wp)     ::  alpha_cc          !:                   
286       REAL(wp)     ::  autocon           !:
287       REAL(wp)     ::  dissipation       !:
288       REAL(wp)     ::  k_au              !:
289       REAL(wp)     ::  l_mix             !:
290       REAL(wp)     ::  nu_c              !:
291       REAL(wp)     ::  phi_au            !:
292       REAL(wp)     ::  r_cc              !:
293       REAL(wp)     ::  rc                !:
294       REAL(wp)     ::  re_lambda         !:
295       REAL(wp)     ::  selfcoll          !:
296       REAL(wp)     ::  sigma_cc          !:
297       REAL(wp)     ::  tau_cloud         !:
298       REAL(wp)     ::  xc                !:
299
300       CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
301
302       DO  i = nxl, nxr
303          DO  j = nys, nyn
304             DO  k = nzb_s_inner(j,i)+1, nzt
305
306                IF ( qc(k,j,i) > eps_sb )  THEN
307
308                   k_au = k_cc / ( 20.0_wp * x0 )
309!
310!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
311!--                (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
312                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
313!
314!--                Universal function for autoconversion process
315!--                (Seifert and Beheng, 2006):
316                   phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
317                            ( 1.0_wp - tau_cloud**0.68_wp )**3
318!
319!--                Shape parameter of gamma distribution (Geoffroy et al., 2010):
320!--                (Use constant nu_c = 1.0_wp instead?)
321                   nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
322!
323!--                Mean weight of cloud droplets:
324                   xc = hyrho(k) * qc(k,j,i) / nc_const
325!
326!--                Parameterized turbulence effects on autoconversion (Seifert,
327!--                Nuijens and Stevens, 2010)
328                   IF ( turbulence )  THEN
329!
330!--                   Weight averaged radius of cloud droplets:
331                      rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
332
333                      alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
334                      r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
335                      sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
336!
337!--                   Mixing length (neglecting distance to ground and
338!--                   stratification)
339                      l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
340!
341!--                   Limit dissipation rate according to Seifert, Nuijens and
342!--                   Stevens (2010)
343                      dissipation = MIN( 0.06_wp, diss(k,j,i) )
344!
345!--                   Compute Taylor-microscale Reynolds number:
346                      re_lambda = 6.0_wp / 11.0_wp *                           &
347                                  ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *   &
348                                  SQRT( 15.0_wp / kin_vis_air ) *              &
349                                  dissipation**( 1.0_wp / 6.0_wp )
350!
351!--                   The factor of 1.0E4 is needed to convert the dissipation
352!--                   rate from m2 s-3 to cm2 s-3.
353                      k_au = k_au * ( 1.0_wp +                                 &
354                                      dissipation * 1.0E4_wp *                 &
355                                      ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
356                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &     
357                                                                      r_cc ) / &
358                                                        sigma_cc )**2          &
359                                                      ) + beta_cc              &
360                                      )                                        &
361                                    )
362                   ENDIF
363!
364!--                Autoconversion rate (Seifert and Beheng, 2006):
365                   autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /    &
366                             ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *     &
367                             ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * &
368                             rho_surface
369                   autocon = MIN( autocon, qc(k,j,i) / dt_micro )
370
371                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro
372                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro 
373                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro
374
375                ENDIF
376
377             ENDDO
378          ENDDO
379       ENDDO
380
381       CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
382
383    END SUBROUTINE autoconversion
384
385
386    SUBROUTINE accretion
387
388       USE arrays_3d,                                                          &
389           ONLY:  diss, qc, qr
390
391       USE cloud_parameters,                                                   &
392           ONLY:  eps_sb, hyrho, k_cr0
393
394       USE control_parameters,                                                 &
395           ONLY:  dt_micro, rho_surface, turbulence
396
397       USE cpulog,                                                             &
398           ONLY:  cpu_log, log_point_s
399
400       USE indices,                                                            &
401           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
402
403       USE kinds
404
405       IMPLICIT NONE
406
407       INTEGER(iwp) ::  i                 !:
408       INTEGER(iwp) ::  j                 !:
409       INTEGER(iwp) ::  k                 !:
410
411       REAL(wp)     ::  accr              !:
412       REAL(wp)     ::  k_cr              !:
413       REAL(wp)     ::  phi_ac            !:
414       REAL(wp)     ::  tau_cloud         !:
415       REAL(wp)     ::  xc                !:
416
417       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
418
419       DO  i = nxl, nxr
420          DO  j = nys, nyn
421             DO  k = nzb_s_inner(j,i)+1, nzt
422
423                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb ) )  THEN
424!
425!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
426                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ) 
427!
428!--                Universal function for accretion process (Seifert and
429!--                Beheng, 2001):
430                   phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
431!
432!--                Parameterized turbulence effects on autoconversion (Seifert,
433!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
434!--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
435                   IF ( turbulence )  THEN
436                      k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                      &
437                                       MIN( 600.0_wp,                          &
438                                            diss(k,j,i) * 1.0E4_wp )**0.25_wp  &
439                                     )
440                   ELSE
441                      k_cr = k_cr0                       
442                   ENDIF
443!
444!--                Accretion rate (Seifert and Beheng, 2006):
445                   accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *              &
446                          SQRT( rho_surface * hyrho(k) )
447                   accr = MIN( accr, qc(k,j,i) / dt_micro )
448
449                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro 
450                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro
451
452                ENDIF
453
454             ENDDO
455          ENDDO
456       ENDDO
457
458       CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
459
460    END SUBROUTINE accretion
461
462
463    SUBROUTINE selfcollection_breakup
464
465       USE arrays_3d,                                                          &
466           ONLY:  nr, qr
467
468       USE cloud_parameters,                                                   &
469           ONLY:  dpirho_l, eps_sb, hyrho, k_br, k_rr
470
471       USE control_parameters,                                                 &
472           ONLY:  dt_micro, rho_surface
473
474       USE cpulog,                                                             &
475           ONLY:  cpu_log, log_point_s
476
477       USE indices,                                                            &
478           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
479
480       USE kinds
481   
482       IMPLICIT NONE
483
484       INTEGER(iwp) ::  i                 !:
485       INTEGER(iwp) ::  j                 !:
486       INTEGER(iwp) ::  k                 !:
487
488       REAL(wp)     ::  breakup           !:
489       REAL(wp)     ::  dr                !:
490       REAL(wp)     ::  phi_br            !:
491       REAL(wp)     ::  selfcoll          !:
492
493       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
494
495       DO  i = nxl, nxr
496          DO  j = nys, nyn
497             DO  k = nzb_s_inner(j,i)+1, nzt
498                IF ( qr(k,j,i) > eps_sb )  THEN
499!
500!--                Selfcollection rate (Seifert and Beheng, 2001):
501                   selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) *                   &
502                              SQRT( hyrho(k) * rho_surface )
503!
504!--                Weight averaged diameter of rain drops:
505                   dr = ( hyrho(k) * qr(k,j,i) /                               &
506                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
507!
508!--                Collisional breakup rate (Seifert, 2008):
509                   IF ( dr >= 0.3E-3_wp )  THEN
510                      phi_br  = k_br * ( dr - 1.1E-3_wp )
511                      breakup = selfcoll * ( phi_br + 1.0_wp )
512                   ELSE
513                      breakup = 0.0_wp
514                   ENDIF
515
516                   selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
517                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro
518
519                ENDIF         
520             ENDDO
521          ENDDO
522       ENDDO
523
524       CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
525
526    END SUBROUTINE selfcollection_breakup
527
528
529    SUBROUTINE evaporation_rain
530
531!
532!--    Evaporation of precipitable water. Condensation is neglected for
533!--    precipitable water.
534
535       USE arrays_3d,                                                          &
536           ONLY:  hyp, nr, pt, q,  qc, qr
537
538       USE cloud_parameters,                                                   &
539           ONLY:  a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
540                  dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r,   &
541                  l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l,      &
542                  t_d_pt, ventilation_effect
543
544       USE constants,                                                          &
545           ONLY:  pi
546
547       USE control_parameters,                                                 &
548           ONLY:  dt_micro
549
550       USE cpulog,                                                             &
551           ONLY:  cpu_log, log_point_s
552
553       USE indices,                                                            &
554           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
555
556       USE kinds
557
558       IMPLICIT NONE
559
560       INTEGER(iwp) ::  i                 !:
561       INTEGER(iwp) ::  j                 !:
562       INTEGER(iwp) ::  k                 !:
563
564       REAL(wp)     ::  alpha             !:
565       REAL(wp)     ::  dr                !:
566       REAL(wp)     ::  e_s               !:
567       REAL(wp)     ::  evap              !:
568       REAL(wp)     ::  evap_nr           !:
569       REAL(wp)     ::  f_vent            !:
570       REAL(wp)     ::  g_evap            !:
571       REAL(wp)     ::  lambda_r          !:
572       REAL(wp)     ::  mu_r              !:
573       REAL(wp)     ::  mu_r_2            !:
574       REAL(wp)     ::  mu_r_5d2          !:
575       REAL(wp)     ::  nr_0              !:
576       REAL(wp)     ::  q_s               !:
577       REAL(wp)     ::  sat               !:
578       REAL(wp)     ::  t_l               !:
579       REAL(wp)     ::  temp              !:
580       REAL(wp)     ::  xr                !:
581
582       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
583
584       DO  i = nxl, nxr
585          DO  j = nys, nyn
586             DO  k = nzb_s_inner(j,i)+1, nzt
587                IF ( qr(k,j,i) > eps_sb )  THEN
588!
589!--                Actual liquid water temperature:
590                   t_l = t_d_pt(k) * pt(k,j,i)
591!
592!--                Saturation vapor pressure at t_l:
593                   e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
594                                          ( t_l - 35.86_wp )                   &
595                                        )
596!
597!--                Computation of saturation humidity:
598                   q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
599                   alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
600                   q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /               &
601                           ( 1.0_wp + alpha * q_s )
602!
603!--                Supersaturation:
604                   sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
605!
606!--                Evaporation needs only to be calculated in subsaturated regions
607                   IF ( sat < 0.0_wp )  THEN
608!
609!--                   Actual temperature:
610                      temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
611             
612                      g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *     &
613                                          l_v / ( thermal_conductivity_l * temp ) &
614                                          + r_v * temp / ( diff_coeff_l * e_s )   &
615                                        )
616!
617!--                   Mean weight of rain drops
618                      xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
619!
620!--                   Weight averaged diameter of rain drops:
621                      dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
622!
623!--                   Compute ventilation factor and intercept parameter
624!--                   (Seifert and Beheng, 2006; Seifert, 2008):
625                      IF ( ventilation_effect )  THEN
626!
627!--                      Shape parameter of gamma distribution (Milbrandt and Yau,
628!--                      2005; Stevens and Seifert, 2008):
629                         mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *          &
630                                                          ( dr - 1.4E-3_wp ) ) )
631!
632!--                      Slope parameter of gamma distribution (Seifert, 2008):
633                         lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *  &
634                                      ( mu_r + 1.0_wp )                        &
635                                    )**( 1.0_wp / 3.0_wp ) / dr
636
637                         mu_r_2   = mu_r + 2.0_wp
638                         mu_r_5d2 = mu_r + 2.5_wp 
639
640                         f_vent = a_vent * gamm( mu_r_2 ) *                     &
641                                  lambda_r**( -mu_r_2 ) + b_vent *              &
642                                  schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
643                                  gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) *  &
644                                  ( 1.0_wp -                                    &
645                                    0.5_wp * ( b_term / a_term ) *              &
646                                    ( lambda_r / ( c_term + lambda_r )          &
647                                    )**mu_r_5d2 -                               &
648                                    0.125_wp * ( b_term / a_term )**2 *         &
649                                    ( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
650                                    )**mu_r_5d2 -                               &
651                                    0.0625_wp * ( b_term / a_term )**3 *        &
652                                    ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
653                                    )**mu_r_5d2 -                               &
654                                    0.0390625_wp * ( b_term / a_term )**4 *     & 
655                                    ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
656                                    )**mu_r_5d2                                 &
657                                  )
658
659                         nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
660                                  gamm( mu_r + 1.0_wp ) 
661                      ELSE
662                         f_vent = 1.0_wp
663                         nr_0   = nr(k,j,i) * dr
664                      ENDIF
665   !
666   !--                Evaporation rate of rain water content (Seifert and
667   !--                Beheng, 2006):
668                      evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
669                                hyrho(k)
670                      evap    = MAX( evap, -qr(k,j,i) / dt_micro )
671                      evap_nr = MAX( c_evap * evap / xr * hyrho(k),            &
672                                     -nr(k,j,i) / dt_micro )
673
674                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro
675                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro
676
677                   ENDIF
678                ENDIF         
679
680             ENDDO
681          ENDDO
682       ENDDO
683
684       CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
685
686    END SUBROUTINE evaporation_rain
687
688
689    SUBROUTINE sedimentation_cloud
690
691       USE arrays_3d,                                                          &
692           ONLY:  ddzu, dzu, pt, q, qc
693
694       USE cloud_parameters,                                                   &
695           ONLY:  eps_sb, hyrho, l_d_cp, nc_const, pt_d_t, sed_qc_const
696
697       USE constants,                                                          &
698           ONLY:  pi
699
700       USE control_parameters,                                                 &
701           ONLY:  dt_do2d_xy, dt_micro, intermediate_timestep_count
702
703       USE cpulog,                                                             &
704           ONLY:  cpu_log, log_point_s
705
706       USE indices,                                                            &
707           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
708
709       USE kinds
710       
711       IMPLICIT NONE
712
713       INTEGER(iwp) ::  i                 !:
714       INTEGER(iwp) ::  j                 !:
715       INTEGER(iwp) ::  k                 !:
716
717       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !:
718
719       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
720
721!
722!--    Sedimentation of cloud droplets (Ackermann et al., 2009, MWR):
723       sed_qc(nzt+1) = 0.0_wp
724
725       DO  i = nxl, nxr
726          DO  j = nys, nyn
727             DO  k = nzt, nzb_s_inner(j,i)+1, -1
728
729                IF ( qc(k,j,i) > eps_sb )  THEN
730                   sed_qc(k) = sed_qc_const * nc_const**( -2.0_wp / 3.0_wp ) * &
731                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )
732                ELSE
733                   sed_qc(k) = 0.0_wp
734                ENDIF
735
736                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
737                                            dt_micro + sed_qc(k+1)             &
738                               )
739
740                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
741                                        ddzu(k+1) / hyrho(k) * dt_micro
742                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
743                                        ddzu(k+1) / hyrho(k) * dt_micro
744                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
745                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
746                                        pt_d_t(k) * dt_micro
747
748             ENDDO
749          ENDDO
750       ENDDO
751
752       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
753
754    END SUBROUTINE sedimentation_cloud
755
756
757    SUBROUTINE sedimentation_rain
758
759       USE arrays_3d,                                                          &
760           ONLY:  ddzu, dzu, nr, pt, q, qr
761
762       USE cloud_parameters,                                                   &
763           ONLY:  a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho,        &
764                  limiter_sedimentation, l_d_cp, precipitation_amount, prr,    &
765                  pt_d_t, stp
766
767       USE control_parameters,                                                 &
768           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro,     &
769                  dt_3d, intermediate_timestep_count,                          &
770                  intermediate_timestep_count_max,                             &
771                  precipitation_amount_interval, time_do2d_xy
772
773       USE cpulog,                                                             &
774           ONLY:  cpu_log, log_point_s
775
776       USE indices,                                                            &
777           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
778
779       USE kinds
780
781       USE statistics,                                                         &
782           ONLY:  weight_substep
783       
784       IMPLICIT NONE
785
786       INTEGER(iwp) ::  i                          !:
787       INTEGER(iwp) ::  j                          !:
788       INTEGER(iwp) ::  k                          !:
789       INTEGER(iwp) ::  k_run                      !:
790
791       REAL(wp)     ::  c_run                      !:
792       REAL(wp)     ::  d_max                      !:
793       REAL(wp)     ::  d_mean                     !:
794       REAL(wp)     ::  d_min                      !:
795       REAL(wp)     ::  dr                         !:
796       REAL(wp)     ::  dt_sedi                    !:
797       REAL(wp)     ::  flux                       !:
798       REAL(wp)     ::  lambda_r                   !:
799       REAL(wp)     ::  mu_r                       !:
800       REAL(wp)     ::  z_run                      !:
801
802       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !:
803       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !:
804       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_nr     !:
805       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_qr     !:
806       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !:
807       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !:
808       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !:
809       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !:
810       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !:
811       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !:
812
813       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
814!
815!--    Computation of sedimentation flux. Implementation according to Stevens
816!--    and Seifert (2008). Code is based on UCLA-LES.
817       IF ( intermediate_timestep_count == 1 )  prr(:,:,:) = 0.0_wp
818!
819!--    Compute velocities
820       DO  i = nxl, nxr
821          DO  j = nys, nyn
822             DO  k = nzb_s_inner(j,i)+1, nzt
823                IF ( qr(k,j,i) > eps_sb )  THEN
824!
825!--                Weight averaged diameter of rain drops:
826                   dr = ( hyrho(k) * qr(k,j,i) /                               &
827                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
828!
829!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
830!--                Stevens and Seifert, 2008):
831                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
832                                                     ( dr - 1.4E-3_wp ) ) )
833!
834!--                Slope parameter of gamma distribution (Seifert, 2008):
835                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
836                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
837
838                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
839                                               a_term - b_term * ( 1.0_wp +    &
840                                                  c_term /                     &
841                                                  lambda_r )**( -1.0_wp *      &
842                                                  ( mu_r + 1.0_wp ) )          &
843                                              )                                &
844                                )
845
846                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
847                                               a_term - b_term * ( 1.0_wp +    &
848                                                  c_term /                     &
849                                                  lambda_r )**( -1.0_wp *      &
850                                                  ( mu_r + 4.0_wp ) )          &
851                                             )                                 &
852                                )
853                ELSE
854                   w_nr(k) = 0.0_wp
855                   w_qr(k) = 0.0_wp
856                ENDIF
857             ENDDO
858!
859!--          Adjust boundary values
860             w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
861             w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
862             w_nr(nzt+1) = 0.0_wp
863             w_qr(nzt+1) = 0.0_wp
864!
865!--          Compute Courant number
866             DO  k = nzb_s_inner(j,i)+1, nzt
867                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
868                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
869                          dt_micro * ddzu(k)
870                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
871                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
872                          dt_micro * ddzu(k)
873             ENDDO     
874!
875!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
876             IF ( limiter_sedimentation )  THEN
877
878                DO k = nzb_s_inner(j,i)+1, nzt
879                   d_mean = 0.5_wp * ( qr(k+1,j,i) + qr(k-1,j,i) )
880                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
881                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
882
883                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
884                                                              2.0_wp * d_max,  &
885                                                              ABS( d_mean ) )
886
887                   d_mean = 0.5_wp * ( nr(k+1,j,i) + nr(k-1,j,i) )
888                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
889                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
890
891                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
892                                                              2.0_wp * d_max,  &
893                                                              ABS( d_mean ) )
894                ENDDO
895
896             ELSE
897
898                nr_slope = 0.0_wp
899                qr_slope = 0.0_wp
900
901             ENDIF
902
903             sed_nr(nzt+1) = 0.0_wp
904             sed_qr(nzt+1) = 0.0_wp
905!
906!--          Compute sedimentation flux
907             DO  k = nzt, nzb_s_inner(j,i)+1, -1
908!
909!--             Sum up all rain drop number densities which contribute to the flux
910!--             through k-1/2
911                flux  = 0.0_wp
912                z_run = 0.0_wp ! height above z(k)
913                k_run = k
914                c_run = MIN( 1.0_wp, c_nr(k) )
915                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
916                   flux  = flux + hyrho(k_run) *                               &
917                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
918                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)
919                   z_run = z_run + dzu(k_run)
920                   k_run = k_run + 1
921                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
922                ENDDO
923!
924!--             It is not allowed to sediment more rain drop number density than
925!--             available
926                flux = MIN( flux,                                              &
927                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
928                            dt_micro                                           &
929                          )
930
931                sed_nr(k) = flux / dt_micro
932                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
933                                        ddzu(k+1) / hyrho(k) * dt_micro
934!
935!--             Sum up all rain water content which contributes to the flux
936!--             through k-1/2
937                flux  = 0.0_wp
938                z_run = 0.0_wp ! height above z(k)
939                k_run = k
940                c_run = MIN( 1.0_wp, c_qr(k) )
941
942                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
943
944                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
945                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
946                                  0.5_wp ) * c_run * dzu(k_run)
947                   z_run = z_run + dzu(k_run)
948                   k_run = k_run + 1
949                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
950
951                ENDDO
952!
953!--             It is not allowed to sediment more rain water content than
954!--             available
955                flux = MIN( flux,                                              &
956                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
957                            dt_micro                                           &
958                          )
959
960                sed_qr(k) = flux / dt_micro
961
962                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
963                                        ddzu(k+1) / hyrho(k) * dt_micro
964                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
965                                        ddzu(k+1) / hyrho(k) * dt_micro 
966                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
967                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
968                                        pt_d_t(k) * dt_micro
969!
970!--             Compute the rain rate
971                IF ( call_microphysics_at_all_substeps )  THEN
972                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *            &
973                                weight_substep(intermediate_timestep_count)
974                ELSE
975                   prr(k,j,i) = sed_qr(k) / hyrho(k)
976                ENDIF
977
978             ENDDO
979          ENDDO
980       ENDDO
981
982!
983!--    Precipitation amount
984       IF ( intermediate_timestep_count == intermediate_timestep_count_max     &
985            .AND.  ( dt_do2d_xy - time_do2d_xy ) <                             &
986                   precipitation_amount_interval )  THEN
987          DO  i = nxl, nxr
988             DO  j = nys, nyn
989                precipitation_amount(j,i) = precipitation_amount(j,i) +        &
990                                            prr(nzb_s_inner(j,i)+1,j,i) *      &
991                                            hyrho(nzb_s_inner(j,i)+1) * dt_3d
992             ENDDO
993          ENDDO
994       ENDIF
995
996       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
997
998    END SUBROUTINE sedimentation_rain
999
1000
1001!------------------------------------------------------------------------------!
1002! Call for grid point i,j
1003!------------------------------------------------------------------------------!
1004
1005    SUBROUTINE microphysics_control_ij( i, j )
1006
1007       USE arrays_3d,                                                          &
1008           ONLY:  hyp, nc_1d, nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc,      &
1009                  qc_1d, qr, qr_1d, zu
1010
1011       USE cloud_parameters,                                                   &
1012           ONLY:  cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
1013
1014       USE control_parameters,                                                 &
1015           ONLY:  call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
1016                  g, intermediate_timestep_count, large_scale_forcing,         &
1017                  lsf_surf, precipitation, pt_surface,                         &
1018                  rho_surface,surface_pressure
1019
1020       USE indices,                                                            &
1021           ONLY:  nzb, nzt
1022
1023       USE kinds
1024
1025       USE statistics,                                                         &
1026           ONLY:  weight_pres
1027
1028       IMPLICIT NONE
1029
1030       INTEGER(iwp) ::  i                 !:
1031       INTEGER(iwp) ::  j                 !:
1032       INTEGER(iwp) ::  k                 !:
1033
1034       REAL(wp)     ::  t_surface         !:
1035
1036       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
1037!
1038!--       Calculate:
1039!--       pt / t : ratio of potential and actual temperature (pt_d_t)
1040!--       t / pt : ratio of actual and potential temperature (t_d_pt)
1041!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
1042          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1043          DO  k = nzb, nzt+1
1044             hyp(k)    = surface_pressure * 100.0_wp * &
1045                         ( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
1046             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
1047             t_d_pt(k) = 1.0_wp / pt_d_t(k)
1048             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
1049          ENDDO
1050!
1051!--       Compute reference density
1052          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
1053       ENDIF
1054
1055!
1056!--    Compute length of time step
1057       IF ( call_microphysics_at_all_substeps )  THEN
1058          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
1059       ELSE
1060          dt_micro = dt_3d
1061       ENDIF
1062
1063!
1064!--    Use 1d arrays
1065       q_1d(:)  = q(:,j,i)
1066       pt_1d(:) = pt(:,j,i)
1067       qc_1d(:) = qc(:,j,i)
1068       nc_1d(:) = nc_const
1069       IF ( precipitation )  THEN
1070          qr_1d(:) = qr(:,j,i)
1071          nr_1d(:) = nr(:,j,i)
1072       ENDIF
1073
1074!
1075!--    Compute cloud physics
1076       IF ( precipitation )  THEN
1077          CALL adjust_cloud( i,j ) 
1078          CALL autoconversion( i,j )
1079          CALL accretion( i,j )
1080          CALL selfcollection_breakup( i,j )
1081          CALL evaporation_rain( i,j )
1082          CALL sedimentation_rain( i,j )
1083       ENDIF
1084
1085       IF ( drizzle )  CALL sedimentation_cloud( i,j )
1086
1087!
1088!--    Store results on the 3d arrays
1089       q(:,j,i)  = q_1d(:)
1090       pt(:,j,i) = pt_1d(:)
1091       IF ( precipitation )  THEN
1092          qr(:,j,i) = qr_1d(:)
1093          nr(:,j,i) = nr_1d(:)
1094       ENDIF
1095
1096    END SUBROUTINE microphysics_control_ij
1097
1098    SUBROUTINE adjust_cloud_ij( i, j )
1099
1100       USE arrays_3d,                                                          &
1101           ONLY:  qr_1d, nr_1d
1102
1103       USE cloud_parameters,                                                   &
1104           ONLY:  eps_sb, xrmin, xrmax, hyrho, k_cc, x0
1105
1106       USE indices,                                                            &
1107           ONLY:  nzb, nzb_s_inner, nzt
1108
1109       USE kinds
1110
1111       IMPLICIT NONE
1112
1113       INTEGER(iwp) ::  i                 !:
1114       INTEGER(iwp) ::  j                 !:
1115       INTEGER(iwp) ::  k                 !:
1116!
1117!--    Adjust number of raindrops to avoid nonlinear effects in
1118!--    sedimentation and evaporation of rain drops due to too small or
1119!--    too big weights of rain drops (Stevens and Seifert, 2008).
1120!--    The same procedure is applied to cloud droplets if they are determined
1121!--    prognostically. 
1122       DO  k = nzb_s_inner(j,i)+1, nzt
1123
1124          IF ( qr_1d(k) <= eps_sb )  THEN
1125             qr_1d(k) = 0.0_wp
1126             nr_1d(k) = 0.0_wp
1127          ELSE
1128!
1129!--          Adjust number of raindrops to avoid nonlinear effects in
1130!--          sedimentation and evaporation of rain drops due to too small or
1131!--          too big weights of rain drops (Stevens and Seifert, 2008).
1132             IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) )  THEN
1133                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin
1134             ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) )  THEN
1135                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax
1136             ENDIF
1137
1138          ENDIF
1139
1140       ENDDO
1141
1142    END SUBROUTINE adjust_cloud_ij
1143
1144
1145    SUBROUTINE autoconversion_ij( i, j )
1146
1147       USE arrays_3d,                                                          &
1148           ONLY:  diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d
1149
1150       USE cloud_parameters,                                                   &
1151           ONLY:  a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3,        &
1152                  c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, x0
1153
1154       USE control_parameters,                                                 &
1155           ONLY:  dt_micro, rho_surface, turbulence
1156
1157       USE grid_variables,                                                     &
1158           ONLY:  dx, dy
1159
1160       USE indices,                                                            &
1161           ONLY:  nzb, nzb_s_inner, nzt
1162
1163       USE kinds
1164
1165       IMPLICIT NONE
1166
1167       INTEGER(iwp) ::  i                 !:
1168       INTEGER(iwp) ::  j                 !:
1169       INTEGER(iwp) ::  k                 !:
1170
1171       REAL(wp)     ::  alpha_cc          !:                   
1172       REAL(wp)     ::  autocon           !:
1173       REAL(wp)     ::  dissipation       !:
1174       REAL(wp)     ::  k_au              !:
1175       REAL(wp)     ::  l_mix             !:
1176       REAL(wp)     ::  nu_c              !:
1177       REAL(wp)     ::  phi_au            !:
1178       REAL(wp)     ::  r_cc              !:
1179       REAL(wp)     ::  rc                !:
1180       REAL(wp)     ::  re_lambda         !:
1181       REAL(wp)     ::  selfcoll          !:
1182       REAL(wp)     ::  sigma_cc          !:
1183       REAL(wp)     ::  tau_cloud         !:
1184       REAL(wp)     ::  xc                !:
1185
1186       DO  k = nzb_s_inner(j,i)+1, nzt
1187
1188          IF ( qc_1d(k) > eps_sb )  THEN
1189
1190             k_au = k_cc / ( 20.0_wp * x0 )
1191!
1192!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
1193!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
1194             tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
1195!
1196!--          Universal function for autoconversion process
1197!--          (Seifert and Beheng, 2006):
1198             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
1199!
1200!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
1201!--          (Use constant nu_c = 1.0_wp instead?)
1202             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
1203!
1204!--          Mean weight of cloud droplets:
1205             xc = hyrho(k) * qc_1d(k) / nc_1d(k)
1206!
1207!--          Parameterized turbulence effects on autoconversion (Seifert,
1208!--          Nuijens and Stevens, 2010)
1209             IF ( turbulence )  THEN
1210!
1211!--             Weight averaged radius of cloud droplets:
1212                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
1213
1214                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
1215                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
1216                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
1217!
1218!--             Mixing length (neglecting distance to ground and stratification)
1219                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
1220!
1221!--             Limit dissipation rate according to Seifert, Nuijens and
1222!--             Stevens (2010)
1223                dissipation = MIN( 0.06_wp, diss(k,j,i) )
1224!
1225!--             Compute Taylor-microscale Reynolds number:
1226                re_lambda = 6.0_wp / 11.0_wp *                                 &
1227                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
1228                            SQRT( 15.0_wp / kin_vis_air ) *                    &
1229                            dissipation**( 1.0_wp / 6.0_wp )
1230!
1231!--             The factor of 1.0E4 is needed to convert the dissipation rate
1232!--             from m2 s-3 to cm2 s-3.
1233                k_au = k_au * ( 1.0_wp +                                       &
1234                                dissipation * 1.0E4_wp *                       &
1235                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
1236                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
1237                                                  sigma_cc )**2                &
1238                                                ) + beta_cc                    &
1239                                )                                              &
1240                              )
1241             ENDIF
1242!
1243!--          Autoconversion rate (Seifert and Beheng, 2006):
1244             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
1245                       ( nu_c + 1.0_wp )**2 * qc_1d(k)**2 * xc**2 *            &
1246                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
1247                       rho_surface
1248             autocon = MIN( autocon, qc_1d(k) / dt_micro )
1249
1250             qr_1d(k) = qr_1d(k) + autocon * dt_micro
1251             qc_1d(k) = qc_1d(k) - autocon * dt_micro 
1252             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
1253
1254          ENDIF
1255
1256       ENDDO
1257
1258    END SUBROUTINE autoconversion_ij
1259
1260
1261    SUBROUTINE accretion_ij( i, j )
1262
1263       USE arrays_3d,                                                          &
1264           ONLY:  diss, qc_1d, qr_1d
1265
1266       USE cloud_parameters,                                                   &
1267           ONLY:  eps_sb, hyrho, k_cr0
1268
1269       USE control_parameters,                                                 &
1270           ONLY:  dt_micro, rho_surface, turbulence
1271
1272       USE indices,                                                            &
1273           ONLY:  nzb, nzb_s_inner, nzt
1274
1275       USE kinds
1276
1277       IMPLICIT NONE
1278
1279       INTEGER(iwp) ::  i                 !:
1280       INTEGER(iwp) ::  j                 !:
1281       INTEGER(iwp) ::  k                 !:
1282
1283       REAL(wp)     ::  accr              !:
1284       REAL(wp)     ::  k_cr              !:
1285       REAL(wp)     ::  phi_ac            !:
1286       REAL(wp)     ::  tau_cloud         !:
1287       REAL(wp)     ::  xc                !:
1288
1289       DO  k = nzb_s_inner(j,i)+1, nzt
1290          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
1291!
1292!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
1293             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 
1294!
1295!--          Universal function for accretion process
1296!--          (Seifert and Beheng, 2001):
1297             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
1298!
1299!--          Parameterized turbulence effects on autoconversion (Seifert,
1300!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
1301!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
1302             IF ( turbulence )  THEN
1303                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
1304                                 MIN( 600.0_wp,                                &
1305                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
1306                               )
1307             ELSE
1308                k_cr = k_cr0                       
1309             ENDIF
1310!
1311!--          Accretion rate (Seifert and Beheng, 2006):
1312             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
1313             accr = MIN( accr, qc_1d(k) / dt_micro )
1314
1315             qr_1d(k) = qr_1d(k) + accr * dt_micro 
1316             qc_1d(k) = qc_1d(k) - accr * dt_micro
1317
1318          ENDIF
1319
1320       ENDDO
1321
1322    END SUBROUTINE accretion_ij
1323
1324
1325    SUBROUTINE selfcollection_breakup_ij( i, j )
1326
1327       USE arrays_3d,                                                          &
1328           ONLY:  nr_1d, qr_1d
1329
1330       USE cloud_parameters,                                                   &
1331           ONLY:  dpirho_l, eps_sb, hyrho, k_br, k_rr
1332
1333       USE control_parameters,                                                 &
1334           ONLY:  dt_micro, rho_surface
1335
1336       USE indices,                                                            &
1337           ONLY:  nzb, nzb_s_inner, nzt
1338
1339       USE kinds
1340   
1341       IMPLICIT NONE
1342
1343       INTEGER(iwp) ::  i                 !:
1344       INTEGER(iwp) ::  j                 !:
1345       INTEGER(iwp) ::  k                 !:
1346
1347       REAL(wp)     ::  breakup           !:
1348       REAL(wp)     ::  dr                !:
1349       REAL(wp)     ::  phi_br            !:
1350       REAL(wp)     ::  selfcoll          !:
1351
1352       DO  k = nzb_s_inner(j,i)+1, nzt
1353          IF ( qr_1d(k) > eps_sb )  THEN
1354!
1355!--          Selfcollection rate (Seifert and Beheng, 2001):
1356             selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
1357!
1358!--          Weight averaged diameter of rain drops:
1359             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
1360!
1361!--          Collisional breakup rate (Seifert, 2008):
1362             IF ( dr >= 0.3E-3_wp )  THEN
1363                phi_br  = k_br * ( dr - 1.1E-3_wp )
1364                breakup = selfcoll * ( phi_br + 1.0_wp )
1365             ELSE
1366                breakup = 0.0_wp
1367             ENDIF
1368
1369             selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
1370             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
1371
1372          ENDIF         
1373       ENDDO
1374
1375    END SUBROUTINE selfcollection_breakup_ij
1376
1377
1378    SUBROUTINE evaporation_rain_ij( i, j )
1379!
1380!--    Evaporation of precipitable water. Condensation is neglected for
1381!--    precipitable water.
1382
1383       USE arrays_3d,                                                          &
1384           ONLY:  hyp, nr_1d, pt_1d, q_1d,  qc_1d, qr_1d
1385
1386       USE cloud_parameters,                                                   &
1387           ONLY:  a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
1388                  dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r,   &
1389                  l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l,      &
1390                  t_d_pt, ventilation_effect
1391
1392       USE constants,                                                          &
1393           ONLY:  pi
1394
1395       USE control_parameters,                                                 &
1396           ONLY:  dt_micro
1397
1398       USE indices,                                                            &
1399           ONLY:  nzb, nzb_s_inner, nzt
1400
1401       USE kinds
1402
1403       IMPLICIT NONE
1404
1405       INTEGER(iwp) ::  i                 !:
1406       INTEGER(iwp) ::  j                 !:
1407       INTEGER(iwp) ::  k                 !:
1408
1409       REAL(wp)     ::  alpha             !:
1410       REAL(wp)     ::  dr                !:
1411       REAL(wp)     ::  e_s               !:
1412       REAL(wp)     ::  evap              !:
1413       REAL(wp)     ::  evap_nr           !:
1414       REAL(wp)     ::  f_vent            !:
1415       REAL(wp)     ::  g_evap            !:
1416       REAL(wp)     ::  lambda_r          !:
1417       REAL(wp)     ::  mu_r              !:
1418       REAL(wp)     ::  mu_r_2            !:
1419       REAL(wp)     ::  mu_r_5d2          !:
1420       REAL(wp)     ::  nr_0              !:
1421       REAL(wp)     ::  q_s               !:
1422       REAL(wp)     ::  sat               !:
1423       REAL(wp)     ::  t_l               !:
1424       REAL(wp)     ::  temp              !:
1425       REAL(wp)     ::  xr                !:
1426
1427       DO  k = nzb_s_inner(j,i)+1, nzt
1428          IF ( qr_1d(k) > eps_sb )  THEN
1429!
1430!--          Actual liquid water temperature:
1431             t_l = t_d_pt(k) * pt_1d(k)
1432!
1433!--          Saturation vapor pressure at t_l:
1434             e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /          &
1435                                    ( t_l - 35.86_wp )                         &
1436                                  )
1437!
1438!--          Computation of saturation humidity:
1439             q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
1440             alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
1441             q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
1442!
1443!--          Supersaturation:
1444             sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
1445!
1446!--          Evaporation needs only to be calculated in subsaturated regions
1447             IF ( sat < 0.0_wp )  THEN
1448!
1449!--             Actual temperature:
1450                temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
1451       
1452                g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /  &
1453                                    ( thermal_conductivity_l * temp ) +        &
1454                                    r_v * temp / ( diff_coeff_l * e_s )        &
1455                                  )
1456!
1457!--             Mean weight of rain drops
1458                xr = hyrho(k) * qr_1d(k) / nr_1d(k)
1459!
1460!--             Weight averaged diameter of rain drops:
1461                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
1462!
1463!--             Compute ventilation factor and intercept parameter
1464!--             (Seifert and Beheng, 2006; Seifert, 2008):
1465                IF ( ventilation_effect )  THEN
1466!
1467!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1468!--                Stevens and Seifert, 2008):
1469                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
1470!
1471!--                Slope parameter of gamma distribution (Seifert, 2008):
1472                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
1473                                ( mu_r + 1.0_wp )                              &
1474                              )**( 1.0_wp / 3.0_wp ) / dr
1475
1476                   mu_r_2   = mu_r + 2.0_wp
1477                   mu_r_5d2 = mu_r + 2.5_wp 
1478
1479                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  & 
1480                            b_vent * schmidt_p_1d3 *                           &
1481                            SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
1482                            lambda_r**( -mu_r_5d2 ) *                          &
1483                            ( 1.0_wp -                                         &
1484                              0.5_wp * ( b_term / a_term ) *                   &
1485                              ( lambda_r / ( c_term + lambda_r )               &
1486                              )**mu_r_5d2 -                                    &
1487                              0.125_wp * ( b_term / a_term )**2 *              &
1488                              ( lambda_r / ( 2.0_wp * c_term + lambda_r )      &
1489                              )**mu_r_5d2 -                                    &
1490                              0.0625_wp * ( b_term / a_term )**3 *             &
1491                              ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
1492                              )**mu_r_5d2 -                                    &
1493                              0.0390625_wp * ( b_term / a_term )**4 *          & 
1494                              ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
1495                              )**mu_r_5d2                                      &
1496                            )
1497
1498                   nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) /           &
1499                            gamm( mu_r + 1.0_wp ) 
1500                ELSE
1501                   f_vent = 1.0_wp
1502                   nr_0   = nr_1d(k) * dr
1503                ENDIF
1504!
1505!--             Evaporation rate of rain water content (Seifert and Beheng, 2006):
1506                evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
1507                evap    = MAX( evap, -qr_1d(k) / dt_micro )
1508                evap_nr = MAX( c_evap * evap / xr * hyrho(k),                  &
1509                               -nr_1d(k) / dt_micro )
1510
1511                qr_1d(k) = qr_1d(k) + evap    * dt_micro
1512                nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
1513
1514             ENDIF
1515          ENDIF         
1516
1517       ENDDO
1518
1519    END SUBROUTINE evaporation_rain_ij
1520
1521
1522    SUBROUTINE sedimentation_cloud_ij( i, j )
1523
1524       USE arrays_3d,                                                          &
1525           ONLY:  ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d
1526
1527       USE cloud_parameters,                                                   &
1528           ONLY:  eps_sb, hyrho, l_d_cp, pt_d_t, sed_qc_const
1529
1530       USE constants,                                                          &
1531           ONLY:  pi
1532
1533       USE control_parameters,                                                 &
1534           ONLY:  dt_do2d_xy, dt_micro, intermediate_timestep_count
1535
1536       USE indices,                                                            &
1537           ONLY:  nzb, nzb_s_inner, nzt
1538
1539       USE kinds
1540       
1541       IMPLICIT NONE
1542
1543       INTEGER(iwp) ::  i                 !:
1544       INTEGER(iwp) ::  j                 !:
1545       INTEGER(iwp) ::  k                 !:
1546
1547       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc  !:
1548
1549!
1550!--    Sedimentation of cloud droplets (Ackermann et al., 2009, MWR):
1551       sed_qc(nzt+1) = 0.0_wp
1552
1553       DO  k = nzt, nzb_s_inner(j,i)+1, -1
1554          IF ( qc_1d(k) > eps_sb )  THEN
1555             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) *       &
1556                         ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
1557          ELSE
1558             sed_qc(k) = 0.0_wp
1559          ENDIF
1560
1561          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /          &
1562                                      dt_micro + sed_qc(k+1)                   &
1563                         )
1564
1565          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
1566                                hyrho(k) * dt_micro
1567          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      & 
1568                                hyrho(k) * dt_micro
1569          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
1570                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
1571
1572       ENDDO
1573
1574    END SUBROUTINE sedimentation_cloud_ij
1575
1576
1577    SUBROUTINE sedimentation_rain_ij( i, j )
1578
1579       USE arrays_3d,                                                          &
1580           ONLY:  ddzu, dzu, nr_1d, pt_1d, q_1d, qr_1d
1581
1582       USE cloud_parameters,                                                   &
1583           ONLY:  a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho,        &
1584                  limiter_sedimentation, l_d_cp, precipitation_amount, prr,    &
1585                  pt_d_t, stp
1586
1587       USE control_parameters,                                                 &
1588           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro,     &
1589                  dt_3d, intermediate_timestep_count,                          &
1590                  intermediate_timestep_count_max,                             &
1591                  precipitation_amount_interval, time_do2d_xy
1592
1593       USE indices,                                                            &
1594           ONLY:  nzb, nzb_s_inner, nzt
1595
1596       USE kinds
1597
1598       USE statistics,                                                         &
1599           ONLY:  weight_substep
1600       
1601       IMPLICIT NONE
1602
1603       INTEGER(iwp) ::  i                          !:
1604       INTEGER(iwp) ::  j                          !:
1605       INTEGER(iwp) ::  k                          !:
1606       INTEGER(iwp) ::  k_run                      !:
1607
1608       REAL(wp)     ::  c_run                      !:
1609       REAL(wp)     ::  d_max                      !:
1610       REAL(wp)     ::  d_mean                     !:
1611       REAL(wp)     ::  d_min                      !:
1612       REAL(wp)     ::  dr                         !:
1613       REAL(wp)     ::  dt_sedi                    !:
1614       REAL(wp)     ::  flux                       !:
1615       REAL(wp)     ::  lambda_r                   !:
1616       REAL(wp)     ::  mu_r                       !:
1617       REAL(wp)     ::  z_run                      !:
1618
1619       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !:
1620       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !:
1621       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_nr     !:
1622       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_qr     !:
1623       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !:
1624       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !:
1625       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !:
1626       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !:
1627       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !:
1628       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !:
1629
1630
1631!
1632!--    Computation of sedimentation flux. Implementation according to Stevens
1633!--    and Seifert (2008). Code is based on UCLA-LES.
1634       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
1635!
1636!--    Compute velocities
1637       DO  k = nzb_s_inner(j,i)+1, nzt
1638          IF ( qr_1d(k) > eps_sb )  THEN
1639!
1640!--          Weight averaged diameter of rain drops:
1641             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
1642!
1643!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1644!--          Stevens and Seifert, 2008):
1645             mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
1646!
1647!--          Slope parameter of gamma distribution (Seifert, 2008):
1648             lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *              &
1649                          ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
1650
1651             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
1652                                         a_term - b_term * ( 1.0_wp +          &
1653                                            c_term / lambda_r )**( -1.0_wp *   &
1654                                            ( mu_r + 1.0_wp ) )                &
1655                                        )                                      &
1656                          )
1657             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
1658                                         a_term - b_term * ( 1.0_wp +          &
1659                                            c_term / lambda_r )**( -1.0_wp *   &
1660                                            ( mu_r + 4.0_wp ) )                &
1661                                       )                                       &
1662                          )
1663          ELSE
1664             w_nr(k) = 0.0_wp
1665             w_qr(k) = 0.0_wp
1666          ENDIF
1667       ENDDO
1668!
1669!--    Adjust boundary values
1670       w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
1671       w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
1672       w_nr(nzt+1) = 0.0_wp
1673       w_qr(nzt+1) = 0.0_wp
1674!
1675!--    Compute Courant number
1676       DO  k = nzb_s_inner(j,i)+1, nzt
1677          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
1678                    dt_micro * ddzu(k)
1679          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
1680                    dt_micro * ddzu(k)
1681       ENDDO     
1682!
1683!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
1684       IF ( limiter_sedimentation )  THEN
1685
1686          DO k = nzb_s_inner(j,i)+1, nzt
1687             d_mean = 0.5_wp * ( qr_1d(k+1) + qr_1d(k-1) )
1688             d_min  = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
1689             d_max  = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
1690
1691             qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
1692                                                        2.0_wp * d_max,        &
1693                                                        ABS( d_mean ) )
1694
1695             d_mean = 0.5_wp * ( nr_1d(k+1) + nr_1d(k-1) )
1696             d_min  = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
1697             d_max  = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
1698
1699             nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
1700                                                        2.0_wp * d_max,        &
1701                                                        ABS( d_mean ) )
1702          ENDDO
1703
1704       ELSE
1705
1706          nr_slope = 0.0_wp
1707          qr_slope = 0.0_wp
1708
1709       ENDIF
1710
1711       sed_nr(nzt+1) = 0.0_wp
1712       sed_qr(nzt+1) = 0.0_wp
1713!
1714!--    Compute sedimentation flux
1715       DO  k = nzt, nzb_s_inner(j,i)+1, -1
1716!
1717!--       Sum up all rain drop number densities which contribute to the flux
1718!--       through k-1/2
1719          flux  = 0.0_wp
1720          z_run = 0.0_wp ! height above z(k)
1721          k_run = k
1722          c_run = MIN( 1.0_wp, c_nr(k) )
1723          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
1724             flux  = flux + hyrho(k_run) *                                     &
1725                     ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) *   &
1726                     0.5_wp ) * c_run * dzu(k_run)
1727             z_run = z_run + dzu(k_run)
1728             k_run = k_run + 1
1729             c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
1730          ENDDO
1731!
1732!--       It is not allowed to sediment more rain drop number density than
1733!--       available
1734          flux = MIN( flux,                                                    &
1735                      hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
1736
1737          sed_nr(k) = flux / dt_micro
1738          nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
1739                                    hyrho(k) * dt_micro
1740!
1741!--       Sum up all rain water content which contributes to the flux
1742!--       through k-1/2
1743          flux  = 0.0_wp
1744          z_run = 0.0_wp ! height above z(k)
1745          k_run = k
1746          c_run = MIN( 1.0_wp, c_qr(k) )
1747
1748          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
1749
1750             flux  = flux + hyrho(k_run) *                                     &
1751                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) *   &
1752                     0.5_wp ) * c_run * dzu(k_run)
1753             z_run = z_run + dzu(k_run)
1754             k_run = k_run + 1
1755             c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
1756
1757          ENDDO
1758!
1759!--       It is not allowed to sediment more rain water content than available
1760          flux = MIN( flux,                                                    &
1761                      hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
1762
1763          sed_qr(k) = flux / dt_micro
1764
1765          qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
1766                                hyrho(k) * dt_micro
1767          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
1768                                hyrho(k) * dt_micro 
1769          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
1770                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
1771!
1772!--       Compute the rain rate
1773          IF ( call_microphysics_at_all_substeps )  THEN
1774             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *                  &
1775                          weight_substep(intermediate_timestep_count)
1776          ELSE
1777             prr(k,j,i) = sed_qr(k) / hyrho(k)
1778          ENDIF
1779
1780       ENDDO
1781
1782!
1783!--    Precipitation amount
1784       IF ( intermediate_timestep_count == intermediate_timestep_count_max     &
1785            .AND.  ( dt_do2d_xy - time_do2d_xy ) <                             &
1786                   precipitation_amount_interval )  THEN
1787
1788          precipitation_amount(j,i) = precipitation_amount(j,i) +              &
1789                                      prr(nzb_s_inner(j,i)+1,j,i) *            &
1790                                      hyrho(nzb_s_inner(j,i)+1) * dt_3d
1791       ENDIF
1792
1793    END SUBROUTINE sedimentation_rain_ij
1794
1795!------------------------------------------------------------------------------!
1796! Call for all optimizations
1797!------------------------------------------------------------------------------!
1798!
1799!-- This function computes the gamma function (Press et al., 1992).
1800!-- The gamma function is needed for the calculation of the evaporation
1801!-- of rain drops.
1802    FUNCTION gamm( xx ) 
1803       
1804       USE cloud_parameters,                                                   &
1805           ONLY:  cof, stp
1806
1807       USE kinds
1808
1809       IMPLICIT NONE
1810
1811       INTEGER(iwp) ::  j            !:
1812
1813       REAL(wp)     ::  gamm         !:
1814       REAL(wp)     ::  ser          !:
1815       REAL(wp)     ::  tmp          !:
1816       REAL(wp)     ::  x_gamm       !:
1817       REAL(wp)     ::  xx           !:
1818       REAL(wp)     ::  y_gamm       !:
1819
1820       x_gamm = xx 
1821       y_gamm = x_gamm 
1822       tmp = x_gamm + 5.5_wp
1823       tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp 
1824       ser = 1.000000000190015_wp
1825
1826       DO  j = 1, 6 
1827          y_gamm = y_gamm + 1.0_wp 
1828          ser    = ser + cof( j ) / y_gamm 
1829       ENDDO
1830
1831!
1832!--    Until this point the algorithm computes the logarithm of the gamma
1833!--    function. Hence, the exponential function is used. 
1834!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
1835       gamm = EXP( tmp ) * stp * ser / x_gamm 
1836
1837       RETURN
1838 
1839    END FUNCTION gamm 
1840
1841 END MODULE microphysics_mod
Note: See TracBrowser for help on using the repository browser.