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

Last change on this file since 2416 was 2375, checked in by schwenkel, 7 years ago

improved aerosol initialization for bulk microphysics

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