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

Last change on this file since 1335 was 1335, checked in by heinze, 10 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 38.0 KB
RevLine 
[1000]1 MODULE microphysics_mod
2
[1093]3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
[1310]17! Copyright 1997-2014 Leibniz Universitaet Hannover
[1093]18!--------------------------------------------------------------------------------!
19!
[1000]20! Current revisions:
[1092]21! ------------------
[1335]22!
23!
[1321]24! Former revisions:
25! -----------------
26! $Id: microphysics.f90 1335 2014-03-25 12:23:36Z heinze $
27!
[1335]28! 1334 2014-03-25 12:21:40Z heinze
29! Bugfix: REAL constants provided with KIND-attribute
30!
[1323]31! 1322 2014-03-20 16:38:49Z raasch
32! REAL constants defined as wp-kind
33!
[1321]34! 1320 2014-03-20 08:40:49Z raasch
[1320]35! ONLY-attribute added to USE-statements,
36! kind-parameters added to all INTEGER and REAL declaration statements,
37! kinds are defined in new module kinds,
38! comment fields (!:) to be used for variable explanations added to
39! all variable declaration statements
[1000]40!
[1242]41! 1241 2013-10-30 11:36:58Z heinze
42! hyp and rho have to be calculated at each time step if data from external
43! file LSF_DATA are used
44!
[1116]45! 1115 2013-03-26 18:16:16Z hoffmann
46! microphyical tendencies are calculated in microphysics_control in an optimized
47! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
48!
[1107]49! 1106 2013-03-04 05:31:38Z raasch
50! small changes in code formatting
51!
[1093]52! 1092 2013-02-02 11:24:22Z raasch
53! unused variables removed
54! file put under GPL
55!
[1066]56! 1065 2012-11-22 17:42:36Z hoffmann
57! Sedimentation process implemented according to Stevens and Seifert (2008).
[1115]58! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
[1066]59! and Stevens, 2010).
60!
[1054]61! 1053 2012-11-13 17:11:03Z hoffmann
62! initial revision
[1000]63!
64! Description:
65! ------------
66! Calculate cloud microphysics according to the two moment bulk
67! scheme by Seifert and Beheng (2006).
68!------------------------------------------------------------------------------!
69
70    PRIVATE
[1115]71    PUBLIC microphysics_control
[1000]72
[1115]73    INTERFACE microphysics_control
74       MODULE PROCEDURE microphysics_control
75       MODULE PROCEDURE microphysics_control_ij
76    END INTERFACE microphysics_control
[1022]77
[1115]78    INTERFACE adjust_cloud
79       MODULE PROCEDURE adjust_cloud
80       MODULE PROCEDURE adjust_cloud_ij
81    END INTERFACE adjust_cloud
82
[1000]83    INTERFACE autoconversion
84       MODULE PROCEDURE autoconversion
85       MODULE PROCEDURE autoconversion_ij
86    END INTERFACE autoconversion
87
88    INTERFACE accretion
89       MODULE PROCEDURE accretion
90       MODULE PROCEDURE accretion_ij
91    END INTERFACE accretion
[1005]92
93    INTERFACE selfcollection_breakup
94       MODULE PROCEDURE selfcollection_breakup
95       MODULE PROCEDURE selfcollection_breakup_ij
96    END INTERFACE selfcollection_breakup
[1012]97
98    INTERFACE evaporation_rain
99       MODULE PROCEDURE evaporation_rain
100       MODULE PROCEDURE evaporation_rain_ij
101    END INTERFACE evaporation_rain
102
103    INTERFACE sedimentation_cloud
104       MODULE PROCEDURE sedimentation_cloud
105       MODULE PROCEDURE sedimentation_cloud_ij
106    END INTERFACE sedimentation_cloud
[1000]107 
[1012]108    INTERFACE sedimentation_rain
109       MODULE PROCEDURE sedimentation_rain
110       MODULE PROCEDURE sedimentation_rain_ij
111    END INTERFACE sedimentation_rain
112
[1000]113 CONTAINS
114
115
116!------------------------------------------------------------------------------!
117! Call for all grid points
118!------------------------------------------------------------------------------!
[1115]119    SUBROUTINE microphysics_control
[1022]120
121       USE arrays_3d
[1241]122       USE cloud_parameters
[1115]123       USE control_parameters
[1241]124       USE grid_variables
[1115]125       USE indices
[1320]126       USE kinds
[1115]127       USE statistics
128
129       IMPLICIT NONE
130
[1320]131       INTEGER(iwp) ::  i                 !:
132       INTEGER(iwp) ::  j                 !:
133       INTEGER(iwp) ::  k                 !:
[1115]134 
135       DO  i = nxl, nxr
136          DO  j = nys, nyn
137             DO  k = nzb_s_inner(j,i)+1, nzt
138
139             ENDDO
140          ENDDO
141       ENDDO
142
143    END SUBROUTINE microphysics_control
144
145    SUBROUTINE adjust_cloud
146
147       USE arrays_3d
[1022]148       USE cloud_parameters
149       USE indices
[1320]150       USE kinds
[1022]151
152       IMPLICIT NONE
153
[1320]154       INTEGER(iwp) ::  i                 !:
155       INTEGER(iwp) ::  j                 !:
156       INTEGER(iwp) ::  k                 !:
[1022]157
158 
159       DO  i = nxl, nxr
160          DO  j = nys, nyn
[1115]161             DO  k = nzb_s_inner(j,i)+1, nzt
[1022]162
163             ENDDO
164          ENDDO
165       ENDDO
166
[1115]167    END SUBROUTINE adjust_cloud
[1022]168
[1106]169
[1000]170    SUBROUTINE autoconversion
171
172       USE arrays_3d
173       USE cloud_parameters
[1115]174       USE control_parameters
175       USE grid_variables
[1000]176       USE indices
[1320]177       USE kinds
[1000]178
179       IMPLICIT NONE
180
[1320]181       INTEGER(iwp) ::  i                 !:
182       INTEGER(iwp) ::  j                 !:
183       INTEGER(iwp) ::  k                 !:
[1000]184
185       DO  i = nxl, nxr
186          DO  j = nys, nyn
[1115]187             DO  k = nzb_s_inner(j,i)+1, nzt
[1000]188
189             ENDDO
190          ENDDO
191       ENDDO
192
193    END SUBROUTINE autoconversion
194
[1106]195
[1005]196    SUBROUTINE accretion
[1000]197
198       USE arrays_3d
199       USE cloud_parameters
[1115]200       USE control_parameters
[1000]201       USE indices
[1320]202       USE kinds
[1005]203
[1000]204       IMPLICIT NONE
205
[1320]206       INTEGER(iwp) ::  i                 !:
207       INTEGER(iwp) ::  j                 !:
208       INTEGER(iwp) ::  k                 !:
[1000]209
[1005]210       DO  i = nxl, nxr
211          DO  j = nys, nyn
[1115]212             DO  k = nzb_s_inner(j,i)+1, nzt
[1000]213
[1005]214             ENDDO
215          ENDDO
[1000]216       ENDDO
217
[1005]218    END SUBROUTINE accretion
[1000]219
[1106]220
[1005]221    SUBROUTINE selfcollection_breakup
[1000]222
223       USE arrays_3d
224       USE cloud_parameters
[1115]225       USE control_parameters
[1000]226       USE indices
[1320]227       USE kinds
[1000]228
229       IMPLICIT NONE
230
[1320]231       INTEGER(iwp) ::  i                 !:
232       INTEGER(iwp) ::  j                 !:
233       INTEGER(iwp) ::  k                 !:
[1000]234
235 
236       DO  i = nxl, nxr
237          DO  j = nys, nyn
[1115]238             DO  k = nzb_s_inner(j,i)+1, nzt
[1000]239
240             ENDDO
241          ENDDO
242       ENDDO
243
[1005]244    END SUBROUTINE selfcollection_breakup
[1000]245
[1106]246
[1012]247    SUBROUTINE evaporation_rain
[1000]248
[1012]249       USE arrays_3d
250       USE cloud_parameters
251       USE constants
[1115]252       USE control_parameters
[1012]253       USE indices
[1320]254       USE kinds
[1012]255
256       IMPLICIT NONE
257
[1320]258       INTEGER(iwp) ::  i                 !:
259       INTEGER(iwp) ::  j                 !:
260       INTEGER(iwp) ::  k                 !:
[1012]261 
262       DO  i = nxl, nxr
263          DO  j = nys, nyn
[1115]264             DO  k = nzb_s_inner(j,i)+1, nzt
[1012]265
266             ENDDO
267          ENDDO
268       ENDDO
269
270    END SUBROUTINE evaporation_rain
271
[1106]272
[1012]273    SUBROUTINE sedimentation_cloud
274
275       USE arrays_3d
276       USE cloud_parameters
277       USE constants
[1115]278       USE control_parameters
[1012]279       USE indices
[1320]280       USE kinds
[1012]281
282       IMPLICIT NONE
283
[1320]284       INTEGER(iwp) ::  i                 !:
285       INTEGER(iwp) ::  j                 !:
286       INTEGER(iwp) ::  k                 !:
[1012]287 
288       DO  i = nxl, nxr
289          DO  j = nys, nyn
[1115]290             DO  k = nzb_s_inner(j,i)+1, nzt
[1012]291
292             ENDDO
293          ENDDO
294       ENDDO
295
296    END SUBROUTINE sedimentation_cloud
297
[1106]298
[1012]299    SUBROUTINE sedimentation_rain
300
301       USE arrays_3d
302       USE cloud_parameters
303       USE constants
[1115]304       USE control_parameters
[1012]305       USE indices
[1320]306       USE kinds
[1115]307       USE statistics
[1012]308
309       IMPLICIT NONE
310
[1320]311       INTEGER(iwp) ::  i                 !:
312       INTEGER(iwp) ::  j                 !:
313       INTEGER(iwp) ::  k                 !:
[1012]314 
315       DO  i = nxl, nxr
316          DO  j = nys, nyn
[1115]317             DO  k = nzb_s_inner(j,i)+1, nzt
[1012]318
319             ENDDO
320          ENDDO
321       ENDDO
322
323    END SUBROUTINE sedimentation_rain
324
325
[1000]326!------------------------------------------------------------------------------!
327! Call for grid point i,j
328!------------------------------------------------------------------------------!
[1022]329
[1115]330    SUBROUTINE microphysics_control_ij( i, j )
331
[1320]332       USE arrays_3d,                                                          &
333           ONLY:  hyp, nc_1d,  nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc,     &
334                  qc_1d, qr, qr_1d, tend_nr, tend_pt, tend_q, tend_qr, zu
[1115]335
[1320]336       USE cloud_parameters,                                                   &
337           ONLY:  cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
338
339       USE control_parameters,                                                 &
340           ONLY:  drizzle, dt_3d, dt_micro, g, intermediate_timestep_count,    &
341                  large_scale_forcing, lsf_surf, precipitation, pt_surface,    &
342                  rho_surface,surface_pressure
343
344       USE indices,                                                            &
345           ONLY:  nzb, nzt
346
347       USE kinds
348
349       USE statistics,                                                         &
350           ONLY:  weight_pres
351
[1022]352       IMPLICIT NONE
353
[1320]354       INTEGER(iwp) ::  i                 !:
355       INTEGER(iwp) ::  j                 !:
356       INTEGER(iwp) ::  k                 !:
[1115]357
[1320]358       REAL(wp)     ::  t_surface         !:
359
[1241]360       IF ( large_scale_forcing .AND. lsf_surf ) THEN
361!
362!--       Calculate:
363!--       pt / t : ratio of potential and actual temperature (pt_d_t)
364!--       t / pt : ratio of actual and potential temperature (t_d_pt)
365!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
[1334]366          t_surface = pt_surface * ( surface_pressure / 1000.0 )**0.286_wp
[1241]367          DO  k = nzb, nzt+1
368             hyp(k)    = surface_pressure * 100.0 * &
[1334]369                         ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0_wp/0.286_wp)
370             pt_d_t(k) = ( 100000.0 / hyp(k) )**0.286_wp
[1241]371             t_d_pt(k) = 1.0 / pt_d_t(k)
372             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
373          ENDDO
374!
375!--       Compute reference density
376          rho_surface = surface_pressure * 100.0 / ( r_d * t_surface )
377       ENDIF
378
379
[1115]380       dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
381!
382!--    Adjust unrealistic values
383       IF ( precipitation )  CALL adjust_cloud( i,j ) 
384!
385!--    Use 1-d arrays
386       q_1d(:)  = q(:,j,i)
387       pt_1d(:) = pt(:,j,i)
388       qc_1d(:) = qc(:,j,i)
389       nc_1d(:) = nc_const
390       IF ( precipitation )  THEN
391          qr_1d(:) = qr(:,j,i)
392          nr_1d(:) = nr(:,j,i)
393       ENDIF
394!
395!--    Compute cloud physics
396       IF ( precipitation )  THEN
397          CALL autoconversion( i,j )
398          CALL accretion( i,j )
399          CALL selfcollection_breakup( i,j )
400          CALL evaporation_rain( i,j )
401          CALL sedimentation_rain( i,j )
402       ENDIF
403
404       IF ( drizzle )  CALL sedimentation_cloud( i,j )
405!
406!--    Derive tendencies
407       tend_q(:,j,i)  = ( q_1d(:) - q(:,j,i) ) / dt_micro
408       tend_pt(:,j,i) = ( pt_1d(:) - pt(:,j,i) ) / dt_micro
409       IF ( precipitation )  THEN
410          tend_qr(:,j,i) = ( qr_1d(:) - qr(:,j,i) ) / dt_micro
411          tend_nr(:,j,i) = ( nr_1d(:) - nr(:,j,i) ) / dt_micro
412       ENDIF
413
414    END SUBROUTINE microphysics_control_ij
415
416    SUBROUTINE adjust_cloud_ij( i, j )
417
[1320]418       USE arrays_3d,                                                          &
419           ONLY:  qr, nr
[1115]420
[1320]421       USE cloud_parameters,                                                   &
422           ONLY:  eps_sb, xrmin, xrmax, hyrho, k_cc, x0
423
424       USE indices,                                                            &
425           ONLY:  nzb, nzb_s_inner, nzt
426
427       USE kinds
428
[1115]429       IMPLICIT NONE
430
[1320]431       INTEGER(iwp) ::  i                 !:
432       INTEGER(iwp) ::  j                 !:
433       INTEGER(iwp) ::  k                 !:
[1115]434!
435!--    Adjust number of raindrops to avoid nonlinear effects in
436!--    sedimentation and evaporation of rain drops due to too small or
437!--    too big weights of rain drops (Stevens and Seifert, 2008).
438!--    The same procedure is applied to cloud droplets if they are determined
439!--    prognostically. 
440       DO  k = nzb_s_inner(j,i)+1, nzt
[1022]441
[1065]442          IF ( qr(k,j,i) <= eps_sb )  THEN
443             qr(k,j,i) = 0.0
[1115]444             nr(k,j,i) = 0.0
[1065]445          ELSE
[1022]446!
[1048]447!--          Adjust number of raindrops to avoid nonlinear effects in
448!--          sedimentation and evaporation of rain drops due to too small or
[1065]449!--          too big weights of rain drops (Stevens and Seifert, 2008).
450             IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
451                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
452             ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
453                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
[1048]454             ENDIF
[1115]455
[1022]456          ENDIF
[1115]457
[1022]458       ENDDO
459
[1115]460    END SUBROUTINE adjust_cloud_ij
[1022]461
[1106]462
[1005]463    SUBROUTINE autoconversion_ij( i, j )
[1000]464
[1320]465       USE arrays_3d,                                                          &
466           ONLY:  diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d
[1115]467
[1320]468       USE cloud_parameters,                                                   &
469           ONLY:  a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3,        &
470                  c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, x0
471
472       USE control_parameters,                                                 &
473           ONLY:  dt_micro, rho_surface, turbulence
474
475       USE grid_variables,                                                     &
476           ONLY:  dx, dy
477
478       USE indices,                                                            &
479           ONLY:  nzb, nzb_s_inner, nzt
480
481       USE kinds
482
[1000]483       IMPLICIT NONE
484
[1320]485       INTEGER(iwp) ::  i                 !:
486       INTEGER(iwp) ::  j                 !:
487       INTEGER(iwp) ::  k                 !:
[1000]488
[1320]489       REAL(wp)     ::  alpha_cc          !:                   
490       REAL(wp)     ::  autocon           !:
491       REAL(wp)     ::  epsilon           !:
492       REAL(wp)     ::  k_au              !:
493       REAL(wp)     ::  l_mix             !:
494       REAL(wp)     ::  nu_c              !:
495       REAL(wp)     ::  phi_au            !:
496       REAL(wp)     ::  r_cc              !:
497       REAL(wp)     ::  rc                !:
498       REAL(wp)     ::  re_lambda         !:
499       REAL(wp)     ::  selfcoll          !:
500       REAL(wp)     ::  sigma_cc          !:
501       REAL(wp)     ::  tau_cloud         !:
502       REAL(wp)     ::  xc                !:
[1106]503
[1005]504       k_au = k_cc / ( 20.0 * x0 )
505
[1115]506       DO  k = nzb_s_inner(j,i)+1, nzt
[1000]507
[1115]508          IF ( qc_1d(k) > eps_sb )  THEN
[1012]509!
[1048]510!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
[1115]511!--          (1.0 - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
512             tau_cloud = 1.0 - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
[1012]513!
514!--          Universal function for autoconversion process
515!--          (Seifert and Beheng, 2006):
[1334]516             phi_au    = 600.0 * tau_cloud**0.68_wp * ( 1.0 - tau_cloud**0.68_wp )**3
[1012]517!
518!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
519!--          (Use constant nu_c = 1.0 instead?)
[1334]520             nu_c      = 1.0 !MAX( 0.0_wp, 1580.0 * hyrho(k) * qc(k,j,i) - 0.28_wp )
[1012]521!
522!--          Mean weight of cloud droplets:
[1115]523             xc = hyrho(k) * qc_1d(k) / nc_1d(k)
[1012]524!
[1065]525!--          Parameterized turbulence effects on autoconversion (Seifert,
526!--          Nuijens and Stevens, 2010)
527             IF ( turbulence )  THEN
528!
529!--             Weight averaged radius of cloud droplets:
[1334]530                rc = 0.5 * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
[1065]531
532                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0 + a_3 * nu_c )
533                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0 + b_3 * nu_c )
534                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0 + c_3 * nu_c )
535!
536!--             Mixing length (neglecting distance to ground and stratification)
[1334]537                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
[1065]538!
539!--             Limit dissipation rate according to Seifert, Nuijens and
540!--             Stevens (2010)
[1334]541                epsilon = MIN( 0.06_wp, diss(k,j,i) )
[1065]542!
543!--             Compute Taylor-microscale Reynolds number:
[1334]544                re_lambda = 6.0 / 11.0 * ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *  &
545                            SQRT( 15.0 / kin_vis_air ) * epsilon**( 1.0_wp / 6.0_wp )
[1065]546!
547!--             The factor of 1.0E4 is needed to convert the dissipation rate
548!--             from m2 s-3 to cm2 s-3.
[1334]549                k_au = k_au * ( 1.0_wp +                                       &
550                       epsilon * 1.0E4 * ( re_lambda * 1.0E-3 )**0.25_wp *     &
[1065]551                       ( alpha_cc * EXP( -1.0 * ( ( rc - r_cc ) /              &
552                       sigma_cc )**2 ) + beta_cc ) )
553             ENDIF
554!
[1012]555!--          Autoconversion rate (Seifert and Beheng, 2006):
[1115]556             autocon = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) /    &
[1334]557                       ( nu_c + 1.0 )**2.0_wp * qc_1d(k)**2.0_wp * xc**2.0_wp *   &
558                       ( 1.0 + phi_au / ( 1.0 - tau_cloud )**2.0_wp ) * &
[1115]559                       rho_surface
560             autocon = MIN( autocon, qc_1d(k) / dt_micro )
[1106]561
[1115]562             qr_1d(k) = qr_1d(k) + autocon * dt_micro
563             qc_1d(k) = qc_1d(k) - autocon * dt_micro 
564             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
565
[1005]566          ENDIF
[1000]567
568       ENDDO
569
[1005]570    END SUBROUTINE autoconversion_ij
571
[1106]572
[1005]573    SUBROUTINE accretion_ij( i, j )
574
[1320]575       USE arrays_3d,                                                          &
576           ONLY:  diss, qc_1d, qr_1d
[1115]577
[1320]578       USE cloud_parameters,                                                   &
579           ONLY:  eps_sb, hyrho, k_cr0
580
581       USE control_parameters,                                                 &
582           ONLY:  dt_micro, rho_surface, turbulence
583
584       USE indices,                                                            &
585           ONLY:  nzb, nzb_s_inner, nzt
586
587       USE kinds
588
[1005]589       IMPLICIT NONE
590
[1320]591       INTEGER(iwp) ::  i                 !:
592       INTEGER(iwp) ::  j                 !:
593       INTEGER(iwp) ::  k                 !:
[1005]594
[1320]595       REAL(wp)     ::  accr              !:
596       REAL(wp)     ::  k_cr              !:
597       REAL(wp)     ::  phi_ac            !:
598       REAL(wp)     ::  tau_cloud         !:
599       REAL(wp)     ::  xc                !:
600
[1115]601       DO  k = nzb_s_inner(j,i)+1, nzt
602          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
[1012]603!
[1048]604!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
[1115]605             tau_cloud = 1.0 - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 
[1012]606!
607!--          Universal function for accretion process
[1048]608!--          (Seifert and Beheng, 2001):
[1065]609             phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 ) 
[1334]610             phi_ac = ( phi_ac**2.0_wp )**2.0_wp
[1012]611!
[1065]612!--          Parameterized turbulence effects on autoconversion (Seifert,
613!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
614!--          convert the dissipation (diss) from m2 s-3 to cm2 s-3.
615             IF ( turbulence )  THEN
[1115]616                k_cr = k_cr0 * ( 1.0 + 0.05 *                             &
[1334]617                                 MIN( 600.0_wp, diss(k,j,i) * 1.0E4 )**0.25_wp )
[1065]618             ELSE
619                k_cr = k_cr0                       
620             ENDIF
621!
[1012]622!--          Accretion rate (Seifert and Beheng, 2006):
[1115]623             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac *                 &
[1065]624                    SQRT( rho_surface * hyrho(k) )
[1115]625             accr = MIN( accr, qc_1d(k) / dt_micro )
[1106]626
[1115]627             qr_1d(k) = qr_1d(k) + accr * dt_micro 
628             qc_1d(k) = qc_1d(k) - accr * dt_micro
629
[1005]630          ENDIF
[1106]631
[1005]632       ENDDO
633
[1000]634    END SUBROUTINE accretion_ij
635
[1005]636
637    SUBROUTINE selfcollection_breakup_ij( i, j )
638
[1320]639       USE arrays_3d,                                                          &
640           ONLY:  nr_1d, qr_1d
641
642       USE cloud_parameters,                                                   &
643           ONLY:  dpirho_l, eps_sb, hyrho, k_br, k_rr
644
645       USE control_parameters,                                                 &
646           ONLY:  dt_micro, rho_surface
647
648       USE indices,                                                            &
649           ONLY:  nzb, nzb_s_inner, nzt
650
651       USE kinds
[1005]652   
653       IMPLICIT NONE
654
[1320]655       INTEGER(iwp) ::  i                 !:
656       INTEGER(iwp) ::  j                 !:
657       INTEGER(iwp) ::  k                 !:
[1005]658
[1320]659       REAL(wp)     ::  breakup           !:
660       REAL(wp)     ::  dr                !:
661       REAL(wp)     ::  phi_br            !:
662       REAL(wp)     ::  selfcoll          !:
663
[1115]664       DO  k = nzb_s_inner(j,i)+1, nzt
665          IF ( qr_1d(k) > eps_sb )  THEN
[1012]666!
[1115]667!--          Selfcollection rate (Seifert and Beheng, 2001):
668             selfcoll = k_rr * nr_1d(k) * qr_1d(k) *         &
[1005]669                        SQRT( hyrho(k) * rho_surface )
[1012]670!
[1115]671!--          Weight averaged diameter of rain drops:
[1334]672             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]673!
[1048]674!--          Collisional breakup rate (Seifert, 2008):
[1115]675             IF ( dr >= 0.3E-3 )  THEN
676                phi_br  = k_br * ( dr - 1.1E-3 )
[1005]677                breakup = selfcoll * ( phi_br + 1.0 )
678             ELSE
679                breakup = 0.0
680             ENDIF
[1048]681
[1115]682             selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
683             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
[1106]684
[1005]685          ENDIF         
686       ENDDO
687
688    END SUBROUTINE selfcollection_breakup_ij
689
[1106]690
[1012]691    SUBROUTINE evaporation_rain_ij( i, j )
[1022]692!
693!--    Evaporation of precipitable water. Condensation is neglected for
694!--    precipitable water.
[1012]695
[1320]696       USE arrays_3d,                                                          &
697           ONLY:  hyp, nr_1d, pt_1d, q_1d,  qc_1d, qr_1d
[1048]698
[1320]699       USE cloud_parameters,                                                   &
700           ONLY:  a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
701                  dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r,   &
702                  l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l,      &
703                  t_d_pt, ventilation_effect
704
705       USE constants,                                                          &
706           ONLY:  pi
707
708       USE control_parameters,                                                 &
709           ONLY:  dt_micro
710
711       USE indices,                                                            &
712           ONLY:  nzb, nzb_s_inner, nzt
713
714       USE kinds
715
[1012]716       IMPLICIT NONE
717
[1320]718       INTEGER(iwp) ::  i                 !:
719       INTEGER(iwp) ::  j                 !:
720       INTEGER(iwp) ::  k                 !:
[1012]721
[1320]722       REAL(wp)     ::  alpha             !:
723       REAL(wp)     ::  dr                !:
724       REAL(wp)     ::  e_s               !:
725       REAL(wp)     ::  evap              !:
726       REAL(wp)     ::  evap_nr           !:
727       REAL(wp)     ::  f_vent            !:
728       REAL(wp)     ::  g_evap            !:
729       REAL(wp)     ::  lambda_r          !:
730       REAL(wp)     ::  mu_r              !:
731       REAL(wp)     ::  mu_r_2            !:
732       REAL(wp)     ::  mu_r_5d2          !:
733       REAL(wp)     ::  nr_0              !:
734       REAL(wp)     ::  q_s               !:
735       REAL(wp)     ::  sat               !:
736       REAL(wp)     ::  t_l               !:
737       REAL(wp)     ::  temp              !:
738       REAL(wp)     ::  xr                !:
739
[1115]740       DO  k = nzb_s_inner(j,i)+1, nzt
741          IF ( qr_1d(k) > eps_sb )  THEN
[1012]742!
743!--          Actual liquid water temperature:
[1115]744             t_l = t_d_pt(k) * pt_1d(k)
[1012]745!
746!--          Saturation vapor pressure at t_l:
747             e_s = 610.78 * EXP( 17.269 * ( t_l - 273.16 ) / ( t_l - 35.86 ) )
748!
749!--          Computation of saturation humidity:
750             q_s = 0.622 * e_s / ( hyp(k) - 0.378 * e_s )
751             alpha = 0.622 * l_d_r * l_d_cp / ( t_l * t_l )
[1115]752             q_s = q_s * ( 1.0 + alpha * q_1d(k) ) / ( 1.0 + alpha * q_s )
[1012]753!
[1106]754!--          Supersaturation:
[1334]755             sat = MIN( 0.0_wp, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0 )
[1012]756!
757!--          Actual temperature:
[1115]758             temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
759   
760             g_evap = 1.0 / ( ( l_v / ( r_v * temp ) - 1.0 ) * l_v /   &
761                      ( thermal_conductivity_l * temp ) + r_v * temp / &
762                      ( diff_coeff_l * e_s ) )
[1012]763!
[1115]764!--          Mean weight of rain drops
765             xr = hyrho(k) * qr_1d(k) / nr_1d(k)
[1012]766!
[1115]767!--          Weight averaged diameter of rain drops:
[1334]768             dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]769!
[1049]770!--          Compute ventilation factor and intercept parameter
771!--          (Seifert and Beheng, 2006; Seifert, 2008):
[1048]772             IF ( ventilation_effect )  THEN
[1115]773!
774!--             Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
775!--             Stevens and Seifert, 2008):
776                mu_r = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr - 1.4E-3 ) ) )
777!
778!--             Slope parameter of gamma distribution (Seifert, 2008):
779                lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0 ) *       &
[1334]780                             ( mu_r + 1.0 ) )**( 1.0_wp / 3.0_wp ) / dr
[1115]781
782                mu_r_2   = mu_r + 2.0
783                mu_r_5d2 = mu_r + 2.5 
[1048]784                f_vent = a_vent * gamm( mu_r_2 ) *                            &
[1115]785                         lambda_r**( -mu_r_2 ) +                              &
[1048]786                         b_vent * schmidt_p_1d3 *                             &
787                         SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *    &
[1115]788                         lambda_r**( -mu_r_5d2 ) *                            &
[1048]789                         ( 1.0 - 0.5 * ( b_term / a_term ) *                  &
[1115]790                         ( lambda_r /                                         &
791                         (       c_term + lambda_r ) )**mu_r_5d2 -            &
[1334]792                                 0.125 * ( b_term / a_term )**2.0_wp *        &
[1115]793                         ( lambda_r /                                         &
794                         ( 2.0 * c_term + lambda_r ) )**mu_r_5d2 -            &
[1334]795                                 0.0625 * ( b_term / a_term )**3.0_wp *       &
[1115]796                         ( lambda_r /                                         &
797                         ( 3.0 * c_term + lambda_r ) )**mu_r_5d2 -            &
[1334]798                                 0.0390625 * ( b_term / a_term )**4.0_wp *    &
[1115]799                         ( lambda_r /                                         &
800                         ( 4.0 * c_term + lambda_r ) )**mu_r_5d2 )
801                nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0 ) /                &
802                         gamm( mu_r + 1.0 ) 
[1048]803             ELSE
804                f_vent = 1.0
[1115]805                nr_0   = nr_1d(k) * dr
[1048]806             ENDIF
[1012]807!
[1048]808!--          Evaporation rate of rain water content (Seifert and Beheng, 2006):
[1049]809             evap = 2.0 * pi * nr_0 * g_evap * f_vent * sat /    &
[1048]810                    hyrho(k)
[1106]811
[1115]812             evap    = MAX( evap, -qr_1d(k) / dt_micro )
813             evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
814                            -nr_1d(k) / dt_micro )
815
816             qr_1d(k) = qr_1d(k) + evap * dt_micro
817             nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
[1012]818          ENDIF         
[1106]819
[1012]820       ENDDO
821
822    END SUBROUTINE evaporation_rain_ij
823
[1106]824
[1012]825    SUBROUTINE sedimentation_cloud_ij( i, j )
826
[1320]827       USE arrays_3d,                                                          &
828           ONLY:  ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d
829
830       USE cloud_parameters,                                                   &
831           ONLY:  eps_sb, hyrho, k_st, l_d_cp, prr, pt_d_t, rho_l, sigma_gc
832
833       USE constants,                                                          &
834           ONLY:  pi
835
836       USE control_parameters,                                                 &
837           ONLY:  dt_do2d_xy, dt_micro, intermediate_timestep_count
838
839       USE indices,                                                            &
840           ONLY:  nzb, nzb_s_inner, nzt
841
842       USE kinds
[1012]843       
844       IMPLICIT NONE
845
[1320]846       INTEGER(iwp) ::  i                 !:
847       INTEGER(iwp) ::  j                 !:
848       INTEGER(iwp) ::  k                 !:
[1106]849
[1320]850       REAL(wp)     ::  sed_qc_const      !:
[1115]851
[1320]852
853       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc
854
[1012]855!
856!--    Sedimentation of cloud droplets (Heus et al., 2010):
[1334]857       sed_qc_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))**( 2.0_wp / 3.0_wp ) *   &
[1048]858                     EXP( 5.0 * LOG( sigma_gc )**2 )
[1012]859
[1115]860       sed_qc(nzt+1) = 0.0
[1012]861
[1115]862       DO  k = nzt, nzb_s_inner(j,i)+1, -1
863          IF ( qc_1d(k) > eps_sb )  THEN
[1334]864             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) * &
865                        ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
[1115]866          ELSE
867             sed_qc(k) = 0.0
[1012]868          ENDIF
[1115]869
870          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /     &
871                                      dt_micro + sed_qc(k+1) )
872
873          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /  &
874                                hyrho(k) * dt_micro
875          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 
876                                hyrho(k) * dt_micro
877          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
878                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
879
[1012]880       ENDDO
881
882    END SUBROUTINE sedimentation_cloud_ij
883
[1106]884
[1012]885    SUBROUTINE sedimentation_rain_ij( i, j )
886
[1320]887       USE arrays_3d,                                                          &
888           ONLY:  ddzu, dzu, nr_1d, pt_1d, q_1d, qr_1d
889
890       USE cloud_parameters,                                                   &
891           ONLY:  a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho,        &
892                  limiter_sedimentation, l_d_cp, precipitation_amount, prr,    &
893                  pt_d_t, stp
894
895       USE control_parameters,                                                 &
896           ONLY:  dt_do2d_xy, dt_micro, dt_3d, intermediate_timestep_count,    &
897                  intermediate_timestep_count_max,                             &
898                  precipitation_amount_interval, time_do2d_xy
899
900       USE indices,                                                            &
901           ONLY:  nzb, nzb_s_inner, nzt
902
903       USE kinds
904
905       USE statistics,                                                         &
906           ONLY:  weight_substep
[1012]907       
908       IMPLICIT NONE
909
[1320]910       INTEGER(iwp) ::  i                          !:
911       INTEGER(iwp) ::  j                          !:
912       INTEGER(iwp) ::  k                          !:
913       INTEGER(iwp) ::  k_run                      !:
[1012]914
[1320]915       REAL(wp)     ::  c_run                      !:
916       REAL(wp)     ::  d_max                      !:
917       REAL(wp)     ::  d_mean                     !:
918       REAL(wp)     ::  d_min                      !:
919       REAL(wp)     ::  dr                         !:
920       REAL(wp)     ::  dt_sedi                    !:
921       REAL(wp)     ::  flux                       !:
922       REAL(wp)     ::  lambda_r                   !:
923       REAL(wp)     ::  mu_r                       !:
924       REAL(wp)     ::  z_run                      !:
925
926      REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr      !:
927      REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr      !:
928      REAL(wp), DIMENSION(nzb:nzt+1) ::  d_nr      !:
929      REAL(wp), DIMENSION(nzb:nzt+1) ::  d_qr      !:
930      REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope  !:
931      REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope  !:
932      REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr    !:
933      REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr    !:
934      REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr      !:
935      REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr      !:
936
937
[1065]938!
939!--    Computation of sedimentation flux. Implementation according to Stevens
940!--    and Seifert (2008).
[1048]941       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0
[1012]942!
[1065]943!--    Compute velocities
944       DO  k = nzb_s_inner(j,i)+1, nzt
[1115]945          IF ( qr_1d(k) > eps_sb )  THEN
946!
947!--          Weight averaged diameter of rain drops:
[1334]948             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]949!
950!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
951!--          Stevens and Seifert, 2008):
952             mu_r = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr - 1.4E-3 ) ) )
953!
954!--          Slope parameter of gamma distribution (Seifert, 2008):
955             lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0 ) *          &
[1334]956                        ( mu_r + 1.0 ) )**( 1.0_wp / 3.0_wp ) / dr
[1115]957
[1334]958             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0 +    &
[1115]959                       c_term / lambda_r )**( -1.0 * ( mu_r + 1.0 ) ) ) )
[1334]960             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0 +    &
[1115]961                       c_term / lambda_r )**( -1.0 * ( mu_r + 4.0 ) ) ) )
[1065]962          ELSE
963             w_nr(k) = 0.0
964             w_qr(k) = 0.0
965          ENDIF
966       ENDDO
[1048]967!
[1065]968!--    Adjust boundary values
[1115]969       w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
970       w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
971       w_nr(nzt+1) = 0.0
972       w_qr(nzt+1) = 0.0
[1065]973!
974!--    Compute Courant number
[1115]975       DO  k = nzb_s_inner(j,i)+1, nzt
[1065]976          c_nr(k) = 0.25 * ( w_nr(k-1) + 2.0 * w_nr(k) + w_nr(k+1) ) * &
[1115]977                    dt_micro * ddzu(k)
[1065]978          c_qr(k) = 0.25 * ( w_qr(k-1) + 2.0 * w_qr(k) + w_qr(k+1) ) * &
[1115]979                    dt_micro * ddzu(k)
980       ENDDO     
[1065]981!
982!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
983       IF ( limiter_sedimentation )  THEN
984
[1115]985          DO k = nzb_s_inner(j,i)+1, nzt
986             d_mean = 0.5 * ( qr_1d(k+1) + qr_1d(k-1) )
987             d_min  = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
988             d_max  = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
[1065]989
990             qr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
991                                                     ABS( d_mean ) )
992
[1115]993             d_mean = 0.5 * ( nr_1d(k+1) + nr_1d(k-1) )
994             d_min  = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
995             d_max  = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
[1065]996
997             nr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
998                                                     ABS( d_mean ) )
[1022]999          ENDDO
[1048]1000
[1065]1001       ELSE
[1106]1002
[1065]1003          nr_slope = 0.0
1004          qr_slope = 0.0
[1106]1005
[1065]1006       ENDIF
[1115]1007
1008       sed_nr(nzt+1) = 0.0
1009       sed_qr(nzt+1) = 0.0
[1065]1010!
1011!--    Compute sedimentation flux
[1115]1012       DO  k = nzt, nzb_s_inner(j,i)+1, -1
[1065]1013!
1014!--       Sum up all rain drop number densities which contribute to the flux
1015!--       through k-1/2
1016          flux  = 0.0
1017          z_run = 0.0 ! height above z(k)
1018          k_run = k
1019          c_run = MIN( 1.0, c_nr(k) )
[1115]1020          DO WHILE ( c_run > 0.0  .AND.  k_run <= nzt )
[1065]1021             flux  = flux + hyrho(k_run) *                                    &
[1115]1022                     ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0 - c_run ) *     &
[1065]1023                     0.5 ) * c_run * dzu(k_run)
1024             z_run = z_run + dzu(k_run)
1025             k_run = k_run + 1
1026             c_run = MIN( 1.0, c_nr(k_run) - z_run * ddzu(k_run) )
[1022]1027          ENDDO
1028!
[1065]1029!--       It is not allowed to sediment more rain drop number density than
1030!--       available
1031          flux = MIN( flux,                                                   &
[1115]1032                      hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
[1065]1033
[1115]1034          sed_nr(k) = flux / dt_micro
1035          nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /    &
1036                                 hyrho(k) * dt_micro
[1065]1037!
1038!--       Sum up all rain water content which contributes to the flux
1039!--       through k-1/2
1040          flux  = 0.0
1041          z_run = 0.0 ! height above z(k)
1042          k_run = k
1043          c_run = MIN( 1.0, c_qr(k) )
[1106]1044
[1065]1045          DO WHILE ( c_run > 0.0  .AND.  k_run <= nzt-1 )
[1106]1046
[1065]1047             flux  = flux + hyrho(k_run) *                                    &
[1115]1048                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0 - c_run ) *    &
[1065]1049                     0.5 ) * c_run * dzu(k_run)
1050             z_run = z_run + dzu(k_run)
1051             k_run = k_run + 1
1052             c_run = MIN( 1.0, c_qr(k_run) - z_run * ddzu(k_run) )
[1106]1053
[1065]1054          ENDDO
1055!
1056!--       It is not allowed to sediment more rain water content than available
1057          flux = MIN( flux,                                                   &
[1115]1058                      hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
[1065]1059
[1115]1060          sed_qr(k) = flux / dt_micro
1061
1062          qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
1063                                hyrho(k) * dt_micro
1064          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
1065                                hyrho(k) * dt_micro 
1066          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
1067                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
[1065]1068!
1069!--       Compute the rain rate
1070          prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *                    &
[1115]1071                       weight_substep(intermediate_timestep_count)
[1065]1072       ENDDO
[1115]1073
[1065]1074!
[1048]1075!--    Precipitation amount
1076       IF ( intermediate_timestep_count == intermediate_timestep_count_max    &
1077            .AND.  ( dt_do2d_xy - time_do2d_xy ) <                            &
1078            precipitation_amount_interval )  THEN
[1012]1079
[1048]1080          precipitation_amount(j,i) = precipitation_amount(j,i) +   &
[1115]1081                                      prr(nzb_s_inner(j,i)+1,j,i) *      &
1082                                      hyrho(nzb_s_inner(j,i)+1) * dt_3d
[1048]1083       ENDIF
1084
[1012]1085    END SUBROUTINE sedimentation_rain_ij
1086
[1106]1087
[1012]1088!
1089!-- This function computes the gamma function (Press et al., 1992).
1090!-- The gamma function is needed for the calculation of the evaporation
1091!-- of rain drops.
1092    FUNCTION gamm( xx ) 
[1048]1093       
[1320]1094       USE cloud_parameters,                                                   &
1095           ONLY:  cof, stp
1096
1097       USE kinds
1098
[1012]1099       IMPLICIT NONE
[1106]1100
[1320]1101       INTEGER(iwp) ::  j            !:
1102
1103       REAL(wp)     ::  gamm         !:
1104       REAL(wp)     ::  ser          !:
1105       REAL(wp)     ::  tmp          !:
1106       REAL(wp)     ::  x_gamm       !:
1107       REAL(wp)     ::  xx           !:
1108       REAL(wp)     ::  y_gamm       !:
1109
[1012]1110       x_gamm = xx 
1111       y_gamm = x_gamm 
1112       tmp = x_gamm + 5.5 
1113       tmp = ( x_gamm + 0.5 ) * LOG( tmp ) - tmp 
[1334]1114       ser = 1.000000000190015_wp
[1106]1115
1116       DO  j = 1, 6 
[1012]1117          y_gamm = y_gamm + 1.0 
1118          ser    = ser + cof( j ) / y_gamm 
[1106]1119       ENDDO
1120
[1012]1121!
1122!--    Until this point the algorithm computes the logarithm of the gamma
1123!--    function. Hence, the exponential function is used. 
1124!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
1125       gamm = EXP( tmp ) * stp * ser / x_gamm 
[1106]1126
[1012]1127       RETURN
1128 
1129    END FUNCTION gamm 
1130
1131 END MODULE microphysics_mod
Note: See TracBrowser for help on using the repository browser.