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

Last change on this file since 3046 was 3046, checked in by Giersch, 3 years ago

Remaining error messages revised, comments extended

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