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

Last change on this file since 1071 was 1071, checked in by franke, 9 years ago

ventilation effect for cloud droplets included, calculation of cloud droplet collisions now uses formulation by Wang, bugfixes in Rosenbrock method and in calculation of collision efficiencies

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