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

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

REAL functions and a lot of REAL constants provided with KIND-attribute,
some small bugfixes

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