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

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

last commit documented

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