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

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

last commit documented

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