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

Last change on this file since 3039 was 3039, checked in by schwenkel, 3 years ago

bugfix for lcm with grid stretching

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