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

Last change on this file since 1267 was 1093, checked in by raasch, 12 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 17.7 KB
RevLine 
[849]1 SUBROUTINE lpm_droplet_condensation
2
[1036]3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
[849]20! Current revisions:
21! ------------------
[1093]22!
[1072]23!
24! Former revisions:
25! -----------------
26! $Id: lpm_droplet_condensation.f90 1093 2013-02-02 12:58:49Z heinze $
27!
[1093]28! 1092 2013-02-02 11:24:22Z raasch
29! unused variables removed
30!
[1072]31! 1071 2012-11-29 16:54:55Z franke
[1071]32! Ventilation effect for evaporation of large droplets included
33! Check for unreasonable results included in calculation of Rosenbrock method
34! since physically unlikely results were observed and for the same
35! reason the first internal time step in Rosenbrock method should be < 1.0E02 in
36! case of evaporation
37! Unnecessary calculation of ql_int removed
38! Unnecessary calculations in Rosenbrock method (d2rdt2, drdt_m, dt_ros_last)
39! removed
40! Bugfix: factor in calculation of surface tension changed from 0.00155 to
41! 0.000155
[849]42!
[1037]43! 1036 2012-10-22 13:43:42Z raasch
44! code put under GPL (PALM 3.9)
45!
[850]46! 849 2012-03-15 10:35:09Z raasch
47! initial revision (former part of advec_particles)
[849]48!
[850]49!
[849]50! Description:
51! ------------
52! Calculates change in droplet radius by condensation/evaporation, using
53! either an analytic formula or by numerically integrating the radius growth
54! equation including curvature and solution effects using Rosenbrocks method
55! (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
56! The analytical formula and growth equation follow those given in
57! Rogers and Yau (A short course in cloud physics, 3rd edition, p. 102/103).
58!------------------------------------------------------------------------------!
59
60    USE arrays_3d
61    USE cloud_parameters
62    USE constants
63    USE control_parameters
64    USE cpulog
65    USE grid_variables
66    USE interfaces
67    USE lpm_collision_kernels_mod
68    USE particle_attributes
69
70    IMPLICIT NONE
71
[1071]72    INTEGER ::  i, internal_timestep_count, j, jtry, k, n, ros_count
[849]73
74    INTEGER, PARAMETER ::  maxtry = 40
75
[1071]76    LOGICAL ::  repeat
77
[849]78    REAL ::  aa, afactor, arg, bb, cc, dd, ddenom, delta_r, drdt, drdt_ini,    &
[1071]79             dt_ros, dt_ros_next, dt_ros_sum, dt_ros_sum_ini, d2rdtdr, errmax, &
80             err_ros, g1, g2, g3, g4, e_a, e_s, gg, new_r, p_int, pt_int,      &
[1092]81             pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, r_ros, r_ros_ini,    &
82             sigma, t_int, x, y, re_p
[849]83
84!
85!-- Parameters for Rosenbrock method
86    REAL, PARAMETER ::  a21 = 2.0, a31 = 48.0/25.0, a32 = 6.0/25.0,        &
[1071]87                        b1 = 19.0/9.0, b2 = 0.5, b3 = 25.0/108.0,          &
88                        b4 = 125.0/108.0, c21 = -8.0, c31 = 372.0/25.0,    &
89                        c32 = 12.0/5.0, c41 = -112.0/125.0,                &
90                        c42 = -54.0/125.0, c43 = -2.0/5.0,                 &
[849]91                        errcon = 0.1296, e1 = 17.0/54.0, e2 = 7.0/36.0,    &
92                        e3 = 0.0, e4 = 125.0/108.0, gam = 0.5, grow = 1.5, &
93                        pgrow = -0.25, pshrnk = -1.0/3.0, shrnk = 0.5
94
95
96    CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'start' )
97
98    DO  n = 1, number_of_particles
99!
100!--    Interpolate temperature and humidity.
101!--    First determine left, south, and bottom index of the arrays.
102       i = particles(n)%x * ddx
103       j = particles(n)%y * ddy
104       k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
105           + offset_ocean_nzt                     ! only exact if equidistant
106
107       x  = particles(n)%x - i * dx
108       y  = particles(n)%y - j * dy
109       aa = x**2          + y**2
110       bb = ( dx - x )**2 + y**2
111       cc = x**2          + ( dy - y )**2
112       dd = ( dx - x )**2 + ( dy - y )**2
113       gg = aa + bb + cc + dd
114
115       pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) * pt(k,j,i+1)   &
116                  + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
117                  ) / ( 3.0 * gg )
118
119       pt_int_u = ( ( gg-aa ) * pt(k+1,j,i)   + ( gg-bb ) * pt(k+1,j,i+1)   &
120                  + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
121                  ) / ( 3.0 * gg )
122
123       pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
124                           ( pt_int_u - pt_int_l )
125
126       q_int_l = ( ( gg - aa ) * q(k,j,i)   + ( gg - bb ) * q(k,j,i+1)   &
127                 + ( gg - cc ) * q(k,j+1,i) + ( gg - dd ) * q(k,j+1,i+1) &
128                 ) / ( 3.0 * gg )
129
130       q_int_u = ( ( gg-aa ) * q(k+1,j,i)   + ( gg-bb ) * q(k+1,j,i+1)   &
131                 + ( gg-cc ) * q(k+1,j+1,i) + ( gg-dd ) * q(k+1,j+1,i+1) &
132                 ) / ( 3.0 * gg )
133
134       q_int = q_int_l + ( particles(n)%z - zu(k) ) / dz * &
135                           ( q_int_u - q_int_l )
136
137!
138!--    Calculate real temperature and saturation vapor pressure
139       p_int = hyp(k) + ( particles(n)%z - zu(k) ) / dz * ( hyp(k+1)-hyp(k) )
140       t_int = pt_int * ( p_int / 100000.0 )**0.286
141
142       e_s = 611.0 * EXP( l_d_rv * ( 3.6609E-3 - 1.0 / t_int ) )
143
144!
145!--    Current vapor pressure
146       e_a = q_int * p_int / ( 0.378 * q_int + 0.622 )
147
148!
149!--    Thermal conductivity for water (from Rogers and Yau, Table 7.1),
150!--    diffusivity for water vapor (after Hall und Pruppacher, 1976)
151       thermal_conductivity_l = 7.94048E-05 * t_int + 0.00227011 
152       diff_coeff_l           = 0.211E-4 * ( t_int / 273.15 )**1.94 * &
153                                ( 101325.0 / p_int)
154!
155!--    Change in radius by condensation/evaporation
[1071]156       IF ( particles(n)%radius >= 4.0E-5  .AND.  e_a/e_s < 1.0 )  THEN
[849]157!
[1071]158!--       Approximation for large radii, where curvature and solution effects
159!--       can be neglected but ventilation effect has to be included in case of
160!--       evaporation.
161!--       First calculate the droplet's Reynolds number
162          re_p = 2.0 * particles(n)%radius * ABS( particles(n)%speed_z ) / &
163                 molecular_viscosity
164!
165!--       Ventilation coefficient after Rogers and Yau, 1989
166          IF ( re_p > 2.5 )  THEN
167             afactor = 0.78 + 0.28 * SQRT( re_p )
168          ELSE
169             afactor = 1.0 + 0.09 * re_p
170          ENDIF
171
172          arg = particles(n)%radius**2 + 2.0 * dt_3d * afactor *              &
173                           ( e_a / e_s - 1.0 ) /                              &
174                           ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / &
175                             thermal_conductivity_l +                         &
176                             rho_l * r_v * t_int / diff_coeff_l / e_s )
177
178          new_r = SQRT( arg )
179
180       ELSEIF ( particles(n)%radius >= 1.0E-6  .OR.  &
181                .NOT. curvature_solution_effects )  THEN
182!
183!--       Approximation for larger radii in case that curvature and solution
184!--       effects are neglected and ventilation effects does not play a role
[849]185          arg = particles(n)%radius**2 + 2.0 * dt_3d *                        &
186                           ( e_a / e_s - 1.0 ) /                              &
187                           ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / &
188                             thermal_conductivity_l +                         &
189                             rho_l * r_v * t_int / diff_coeff_l / e_s )
190          IF ( arg < 1.0E-16 )  THEN
191             new_r = 1.0E-8
192          ELSE
193             new_r = SQRT( arg )
194          ENDIF
195       ENDIF
196
197       IF ( curvature_solution_effects  .AND.                              &
198            ( ( particles(n)%radius < 1.0E-6 ) .OR. ( new_r < 1.0E-6 ) ) ) &
199       THEN
200!
201!--       Curvature and solutions effects are included in growth equation.
202!--       Change in Radius is calculated with 4th-order Rosenbrock method
203!--       for stiff o.d.e's with monitoring local truncation error to adjust
204!--       stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
205!--       For larger radii the simple analytic method (see ELSE) gives
206!--       almost the same results.
[1071]207
208          ros_count = 0
209          repeat = .TRUE.
[849]210!
[1071]211!--       Carry out the Rosenbrock algorithm. In case of unreasonable results
212!--       the switch "repeat" will be set true and the algorithm will be carried
213!--       out again with the internal time step set to its initial (small) value.
214!--       Unreasonable results may occur if the external conditions, especially the
215!--       supersaturation, has significantly changed compared to the last PALM
216!--       timestep.
217          DO WHILE ( repeat )
[849]218
[1071]219             repeat = .FALSE.
220!
221!--          Surface tension after (Straka, 2009)
222             sigma = 0.0761 - 0.000155 * ( t_int - 273.15 )
[849]223
[1071]224             r_ros = particles(n)%radius
225             dt_ros_sum  = 0.0      ! internal integrated time (s)
226             internal_timestep_count = 0
[849]227
[1071]228             ddenom  = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) +  &
229                               ( l_v / ( r_v * t_int ) - 1.0 ) *               &
230                               rho_l * l_v / ( thermal_conductivity_l * t_int )&
231                             )
232
233             afactor = 2.0 * sigma / ( rho_l * r_v * t_int )
234
[849]235!
[1071]236!--          Take internal time step values from the end of last PALM time step
237             dt_ros_next = particles(n)%rvar1
238
[849]239!
[1071]240!--          Internal time step should not be > 1.0E-2 in case of evaporation
241!--          because larger values may lead to secondary solutions which are
242!--          physically unlikely
243             IF ( dt_ros_next > 1.0E-2  .AND.  e_a/e_s < 1.0 )  THEN
244                dt_ros_next = 1.0E-3
245             ENDIF
[849]246!
[1071]247!--          If calculation of Rosenbrock method is repeated due to unreasonalble
248!--          results during previous try the initial internal time step has to be
249!--          reduced
250             IF ( ros_count > 1 )  THEN
251                dt_ros_next = dt_ros_next - ( 0.2 * dt_ros_next )
252             ELSEIF ( ros_count > 5 )  THEN
[849]253!
[1071]254!--             Prevent creation of infinite loop
255                message_string = 'ros_count > 5 in Rosenbrock method'
256                CALL message( 'lpm_droplet_condensation', 'PA0018', 2, 2, &
257                               0, 6, 0 )
258             ENDIF
259
[849]260!
[1071]261!--          Internal time step must not be larger than PALM time step
262             dt_ros = MIN( dt_ros_next, dt_3d )
263!
264!--          Integrate growth equation in time unless PALM time step is reached
265             DO WHILE ( dt_ros_sum < dt_3d )
[849]266
[1071]267                internal_timestep_count = internal_timestep_count + 1
[849]268
269!
[1071]270!--             Derivative at starting value
271                drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + &
272                                          bfactor / r_ros**3 )
273                drdt_ini       = drdt
274                dt_ros_sum_ini = dt_ros_sum
275                r_ros_ini      = r_ros
[849]276
277!
[1071]278!--             Calculate radial derivative of dr/dt
279                d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + &
280                                     2.0 * afactor / r_ros**3 -     &
281                                     4.0 * bfactor / r_ros**5 )
[849]282!
[1071]283!--             Adjust stepsize unless required accuracy is reached
284                DO  jtry = 1, maxtry+1
[849]285
[1071]286                   IF ( jtry == maxtry+1 )  THEN
287                      message_string = 'maxtry > 40 in Rosenbrock method'
288                      CALL message( 'lpm_droplet_condensation', 'PA0347', 2, 2, &
289                                    0, 6, 0 )
290                   ENDIF
[849]291
[1071]292                   aa    = 1.0 / ( gam * dt_ros ) - d2rdtdr
293                   g1    = drdt_ini / aa
294                   r_ros = r_ros_ini + a21 * g1
295                   drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
296                                              afactor / r_ros + &
297                                              bfactor / r_ros**3 )
[849]298
[1071]299                   g2    = ( drdt + c21 * g1 / dt_ros )&
300                           / aa
301                   r_ros = r_ros_ini + a31 * g1 + a32 * g2
302                   drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
303                                              afactor / r_ros + &
304                                              bfactor / r_ros**3 )
[849]305
[1071]306                   g3    = ( drdt +  &
307                             ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
308                   g4    = ( drdt +  &
309                             ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
310                   r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
[849]311
[1071]312                   dt_ros_sum = dt_ros_sum_ini + dt_ros
[849]313
[1071]314                   IF ( dt_ros_sum == dt_ros_sum_ini )  THEN
315                      message_string = 'zero stepsize in Rosenbrock method'
316                      CALL message( 'lpm_droplet_condensation', 'PA0348', 2, 2, &
317                                    0, 6, 0 )
318                   ENDIF
[849]319!
[1071]320!--                Calculate error
321                   err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4
322                   errmax  = 0.0
323                   errmax  = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
[849]324!
[1071]325!--                Leave loop if accuracy is sufficient, otherwise try again
326!--                with a reduced stepsize
327                   IF ( errmax <= 1.0 )  THEN
328                      EXIT
329                   ELSE
330                      dt_ros = SIGN( MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), &
331                                          shrnk * ABS( dt_ros ) ), dt_ros )
332                   ENDIF
333
334                ENDDO  ! loop for stepsize adjustment
335
336!
337!--             Calculate next internal time step
338                IF ( errmax > errcon )  THEN
339                   dt_ros_next = 0.9 * dt_ros * errmax**pgrow
[849]340                ELSE
[1071]341                   dt_ros_next = grow * dt_ros
[849]342                ENDIF
343
[1071]344!
345!--             Estimated time step is reduced if the PALM time step is exceeded
346                IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d )  THEN
347                   dt_ros = dt_3d - dt_ros_sum
348                ELSE
349                   dt_ros = dt_ros_next
350                ENDIF
[849]351
[1071]352             ENDDO
[849]353!
[1071]354!--          Store internal time step value for next PALM step
355             particles(n)%rvar1 = dt_ros_next
[849]356
[1071]357             new_r = r_ros
[849]358!
[1071]359!--          Radius should not fall below 1E-8 because Rosenbrock method may
360!--          lead to errors otherwise
361             new_r = MAX( new_r, 1.0E-8 )
362!
363!--          Check if calculated droplet radius change is reasonable since in
364!--          case of droplet evaporation the Rosenbrock method may lead to
365!--          secondary solutions which are physically unlikely.
366!--          Due to the solution effect the droplets may grow for relative
367!--          humidities below 100%, but change of radius should not be too large.
368!--          In case of unreasonable droplet growth the Rosenbrock method is
369!--          recalculated using a smaller initial time step.
370!--          Limiting values are tested for droplets down to 1.0E-7
371             IF ( new_r - particles(n)%radius >= 3.0E-7  .AND.  &
372                  e_a/e_s < 0.97 )  THEN
373                ros_count = ros_count + 1
374                repeat = .TRUE.
[849]375             ENDIF
376
[1071]377          ENDDO    ! Rosenbrock method
[849]378
379       ENDIF
380
381       delta_r = new_r - particles(n)%radius
382
383!
384!--    Sum up the change in volume of liquid water for the respective grid
385!--    volume (this is needed later in lpm_calc_liquid_water_content for
386!--    calculating the release of latent heat)
387       i = ( particles(n)%x + 0.5 * dx ) * ddx
388       j = ( particles(n)%y + 0.5 * dy ) * ddy
389       k = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
390           ! only exact if equidistant
391
392       ql_c(k,j,i) = ql_c(k,j,i) + particles(n)%weight_factor *            &
393                                   rho_l * 1.33333333 * pi *               &
394                                   ( new_r**3 - particles(n)%radius**3 ) / &
395                                   ( rho_surface * dx * dy * dz )
396       IF ( ql_c(k,j,i) > 100.0 )  THEN
397          WRITE( message_string, * ) 'k=',k,' j=',j,' i=',i,      &
398                       ' ql_c=',ql_c(k,j,i), ' &part(',n,')%wf=', &
399                       particles(n)%weight_factor,' delta_r=',delta_r
400          CALL message( 'lpm_droplet_condensation', 'PA0143', 2, 2, -1, 6, 1 )
401       ENDIF
402
403!
404!--    Change the droplet radius
405       IF ( ( new_r - particles(n)%radius ) < 0.0  .AND.  new_r < 0.0 ) &
406       THEN
407          WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i,   &
408                       ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int,  &
409                       ' &delta_r=',delta_r,                      &
410                       ' particle_radius=',particles(n)%radius
411          CALL message( 'lpm_droplet_condensation', 'PA0144', 2, 2, -1, 6, 1 )
412       ENDIF
413
414!
415!--    Sum up the total volume of liquid water (needed below for
416!--    re-calculating the weighting factors)
417       ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * new_r**3
418
419       particles(n)%radius = new_r
420
421!
422!--    Determine radius class of the particle needed for collision
423       IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  use_kernel_tables ) &
424       THEN
425          particles(n)%class = ( LOG( new_r ) - rclass_lbound ) /  &
426                               ( rclass_ubound - rclass_lbound ) * &
427                               radius_classes
428          particles(n)%class = MIN( particles(n)%class, radius_classes )
429          particles(n)%class = MAX( particles(n)%class, 1 )
430       ENDIF
431
432    ENDDO
433
434    CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'stop' )
435
436
437 END SUBROUTINE lpm_droplet_condensation
Note: See TracBrowser for help on using the repository browser.