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

Last change on this file since 3272 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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