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

Last change on this file since 2716 was 2716, checked in by kanani, 6 years ago

Correction of "Former revisions" section

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