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

Last change on this file since 1834 was 1832, checked in by hoffmann, 5 years ago

last commit documented

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