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

Last change on this file since 1092 was 1092, checked in by raasch, 9 years ago

unused variables remove from several routines

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