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

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

last commit documented

  • Property svn:keywords set to Id
File size: 38.0 KB
Line 
1 MODULE microphysics_mod
2
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!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: microphysics.f90 1335 2014-03-25 12:23:36Z heinze $
27!
28! 1334 2014-03-25 12:21:40Z heinze
29! Bugfix: REAL constants provided with KIND-attribute
30!
31! 1322 2014-03-20 16:38:49Z raasch
32! REAL constants defined as wp-kind
33!
34! 1320 2014-03-20 08:40:49Z raasch
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
40!
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!
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!
49! 1106 2013-03-04 05:31:38Z raasch
50! small changes in code formatting
51!
52! 1092 2013-02-02 11:24:22Z raasch
53! unused variables removed
54! file put under GPL
55!
56! 1065 2012-11-22 17:42:36Z hoffmann
57! Sedimentation process implemented according to Stevens and Seifert (2008).
58! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
59! and Stevens, 2010).
60!
61! 1053 2012-11-13 17:11:03Z hoffmann
62! initial revision
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
71    PUBLIC microphysics_control
72
73    INTERFACE microphysics_control
74       MODULE PROCEDURE microphysics_control
75       MODULE PROCEDURE microphysics_control_ij
76    END INTERFACE microphysics_control
77
78    INTERFACE adjust_cloud
79       MODULE PROCEDURE adjust_cloud
80       MODULE PROCEDURE adjust_cloud_ij
81    END INTERFACE adjust_cloud
82
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
92
93    INTERFACE selfcollection_breakup
94       MODULE PROCEDURE selfcollection_breakup
95       MODULE PROCEDURE selfcollection_breakup_ij
96    END INTERFACE selfcollection_breakup
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
107 
108    INTERFACE sedimentation_rain
109       MODULE PROCEDURE sedimentation_rain
110       MODULE PROCEDURE sedimentation_rain_ij
111    END INTERFACE sedimentation_rain
112
113 CONTAINS
114
115
116!------------------------------------------------------------------------------!
117! Call for all grid points
118!------------------------------------------------------------------------------!
119    SUBROUTINE microphysics_control
120
121       USE arrays_3d
122       USE cloud_parameters
123       USE control_parameters
124       USE grid_variables
125       USE indices
126       USE kinds
127       USE statistics
128
129       IMPLICIT NONE
130
131       INTEGER(iwp) ::  i                 !:
132       INTEGER(iwp) ::  j                 !:
133       INTEGER(iwp) ::  k                 !:
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
148       USE cloud_parameters
149       USE indices
150       USE kinds
151
152       IMPLICIT NONE
153
154       INTEGER(iwp) ::  i                 !:
155       INTEGER(iwp) ::  j                 !:
156       INTEGER(iwp) ::  k                 !:
157
158 
159       DO  i = nxl, nxr
160          DO  j = nys, nyn
161             DO  k = nzb_s_inner(j,i)+1, nzt
162
163             ENDDO
164          ENDDO
165       ENDDO
166
167    END SUBROUTINE adjust_cloud
168
169
170    SUBROUTINE autoconversion
171
172       USE arrays_3d
173       USE cloud_parameters
174       USE control_parameters
175       USE grid_variables
176       USE indices
177       USE kinds
178
179       IMPLICIT NONE
180
181       INTEGER(iwp) ::  i                 !:
182       INTEGER(iwp) ::  j                 !:
183       INTEGER(iwp) ::  k                 !:
184
185       DO  i = nxl, nxr
186          DO  j = nys, nyn
187             DO  k = nzb_s_inner(j,i)+1, nzt
188
189             ENDDO
190          ENDDO
191       ENDDO
192
193    END SUBROUTINE autoconversion
194
195
196    SUBROUTINE accretion
197
198       USE arrays_3d
199       USE cloud_parameters
200       USE control_parameters
201       USE indices
202       USE kinds
203
204       IMPLICIT NONE
205
206       INTEGER(iwp) ::  i                 !:
207       INTEGER(iwp) ::  j                 !:
208       INTEGER(iwp) ::  k                 !:
209
210       DO  i = nxl, nxr
211          DO  j = nys, nyn
212             DO  k = nzb_s_inner(j,i)+1, nzt
213
214             ENDDO
215          ENDDO
216       ENDDO
217
218    END SUBROUTINE accretion
219
220
221    SUBROUTINE selfcollection_breakup
222
223       USE arrays_3d
224       USE cloud_parameters
225       USE control_parameters
226       USE indices
227       USE kinds
228
229       IMPLICIT NONE
230
231       INTEGER(iwp) ::  i                 !:
232       INTEGER(iwp) ::  j                 !:
233       INTEGER(iwp) ::  k                 !:
234
235 
236       DO  i = nxl, nxr
237          DO  j = nys, nyn
238             DO  k = nzb_s_inner(j,i)+1, nzt
239
240             ENDDO
241          ENDDO
242       ENDDO
243
244    END SUBROUTINE selfcollection_breakup
245
246
247    SUBROUTINE evaporation_rain
248
249       USE arrays_3d
250       USE cloud_parameters
251       USE constants
252       USE control_parameters
253       USE indices
254       USE kinds
255
256       IMPLICIT NONE
257
258       INTEGER(iwp) ::  i                 !:
259       INTEGER(iwp) ::  j                 !:
260       INTEGER(iwp) ::  k                 !:
261 
262       DO  i = nxl, nxr
263          DO  j = nys, nyn
264             DO  k = nzb_s_inner(j,i)+1, nzt
265
266             ENDDO
267          ENDDO
268       ENDDO
269
270    END SUBROUTINE evaporation_rain
271
272
273    SUBROUTINE sedimentation_cloud
274
275       USE arrays_3d
276       USE cloud_parameters
277       USE constants
278       USE control_parameters
279       USE indices
280       USE kinds
281
282       IMPLICIT NONE
283
284       INTEGER(iwp) ::  i                 !:
285       INTEGER(iwp) ::  j                 !:
286       INTEGER(iwp) ::  k                 !:
287 
288       DO  i = nxl, nxr
289          DO  j = nys, nyn
290             DO  k = nzb_s_inner(j,i)+1, nzt
291
292             ENDDO
293          ENDDO
294       ENDDO
295
296    END SUBROUTINE sedimentation_cloud
297
298
299    SUBROUTINE sedimentation_rain
300
301       USE arrays_3d
302       USE cloud_parameters
303       USE constants
304       USE control_parameters
305       USE indices
306       USE kinds
307       USE statistics
308
309       IMPLICIT NONE
310
311       INTEGER(iwp) ::  i                 !:
312       INTEGER(iwp) ::  j                 !:
313       INTEGER(iwp) ::  k                 !:
314 
315       DO  i = nxl, nxr
316          DO  j = nys, nyn
317             DO  k = nzb_s_inner(j,i)+1, nzt
318
319             ENDDO
320          ENDDO
321       ENDDO
322
323    END SUBROUTINE sedimentation_rain
324
325
326!------------------------------------------------------------------------------!
327! Call for grid point i,j
328!------------------------------------------------------------------------------!
329
330    SUBROUTINE microphysics_control_ij( i, j )
331
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
335
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
352       IMPLICIT NONE
353
354       INTEGER(iwp) ::  i                 !:
355       INTEGER(iwp) ::  j                 !:
356       INTEGER(iwp) ::  k                 !:
357
358       REAL(wp)     ::  t_surface         !:
359
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)
366          t_surface = pt_surface * ( surface_pressure / 1000.0 )**0.286_wp
367          DO  k = nzb, nzt+1
368             hyp(k)    = surface_pressure * 100.0 * &
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
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
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
418       USE arrays_3d,                                                          &
419           ONLY:  qr, nr
420
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
429       IMPLICIT NONE
430
431       INTEGER(iwp) ::  i                 !:
432       INTEGER(iwp) ::  j                 !:
433       INTEGER(iwp) ::  k                 !:
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
441
442          IF ( qr(k,j,i) <= eps_sb )  THEN
443             qr(k,j,i) = 0.0
444             nr(k,j,i) = 0.0
445          ELSE
446!
447!--          Adjust number of raindrops to avoid nonlinear effects in
448!--          sedimentation and evaporation of rain drops due to too small or
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
454             ENDIF
455
456          ENDIF
457
458       ENDDO
459
460    END SUBROUTINE adjust_cloud_ij
461
462
463    SUBROUTINE autoconversion_ij( i, j )
464
465       USE arrays_3d,                                                          &
466           ONLY:  diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d
467
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
483       IMPLICIT NONE
484
485       INTEGER(iwp) ::  i                 !:
486       INTEGER(iwp) ::  j                 !:
487       INTEGER(iwp) ::  k                 !:
488
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                !:
503
504       k_au = k_cc / ( 20.0 * x0 )
505
506       DO  k = nzb_s_inner(j,i)+1, nzt
507
508          IF ( qc_1d(k) > eps_sb )  THEN
509!
510!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
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) )
513!
514!--          Universal function for autoconversion process
515!--          (Seifert and Beheng, 2006):
516             phi_au    = 600.0 * tau_cloud**0.68_wp * ( 1.0 - tau_cloud**0.68_wp )**3
517!
518!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
519!--          (Use constant nu_c = 1.0 instead?)
520             nu_c      = 1.0 !MAX( 0.0_wp, 1580.0 * hyrho(k) * qc(k,j,i) - 0.28_wp )
521!
522!--          Mean weight of cloud droplets:
523             xc = hyrho(k) * qc_1d(k) / nc_1d(k)
524!
525!--          Parameterized turbulence effects on autoconversion (Seifert,
526!--          Nuijens and Stevens, 2010)
527             IF ( turbulence )  THEN
528!
529!--             Weight averaged radius of cloud droplets:
530                rc = 0.5 * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
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)
537                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
538!
539!--             Limit dissipation rate according to Seifert, Nuijens and
540!--             Stevens (2010)
541                epsilon = MIN( 0.06_wp, diss(k,j,i) )
542!
543!--             Compute Taylor-microscale Reynolds number:
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 )
546!
547!--             The factor of 1.0E4 is needed to convert the dissipation rate
548!--             from m2 s-3 to cm2 s-3.
549                k_au = k_au * ( 1.0_wp +                                       &
550                       epsilon * 1.0E4 * ( re_lambda * 1.0E-3 )**0.25_wp *     &
551                       ( alpha_cc * EXP( -1.0 * ( ( rc - r_cc ) /              &
552                       sigma_cc )**2 ) + beta_cc ) )
553             ENDIF
554!
555!--          Autoconversion rate (Seifert and Beheng, 2006):
556             autocon = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) /    &
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 ) * &
559                       rho_surface
560             autocon = MIN( autocon, qc_1d(k) / dt_micro )
561
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
566          ENDIF
567
568       ENDDO
569
570    END SUBROUTINE autoconversion_ij
571
572
573    SUBROUTINE accretion_ij( i, j )
574
575       USE arrays_3d,                                                          &
576           ONLY:  diss, qc_1d, qr_1d
577
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
589       IMPLICIT NONE
590
591       INTEGER(iwp) ::  i                 !:
592       INTEGER(iwp) ::  j                 !:
593       INTEGER(iwp) ::  k                 !:
594
595       REAL(wp)     ::  accr              !:
596       REAL(wp)     ::  k_cr              !:
597       REAL(wp)     ::  phi_ac            !:
598       REAL(wp)     ::  tau_cloud         !:
599       REAL(wp)     ::  xc                !:
600
601       DO  k = nzb_s_inner(j,i)+1, nzt
602          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
603!
604!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
605             tau_cloud = 1.0 - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 
606!
607!--          Universal function for accretion process
608!--          (Seifert and Beheng, 2001):
609             phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 ) 
610             phi_ac = ( phi_ac**2.0_wp )**2.0_wp
611!
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
616                k_cr = k_cr0 * ( 1.0 + 0.05 *                             &
617                                 MIN( 600.0_wp, diss(k,j,i) * 1.0E4 )**0.25_wp )
618             ELSE
619                k_cr = k_cr0                       
620             ENDIF
621!
622!--          Accretion rate (Seifert and Beheng, 2006):
623             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac *                 &
624                    SQRT( rho_surface * hyrho(k) )
625             accr = MIN( accr, qc_1d(k) / dt_micro )
626
627             qr_1d(k) = qr_1d(k) + accr * dt_micro
628             qc_1d(k) = qc_1d(k) - accr * dt_micro
629
630          ENDIF
631
632       ENDDO
633
634    END SUBROUTINE accretion_ij
635
636
637    SUBROUTINE selfcollection_breakup_ij( i, j )
638
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
652   
653       IMPLICIT NONE
654
655       INTEGER(iwp) ::  i                 !:
656       INTEGER(iwp) ::  j                 !:
657       INTEGER(iwp) ::  k                 !:
658
659       REAL(wp)     ::  breakup           !:
660       REAL(wp)     ::  dr                !:
661       REAL(wp)     ::  phi_br            !:
662       REAL(wp)     ::  selfcoll          !:
663
664       DO  k = nzb_s_inner(j,i)+1, nzt
665          IF ( qr_1d(k) > eps_sb )  THEN
666!
667!--          Selfcollection rate (Seifert and Beheng, 2001):
668             selfcoll = k_rr * nr_1d(k) * qr_1d(k) *         &
669                        SQRT( hyrho(k) * rho_surface )
670!
671!--          Weight averaged diameter of rain drops:
672             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
673!
674!--          Collisional breakup rate (Seifert, 2008):
675             IF ( dr >= 0.3E-3 )  THEN
676                phi_br  = k_br * ( dr - 1.1E-3 )
677                breakup = selfcoll * ( phi_br + 1.0 )
678             ELSE
679                breakup = 0.0
680             ENDIF
681
682             selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
683             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
684
685          ENDIF         
686       ENDDO
687
688    END SUBROUTINE selfcollection_breakup_ij
689
690
691    SUBROUTINE evaporation_rain_ij( i, j )
692!
693!--    Evaporation of precipitable water. Condensation is neglected for
694!--    precipitable water.
695
696       USE arrays_3d,                                                          &
697           ONLY:  hyp, nr_1d, pt_1d, q_1d,  qc_1d, qr_1d
698
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
716       IMPLICIT NONE
717
718       INTEGER(iwp) ::  i                 !:
719       INTEGER(iwp) ::  j                 !:
720       INTEGER(iwp) ::  k                 !:
721
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
740       DO  k = nzb_s_inner(j,i)+1, nzt
741          IF ( qr_1d(k) > eps_sb )  THEN
742!
743!--          Actual liquid water temperature:
744             t_l = t_d_pt(k) * pt_1d(k)
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 )
752             q_s = q_s * ( 1.0 + alpha * q_1d(k) ) / ( 1.0 + alpha * q_s )
753!
754!--          Supersaturation:
755             sat = MIN( 0.0_wp, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0 )
756!
757!--          Actual temperature:
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 ) )
763!
764!--          Mean weight of rain drops
765             xr = hyrho(k) * qr_1d(k) / nr_1d(k)
766!
767!--          Weight averaged diameter of rain drops:
768             dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
769!
770!--          Compute ventilation factor and intercept parameter
771!--          (Seifert and Beheng, 2006; Seifert, 2008):
772             IF ( ventilation_effect )  THEN
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 ) *       &
780                             ( mu_r + 1.0 ) )**( 1.0_wp / 3.0_wp ) / dr
781
782                mu_r_2   = mu_r + 2.0
783                mu_r_5d2 = mu_r + 2.5 
784                f_vent = a_vent * gamm( mu_r_2 ) *                            &
785                         lambda_r**( -mu_r_2 ) +                              &
786                         b_vent * schmidt_p_1d3 *                             &
787                         SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *    &
788                         lambda_r**( -mu_r_5d2 ) *                            &
789                         ( 1.0 - 0.5 * ( b_term / a_term ) *                  &
790                         ( lambda_r /                                         &
791                         (       c_term + lambda_r ) )**mu_r_5d2 -            &
792                                 0.125 * ( b_term / a_term )**2.0_wp *        &
793                         ( lambda_r /                                         &
794                         ( 2.0 * c_term + lambda_r ) )**mu_r_5d2 -            &
795                                 0.0625 * ( b_term / a_term )**3.0_wp *       &
796                         ( lambda_r /                                         &
797                         ( 3.0 * c_term + lambda_r ) )**mu_r_5d2 -            &
798                                 0.0390625 * ( b_term / a_term )**4.0_wp *    &
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 ) 
803             ELSE
804                f_vent = 1.0
805                nr_0   = nr_1d(k) * dr
806             ENDIF
807!
808!--          Evaporation rate of rain water content (Seifert and Beheng, 2006):
809             evap = 2.0 * pi * nr_0 * g_evap * f_vent * sat /    &
810                    hyrho(k)
811
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
818          ENDIF         
819
820       ENDDO
821
822    END SUBROUTINE evaporation_rain_ij
823
824
825    SUBROUTINE sedimentation_cloud_ij( i, j )
826
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
843       
844       IMPLICIT NONE
845
846       INTEGER(iwp) ::  i                 !:
847       INTEGER(iwp) ::  j                 !:
848       INTEGER(iwp) ::  k                 !:
849
850       REAL(wp)     ::  sed_qc_const      !:
851
852
853       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc
854
855!
856!--    Sedimentation of cloud droplets (Heus et al., 2010):
857       sed_qc_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))**( 2.0_wp / 3.0_wp ) *   &
858                     EXP( 5.0 * LOG( sigma_gc )**2 )
859
860       sed_qc(nzt+1) = 0.0
861
862       DO  k = nzt, nzb_s_inner(j,i)+1, -1
863          IF ( qc_1d(k) > eps_sb )  THEN
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 )
866          ELSE
867             sed_qc(k) = 0.0
868          ENDIF
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
880       ENDDO
881
882    END SUBROUTINE sedimentation_cloud_ij
883
884
885    SUBROUTINE sedimentation_rain_ij( i, j )
886
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
907       
908       IMPLICIT NONE
909
910       INTEGER(iwp) ::  i                          !:
911       INTEGER(iwp) ::  j                          !:
912       INTEGER(iwp) ::  k                          !:
913       INTEGER(iwp) ::  k_run                      !:
914
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
938!
939!--    Computation of sedimentation flux. Implementation according to Stevens
940!--    and Seifert (2008).
941       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0
942!
943!--    Compute velocities
944       DO  k = nzb_s_inner(j,i)+1, nzt
945          IF ( qr_1d(k) > eps_sb )  THEN
946!
947!--          Weight averaged diameter of rain drops:
948             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * 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 * ( 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 ) *          &
956                        ( mu_r + 1.0 ) )**( 1.0_wp / 3.0_wp ) / dr
957
958             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0 +    &
959                       c_term / lambda_r )**( -1.0 * ( mu_r + 1.0 ) ) ) )
960             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0 +    &
961                       c_term / lambda_r )**( -1.0 * ( mu_r + 4.0 ) ) ) )
962          ELSE
963             w_nr(k) = 0.0
964             w_qr(k) = 0.0
965          ENDIF
966       ENDDO
967!
968!--    Adjust boundary values
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
973!
974!--    Compute Courant number
975       DO  k = nzb_s_inner(j,i)+1, nzt
976          c_nr(k) = 0.25 * ( w_nr(k-1) + 2.0 * w_nr(k) + w_nr(k+1) ) * &
977                    dt_micro * ddzu(k)
978          c_qr(k) = 0.25 * ( w_qr(k-1) + 2.0 * w_qr(k) + w_qr(k+1) ) * &
979                    dt_micro * ddzu(k)
980       ENDDO     
981!
982!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
983       IF ( limiter_sedimentation )  THEN
984
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)
989
990             qr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
991                                                     ABS( d_mean ) )
992
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)
996
997             nr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
998                                                     ABS( d_mean ) )
999          ENDDO
1000
1001       ELSE
1002
1003          nr_slope = 0.0
1004          qr_slope = 0.0
1005
1006       ENDIF
1007
1008       sed_nr(nzt+1) = 0.0
1009       sed_qr(nzt+1) = 0.0
1010!
1011!--    Compute sedimentation flux
1012       DO  k = nzt, nzb_s_inner(j,i)+1, -1
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) )
1020          DO WHILE ( c_run > 0.0  .AND.  k_run <= nzt )
1021             flux  = flux + hyrho(k_run) *                                    &
1022                     ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0 - c_run ) *     &
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) )
1027          ENDDO
1028!
1029!--       It is not allowed to sediment more rain drop number density than
1030!--       available
1031          flux = MIN( flux,                                                   &
1032                      hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
1033
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
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) )
1044
1045          DO WHILE ( c_run > 0.0  .AND.  k_run <= nzt-1 )
1046
1047             flux  = flux + hyrho(k_run) *                                    &
1048                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0 - c_run ) *    &
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) )
1053
1054          ENDDO
1055!
1056!--       It is not allowed to sediment more rain water content than available
1057          flux = MIN( flux,                                                   &
1058                      hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
1059
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
1068!
1069!--       Compute the rain rate
1070          prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *                    &
1071                       weight_substep(intermediate_timestep_count)
1072       ENDDO
1073
1074!
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
1079
1080          precipitation_amount(j,i) = precipitation_amount(j,i) +   &
1081                                      prr(nzb_s_inner(j,i)+1,j,i) *      &
1082                                      hyrho(nzb_s_inner(j,i)+1) * dt_3d
1083       ENDIF
1084
1085    END SUBROUTINE sedimentation_rain_ij
1086
1087
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 ) 
1093       
1094       USE cloud_parameters,                                                   &
1095           ONLY:  cof, stp
1096
1097       USE kinds
1098
1099       IMPLICIT NONE
1100
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
1110       x_gamm = xx
1111       y_gamm = x_gamm
1112       tmp = x_gamm + 5.5 
1113       tmp = ( x_gamm + 0.5 ) * LOG( tmp ) - tmp
1114       ser = 1.000000000190015_wp
1115
1116       DO  j = 1, 6 
1117          y_gamm = y_gamm + 1.0 
1118          ser    = ser + cof( j ) / y_gamm
1119       ENDDO
1120
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
1126
1127       RETURN
1128 
1129    END FUNCTION gamm
1130
1131 END MODULE microphysics_mod
Note: See TracBrowser for help on using the repository browser.