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

Last change on this file since 1318 was 1318, checked in by raasch, 8 years ago

former files/routines cpu_log and cpu_statistics combined to one module,
which also includes the former data module cpulog from the modules-file,
module interfaces removed

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