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

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

last commit documented

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