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

Last change on this file since 1846 was 1846, checked in by raasch, 7 years ago

last commit documented

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