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

Last change on this file since 1832 was 1832, checked in by hoffmann, 6 years ago

last commit documented

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