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

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