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

Last change on this file since 2101 was 2101, checked in by suehring, 7 years ago

last commit documented

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