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

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

improved version of two-moment cloud physics

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