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

Last change on this file since 3236 was 3049, checked in by Giersch, 6 years ago

Revision history corrected

  • Property svn:keywords set to Id
File size: 17.8 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 3049 2018-05-29 13:52:36Z raasch $
27! Error messages revised
28!
29! 3045 2018-05-28 07:55:41Z Giersch
30! Error messages revised
31!
32! 3039 2018-05-24 13:13:11Z schwenkel
33! bugfix for lcm with grid stretching
34!
35! 2718 2018-01-02 08:49:38Z maronga
36! Corrected "Former revisions" section
37!
38! 2696 2017-12-14 17:12:51Z kanani
39! Change in file header (GPL part)
40!
41! 2608 2017-11-13 14:04:26Z schwenkel
42! Calculation of magnus equation in external module (diagnostic_quantities_mod).
43!
44! 2375 2017-08-29 14:10:28Z schwenkel
45! Changed ONLY-dependencies
46!
47! 2312 2017-07-14 20:26:51Z hoffmann
48! Rosenbrock scheme improved. Gas-kinetic effect added.
49!
50! 2000 2016-08-20 18:09:15Z knoop
51! Forced header and separation lines into 80 columns
52!
53! 1890 2016-04-22 08:52:11Z hoffmann
54! Some improvements of the Rosenbrock method. If the Rosenbrock method needs more
55! than 40 iterations to find a sufficient time setp, the model is not aborted.
56! This might lead to small erros. But they can be assumend as negligible, since
57! the maximum timesetp of the Rosenbrock method after 40 iterations will be
58! smaller than 10^-16 s.
59!
60! 1871 2016-04-15 11:46:09Z hoffmann
61! Initialization of aerosols added.
62!
63! 1849 2016-04-08 11:33:18Z hoffmann
64! Interpolation of supersaturation has been removed because it is not in
65! accordance with the release/depletion of latent heat/water vapor in
66! interaction_droplets_ptq.
67! Calculation of particle Reynolds number has been corrected.
68! eps_ros added from modules.
69!
70! 1831 2016-04-07 13:15:51Z hoffmann
71! curvature_solution_effects moved to particle_attributes
72!
73! 1822 2016-04-07 07:49:42Z hoffmann
74! Unused variables removed.
75!
76! 1682 2015-10-07 23:56:08Z knoop
77! Code annotations made doxygen readable
78!
79! 1359 2014-04-11 17:15:14Z hoffmann
80! New particle structure integrated.
81! Kind definition added to all floating point numbers.
82!
83! 1346 2014-03-27 13:18:20Z heinze
84! Bugfix: REAL constants provided with KIND-attribute especially in call of
85! intrinsic function like MAX, MIN, SIGN
86!
87! 1322 2014-03-20 16:38:49Z raasch
88! REAL constants defined as wp-kind
89!
90! 1320 2014-03-20 08:40:49Z raasch
91! ONLY-attribute added to USE-statements,
92! kind-parameters added to all INTEGER and REAL declaration statements,
93! kinds are defined in new module kinds,
94! comment fields (!:) to be used for variable explanations added to
95! all variable declaration statements
96!
97! 1318 2014-03-17 13:35:16Z raasch
98! module interfaces removed
99!
100! 1092 2013-02-02 11:24:22Z raasch
101! unused variables removed
102!
103! 1071 2012-11-29 16:54:55Z franke
104! Ventilation effect for evaporation of large droplets included
105! Check for unreasonable results included in calculation of Rosenbrock method
106! since physically unlikely results were observed and for the same
107! reason the first internal time step in Rosenbrock method should be < 1.0E02 in
108! case of evaporation
109! Unnecessary calculation of ql_int removed
110! Unnecessary calculations in Rosenbrock method (d2rdt2, drdt_m, dt_ros_last)
111! removed
112! Bugfix: factor in calculation of surface tension changed from 0.00155 to
113! 0.000155
114!
115! 1036 2012-10-22 13:43:42Z raasch
116! code put under GPL (PALM 3.9)
117!
118! 849 2012-03-15 10:35:09Z raasch
119! initial revision (former part of advec_particles)
120!
121!
122! Description:
123! ------------
124!> Calculates change in droplet radius by condensation/evaporation, using
125!> either an analytic formula or by numerically integrating the radius growth
126!> equation including curvature and solution effects using Rosenbrocks method
127!> (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
128!> The analytical formula and growth equation follow those given in
129!> Rogers and Yau (A short course in cloud physics, 3rd edition, p. 102/103).
130!------------------------------------------------------------------------------!
131 SUBROUTINE lpm_droplet_condensation (ip,jp,kp)
132
133
134    USE arrays_3d,                                                             &
135        ONLY:  dzw, hyp, pt, q, ql_c, ql_v
136
137    USE cloud_parameters,                                                      &
138        ONLY:  l_d_rv, l_v, molecular_weight_of_solute,                        &
139               molecular_weight_of_water, rho_l, rho_s, r_v, vanthoff
140
141    USE constants,                                                             &
142        ONLY:  pi
143
144    USE control_parameters,                                                    &
145        ONLY:  dt_3d, dz, message_string, molecular_viscosity, rho_surface
146
147    USE cpulog,                                                                &
148        ONLY:  cpu_log, log_point_s
149
150    USE diagnostic_quantities_mod,                                             &
151        ONLY:  magnus
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, hall_kernel, number_of_particles,   &
163               particles, radius_classes, use_kernel_tables, wang_kernel
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) * ( hyp(kp) / 100000.0_wp )**0.286_wp
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) + 0.622_wp )
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.