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

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

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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