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

Last change on this file since 1106 was 1106, checked in by raasch, 12 years ago

New:
---

Porting of FFT-solver for serial runs to GPU using CUDA FFT,
preprocessor lines in transpose routines rearranged, so that routines can also
be used in serial (non-parallel) mode,
transpositions also carried out in serial mode, routines fftx, fftxp replaced
by calls of fft_x, fft_x replaced by fft_x_1d in the 1D-decomposition routines
(Makefile, Makefile_check, fft_xy, poisfft, poisfft_hybrid, transpose, new: cuda_fft_interfaces)

--stdin argument for mpiexec on lckyuh, -y and -Y settings output to header (mrun)

Changed:


Module array_kind renamed precision_kind
(check_open, data_output_3d, fft_xy, modules, user_data_output_3d)

some format changes for coupled atmosphere-ocean runs (header)
small changes in code formatting (microphysics, prognostic_equations)

Errors:


bugfix: default value (0) assigned to coupling_start_time (modules)
bugfix: initial time for preruns of coupled runs is output as -coupling_start_time (data_output_profiles)

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