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

Last change on this file since 2946 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

  • Property svn:keywords set to Id
File size: 17.5 KB
RevLine 
[1682]1!> @file lpm_droplet_condensation.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]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!
[2718]17! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[849]20! Current revisions:
21! ------------------
[2375]22!
23!
[1891]24! Former revisions:
25! -----------------
26! $Id: lpm_droplet_condensation.f90 2718 2018-01-02 08:49:38Z suehring $
[2716]27! Corrected "Former revisions" section
28!
29! 2696 2017-12-14 17:12:51Z kanani
30! Change in file header (GPL part)
31!
32! 2608 2017-11-13 14:04:26Z schwenkel
[2608]33! Calculation of magnus equation in external module (diagnostic_quantities_mod).
34!
35! 2375 2017-08-29 14:10:28Z schwenkel
[2375]36! Changed ONLY-dependencies
37!
38! 2312 2017-07-14 20:26:51Z hoffmann
[2312]39! Rosenbrock scheme improved. Gas-kinetic effect added.
[1891]40!
[2001]41! 2000 2016-08-20 18:09:15Z knoop
42! Forced header and separation lines into 80 columns
[2312]43!
[1891]44! 1890 2016-04-22 08:52:11Z hoffmann
[1890]45! Some improvements of the Rosenbrock method. If the Rosenbrock method needs more
46! than 40 iterations to find a sufficient time setp, the model is not aborted.
[2312]47! This might lead to small erros. But they can be assumend as negligible, since
48! the maximum timesetp of the Rosenbrock method after 40 iterations will be
49! smaller than 10^-16 s.
50!
[1872]51! 1871 2016-04-15 11:46:09Z hoffmann
52! Initialization of aerosols added.
53!
[1851]54! 1849 2016-04-08 11:33:18Z hoffmann
[1852]55! Interpolation of supersaturation has been removed because it is not in
56! accordance with the release/depletion of latent heat/water vapor in
[1849]57! interaction_droplets_ptq.
58! Calculation of particle Reynolds number has been corrected.
[1852]59! eps_ros added from modules.
[1849]60!
[1832]61! 1831 2016-04-07 13:15:51Z hoffmann
62! curvature_solution_effects moved to particle_attributes
63!
[1823]64! 1822 2016-04-07 07:49:42Z hoffmann
65! Unused variables removed.
66!
[1683]67! 1682 2015-10-07 23:56:08Z knoop
[2312]68! Code annotations made doxygen readable
69!
[1360]70! 1359 2014-04-11 17:15:14Z hoffmann
[2312]71! New particle structure integrated.
[1360]72! Kind definition added to all floating point numbers.
73!
[1347]74! 1346 2014-03-27 13:18:20Z heinze
[2312]75! Bugfix: REAL constants provided with KIND-attribute especially in call of
[1347]76! intrinsic function like MAX, MIN, SIGN
77!
[1323]78! 1322 2014-03-20 16:38:49Z raasch
79! REAL constants defined as wp-kind
80!
[1321]81! 1320 2014-03-20 08:40:49Z raasch
[1320]82! ONLY-attribute added to USE-statements,
83! kind-parameters added to all INTEGER and REAL declaration statements,
84! kinds are defined in new module kinds,
85! comment fields (!:) to be used for variable explanations added to
86! all variable declaration statements
[1072]87!
[1319]88! 1318 2014-03-17 13:35:16Z raasch
89! module interfaces removed
90!
[1093]91! 1092 2013-02-02 11:24:22Z raasch
92! unused variables removed
93!
[1072]94! 1071 2012-11-29 16:54:55Z franke
[1071]95! Ventilation effect for evaporation of large droplets included
96! Check for unreasonable results included in calculation of Rosenbrock method
97! since physically unlikely results were observed and for the same
98! reason the first internal time step in Rosenbrock method should be < 1.0E02 in
99! case of evaporation
100! Unnecessary calculation of ql_int removed
101! Unnecessary calculations in Rosenbrock method (d2rdt2, drdt_m, dt_ros_last)
102! removed
103! Bugfix: factor in calculation of surface tension changed from 0.00155 to
104! 0.000155
[849]105!
[1037]106! 1036 2012-10-22 13:43:42Z raasch
107! code put under GPL (PALM 3.9)
108!
[850]109! 849 2012-03-15 10:35:09Z raasch
110! initial revision (former part of advec_particles)
[849]111!
[850]112!
[849]113! Description:
114! ------------
[1682]115!> Calculates change in droplet radius by condensation/evaporation, using
116!> either an analytic formula or by numerically integrating the radius growth
117!> equation including curvature and solution effects using Rosenbrocks method
118!> (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
119!> The analytical formula and growth equation follow those given in
120!> Rogers and Yau (A short course in cloud physics, 3rd edition, p. 102/103).
[849]121!------------------------------------------------------------------------------!
[1682]122 SUBROUTINE lpm_droplet_condensation (ip,jp,kp)
[849]123
[2312]124
[1320]125    USE arrays_3d,                                                             &
[2312]126        ONLY:  hyp, pt, q, ql_c, ql_v
[849]127
[1320]128    USE cloud_parameters,                                                      &
[2375]129        ONLY:  l_d_rv, l_v, molecular_weight_of_solute,                        &
130               molecular_weight_of_water, rho_l, rho_s, r_v, vanthoff
[849]131
[1320]132    USE constants,                                                             &
133        ONLY:  pi
[849]134
[1320]135    USE control_parameters,                                                    &
[1822]136        ONLY:  dt_3d, dz, message_string, molecular_viscosity, rho_surface
137
[1320]138    USE cpulog,                                                                &
139        ONLY:  cpu_log, log_point_s
[849]140
[2608]141    USE diagnostic_quantities_mod,                                             &
142        ONLY:  magnus
143
[1320]144    USE grid_variables,                                                        &
[1822]145        ONLY:  dx, dy
[1071]146
[1320]147    USE lpm_collision_kernels_mod,                                             &
148        ONLY:  rclass_lbound, rclass_ubound
[849]149
[1320]150    USE kinds
151
152    USE particle_attributes,                                                   &
[2375]153        ONLY:  curvature_solution_effects, hall_kernel, number_of_particles,   &
154               particles, radius_classes, use_kernel_tables, wang_kernel
[1320]155
156    IMPLICIT NONE
157
[1682]158    INTEGER(iwp) :: ip                         !<
159    INTEGER(iwp) :: jp                         !<
160    INTEGER(iwp) :: kp                         !<
161    INTEGER(iwp) :: n                          !<
[1320]162
[1849]163    REAL(wp) ::  afactor                       !< curvature effects
[1682]164    REAL(wp) ::  arg                           !<
[1849]165    REAL(wp) ::  bfactor                       !< solute effects
[1682]166    REAL(wp) ::  ddenom                        !<
167    REAL(wp) ::  delta_r                       !<
[1849]168    REAL(wp) ::  diameter                      !< diameter of cloud droplets
[2312]169    REAL(wp) ::  diff_coeff                    !< diffusivity for water vapor
[1682]170    REAL(wp) ::  drdt                          !<
171    REAL(wp) ::  dt_ros                        !<
172    REAL(wp) ::  dt_ros_sum                    !<
173    REAL(wp) ::  d2rdtdr                       !<
[1849]174    REAL(wp) ::  e_a                           !< current vapor pressure
175    REAL(wp) ::  e_s                           !< current saturation vapor pressure
[2312]176    REAL(wp) ::  error                         !< local truncation error in Rosenbrock
177    REAL(wp) ::  k1                            !<
178    REAL(wp) ::  k2                            !<
179    REAL(wp) ::  r_err                         !< First order estimate of Rosenbrock radius
180    REAL(wp) ::  r_ros                         !< Rosenbrock radius
181    REAL(wp) ::  r_ros_ini                     !< initial Rosenbrock radius
182    REAL(wp) ::  r0                            !< gas-kinetic lengthscale
183    REAL(wp) ::  sigma                         !< surface tension of water
184    REAL(wp) ::  thermal_conductivity          !< thermal conductivity for water
[1849]185    REAL(wp) ::  t_int                         !< temperature
186    REAL(wp) ::  w_s                           !< terminal velocity of droplets
[2312]187    REAL(wp) ::  re_p                          !< particle Reynolds number
[1849]188!
[2312]189!-- Parameters for Rosenbrock method (see Verwer et al., 1999)
190    REAL(wp), PARAMETER :: prec = 1.0E-3_wp     !< precision of Rosenbrock solution
191    REAL(wp), PARAMETER :: q_increase = 1.5_wp  !< increase factor in timestep
192    REAL(wp), PARAMETER :: q_decrease = 0.9_wp  !< decrease factor in timestep
193    REAL(wp), PARAMETER :: gamma = 0.292893218814_wp !< = 1.0 - 1.0 / SQRT(2.0)
[849]194!
[1849]195!-- Parameters for terminal velocity
196    REAL(wp), PARAMETER ::  a_rog = 9.65_wp      !< parameter for fall velocity
197    REAL(wp), PARAMETER ::  b_rog = 10.43_wp     !< parameter for fall velocity
198    REAL(wp), PARAMETER ::  c_rog = 0.6_wp       !< parameter for fall velocity
199    REAL(wp), PARAMETER ::  k_cap_rog = 4.0_wp   !< parameter for fall velocity
200    REAL(wp), PARAMETER ::  k_low_rog = 12.0_wp  !< parameter for fall velocity
201    REAL(wp), PARAMETER ::  d0_rog = 0.745_wp    !< separation diameter
[849]202
[1849]203    REAL(wp), DIMENSION(number_of_particles) ::  ventilation_effect     !<
204    REAL(wp), DIMENSION(number_of_particles) ::  new_r                  !<
[849]205
[1849]206    CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'start' )
[849]207
208!
[2312]209!-- Absolute temperature
[1849]210    t_int = pt(kp,jp,ip) * ( hyp(kp) / 100000.0_wp )**0.286_wp
[849]211!
[2312]212!-- Saturation vapor pressure (Eq. 10 in Bolton, 1980)
[2608]213    e_s = magnus( t_int )
[1849]214!
[2312]215!-- Current vapor pressure
216    e_a = q(kp,jp,ip) * hyp(kp) / ( q(kp,jp,ip) + 0.622_wp )
217!
218!-- Thermal conductivity for water (from Rogers and Yau, Table 7.1)
219    thermal_conductivity = 7.94048E-05_wp * t_int + 0.00227011_wp
220!
221!-- Moldecular diffusivity of water vapor in air (Hall und Pruppacher, 1976)
222    diff_coeff           = 0.211E-4_wp * ( t_int / 273.15_wp )**1.94_wp * &
223                           ( 101325.0_wp / hyp(kp) )
224!
225!-- Lengthscale for gas-kinetic effects (from Mordy, 1959, p. 23):
226    r0 = diff_coeff / 0.036_wp * SQRT( 2.0_wp * pi / ( r_v * t_int ) )
227!
228!-- Calculate effects of heat conductivity and diffusion of water vapor on the
229!-- diffusional growth process (usually known as 1.0 / (F_k + F_d) )
230    ddenom  = 1.0_wp / ( rho_l * r_v * t_int / ( e_s * diff_coeff ) +          &
[1849]231                         ( l_v / ( r_v * t_int ) - 1.0_wp ) * rho_l *          &
[2312]232                         l_v / ( thermal_conductivity * t_int )                &
[1849]233                       )
[1359]234    new_r = 0.0_wp
[1849]235!
236!-- Determine ventilation effect on evaporation of large drops
[1359]237    DO  n = 1, number_of_particles
[1849]238
239       IF ( particles(n)%radius >= 4.0E-5_wp  .AND.  e_a / e_s < 1.0_wp )  THEN
[849]240!
[1849]241!--       Terminal velocity is computed for vertical direction (Rogers et al.,
242!--       1993, J. Appl. Meteorol.)
243          diameter = particles(n)%radius * 2000.0_wp !diameter in mm
244          IF ( diameter <= d0_rog )  THEN
245             w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
246          ELSE
247             w_s = a_rog - b_rog * EXP( -c_rog * diameter )
248          ENDIF
[849]249!
[2312]250!--       Calculate droplet's Reynolds number
[1849]251          re_p = 2.0_wp * particles(n)%radius * w_s / molecular_viscosity
[1071]252!
[1359]253!--       Ventilation coefficient (Rogers and Yau, 1989):
254          IF ( re_p > 2.5_wp )  THEN
[1849]255             ventilation_effect(n) = 0.78_wp + 0.28_wp * SQRT( re_p )
[1071]256          ELSE
[1849]257             ventilation_effect(n) = 1.0_wp + 0.09_wp * re_p
[1071]258          ENDIF
[1849]259       ELSE
[1071]260!
[2312]261!--       For small droplets or in supersaturated environments, the ventilation
[1849]262!--       effect does not play a role
263          ventilation_effect(n) = 1.0_wp
[849]264       ENDIF
[1359]265    ENDDO
[849]266
[2312]267    IF( .NOT. curvature_solution_effects ) then
[849]268!
[2312]269!--    Use analytic model for diffusional growth including gas-kinetic
270!--    effects (Mordy, 1959) but without the impact of aerosols.
[1849]271       DO  n = 1, number_of_particles
[2312]272          arg      = ( particles(n)%radius + r0 )**2 + 2.0_wp * dt_3d * ddenom * &
273                                                       ventilation_effect(n) *   &
274                                                       ( e_a / e_s - 1.0_wp )
275          arg      = MAX( arg, ( 0.01E-6 + r0 )**2 )
276          new_r(n) = SQRT( arg ) - r0
[1849]277       ENDDO
[1359]278
[2312]279    ELSE
[1849]280!
[2312]281!--    Integrate the diffusional growth including gas-kinetic (Mordy, 1959),
282!--    as well as curvature and solute effects (e.g., Köhler, 1936).
[849]283!
[2312]284!--    Curvature effect (afactor) with surface tension (sigma) by Straka (2009)
285       sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
[1071]286!
[2312]287!--    Solute effect (afactor)
288       afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
[849]289
[2312]290       DO  n = 1, number_of_particles
[849]291!
[2312]292!--       Solute effect (bfactor)
293          bfactor = vanthoff * rho_s * particles(n)%aux1**3 *                    &
294                    molecular_weight_of_water / ( rho_l * molecular_weight_of_solute )
[1071]295
[2312]296          dt_ros     = particles(n)%aux2  ! use previously stored Rosenbrock timestep
297          dt_ros_sum = 0.0_wp
[1871]298
[2312]299          r_ros     = particles(n)%radius  ! initialize Rosenbrock particle radius
300          r_ros_ini = r_ros
[849]301!
[2312]302!--       Integrate growth equation using a 2nd-order Rosenbrock method
303!--       (see Verwer et al., 1999, Eq. (3.2)). The Rosenbrock method adjusts
304!--       its with internal timestep to minimize the local truncation error.
305          DO WHILE ( dt_ros_sum < dt_3d )
[1071]306
[2312]307             dt_ros = MIN( dt_ros, dt_3d - dt_ros_sum )
[1871]308
[2312]309             DO
[849]310
[2312]311                drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -    &
[1849]312                                                          afactor / r_ros +    &
313                                                          bfactor / r_ros**3   &
[2312]314                                                        ) / ( r_ros + r0 )
[1849]315
[2312]316                d2rdtdr = -ddenom * ventilation_effect(n) * (                  &
317                                                (e_a / e_s - 1.0) * r_ros**4 - &
318                                                afactor * r0 * r_ros**2 -      &
319                                                2.0 * afactor * r_ros**3 +     &
320                                                3.0 * bfactor * r0 +           &
321                                                4.0 * bfactor * r_ros          &
322                                                            )                  &
323                          / ( r_ros**4 * ( r_ros + r0 )**2 )
[849]324
[2312]325                k1    = drdt / ( 1.0 - gamma * dt_ros * d2rdtdr )
[849]326
[2312]327                r_ros = MAX(r_ros_ini + k1 * dt_ros, particles(n)%aux1)
328                r_err = r_ros
[849]329
[2312]330                drdt  = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -   &
331                                                           afactor / r_ros +   &
332                                                           bfactor / r_ros**3  &
333                                                         ) / ( r_ros + r0 )
[849]334
[2312]335                k2 = ( drdt - dt_ros * 2.0 * gamma * d2rdtdr * k1 ) / &
336                     ( 1.0 - dt_ros * gamma * d2rdtdr )
[849]337
[2312]338                r_ros = MAX(r_ros_ini + dt_ros * ( 1.5 * k1 + 0.5 * k2), particles(n)%aux1)
339   !
340   !--          Check error of the solution, and reduce dt_ros if necessary.
341                error = ABS(r_err - r_ros) / r_ros
342                IF ( error .GT. prec )  THEN
343                   dt_ros = SQRT( q_decrease * prec / error ) * dt_ros
344                   r_ros  = r_ros_ini
[849]345                ELSE
[2312]346                   dt_ros_sum = dt_ros_sum + dt_ros
347                   dt_ros     = q_increase * dt_ros
348                   r_ros_ini  = r_ros
349                   EXIT
[849]350                ENDIF
351
[2312]352             END DO
[849]353
[2312]354          END DO !Rosenbrock loop
[849]355!
[2312]356!--       Store new particle radius
357          new_r(n) = r_ros
[849]358!
[2312]359!--       Store internal time step value for next PALM step
360          particles(n)%aux2 = dt_ros
[849]361
[2312]362       ENDDO !Particle loop
[849]363
[2312]364    ENDIF
[849]365
[2312]366    DO  n = 1, number_of_particles
[849]367!
[2312]368!--    Sum up the change in liquid water for the respective grid
369!--    box for the computation of the release/depletion of water vapor
370!--    and heat.
[1890]371       ql_c(kp,jp,ip) = ql_c(kp,jp,ip) + particles(n)%weight_factor *          &
[1359]372                                   rho_l * 1.33333333_wp * pi *                &
373                                   ( new_r(n)**3 - particles(n)%radius**3 ) /  &
[849]374                                   ( rho_surface * dx * dy * dz )
[2312]375!
376!--    Check if the increase in liqid water is not too big. If this is the case,
377!--    the model timestep might be too long.
[1890]378       IF ( ql_c(kp,jp,ip) > 100.0_wp )  THEN
379          WRITE( message_string, * ) 'k=',kp,' j=',jp,' i=',ip,      &
380                       ' ql_c=',ql_c(kp,jp,ip), ' &part(',n,')%wf=', &
[849]381                       particles(n)%weight_factor,' delta_r=',delta_r
382          CALL message( 'lpm_droplet_condensation', 'PA0143', 2, 2, -1, 6, 1 )
383       ENDIF
384!
[2312]385!--    Check if the change in the droplet radius is not too big. If this is the
386!--    case, the model timestep might be too long.
387       delta_r = new_r(n) - particles(n)%radius
388       IF ( delta_r < 0.0_wp  .AND. new_r(n) < 0.0_wp )  THEN
[1890]389          WRITE( message_string, * ) '#1 k=',kp,' j=',jp,' i=',ip,    &
[1849]390                       ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int,      &
[1890]391                       ' &delta_r=',delta_r,                          &
[849]392                       ' particle_radius=',particles(n)%radius
393          CALL message( 'lpm_droplet_condensation', 'PA0144', 2, 2, -1, 6, 1 )
394       ENDIF
395!
396!--    Sum up the total volume of liquid water (needed below for
397!--    re-calculating the weighting factors)
[1890]398       ql_v(kp,jp,ip) = ql_v(kp,jp,ip) + particles(n)%weight_factor * new_r(n)**3
[849]399!
400!--    Determine radius class of the particle needed for collision
[2312]401       IF ( use_kernel_tables )  THEN
[1359]402          particles(n)%class = ( LOG( new_r(n) ) - rclass_lbound ) /           &
403                               ( rclass_ubound - rclass_lbound ) *             &
[849]404                               radius_classes
405          particles(n)%class = MIN( particles(n)%class, radius_classes )
406          particles(n)%class = MAX( particles(n)%class, 1 )
407       ENDIF
[2312]408 !
409 !--   Store new radius to particle features
410       particles(n)%radius = new_r(n)
[849]411
412    ENDDO
413
414    CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'stop' )
415
416
417 END SUBROUTINE lpm_droplet_condensation
Note: See TracBrowser for help on using the repository browser.