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

Last change on this file since 1323 was 1323, checked in by raasch, 7 years ago

last commit documented

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