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

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

Bugfix: REAL constants provided with KIND-attribute especially in call of intrinsic function like MAX, MIN, SIGN

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