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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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