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

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

changes in LPM and bulk cloud microphysics

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