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

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

REAL functions and a lot of REAL constants provided with KIND-attribute,
some small bugfixes

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