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

Last change on this file since 1823 was 1823, checked in by hoffmann, 5 years ago

last commit documented

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