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

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

last commit documented

  • Property svn:keywords set to Id
File size: 81.6 KB
RevLine 
[1682]1!> @file microphysics.f90
[1093]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!
[1818]16! Copyright 1997-2016 Leibniz Universitaet Hannover
[1093]17!--------------------------------------------------------------------------------!
18!
[1000]19! Current revisions:
[1092]20! ------------------
[1822]21!
[1823]22!
[1321]23! Former revisions:
24! -----------------
25! $Id: microphysics.f90 1823 2016-04-07 08:57:52Z hoffmann $
26!
[1823]27! 1822 2016-04-07 07:49:42Z hoffmann
28! Unused variables removed.
29! Kessler scheme integrated.
30!
[1692]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!
[1683]35! 1682 2015-10-07 23:56:08Z knoop
36! Code annotations made doxygen readable
37!
[1647]38! 1646 2015-09-02 16:00:10Z hoffmann
39! Bugfix: Wrong computation of d_mean.
40!
[1362]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!
[1354]46! 1353 2014-04-08 15:21:23Z heinze
47! REAL constants provided with KIND-attribute
48!
[1347]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!
[1335]53! 1334 2014-03-25 12:21:40Z heinze
54! Bugfix: REAL constants provided with KIND-attribute
55!
[1323]56! 1322 2014-03-20 16:38:49Z raasch
57! REAL constants defined as wp-kind
58!
[1321]59! 1320 2014-03-20 08:40:49Z raasch
[1320]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
[1000]65!
[1242]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!
[1116]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!
[1107]74! 1106 2013-03-04 05:31:38Z raasch
75! small changes in code formatting
76!
[1093]77! 1092 2013-02-02 11:24:22Z raasch
78! unused variables removed
79! file put under GPL
80!
[1066]81! 1065 2012-11-22 17:42:36Z hoffmann
82! Sedimentation process implemented according to Stevens and Seifert (2008).
[1115]83! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
[1066]84! and Stevens, 2010).
85!
[1054]86! 1053 2012-11-13 17:11:03Z hoffmann
87! initial revision
[1000]88!
89! Description:
90! ------------
[1682]91!> Calculate cloud microphysics according to the two moment bulk
92!> scheme by Seifert and Beheng (2006).
[1000]93!------------------------------------------------------------------------------!
[1682]94 MODULE microphysics_mod
95 
[1000]96
97    PRIVATE
[1115]98    PUBLIC microphysics_control
[1000]99
[1115]100    INTERFACE microphysics_control
101       MODULE PROCEDURE microphysics_control
102       MODULE PROCEDURE microphysics_control_ij
103    END INTERFACE microphysics_control
[1022]104
[1115]105    INTERFACE adjust_cloud
106       MODULE PROCEDURE adjust_cloud
107       MODULE PROCEDURE adjust_cloud_ij
108    END INTERFACE adjust_cloud
109
[1000]110    INTERFACE autoconversion
111       MODULE PROCEDURE autoconversion
112       MODULE PROCEDURE autoconversion_ij
113    END INTERFACE autoconversion
114
[1822]115    INTERFACE autoconversion_kessler
116       MODULE PROCEDURE autoconversion_kessler
117       MODULE PROCEDURE autoconversion_kessler_ij
118    END INTERFACE autoconversion_kessler
119
[1000]120    INTERFACE accretion
121       MODULE PROCEDURE accretion
122       MODULE PROCEDURE accretion_ij
123    END INTERFACE accretion
[1005]124
125    INTERFACE selfcollection_breakup
126       MODULE PROCEDURE selfcollection_breakup
127       MODULE PROCEDURE selfcollection_breakup_ij
128    END INTERFACE selfcollection_breakup
[1012]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
[1000]139 
[1012]140    INTERFACE sedimentation_rain
141       MODULE PROCEDURE sedimentation_rain
142       MODULE PROCEDURE sedimentation_rain_ij
143    END INTERFACE sedimentation_rain
144
[1691]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
[1000]151 CONTAINS
152
153
154!------------------------------------------------------------------------------!
[1682]155! Description:
156! ------------
157!> Call for all grid points
[1000]158!------------------------------------------------------------------------------!
[1115]159    SUBROUTINE microphysics_control
[1022]160
[1361]161       USE arrays_3d,                                                          &
[1822]162           ONLY:  hyp, pt_init, zu
[1361]163
164       USE cloud_parameters,                                                   &
[1822]165           ONLY:  cp, hyrho, prr, pt_d_t, r_d, t_d_pt
[1361]166
167       USE control_parameters,                                                 &
168           ONLY:  call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
[1822]169                  g, intermediate_timestep_count, large_scale_forcing,         &
170                  lsf_surf, microphysics_kessler, microphysics_seifert,        &
171                  pt_surface, rho_surface,surface_pressure
[1361]172
173       USE indices,                                                            &
174           ONLY:  nzb, nzt
175
[1320]176       USE kinds
[1115]177
[1361]178       USE statistics,                                                         &
179           ONLY:  weight_pres
180
[1115]181       IMPLICIT NONE
182
[1682]183       INTEGER(iwp) ::  k                 !<
[1115]184
[1682]185       REAL(wp)     ::  t_surface         !<
[1361]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) )       
[1115]201          ENDDO
[1822]202
[1361]203!
204!--       Compute reference density
205          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
206       ENDIF
[1115]207
[1361]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!
[1822]217!--    Reset precipitation rate
218       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
219
220!
[1361]221!--    Compute cloud physics
[1822]222       IF ( microphysics_kessler )  THEN
223
224          CALL autoconversion_kessler
225
226       ELSEIF ( microphysics_seifert )  THEN
227
[1361]228          CALL adjust_cloud
229          CALL autoconversion
230          CALL accretion
231          CALL selfcollection_breakup
232          CALL evaporation_rain
233          CALL sedimentation_rain
234
[1822]235          IF ( drizzle )  CALL sedimentation_cloud
[1361]236
[1691]237       ENDIF
238
[1822]239       CALL calc_precipitation_amount
240
[1115]241    END SUBROUTINE microphysics_control
242
[1682]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!------------------------------------------------------------------------------!
[1115]250    SUBROUTINE adjust_cloud
251
[1361]252       USE arrays_3d,                                                          &
253           ONLY:  qr, nr
254
255       USE cloud_parameters,                                                   &
[1822]256           ONLY:  eps_sb, xrmin, xrmax, hyrho
[1361]257
258       USE cpulog,                                                             &
259           ONLY:  cpu_log, log_point_s
260
261       USE indices,                                                            &
[1822]262           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]263
[1320]264       USE kinds
[1022]265
266       IMPLICIT NONE
267
[1682]268       INTEGER(iwp) ::  i                 !<
269       INTEGER(iwp) ::  j                 !<
270       INTEGER(iwp) ::  k                 !<
[1022]271
[1361]272       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
273
[1022]274       DO  i = nxl, nxr
275          DO  j = nys, nyn
[1115]276             DO  k = nzb_s_inner(j,i)+1, nzt
[1361]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
[1022]287             ENDDO
288          ENDDO
289       ENDDO
290
[1361]291       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
292
[1115]293    END SUBROUTINE adjust_cloud
[1022]294
[1106]295
[1682]296!------------------------------------------------------------------------------!
297! Description:
298! ------------
299!> Autoconversion rate (Seifert and Beheng, 2006).
300!------------------------------------------------------------------------------!
[1000]301    SUBROUTINE autoconversion
302
[1361]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,                                                            &
[1822]321           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]322
[1320]323       USE kinds
[1000]324
325       IMPLICIT NONE
326
[1682]327       INTEGER(iwp) ::  i                 !<
328       INTEGER(iwp) ::  j                 !<
329       INTEGER(iwp) ::  k                 !<
[1000]330
[1682]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                !<
[1361]344
345       CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
346
[1000]347       DO  i = nxl, nxr
348          DO  j = nys, nyn
[1115]349             DO  k = nzb_s_inner(j,i)+1, nzt
[1000]350
[1361]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
[1000]422             ENDDO
423          ENDDO
424       ENDDO
425
[1361]426       CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
427
[1000]428    END SUBROUTINE autoconversion
429
[1106]430
[1682]431!------------------------------------------------------------------------------!
432! Description:
433! ------------
[1822]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! ------------
[1682]489!> Accretion rate (Seifert and Beheng, 2006).
490!------------------------------------------------------------------------------!
[1005]491    SUBROUTINE accretion
[1000]492
[1361]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,                                                            &
[1822]506           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]507
[1320]508       USE kinds
[1005]509
[1000]510       IMPLICIT NONE
511
[1682]512       INTEGER(iwp) ::  i                 !<
513       INTEGER(iwp) ::  j                 !<
514       INTEGER(iwp) ::  k                 !<
[1000]515
[1682]516       REAL(wp)     ::  accr              !<
517       REAL(wp)     ::  k_cr              !<
518       REAL(wp)     ::  phi_ac            !<
519       REAL(wp)     ::  tau_cloud         !<
[1361]520
521       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
522
[1005]523       DO  i = nxl, nxr
524          DO  j = nys, nyn
[1115]525             DO  k = nzb_s_inner(j,i)+1, nzt
[1000]526
[1361]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
[1005]558             ENDDO
559          ENDDO
[1000]560       ENDDO
561
[1361]562       CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
563
[1005]564    END SUBROUTINE accretion
[1000]565
[1106]566
[1682]567!------------------------------------------------------------------------------!
568! Description:
569! ------------
570!> Collisional breakup rate (Seifert, 2008).
571!------------------------------------------------------------------------------!
[1005]572    SUBROUTINE selfcollection_breakup
[1000]573
[1361]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,                                                            &
[1822]587           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]588
[1320]589       USE kinds
[1361]590   
[1000]591       IMPLICIT NONE
592
[1682]593       INTEGER(iwp) ::  i                 !<
594       INTEGER(iwp) ::  j                 !<
595       INTEGER(iwp) ::  k                 !<
[1000]596
[1682]597       REAL(wp)     ::  breakup           !<
598       REAL(wp)     ::  dr                !<
599       REAL(wp)     ::  phi_br            !<
600       REAL(wp)     ::  selfcoll          !<
[1361]601
602       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
603
[1000]604       DO  i = nxl, nxr
605          DO  j = nys, nyn
[1115]606             DO  k = nzb_s_inner(j,i)+1, nzt
[1361]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
[1000]624
[1361]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         
[1000]629             ENDDO
630          ENDDO
631       ENDDO
632
[1361]633       CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
634
[1005]635    END SUBROUTINE selfcollection_breakup
[1000]636
[1106]637
[1682]638!------------------------------------------------------------------------------!
639! Description:
640! ------------
641!> Evaporation of precipitable water. Condensation is neglected for
642!> precipitable water.
643!------------------------------------------------------------------------------!
[1012]644    SUBROUTINE evaporation_rain
[1000]645
[1361]646       USE arrays_3d,                                                          &
647           ONLY:  hyp, nr, pt, q,  qc, qr
648
649       USE cloud_parameters,                                                   &
[1822]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
[1361]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,                                                            &
[1822]665           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]666
[1320]667       USE kinds
[1012]668
669       IMPLICIT NONE
670
[1682]671       INTEGER(iwp) ::  i                 !<
672       INTEGER(iwp) ::  j                 !<
673       INTEGER(iwp) ::  k                 !<
[1361]674
[1682]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                !<
[1361]692
693       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
694
[1012]695       DO  i = nxl, nxr
696          DO  j = nys, nyn
[1115]697             DO  k = nzb_s_inner(j,i)+1, nzt
[1361]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
[1012]747
[1361]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
[1012]791             ENDDO
792          ENDDO
793       ENDDO
794
[1361]795       CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
796
[1012]797    END SUBROUTINE evaporation_rain
798
[1106]799
[1682]800!------------------------------------------------------------------------------!
801! Description:
802! ------------
803!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
804!------------------------------------------------------------------------------!
[1012]805    SUBROUTINE sedimentation_cloud
806
[1361]807       USE arrays_3d,                                                          &
808           ONLY:  ddzu, dzu, pt, q, qc
809
810       USE cloud_parameters,                                                   &
[1691]811           ONLY:  eps_sb, hyrho, l_d_cp, nc_const, prr, pt_d_t, sed_qc_const
[1361]812
813       USE control_parameters,                                                 &
[1691]814           ONLY:  call_microphysics_at_all_substeps, dt_micro,                 &
[1822]815                  intermediate_timestep_count
[1361]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
[1320]823       USE kinds
[1691]824
825       USE statistics,                                                         &
826           ONLY:  weight_substep
827   
828
[1012]829       IMPLICIT NONE
830
[1682]831       INTEGER(iwp) ::  i                 !<
832       INTEGER(iwp) ::  j                 !<
833       INTEGER(iwp) ::  k                 !<
[1361]834
[1682]835       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
[1361]836
837       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
838
839       sed_qc(nzt+1) = 0.0_wp
840
[1012]841       DO  i = nxl, nxr
842          DO  j = nys, nyn
[1361]843             DO  k = nzt, nzb_s_inner(j,i)+1, -1
[1012]844
[1361]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
[1691]864!
865!--             Compute the precipitation rate due to cloud (fog) droplets
[1822]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)
[1691]871                ENDIF
872
[1012]873             ENDDO
874          ENDDO
875       ENDDO
876
[1361]877       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
878
[1012]879    END SUBROUTINE sedimentation_cloud
880
[1106]881
[1682]882!------------------------------------------------------------------------------!
883! Description:
884! ------------
885!> Computation of sedimentation flux. Implementation according to Stevens
886!> and Seifert (2008). Code is based on UCLA-LES.
887!------------------------------------------------------------------------------!
[1012]888    SUBROUTINE sedimentation_rain
889
[1361]890       USE arrays_3d,                                                          &
891           ONLY:  ddzu, dzu, nr, pt, q, qr
892
893       USE cloud_parameters,                                                   &
[1822]894           ONLY:  a_term, b_term, c_term, dpirho_l, eps_sb, hyrho,             &
895                  limiter_sedimentation, l_d_cp, prr, pt_d_t
[1361]896
897       USE control_parameters,                                                 &
[1691]898           ONLY:  call_microphysics_at_all_substeps, dt_micro,                 &
899                  intermediate_timestep_count
[1361]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
[1320]906       USE kinds
[1012]907
[1361]908       USE statistics,                                                         &
909           ONLY:  weight_substep
910       
[1012]911       IMPLICIT NONE
912
[1682]913       INTEGER(iwp) ::  i                          !<
914       INTEGER(iwp) ::  j                          !<
915       INTEGER(iwp) ::  k                          !<
916       INTEGER(iwp) ::  k_run                      !<
[1361]917
[1682]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                      !<
[1361]927
[1682]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     !<
[1361]936
937       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
[1682]938
[1361]939!
940!--    Compute velocities
[1012]941       DO  i = nxl, nxr
942          DO  j = nys, nyn
[1115]943             DO  k = nzb_s_inner(j,i)+1, nzt
[1361]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
[1012]958
[1361]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
[1012]978             ENDDO
[1361]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
[1646]1000                   d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
[1361]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
[1646]1008                   d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
[1361]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
[1691]1093                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
1094                                * weight_substep(intermediate_timestep_count)
[1361]1095                ELSE
[1691]1096                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
[1361]1097                ENDIF
1098
1099             ENDDO
[1012]1100          ENDDO
1101       ENDDO
1102
[1691]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
[1361]1140          DO  i = nxl, nxr
1141             DO  j = nys, nyn
[1691]1142
[1361]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
[1691]1146
[1361]1147             ENDDO
1148          ENDDO
1149       ENDIF
1150
[1691]1151    END SUBROUTINE calc_precipitation_amount
[1361]1152
[1012]1153
[1000]1154!------------------------------------------------------------------------------!
[1682]1155! Description:
1156! ------------
1157!> Call for grid point i,j
[1000]1158!------------------------------------------------------------------------------!
[1022]1159
[1115]1160    SUBROUTINE microphysics_control_ij( i, j )
1161
[1320]1162       USE arrays_3d,                                                          &
[1361]1163           ONLY:  hyp, nc_1d, nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc,      &
1164                  qc_1d, qr, qr_1d, zu
[1115]1165
[1320]1166       USE cloud_parameters,                                                   &
[1822]1167           ONLY:  cp, hyrho, nc_const, prr, pt_d_t, r_d, t_d_pt
[1320]1168
1169       USE control_parameters,                                                 &
[1361]1170           ONLY:  call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
1171                  g, intermediate_timestep_count, large_scale_forcing,         &
[1822]1172                  lsf_surf, microphysics_seifert, microphysics_kessler,        &
1173                  pt_surface, rho_surface, surface_pressure
[1320]1174
1175       USE indices,                                                            &
1176           ONLY:  nzb, nzt
1177
1178       USE kinds
1179
1180       USE statistics,                                                         &
1181           ONLY:  weight_pres
1182
[1022]1183       IMPLICIT NONE
1184
[1682]1185       INTEGER(iwp) ::  i                 !<
1186       INTEGER(iwp) ::  j                 !<
1187       INTEGER(iwp) ::  k                 !<
[1115]1188
[1682]1189       REAL(wp)     ::  t_surface         !<
[1320]1190
[1361]1191       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
[1241]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)
[1353]1197          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
[1241]1198          DO  k = nzb, nzt+1
[1353]1199             hyp(k)    = surface_pressure * 100.0_wp * &
[1361]1200                         ( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
[1353]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)
[1241]1203             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
1204          ENDDO
1205!
1206!--       Compute reference density
[1353]1207          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
[1241]1208       ENDIF
1209
[1361]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
[1241]1217
[1115]1218!
[1361]1219!--    Use 1d arrays
[1115]1220       q_1d(:)  = q(:,j,i)
1221       pt_1d(:) = pt(:,j,i)
1222       qc_1d(:) = qc(:,j,i)
1223       nc_1d(:) = nc_const
[1822]1224       IF ( microphysics_seifert )  THEN
[1115]1225          qr_1d(:) = qr(:,j,i)
1226          nr_1d(:) = nr(:,j,i)
1227       ENDIF
[1361]1228
[1115]1229!
[1822]1230!--    Reset precipitation rate
1231       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
1232
1233!
[1115]1234!--    Compute cloud physics
[1822]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 )
[1115]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
[1822]1248          IF ( drizzle )  CALL sedimentation_cloud( i,j )
[1361]1249
[1691]1250       ENDIF
1251
[1822]1252       CALL calc_precipitation_amount( i,j )
1253
[1115]1254!
[1361]1255!--    Store results on the 3d arrays
1256       q(:,j,i)  = q_1d(:)
1257       pt(:,j,i) = pt_1d(:)
[1822]1258       IF ( microphysics_seifert )  THEN
[1361]1259          qr(:,j,i) = qr_1d(:)
1260          nr(:,j,i) = nr_1d(:)
[1115]1261       ENDIF
1262
1263    END SUBROUTINE microphysics_control_ij
1264
[1682]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!------------------------------------------------------------------------------!
[1115]1274    SUBROUTINE adjust_cloud_ij( i, j )
1275
[1320]1276       USE arrays_3d,                                                          &
[1361]1277           ONLY:  qr_1d, nr_1d
[1115]1278
[1320]1279       USE cloud_parameters,                                                   &
[1822]1280           ONLY:  eps_sb, xrmin, xrmax, hyrho
[1320]1281
1282       USE indices,                                                            &
[1822]1283           ONLY:  nzb_s_inner, nzt
[1320]1284
1285       USE kinds
1286
[1115]1287       IMPLICIT NONE
1288
[1682]1289       INTEGER(iwp) ::  i                 !<
1290       INTEGER(iwp) ::  j                 !<
1291       INTEGER(iwp) ::  k                 !<
1292
[1115]1293       DO  k = nzb_s_inner(j,i)+1, nzt
[1022]1294
[1361]1295          IF ( qr_1d(k) <= eps_sb )  THEN
1296             qr_1d(k) = 0.0_wp
1297             nr_1d(k) = 0.0_wp
[1065]1298          ELSE
[1022]1299!
[1048]1300!--          Adjust number of raindrops to avoid nonlinear effects in
1301!--          sedimentation and evaporation of rain drops due to too small or
[1065]1302!--          too big weights of rain drops (Stevens and Seifert, 2008).
[1361]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
[1048]1307             ENDIF
[1115]1308
[1022]1309          ENDIF
[1115]1310
[1022]1311       ENDDO
1312
[1115]1313    END SUBROUTINE adjust_cloud_ij
[1022]1314
[1106]1315
[1682]1316!------------------------------------------------------------------------------!
1317! Description:
1318! ------------
1319!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
1320!------------------------------------------------------------------------------!
[1005]1321    SUBROUTINE autoconversion_ij( i, j )
[1000]1322
[1320]1323       USE arrays_3d,                                                          &
1324           ONLY:  diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d
[1115]1325
[1320]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,        &
[1822]1328                  c_const, dpirho_l, eps_sb, hyrho, kin_vis_air, k_cc, x0
[1320]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,                                                            &
[1822]1337           ONLY:  nzb_s_inner, nzt
[1320]1338
1339       USE kinds
1340
[1000]1341       IMPLICIT NONE
1342
[1682]1343       INTEGER(iwp) ::  i                 !<
1344       INTEGER(iwp) ::  j                 !<
1345       INTEGER(iwp) ::  k                 !<
[1000]1346
[1682]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                !<
[1106]1360
[1115]1361       DO  k = nzb_s_inner(j,i)+1, nzt
[1000]1362
[1115]1363          IF ( qc_1d(k) > eps_sb )  THEN
[1361]1364
1365             k_au = k_cc / ( 20.0_wp * x0 )
[1012]1366!
[1048]1367!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]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) )
[1012]1370!
1371!--          Universal function for autoconversion process
1372!--          (Seifert and Beheng, 2006):
[1361]1373             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
[1012]1374!
1375!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
[1353]1376!--          (Use constant nu_c = 1.0_wp instead?)
[1361]1377             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
[1012]1378!
1379!--          Mean weight of cloud droplets:
[1115]1380             xc = hyrho(k) * qc_1d(k) / nc_1d(k)
[1012]1381!
[1065]1382!--          Parameterized turbulence effects on autoconversion (Seifert,
1383!--          Nuijens and Stevens, 2010)
1384             IF ( turbulence )  THEN
1385!
1386!--             Weight averaged radius of cloud droplets:
[1353]1387                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
[1065]1388
[1353]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 )
[1065]1392!
1393!--             Mixing length (neglecting distance to ground and stratification)
[1334]1394                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
[1065]1395!
1396!--             Limit dissipation rate according to Seifert, Nuijens and
1397!--             Stevens (2010)
[1361]1398                dissipation = MIN( 0.06_wp, diss(k,j,i) )
[1065]1399!
1400!--             Compute Taylor-microscale Reynolds number:
[1361]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 )
[1065]1405!
1406!--             The factor of 1.0E4 is needed to convert the dissipation rate
1407!--             from m2 s-3 to cm2 s-3.
[1361]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                              )
[1065]1416             ENDIF
1417!
[1012]1418!--          Autoconversion rate (Seifert and Beheng, 2006):
[1361]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 ) *       &
[1115]1422                       rho_surface
1423             autocon = MIN( autocon, qc_1d(k) / dt_micro )
[1106]1424
[1115]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
[1005]1429          ENDIF
[1000]1430
1431       ENDDO
1432
[1005]1433    END SUBROUTINE autoconversion_ij
1434
[1822]1435!------------------------------------------------------------------------------!
1436! Description:
1437! ------------
1438!> Autoconversion process (Kessler, 1969).
1439!------------------------------------------------------------------------------!
1440    SUBROUTINE autoconversion_kessler_ij( i, j )
[1106]1441
[1822]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
[1682]1486!------------------------------------------------------------------------------!
1487! Description:
1488! ------------
1489!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
1490!------------------------------------------------------------------------------!
[1005]1491    SUBROUTINE accretion_ij( i, j )
1492
[1320]1493       USE arrays_3d,                                                          &
1494           ONLY:  diss, qc_1d, qr_1d
[1115]1495
[1320]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,                                                            &
[1822]1503           ONLY:  nzb_s_inner, nzt
[1320]1504
1505       USE kinds
1506
[1005]1507       IMPLICIT NONE
1508
[1682]1509       INTEGER(iwp) ::  i                 !<
1510       INTEGER(iwp) ::  j                 !<
1511       INTEGER(iwp) ::  k                 !<
[1005]1512
[1682]1513       REAL(wp)     ::  accr              !<
1514       REAL(wp)     ::  k_cr              !<
1515       REAL(wp)     ::  phi_ac            !<
1516       REAL(wp)     ::  tau_cloud         !<
[1320]1517
[1115]1518       DO  k = nzb_s_inner(j,i)+1, nzt
1519          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
[1012]1520!
[1048]1521!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]1522             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 
[1012]1523!
1524!--          Universal function for accretion process
[1048]1525!--          (Seifert and Beheng, 2001):
[1361]1526             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
[1012]1527!
[1065]1528!--          Parameterized turbulence effects on autoconversion (Seifert,
1529!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
[1361]1530!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1065]1531             IF ( turbulence )  THEN
[1361]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                               )
[1065]1536             ELSE
1537                k_cr = k_cr0                       
1538             ENDIF
1539!
[1012]1540!--          Accretion rate (Seifert and Beheng, 2006):
[1361]1541             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
[1115]1542             accr = MIN( accr, qc_1d(k) / dt_micro )
[1106]1543
[1115]1544             qr_1d(k) = qr_1d(k) + accr * dt_micro 
1545             qc_1d(k) = qc_1d(k) - accr * dt_micro
1546
[1005]1547          ENDIF
[1106]1548
[1005]1549       ENDDO
1550
[1000]1551    END SUBROUTINE accretion_ij
1552
[1005]1553
[1682]1554!------------------------------------------------------------------------------!
1555! Description:
1556! ------------
1557!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
1558!------------------------------------------------------------------------------!
[1005]1559    SUBROUTINE selfcollection_breakup_ij( i, j )
1560
[1320]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,                                                            &
[1822]1571           ONLY:  nzb_s_inner, nzt
[1320]1572
1573       USE kinds
[1005]1574   
1575       IMPLICIT NONE
1576
[1682]1577       INTEGER(iwp) ::  i                 !<
1578       INTEGER(iwp) ::  j                 !<
1579       INTEGER(iwp) ::  k                 !<
[1005]1580
[1682]1581       REAL(wp)     ::  breakup           !<
1582       REAL(wp)     ::  dr                !<
1583       REAL(wp)     ::  phi_br            !<
1584       REAL(wp)     ::  selfcoll          !<
[1320]1585
[1115]1586       DO  k = nzb_s_inner(j,i)+1, nzt
1587          IF ( qr_1d(k) > eps_sb )  THEN
[1012]1588!
[1115]1589!--          Selfcollection rate (Seifert and Beheng, 2001):
[1361]1590             selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
[1012]1591!
[1115]1592!--          Weight averaged diameter of rain drops:
[1334]1593             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]1594!
[1048]1595!--          Collisional breakup rate (Seifert, 2008):
[1353]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 )
[1005]1599             ELSE
[1353]1600                breakup = 0.0_wp
[1005]1601             ENDIF
[1048]1602
[1115]1603             selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
1604             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
[1106]1605
[1005]1606          ENDIF         
1607       ENDDO
1608
1609    END SUBROUTINE selfcollection_breakup_ij
1610
[1106]1611
[1682]1612!------------------------------------------------------------------------------!
1613! Description:
1614! ------------
1615!> Evaporation of precipitable water. Condensation is neglected for
1616!> precipitable water. Call for grid point i,j
1617!------------------------------------------------------------------------------!
[1012]1618    SUBROUTINE evaporation_rain_ij( i, j )
1619
[1320]1620       USE arrays_3d,                                                          &
1621           ONLY:  hyp, nr_1d, pt_1d, q_1d,  qc_1d, qr_1d
[1048]1622
[1320]1623       USE cloud_parameters,                                                   &
1624           ONLY:  a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
[1822]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,             &
[1320]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,                                                            &
[1822]1636           ONLY:  nzb_s_inner, nzt
[1320]1637
1638       USE kinds
1639
[1012]1640       IMPLICIT NONE
1641
[1682]1642       INTEGER(iwp) ::  i                 !<
1643       INTEGER(iwp) ::  j                 !<
1644       INTEGER(iwp) ::  k                 !<
[1012]1645
[1682]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                !<
[1320]1663
[1115]1664       DO  k = nzb_s_inner(j,i)+1, nzt
1665          IF ( qr_1d(k) > eps_sb )  THEN
[1012]1666!
1667!--          Actual liquid water temperature:
[1115]1668             t_l = t_d_pt(k) * pt_1d(k)
[1012]1669!
1670!--          Saturation vapor pressure at t_l:
[1361]1671             e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /          &
1672                                    ( t_l - 35.86_wp )                         &
1673                                  )
[1012]1674!
1675!--          Computation of saturation humidity:
[1361]1676             q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
[1353]1677             alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
[1361]1678             q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
[1012]1679!
[1106]1680!--          Supersaturation:
[1361]1681             sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
[1012]1682!
[1361]1683!--          Evaporation needs only to be calculated in subsaturated regions
1684             IF ( sat < 0.0_wp )  THEN
[1012]1685!
[1361]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                                  )
[1012]1693!
[1361]1694!--             Mean weight of rain drops
1695                xr = hyrho(k) * qr_1d(k) / nr_1d(k)
[1115]1696!
[1361]1697!--             Weight averaged diameter of rain drops:
1698                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]1699!
[1361]1700!--             Compute ventilation factor and intercept parameter
1701!--             (Seifert and Beheng, 2006; Seifert, 2008):
1702                IF ( ventilation_effect )  THEN
[1115]1703!
[1361]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
[1115]1712
[1361]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
[1012]1741!
[1361]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 )
[1106]1747
[1361]1748                qr_1d(k) = qr_1d(k) + evap    * dt_micro
1749                nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
[1115]1750
[1361]1751             ENDIF
[1012]1752          ENDIF         
[1106]1753
[1012]1754       ENDDO
1755
1756    END SUBROUTINE evaporation_rain_ij
1757
[1106]1758
[1682]1759!------------------------------------------------------------------------------!
1760! Description:
1761! ------------
1762!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
1763!> Call for grid point i,j
1764!------------------------------------------------------------------------------!
[1012]1765    SUBROUTINE sedimentation_cloud_ij( i, j )
1766
[1320]1767       USE arrays_3d,                                                          &
1768           ONLY:  ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d
1769
1770       USE cloud_parameters,                                                   &
[1691]1771           ONLY:  eps_sb, hyrho, l_d_cp, prr, pt_d_t, sed_qc_const
[1320]1772
1773       USE control_parameters,                                                 &
[1822]1774           ONLY:  call_microphysics_at_all_substeps, dt_micro,                 &
1775                  intermediate_timestep_count
[1320]1776
1777       USE indices,                                                            &
1778           ONLY:  nzb, nzb_s_inner, nzt
1779
1780       USE kinds
[1012]1781       
[1691]1782       USE statistics,                                                         &
1783           ONLY:  weight_substep
1784
[1012]1785       IMPLICIT NONE
1786
[1682]1787       INTEGER(iwp) ::  i                 !<
1788       INTEGER(iwp) ::  j                 !<
1789       INTEGER(iwp) ::  k                 !<
[1106]1790
[1682]1791       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc  !<
[1115]1792
[1353]1793       sed_qc(nzt+1) = 0.0_wp
[1012]1794
[1115]1795       DO  k = nzt, nzb_s_inner(j,i)+1, -1
1796          IF ( qc_1d(k) > eps_sb )  THEN
[1361]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 )
[1115]1799          ELSE
[1353]1800             sed_qc(k) = 0.0_wp
[1012]1801          ENDIF
[1115]1802
[1361]1803          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /          &
1804                                      dt_micro + sed_qc(k+1)                   &
1805                         )
[1115]1806
[1361]1807          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
[1115]1808                                hyrho(k) * dt_micro
[1361]1809          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      & 
[1115]1810                                hyrho(k) * dt_micro
[1361]1811          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
[1115]1812                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
1813
[1691]1814!
1815!--       Compute the precipitation rate of cloud (fog) droplets
[1822]1816          IF ( call_microphysics_at_all_substeps )  THEN
1817             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *               &
[1691]1818                             weight_substep(intermediate_timestep_count)
[1822]1819          ELSE
1820             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
[1691]1821          ENDIF
1822
[1012]1823       ENDDO
1824
1825    END SUBROUTINE sedimentation_cloud_ij
1826
[1106]1827
[1682]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!------------------------------------------------------------------------------!
[1012]1834    SUBROUTINE sedimentation_rain_ij( i, j )
1835
[1320]1836       USE arrays_3d,                                                          &
1837           ONLY:  ddzu, dzu, nr_1d, pt_1d, q_1d, qr_1d
1838
1839       USE cloud_parameters,                                                   &
[1822]1840           ONLY:  a_term, b_term, c_term, dpirho_l, eps_sb, hyrho,             &
1841                  limiter_sedimentation, l_d_cp, prr, pt_d_t
[1320]1842
1843       USE control_parameters,                                                 &
[1822]1844           ONLY:  call_microphysics_at_all_substeps, dt_micro,                 &
1845                  intermediate_timestep_count
[1320]1846
1847       USE indices,                                                            &
1848           ONLY:  nzb, nzb_s_inner, nzt
1849
1850       USE kinds
1851
1852       USE statistics,                                                         &
1853           ONLY:  weight_substep
[1012]1854       
1855       IMPLICIT NONE
1856
[1682]1857       INTEGER(iwp) ::  i                          !<
1858       INTEGER(iwp) ::  j                          !<
1859       INTEGER(iwp) ::  k                          !<
1860       INTEGER(iwp) ::  k_run                      !<
[1012]1861
[1682]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                      !<
[1320]1871
[1682]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     !<
[1320]1880
[1012]1881!
[1065]1882!--    Compute velocities
1883       DO  k = nzb_s_inner(j,i)+1, nzt
[1115]1884          IF ( qr_1d(k) > eps_sb )  THEN
1885!
1886!--          Weight averaged diameter of rain drops:
[1334]1887             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]1888!
1889!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1890!--          Stevens and Seifert, 2008):
[1353]1891             mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
[1115]1892!
1893!--          Slope parameter of gamma distribution (Seifert, 2008):
[1361]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
[1115]1896
[1361]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                          )
[1065]1909          ELSE
[1353]1910             w_nr(k) = 0.0_wp
1911             w_qr(k) = 0.0_wp
[1065]1912          ENDIF
1913       ENDDO
[1048]1914!
[1065]1915!--    Adjust boundary values
[1115]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)
[1353]1918       w_nr(nzt+1) = 0.0_wp
1919       w_qr(nzt+1) = 0.0_wp
[1065]1920!
1921!--    Compute Courant number
[1115]1922       DO  k = nzb_s_inner(j,i)+1, nzt
[1361]1923          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
[1115]1924                    dt_micro * ddzu(k)
[1361]1925          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
[1115]1926                    dt_micro * ddzu(k)
1927       ENDDO     
[1065]1928!
1929!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
1930       IF ( limiter_sedimentation )  THEN
1931
[1115]1932          DO k = nzb_s_inner(j,i)+1, nzt
[1646]1933             d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
[1115]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)
[1065]1936
[1361]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 ) )
[1065]1940
[1646]1941             d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
[1115]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)
[1065]1944
[1361]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 ) )
[1022]1948          ENDDO
[1048]1949
[1065]1950       ELSE
[1106]1951
[1353]1952          nr_slope = 0.0_wp
1953          qr_slope = 0.0_wp
[1106]1954
[1065]1955       ENDIF
[1115]1956
[1353]1957       sed_nr(nzt+1) = 0.0_wp
1958       sed_qr(nzt+1) = 0.0_wp
[1065]1959!
1960!--    Compute sedimentation flux
[1115]1961       DO  k = nzt, nzb_s_inner(j,i)+1, -1
[1065]1962!
1963!--       Sum up all rain drop number densities which contribute to the flux
1964!--       through k-1/2
[1353]1965          flux  = 0.0_wp
1966          z_run = 0.0_wp ! height above z(k)
[1065]1967          k_run = k
[1346]1968          c_run = MIN( 1.0_wp, c_nr(k) )
[1353]1969          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
[1361]1970             flux  = flux + hyrho(k_run) *                                     &
1971                     ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) *   &
[1353]1972                     0.5_wp ) * c_run * dzu(k_run)
[1065]1973             z_run = z_run + dzu(k_run)
1974             k_run = k_run + 1
[1346]1975             c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
[1022]1976          ENDDO
1977!
[1065]1978!--       It is not allowed to sediment more rain drop number density than
1979!--       available
[1361]1980          flux = MIN( flux,                                                    &
[1115]1981                      hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
[1065]1982
[1115]1983          sed_nr(k) = flux / dt_micro
[1361]1984          nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
1985                                    hyrho(k) * dt_micro
[1065]1986!
1987!--       Sum up all rain water content which contributes to the flux
1988!--       through k-1/2
[1353]1989          flux  = 0.0_wp
1990          z_run = 0.0_wp ! height above z(k)
[1065]1991          k_run = k
[1346]1992          c_run = MIN( 1.0_wp, c_qr(k) )
[1106]1993
[1361]1994          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
[1106]1995
[1361]1996             flux  = flux + hyrho(k_run) *                                     &
1997                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) *   &
[1353]1998                     0.5_wp ) * c_run * dzu(k_run)
[1065]1999             z_run = z_run + dzu(k_run)
2000             k_run = k_run + 1
[1346]2001             c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
[1106]2002
[1065]2003          ENDDO
2004!
2005!--       It is not allowed to sediment more rain water content than available
[1361]2006          flux = MIN( flux,                                                    &
[1115]2007                      hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
[1065]2008
[1115]2009          sed_qr(k) = flux / dt_micro
2010
[1361]2011          qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
[1115]2012                                hyrho(k) * dt_micro
[1361]2013          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
[1115]2014                                hyrho(k) * dt_micro 
[1361]2015          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
[1115]2016                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
[1065]2017!
2018!--       Compute the rain rate
[1361]2019          IF ( call_microphysics_at_all_substeps )  THEN
[1691]2020             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)                    &
2021                          * weight_substep(intermediate_timestep_count)
[1361]2022          ELSE
[1691]2023             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
[1361]2024          ENDIF
2025
[1065]2026       ENDDO
[1115]2027
[1691]2028    END SUBROUTINE sedimentation_rain_ij
[1012]2029
[1691]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,&
[1822]2045                  precipitation_amount_interval, time_do2d_xy
[1691]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
[1361]2063          precipitation_amount(j,i) = precipitation_amount(j,i) +              &
2064                                      prr(nzb_s_inner(j,i)+1,j,i) *            &
[1115]2065                                      hyrho(nzb_s_inner(j,i)+1) * dt_3d
[1048]2066       ENDIF
2067
[1691]2068    END SUBROUTINE calc_precipitation_amount_ij
[1012]2069
[1361]2070!------------------------------------------------------------------------------!
[1682]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.
[1361]2076!------------------------------------------------------------------------------!
[1012]2077    FUNCTION gamm( xx ) 
[1048]2078       
[1320]2079       USE cloud_parameters,                                                   &
2080           ONLY:  cof, stp
2081
2082       USE kinds
2083
[1012]2084       IMPLICIT NONE
[1106]2085
[1682]2086       INTEGER(iwp) ::  j            !<
[1320]2087
[1682]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       !<
[1320]2094
[1012]2095       x_gamm = xx 
2096       y_gamm = x_gamm 
[1353]2097       tmp = x_gamm + 5.5_wp
2098       tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp 
[1334]2099       ser = 1.000000000190015_wp
[1106]2100
2101       DO  j = 1, 6 
[1353]2102          y_gamm = y_gamm + 1.0_wp 
[1012]2103          ser    = ser + cof( j ) / y_gamm 
[1106]2104       ENDDO
2105
[1012]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 
[1106]2111
[1012]2112       RETURN
2113 
2114    END FUNCTION gamm 
2115
2116 END MODULE microphysics_mod
Note: See TracBrowser for help on using the repository browser.