source: palm/trunk/SOURCE/lpm_droplet_condensation.f90 @ 3655

Last change on this file since 3655 was 3655, checked in by knoop, 5 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

  • Property svn:keywords set to Id
File size: 17.7 KB
Line 
1!> @file lpm_droplet_condensation.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! 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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_droplet_condensation.f90 3655 2019-01-07 16:51:22Z knoop $
27! Modularization of all bulk cloud physics code components
28!
29! 3241 2018-09-12 15:02:00Z raasch
30! unused variables removed
31!
32! 3049 2018-05-29 13:52:36Z Giersch
33! Error messages revised
34!
35! 3045 2018-05-28 07:55:41Z Giersch
36! Error messages revised
37!
38! 3039 2018-05-24 13:13:11Z schwenkel
39! bugfix for lcm with grid stretching
40!
41! 2718 2018-01-02 08:49:38Z maronga
42! Corrected "Former revisions" section
43!
44! 2696 2017-12-14 17:12:51Z kanani
45! Change in file header (GPL part)
46!
47! 2608 2017-11-13 14:04:26Z schwenkel
48! Calculation of magnus equation in external module (diagnostic_quantities_mod).
49!
50! 2375 2017-08-29 14:10:28Z schwenkel
51! Changed ONLY-dependencies
52!
53! 2312 2017-07-14 20:26:51Z hoffmann
54! Rosenbrock scheme improved. Gas-kinetic effect added.
55!
56! 2000 2016-08-20 18:09:15Z knoop
57! Forced header and separation lines into 80 columns
58!
59! 1890 2016-04-22 08:52:11Z hoffmann
60! Some improvements of the Rosenbrock method. If the Rosenbrock method needs more
61! than 40 iterations to find a sufficient time setp, the model is not aborted.
62! This might lead to small erros. But they can be assumend as negligible, since
63! the maximum timesetp of the Rosenbrock method after 40 iterations will be
64! smaller than 10^-16 s.
65!
66! 1871 2016-04-15 11:46:09Z hoffmann
67! Initialization of aerosols added.
68!
69! 1849 2016-04-08 11:33:18Z hoffmann
70! Interpolation of supersaturation has been removed because it is not in
71! accordance with the release/depletion of latent heat/water vapor in
72! interaction_droplets_ptq.
73! Calculation of particle Reynolds number has been corrected.
74! eps_ros added from modules.
75!
76! 1831 2016-04-07 13:15:51Z hoffmann
77! curvature_solution_effects moved to particle_attributes
78!
79! 1822 2016-04-07 07:49:42Z hoffmann
80! Unused variables removed.
81!
82! 1682 2015-10-07 23:56:08Z knoop
83! Code annotations made doxygen readable
84!
85! 1359 2014-04-11 17:15:14Z hoffmann
86! New particle structure integrated.
87! Kind definition added to all floating point numbers.
88!
89! 1346 2014-03-27 13:18:20Z heinze
90! Bugfix: REAL constants provided with KIND-attribute especially in call of
91! intrinsic function like MAX, MIN, SIGN
92!
93! 1322 2014-03-20 16:38:49Z raasch
94! REAL constants defined as wp-kind
95!
96! 1320 2014-03-20 08:40:49Z raasch
97! ONLY-attribute added to USE-statements,
98! kind-parameters added to all INTEGER and REAL declaration statements,
99! kinds are defined in new module kinds,
100! comment fields (!:) to be used for variable explanations added to
101! all variable declaration statements
102!
103! 1318 2014-03-17 13:35:16Z raasch
104! module interfaces removed
105!
106! 1092 2013-02-02 11:24:22Z raasch
107! unused variables removed
108!
109! 1071 2012-11-29 16:54:55Z franke
110! Ventilation effect for evaporation of large droplets included
111! Check for unreasonable results included in calculation of Rosenbrock method
112! since physically unlikely results were observed and for the same
113! reason the first internal time step in Rosenbrock method should be < 1.0E02 in
114! case of evaporation
115! Unnecessary calculation of ql_int removed
116! Unnecessary calculations in Rosenbrock method (d2rdt2, drdt_m, dt_ros_last)
117! removed
118! Bugfix: factor in calculation of surface tension changed from 0.00155 to
119! 0.000155
120!
121! 1036 2012-10-22 13:43:42Z raasch
122! code put under GPL (PALM 3.9)
123!
124! 849 2012-03-15 10:35:09Z raasch
125! initial revision (former part of advec_particles)
126!
127!
128! Description:
129! ------------
130!> Calculates change in droplet radius by condensation/evaporation, using
131!> either an analytic formula or by numerically integrating the radius growth
132!> equation including curvature and solution effects using Rosenbrocks method
133!> (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
134!> The analytical formula and growth equation follow those given in
135!> Rogers and Yau (A short course in cloud physics, 3rd edition, p. 102/103).
136!------------------------------------------------------------------------------!
137 SUBROUTINE lpm_droplet_condensation (ip,jp,kp)
138
139
140    USE arrays_3d,                                                             &
141        ONLY:  dzw, hyp, pt, q, ql_c, ql_v, exner
142
143    USE basic_constants_and_equations_mod,                                     &
144        ONLY:  l_v, molecular_weight_of_solute, molecular_weight_of_water,     &
145               magnus, pi, rho_l, rho_s, rd_d_rv, r_v, vanthoff
146
147    USE control_parameters,                                                    &
148        ONLY:  dt_3d, message_string, molecular_viscosity, rho_surface
149
150    USE cpulog,                                                                &
151        ONLY:  cpu_log, log_point_s
152
153    USE grid_variables,                                                        &
154        ONLY:  dx, dy
155
156    USE lpm_collision_kernels_mod,                                             &
157        ONLY:  rclass_lbound, rclass_ubound
158
159    USE kinds
160
161    USE particle_attributes,                                                   &
162        ONLY:  curvature_solution_effects, number_of_particles,   &
163               particles, radius_classes, use_kernel_tables
164
165    IMPLICIT NONE
166
167    INTEGER(iwp) :: ip                         !<
168    INTEGER(iwp) :: jp                         !<
169    INTEGER(iwp) :: kp                         !<
170    INTEGER(iwp) :: n                          !<
171
172    REAL(wp) ::  afactor                       !< curvature effects
173    REAL(wp) ::  arg                           !<
174    REAL(wp) ::  bfactor                       !< solute effects
175    REAL(wp) ::  ddenom                        !<
176    REAL(wp) ::  delta_r                       !<
177    REAL(wp) ::  diameter                      !< diameter of cloud droplets
178    REAL(wp) ::  diff_coeff                    !< diffusivity for water vapor
179    REAL(wp) ::  drdt                          !<
180    REAL(wp) ::  dt_ros                        !<
181    REAL(wp) ::  dt_ros_sum                    !<
182    REAL(wp) ::  d2rdtdr                       !<
183    REAL(wp) ::  e_a                           !< current vapor pressure
184    REAL(wp) ::  e_s                           !< current saturation vapor pressure
185    REAL(wp) ::  error                         !< local truncation error in Rosenbrock
186    REAL(wp) ::  k1                            !<
187    REAL(wp) ::  k2                            !<
188    REAL(wp) ::  r_err                         !< First order estimate of Rosenbrock radius
189    REAL(wp) ::  r_ros                         !< Rosenbrock radius
190    REAL(wp) ::  r_ros_ini                     !< initial Rosenbrock radius
191    REAL(wp) ::  r0                            !< gas-kinetic lengthscale
192    REAL(wp) ::  sigma                         !< surface tension of water
193    REAL(wp) ::  thermal_conductivity          !< thermal conductivity for water
194    REAL(wp) ::  t_int                         !< temperature
195    REAL(wp) ::  w_s                           !< terminal velocity of droplets
196    REAL(wp) ::  re_p                          !< particle Reynolds number
197!
198!-- Parameters for Rosenbrock method (see Verwer et al., 1999)
199    REAL(wp), PARAMETER :: prec = 1.0E-3_wp     !< precision of Rosenbrock solution
200    REAL(wp), PARAMETER :: q_increase = 1.5_wp  !< increase factor in timestep
201    REAL(wp), PARAMETER :: q_decrease = 0.9_wp  !< decrease factor in timestep
202    REAL(wp), PARAMETER :: gamma = 0.292893218814_wp !< = 1.0 - 1.0 / SQRT(2.0)
203!
204!-- Parameters for terminal velocity
205    REAL(wp), PARAMETER ::  a_rog = 9.65_wp      !< parameter for fall velocity
206    REAL(wp), PARAMETER ::  b_rog = 10.43_wp     !< parameter for fall velocity
207    REAL(wp), PARAMETER ::  c_rog = 0.6_wp       !< parameter for fall velocity
208    REAL(wp), PARAMETER ::  k_cap_rog = 4.0_wp   !< parameter for fall velocity
209    REAL(wp), PARAMETER ::  k_low_rog = 12.0_wp  !< parameter for fall velocity
210    REAL(wp), PARAMETER ::  d0_rog = 0.745_wp    !< separation diameter
211
212    REAL(wp), DIMENSION(number_of_particles) ::  ventilation_effect     !<
213    REAL(wp), DIMENSION(number_of_particles) ::  new_r                  !<
214
215    CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'start' )
216
217!
218!-- Absolute temperature
219    t_int = pt(kp,jp,ip) * exner(kp)
220!
221!-- Saturation vapor pressure (Eq. 10 in Bolton, 1980)
222    e_s = magnus( t_int )
223!
224!-- Current vapor pressure
225    e_a = q(kp,jp,ip) * hyp(kp) / ( q(kp,jp,ip) + rd_d_rv )
226!
227!-- Thermal conductivity for water (from Rogers and Yau, Table 7.1)
228    thermal_conductivity = 7.94048E-05_wp * t_int + 0.00227011_wp
229!
230!-- Moldecular diffusivity of water vapor in air (Hall und Pruppacher, 1976)
231    diff_coeff           = 0.211E-4_wp * ( t_int / 273.15_wp )**1.94_wp * &
232                           ( 101325.0_wp / hyp(kp) )
233!
234!-- Lengthscale for gas-kinetic effects (from Mordy, 1959, p. 23):
235    r0 = diff_coeff / 0.036_wp * SQRT( 2.0_wp * pi / ( r_v * t_int ) )
236!
237!-- Calculate effects of heat conductivity and diffusion of water vapor on the
238!-- diffusional growth process (usually known as 1.0 / (F_k + F_d) )
239    ddenom  = 1.0_wp / ( rho_l * r_v * t_int / ( e_s * diff_coeff ) +          &
240                         ( l_v / ( r_v * t_int ) - 1.0_wp ) * rho_l *          &
241                         l_v / ( thermal_conductivity * t_int )                &
242                       )
243    new_r = 0.0_wp
244!
245!-- Determine ventilation effect on evaporation of large drops
246    DO  n = 1, number_of_particles
247
248       IF ( particles(n)%radius >= 4.0E-5_wp  .AND.  e_a / e_s < 1.0_wp )  THEN
249!
250!--       Terminal velocity is computed for vertical direction (Rogers et al.,
251!--       1993, J. Appl. Meteorol.)
252          diameter = particles(n)%radius * 2000.0_wp !diameter in mm
253          IF ( diameter <= d0_rog )  THEN
254             w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
255          ELSE
256             w_s = a_rog - b_rog * EXP( -c_rog * diameter )
257          ENDIF
258!
259!--       Calculate droplet's Reynolds number
260          re_p = 2.0_wp * particles(n)%radius * w_s / molecular_viscosity
261!
262!--       Ventilation coefficient (Rogers and Yau, 1989):
263          IF ( re_p > 2.5_wp )  THEN
264             ventilation_effect(n) = 0.78_wp + 0.28_wp * SQRT( re_p )
265          ELSE
266             ventilation_effect(n) = 1.0_wp + 0.09_wp * re_p
267          ENDIF
268       ELSE
269!
270!--       For small droplets or in supersaturated environments, the ventilation
271!--       effect does not play a role
272          ventilation_effect(n) = 1.0_wp
273       ENDIF
274    ENDDO
275
276    IF( .NOT. curvature_solution_effects ) then
277!
278!--    Use analytic model for diffusional growth including gas-kinetic
279!--    effects (Mordy, 1959) but without the impact of aerosols.
280       DO  n = 1, number_of_particles
281          arg      = ( particles(n)%radius + r0 )**2 + 2.0_wp * dt_3d * ddenom * &
282                                                       ventilation_effect(n) *   &
283                                                       ( e_a / e_s - 1.0_wp )
284          arg      = MAX( arg, ( 0.01E-6 + r0 )**2 )
285          new_r(n) = SQRT( arg ) - r0
286       ENDDO
287
288    ELSE
289!
290!--    Integrate the diffusional growth including gas-kinetic (Mordy, 1959),
291!--    as well as curvature and solute effects (e.g., Köhler, 1936).
292!
293!--    Curvature effect (afactor) with surface tension (sigma) by Straka (2009)
294       sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
295!
296!--    Solute effect (afactor)
297       afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
298
299       DO  n = 1, number_of_particles
300!
301!--       Solute effect (bfactor)
302          bfactor = vanthoff * rho_s * particles(n)%aux1**3 *                    &
303                    molecular_weight_of_water / ( rho_l * molecular_weight_of_solute )
304
305          dt_ros     = particles(n)%aux2  ! use previously stored Rosenbrock timestep
306          dt_ros_sum = 0.0_wp
307
308          r_ros     = particles(n)%radius  ! initialize Rosenbrock particle radius
309          r_ros_ini = r_ros
310!
311!--       Integrate growth equation using a 2nd-order Rosenbrock method
312!--       (see Verwer et al., 1999, Eq. (3.2)). The Rosenbrock method adjusts
313!--       its with internal timestep to minimize the local truncation error.
314          DO WHILE ( dt_ros_sum < dt_3d )
315
316             dt_ros = MIN( dt_ros, dt_3d - dt_ros_sum )
317
318             DO
319
320                drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -    &
321                                                          afactor / r_ros +    &
322                                                          bfactor / r_ros**3   &
323                                                        ) / ( r_ros + r0 )
324
325                d2rdtdr = -ddenom * ventilation_effect(n) * (                  &
326                                                (e_a / e_s - 1.0) * r_ros**4 - &
327                                                afactor * r0 * r_ros**2 -      &
328                                                2.0 * afactor * r_ros**3 +     &
329                                                3.0 * bfactor * r0 +           &
330                                                4.0 * bfactor * r_ros          &
331                                                            )                  &
332                          / ( r_ros**4 * ( r_ros + r0 )**2 )
333
334                k1    = drdt / ( 1.0 - gamma * dt_ros * d2rdtdr )
335
336                r_ros = MAX(r_ros_ini + k1 * dt_ros, particles(n)%aux1)
337                r_err = r_ros
338
339                drdt  = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -   &
340                                                           afactor / r_ros +   &
341                                                           bfactor / r_ros**3  &
342                                                         ) / ( r_ros + r0 )
343
344                k2 = ( drdt - dt_ros * 2.0 * gamma * d2rdtdr * k1 ) / &
345                     ( 1.0 - dt_ros * gamma * d2rdtdr )
346
347                r_ros = MAX(r_ros_ini + dt_ros * ( 1.5 * k1 + 0.5 * k2), particles(n)%aux1)
348   !
349   !--          Check error of the solution, and reduce dt_ros if necessary.
350                error = ABS(r_err - r_ros) / r_ros
351                IF ( error .GT. prec )  THEN
352                   dt_ros = SQRT( q_decrease * prec / error ) * dt_ros
353                   r_ros  = r_ros_ini
354                ELSE
355                   dt_ros_sum = dt_ros_sum + dt_ros
356                   dt_ros     = q_increase * dt_ros
357                   r_ros_ini  = r_ros
358                   EXIT
359                ENDIF
360
361             END DO
362
363          END DO !Rosenbrock loop
364!
365!--       Store new particle radius
366          new_r(n) = r_ros
367!
368!--       Store internal time step value for next PALM step
369          particles(n)%aux2 = dt_ros
370
371       ENDDO !Particle loop
372
373    ENDIF
374
375    DO  n = 1, number_of_particles
376!
377!--    Sum up the change in liquid water for the respective grid
378!--    box for the computation of the release/depletion of water vapor
379!--    and heat.
380       ql_c(kp,jp,ip) = ql_c(kp,jp,ip) + particles(n)%weight_factor *          &
381                                   rho_l * 1.33333333_wp * pi *                &
382                                   ( new_r(n)**3 - particles(n)%radius**3 ) /  &
383                                   ( rho_surface * dx * dy * dzw(kp) )
384!
385!--    Check if the increase in liqid water is not too big. If this is the case,
386!--    the model timestep might be too long.
387       IF ( ql_c(kp,jp,ip) > 100.0_wp )  THEN
388          WRITE( message_string, * ) 'k=',kp,' j=',jp,' i=',ip,                &
389                       ' ql_c=',ql_c(kp,jp,ip), '&part(',n,')%wf=',            &
390                       particles(n)%weight_factor,' delta_r=',delta_r
391          CALL message( 'lpm_droplet_condensation', 'PA0143', 2, 2, -1, 6, 1 )
392       ENDIF
393!
394!--    Check if the change in the droplet radius is not too big. If this is the
395!--    case, the model timestep might be too long.
396       delta_r = new_r(n) - particles(n)%radius
397       IF ( delta_r < 0.0_wp  .AND. new_r(n) < 0.0_wp )  THEN
398          WRITE( message_string, * ) '#1 k=',kp,' j=',jp,' i=',ip,             &
399                       ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int,               &
400                       '&delta_r=',delta_r,                                    &
401                       ' particle_radius=',particles(n)%radius
402          CALL message( 'lpm_droplet_condensation', 'PA0144', 2, 2, -1, 6, 1 )
403       ENDIF
404!
405!--    Sum up the total volume of liquid water (needed below for
406!--    re-calculating the weighting factors)
407       ql_v(kp,jp,ip) = ql_v(kp,jp,ip) + particles(n)%weight_factor * new_r(n)**3
408!
409!--    Determine radius class of the particle needed for collision
410       IF ( use_kernel_tables )  THEN
411          particles(n)%class = ( LOG( new_r(n) ) - rclass_lbound ) /           &
412                               ( rclass_ubound - rclass_lbound ) *             &
413                               radius_classes
414          particles(n)%class = MIN( particles(n)%class, radius_classes )
415          particles(n)%class = MAX( particles(n)%class, 1 )
416       ENDIF
417 !
418 !--   Store new radius to particle features
419       particles(n)%radius = new_r(n)
420
421    ENDDO
422
423    CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'stop' )
424
425
426 END SUBROUTINE lpm_droplet_condensation
Note: See TracBrowser for help on using the repository browser.