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

Last change on this file since 1574 was 1360, checked in by hoffmann, 11 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 23.8 KB
RevLine 
[1359]1 SUBROUTINE lpm_droplet_condensation (ip,jp,kp)
[849]2
[1036]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!
[1310]17! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]18!--------------------------------------------------------------------------------!
19!
[849]20! Current revisions:
21! ------------------
[1360]22!
23!
[1321]24! Former revisions:
25! -----------------
26! $Id: lpm_droplet_condensation.f90 1360 2014-04-11 17:20:32Z suehring $
27!
[1360]28! 1359 2014-04-11 17:15:14Z hoffmann
29! New particle structure integrated.
30! Kind definition added to all floating point numbers.
31!
[1347]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!
[1323]36! 1322 2014-03-20 16:38:49Z raasch
37! REAL constants defined as wp-kind
38!
[1321]39! 1320 2014-03-20 08:40:49Z raasch
[1320]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
[1072]45!
[1319]46! 1318 2014-03-17 13:35:16Z raasch
47! module interfaces removed
48!
[1093]49! 1092 2013-02-02 11:24:22Z raasch
50! unused variables removed
51!
[1072]52! 1071 2012-11-29 16:54:55Z franke
[1071]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
[849]63!
[1037]64! 1036 2012-10-22 13:43:42Z raasch
65! code put under GPL (PALM 3.9)
66!
[850]67! 849 2012-03-15 10:35:09Z raasch
68! initial revision (former part of advec_particles)
[849]69!
[850]70!
[849]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
[1320]81    USE arrays_3d,                                                             &
82        ONLY:  hyp, pt, q,  ql_c, ql_v, zu
[849]83
[1320]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
[849]87
[1320]88    USE constants,                                                             &
89        ONLY:  pi
[849]90
[1320]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
[849]96
[1320]97    USE grid_variables,                                                        &
98        ONLY:  dx, ddx, dy, ddy
[1071]99
[1320]100    USE lpm_collision_kernels_mod,                                             &
101        ONLY:  rclass_lbound, rclass_ubound
[849]102
[1320]103    USE kinds
104
105    USE particle_attributes,                                                   &
[1359]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
[1320]109
110
111    IMPLICIT NONE
112
113    INTEGER(iwp) :: i                          !:
[1359]114    INTEGER(iwp) :: ip                         !:
[1320]115    INTEGER(iwp) :: internal_timestep_count    !:
116    INTEGER(iwp) :: j                          !:
[1359]117    INTEGER(iwp) :: jp                         !:
[1320]118    INTEGER(iwp) :: jtry                       !:
119    INTEGER(iwp) :: k                          !:
[1359]120    INTEGER(iwp) :: kp                         !:
[1320]121    INTEGER(iwp) :: n                          !:
[1359]122    INTEGER(iwp) :: nb                         !:
[1320]123    INTEGER(iwp) :: ros_count                  !:
124 
[1359]125    INTEGER(iwp), PARAMETER ::  maxtry = 40      !:
[1320]126
[1359]127    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !:
128    INTEGER(iwp), DIMENSION(0:7) ::  start_index !:
[1320]129
[1359]130    LOGICAL ::  repeat                           !:
131
132    LOGICAL, DIMENSION(number_of_particles) ::  flag_1 !:
133
[1320]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 
[849]169!-- Parameters for Rosenbrock method
[1359]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             !:
[849]193
[1359]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                     !:
[849]205
[1320]206
[849]207    CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'start' )
208
[1359]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)
[849]223!
224!--    Interpolate temperature and humidity.
[1359]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
[849]232
[1359]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 )
[849]236
[1359]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 )
[849]240
[1359]241          pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz *                &
242                              ( pt_int_u - pt_int_l )
[849]243
[1359]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 )
[849]247
[1359]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 )
[849]251
[1359]252          q_int = q_int_l + ( zv(n) - zu(k) ) / dz *                           &
253                            ( q_int_u - q_int_l )
[849]254
255!
[1359]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
[849]260
[1359]261          e_s(n) = 611.0_wp * EXP( l_d_rv * ( 3.6609E-3_wp - 1.0_wp /          &
262                                   t_int(n) ) )
[849]263
264!
[1359]265!--       Current vapor pressure
266          e_a(n) = q_int * p_int(n) / ( 0.378_wp * q_int + 0.622_wp )
[849]267
[1359]268       ENDDO
269    ENDDO
270
271    new_r = 0.0_wp
272    flag_1 = .false.
273
274    DO  n = 1, number_of_particles
[849]275!
276!--    Change in radius by condensation/evaporation
[1359]277       IF ( particles(n)%radius >= 4.0E-5_wp  .AND.                            &
278            e_a(n)/e_s(n) < 1.0_wp )  THEN
[849]279!
[1359]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
[1071]286!
[1359]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 )
[1071]290          ELSE
[1359]291             afactor_v(n) = 1.0_wp + 0.09_wp * re_p
[1071]292          ENDIF
[1359]293          flag_1(n) = .TRUE.
294       ELSEIF ( particles(n)%radius >= 1.0E-6_wp  .OR.                         &
[1071]295                .NOT. curvature_solution_effects )  THEN
296!
[1359]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.
[849]301       ENDIF
[1359]302    ENDDO
[849]303
[1359]304    DO  n = 1, number_of_particles
[849]305!
[1359]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!
[849]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.
[1071]335
336          ros_count = 0
337          repeat = .TRUE.
[849]338!
[1071]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.
[1359]342!--       Unreasonable results may occur if the external conditions, especially
343!--       the supersaturation, has significantly changed compared to the last
344!--       PALM timestep.
[1071]345          DO WHILE ( repeat )
[849]346
[1071]347             repeat = .FALSE.
348!
[1359]349!--          Surface tension (Straka, 2009):
350             sigma = 0.0761_wp - 0.000155_wp * ( t_int(n) - 273.15_wp )
[849]351
[1071]352             r_ros = particles(n)%radius
[1359]353             dt_ros_sum  = 0.0_wp      ! internal integrated time (s)
[1071]354             internal_timestep_count = 0
[849]355
[1359]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                                )
[1071]362
[1359]363             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int(n) )
[1071]364
[849]365!
[1071]366!--          Take internal time step values from the end of last PALM time step
367             dt_ros_next = particles(n)%rvar1
368
[849]369!
[1071]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
[1359]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
[1071]375             ENDIF
[849]376!
[1071]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
[1359]381                dt_ros_next = dt_ros_next - ( 0.2_wp * dt_ros_next )
[1071]382             ELSEIF ( ros_count > 5 )  THEN
[849]383!
[1071]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
[849]390!
[1071]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 )
[849]396
[1071]397                internal_timestep_count = internal_timestep_count + 1
[849]398
399!
[1071]400!--             Derivative at starting value
[1359]401                drdt = ddenom / r_ros * ( e_a(n) / e_s(n) - 1.0_wp - afactor / &
402                                          r_ros + bfactor / r_ros**3 )
[1071]403                drdt_ini       = drdt
404                dt_ros_sum_ini = dt_ros_sum
405                r_ros_ini      = r_ros
[849]406
407!
[1071]408!--             Calculate radial derivative of dr/dt
[1359]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 )
[849]412!
[1071]413!--             Adjust stepsize unless required accuracy is reached
414                DO  jtry = 1, maxtry+1
[849]415
[1071]416                   IF ( jtry == maxtry+1 )  THEN
417                      message_string = 'maxtry > 40 in Rosenbrock method'
[1359]418                      CALL message( 'lpm_droplet_condensation', 'PA0347', 2,   &
419                                    2, 0, 6, 0 )
[1071]420                   ENDIF
[849]421
[1359]422                   aa    = 1.0_wp / ( gam * dt_ros ) - d2rdtdr
[1071]423                   g1    = drdt_ini / aa
424                   r_ros = r_ros_ini + a21 * g1
[1359]425                   drdt  = ddenom / r_ros * ( e_a(n) / e_s(n) - 1.0_wp -       &
426                                              afactor / r_ros +                &
[1071]427                                              bfactor / r_ros**3 )
[849]428
[1071]429                   g2    = ( drdt + c21 * g1 / dt_ros )&
430                           / aa
431                   r_ros = r_ros_ini + a31 * g1 + a32 * g2
[1359]432                   drdt  = ddenom / r_ros * ( e_a(n) / e_s(n) - 1.0_wp -       &
433                                              afactor / r_ros +                &
[1071]434                                              bfactor / r_ros**3 )
[849]435
[1071]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
[849]441
[1071]442                   dt_ros_sum = dt_ros_sum_ini + dt_ros
[849]443
[1071]444                   IF ( dt_ros_sum == dt_ros_sum_ini )  THEN
445                      message_string = 'zero stepsize in Rosenbrock method'
[1359]446                      CALL message( 'lpm_droplet_condensation', 'PA0348', 2,   &
447                                    2, 0, 6, 0 )
[1071]448                   ENDIF
[849]449!
[1071]450!--                Calculate error
[1359]451                   err_ros = e1 * g1 + e2 * g2 + e3 * g3 + e4 * g4
452                   errmax  = 0.0_wp
[1071]453                   errmax  = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
[849]454!
[1071]455!--                Leave loop if accuracy is sufficient, otherwise try again
456!--                with a reduced stepsize
[1359]457                   IF ( errmax <= 1.0_wp )  THEN
[1071]458                      EXIT
459                   ELSE
[1359]460                      dt_ros = SIGN( MAX( ABS( 0.9_wp * dt_ros *               &
461                                               errmax**pshrnk ),               &
462                                               shrnk * ABS( dt_ros ) ), dt_ros )
[1071]463                   ENDIF
464
465                ENDDO  ! loop for stepsize adjustment
466
467!
468!--             Calculate next internal time step
469                IF ( errmax > errcon )  THEN
[1359]470                   dt_ros_next = 0.9_wp * dt_ros * errmax**pgrow
[849]471                ELSE
[1071]472                   dt_ros_next = grow * dt_ros
[849]473                ENDIF
474
[1071]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
[849]482
[1071]483             ENDDO
[849]484!
[1071]485!--          Store internal time step value for next PALM step
486             particles(n)%rvar1 = dt_ros_next
[849]487
[1359]488             new_r(n) = r_ros
[849]489!
[1071]490!--          Radius should not fall below 1E-8 because Rosenbrock method may
491!--          lead to errors otherwise
[1359]492             new_r(n) = MAX( new_r(n), 1.0E-8_wp )
[1071]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
[1359]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.
[1071]501!--          Limiting values are tested for droplets down to 1.0E-7
[1359]502             IF ( new_r(n) - particles(n)%radius >= 3.0E-7_wp  .AND.  &
503                  e_a(n)/e_s(n) < 0.97_wp )  THEN
[1071]504                ros_count = ros_count + 1
505                repeat = .TRUE.
[849]506             ENDIF
507
[1071]508          ENDDO    ! Rosenbrock method
[849]509
510       ENDIF
511
[1359]512       delta_r = new_r(n) - particles(n)%radius
[849]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)
[1359]518       i = ip
519       j = jp
520       k = kp
[849]521           ! only exact if equidistant
522
[1359]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 ) /  &
[849]526                                   ( rho_surface * dx * dy * dz )
[1359]527       IF ( ql_c(k,j,i) > 100.0_wp )  THEN
[849]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
[1359]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,                                   &
[849]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)
[1359]548       ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * new_r(n)**3
[849]549
[1359]550       particles(n)%radius = new_r(n)
[849]551
552!
553!--    Determine radius class of the particle needed for collision
[1359]554       IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  use_kernel_tables )     &
[849]555       THEN
[1359]556          particles(n)%class = ( LOG( new_r(n) ) - rclass_lbound ) /           &
557                               ( rclass_ubound - rclass_lbound ) *             &
[849]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.