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

Last change on this file since 1279 was 1242, checked in by heinze, 11 years ago

Last commmit documented

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