source: palm/tags/release-5.0/SOURCE/lpm_droplet_condensation.f90 @ 4418

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

Merge of branch palm4u into trunk

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