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

Last change on this file since 2316 was 2312, checked in by hoffmann, 7 years ago

various improvements of the LCM

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