source: palm/tags/release-4.0/SOURCE/lpm_droplet_condensation.f90 @ 4097

Last change on this file since 4097 was 1360, checked in by hoffmann, 10 years ago

last commit documented

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