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

Last change on this file since 1062 was 1037, checked in by raasch, 12 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 15.0 KB
Line 
1 SUBROUTINE lpm_droplet_condensation
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_droplet_condensation.f90 1037 2012-10-22 14:10:22Z maronga $
27!
28! 1036 2012-10-22 13:43:42Z raasch
29! code put under GPL (PALM 3.9)
30!
31! 849 2012-03-15 10:35:09Z raasch
32! initial revision (former part of advec_particles)
33!
34!
35! Description:
36! ------------
37! Calculates change in droplet radius by condensation/evaporation, using
38! either an analytic formula or by numerically integrating the radius growth
39! equation including curvature and solution effects using Rosenbrocks method
40! (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
41! The analytical formula and growth equation follow those given in
42! Rogers and Yau (A short course in cloud physics, 3rd edition, p. 102/103).
43!------------------------------------------------------------------------------!
44
45    USE arrays_3d
46    USE cloud_parameters
47    USE constants
48    USE control_parameters
49    USE cpulog
50    USE grid_variables
51    USE interfaces
52    USE lpm_collision_kernels_mod
53    USE particle_attributes
54
55    IMPLICIT NONE
56
57    INTEGER ::  i, internal_timestep_count, j, jtry, k, n
58
59    INTEGER, PARAMETER ::  maxtry = 40
60
61    REAL ::  aa, afactor, arg, bb, cc, dd, ddenom, delta_r, drdt, drdt_ini,    &
62             drdt_m, dt_ros, dt_ros_last, dt_ros_next, dt_ros_sum,             &
63             dt_ros_sum_ini, d2rdtdr, d2rdt2, errmax, err_ros, g1, g2, g3, g4, &
64             e_a, e_s, gg, new_r, p_int, pt_int, pt_int_l, pt_int_u, q_int,    &
65             q_int_l, q_int_u, ql_int, ql_int_l, ql_int_u, r_ros, r_ros_ini,   &
66             sigma, t_int, x, y
67
68!
69!-- Parameters for Rosenbrock method
70    REAL, PARAMETER ::  a21 = 2.0, a31 = 48.0/25.0, a32 = 6.0/25.0,        &
71                        a2x = 1.0, a3x = 3.0/5.0, b1 = 19.0/9.0, b2 = 0.5, &
72                        b3 = 25.0/108.0, b4 = 125.0/108.0, c21 = -8.0,     &
73                        c31 = 372.0/25.0, c32 = 12.0/5.0,                  &
74                        c41 = -112.0/125.0, c42 = -54.0/125.0,             &
75                        c43 = -2.0/5.0, c1x = 0.5, c2x= -3.0/2.0,          &
76                        c3x = 121.0/50.0, c4x = 29.0/250.0,                &
77                        errcon = 0.1296, e1 = 17.0/54.0, e2 = 7.0/36.0,    &
78                        e3 = 0.0, e4 = 125.0/108.0, gam = 0.5, grow = 1.5, &
79                        pgrow = -0.25, pshrnk = -1.0/3.0, shrnk = 0.5
80
81
82    CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'start' )
83
84    DO  n = 1, number_of_particles
85!
86!--    Interpolate temperature and humidity.
87!--    First determine left, south, and bottom index of the arrays.
88       i = particles(n)%x * ddx
89       j = particles(n)%y * ddy
90       k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
91           + offset_ocean_nzt                     ! only exact if equidistant
92
93       x  = particles(n)%x - i * dx
94       y  = particles(n)%y - j * dy
95       aa = x**2          + y**2
96       bb = ( dx - x )**2 + y**2
97       cc = x**2          + ( dy - y )**2
98       dd = ( dx - x )**2 + ( dy - y )**2
99       gg = aa + bb + cc + dd
100
101       pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) * pt(k,j,i+1)   &
102                  + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
103                  ) / ( 3.0 * gg )
104
105       pt_int_u = ( ( gg-aa ) * pt(k+1,j,i)   + ( gg-bb ) * pt(k+1,j,i+1)   &
106                  + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
107                  ) / ( 3.0 * gg )
108
109       pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
110                           ( pt_int_u - pt_int_l )
111
112       q_int_l = ( ( gg - aa ) * q(k,j,i)   + ( gg - bb ) * q(k,j,i+1)   &
113                 + ( gg - cc ) * q(k,j+1,i) + ( gg - dd ) * q(k,j+1,i+1) &
114                 ) / ( 3.0 * gg )
115
116       q_int_u = ( ( gg-aa ) * q(k+1,j,i)   + ( gg-bb ) * q(k+1,j,i+1)   &
117                 + ( gg-cc ) * q(k+1,j+1,i) + ( gg-dd ) * q(k+1,j+1,i+1) &
118                 ) / ( 3.0 * gg )
119
120       q_int = q_int_l + ( particles(n)%z - zu(k) ) / dz * &
121                           ( q_int_u - q_int_l )
122
123       ql_int_l = ( ( gg - aa ) * ql(k,j,i)   + ( gg - bb ) * ql(k,j,i+1)   &
124                  + ( gg - cc ) * ql(k,j+1,i) + ( gg - dd ) * ql(k,j+1,i+1) &
125                  ) / ( 3.0 * gg )
126
127       ql_int_u = ( ( gg-aa ) * ql(k+1,j,i)   + ( gg-bb ) * ql(k+1,j,i+1)   &
128                  + ( gg-cc ) * ql(k+1,j+1,i) + ( gg-dd ) * ql(k+1,j+1,i+1) &
129                  ) / ( 3.0 * gg )
130
131       ql_int = ql_int_l + ( particles(n)%z - zu(k) ) / dz * &
132                             ( ql_int_u - ql_int_l )
133
134!
135!--    Calculate real temperature and saturation vapor pressure
136       p_int = hyp(k) + ( particles(n)%z - zu(k) ) / dz * ( hyp(k+1)-hyp(k) )
137       t_int = pt_int * ( p_int / 100000.0 )**0.286
138
139       e_s = 611.0 * EXP( l_d_rv * ( 3.6609E-3 - 1.0 / t_int ) )
140
141!
142!--    Current vapor pressure
143       e_a = q_int * p_int / ( 0.378 * q_int + 0.622 )
144
145!
146!--    Thermal conductivity for water (from Rogers and Yau, Table 7.1),
147!--    diffusivity for water vapor (after Hall und Pruppacher, 1976)
148       thermal_conductivity_l = 7.94048E-05 * t_int + 0.00227011 
149       diff_coeff_l           = 0.211E-4 * ( t_int / 273.15 )**1.94 * &
150                                ( 101325.0 / p_int)
151!
152!--    Change in radius by condensation/evaporation
153       IF ( particles(n)%radius >= 1.0E-6  .OR.  &
154            .NOT. curvature_solution_effects )  THEN
155!
156!--       Approximation for large radii, where curvature and solution
157!--       effects can be neglected
158          arg = particles(n)%radius**2 + 2.0 * dt_3d *                        &
159                           ( e_a / e_s - 1.0 ) /                              &
160                           ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / &
161                             thermal_conductivity_l +                         &
162                             rho_l * r_v * t_int / diff_coeff_l / e_s )
163          IF ( arg < 1.0E-16 )  THEN
164             new_r = 1.0E-8
165          ELSE
166             new_r = SQRT( arg )
167          ENDIF
168       ENDIF
169
170       IF ( curvature_solution_effects  .AND.                              &
171            ( ( particles(n)%radius < 1.0E-6 ) .OR. ( new_r < 1.0E-6 ) ) ) &
172       THEN
173!
174!--       Curvature and solutions effects are included in growth equation.
175!--       Change in Radius is calculated with 4th-order Rosenbrock method
176!--       for stiff o.d.e's with monitoring local truncation error to adjust
177!--       stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
178!--       For larger radii the simple analytic method (see ELSE) gives
179!--       almost the same results.
180!
181!--       Surface tension after (Straka, 2009)
182          sigma = 0.0761 - 0.00155 * ( t_int - 273.15 )
183
184          r_ros = particles(n)%radius
185          dt_ros_sum  = 0.0      ! internal integrated time (s)
186          internal_timestep_count = 0
187
188          ddenom  = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) +  &
189                            ( l_v / ( r_v * t_int ) - 1.0 ) *               &
190                            rho_l * l_v / ( thermal_conductivity_l * t_int )&
191                          )
192             
193          afactor = 2.0 * sigma / ( rho_l * r_v * t_int )
194
195          IF ( particles(n)%rvar3 == -9999999.9 )  THEN
196!
197!--          First particle timestep. Derivative has to be calculated.
198             drdt_m  = ddenom / r_ros * ( e_a / e_s - 1.0 -    &
199                                          afactor / r_ros +    &
200                                          bfactor / r_ros**3 )
201          ELSE
202!
203!--          Take value from last PALM timestep
204             drdt_m = particles(n)%rvar3
205          ENDIF
206!
207!--       Take internal timestep values from the end of last PALM timestep
208          dt_ros_last = particles(n)%rvar1
209          dt_ros_next = particles(n)%rvar2
210!
211!--       Internal timestep must not be larger than PALM timestep
212          dt_ros = MIN( dt_ros_next, dt_3d )
213!
214!--       Integrate growth equation in time unless PALM timestep is reached
215          DO WHILE ( dt_ros_sum < dt_3d )
216
217             internal_timestep_count = internal_timestep_count + 1
218
219!
220!--          Derivative at starting value
221             drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + &
222                                       bfactor / r_ros**3 )
223             drdt_ini       = drdt
224             dt_ros_sum_ini = dt_ros_sum
225             r_ros_ini      = r_ros
226
227!
228!--          Calculate time derivative of dr/dt
229             d2rdt2 = ( drdt - drdt_m ) / dt_ros_last 
230
231!
232!--          Calculate radial derivative of dr/dt
233             d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + &
234                                  2.0 * afactor / r_ros**3 -     &
235                                  4.0 * bfactor / r_ros**5 )
236!
237!--          Adjust stepsize unless required accuracy is reached
238             DO  jtry = 1, maxtry+1
239
240                IF ( jtry == maxtry+1 )  THEN
241                   message_string = 'maxtry > 40 in Rosenbrock method'
242                   CALL message( 'lpm_droplet_condensation', 'PA0347', 2, 2, &
243                                 0, 6, 0 )
244                ENDIF
245
246                aa    = 1.0 / ( gam * dt_ros ) - d2rdtdr
247                g1    = ( drdt_ini + dt_ros * c1x * d2rdt2 ) / aa
248                r_ros = r_ros_ini + a21 * g1
249                drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
250                                           afactor / r_ros + &
251                                           bfactor / r_ros**3 )
252
253                g2    = ( drdt + dt_ros * c2x * d2rdt2 + c21 * g1 / dt_ros )&
254                        / aa
255                r_ros = r_ros_ini + a31 * g1 + a32 * g2
256                drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
257                                           afactor / r_ros + &
258                                           bfactor / r_ros**3 )
259
260                g3    = ( drdt + dt_ros * c3x * d2rdt2 + &
261                          ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
262                g4    = ( drdt + dt_ros * c4x * d2rdt2 + &
263                          ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
264                r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
265
266                dt_ros_sum = dt_ros_sum_ini + dt_ros
267
268                IF ( dt_ros_sum == dt_ros_sum_ini )  THEN
269                   message_string = 'zero stepsize in Rosenbrock method'
270                   CALL message( 'lpm_droplet_condensation', 'PA0348', 2, 2, &
271                                 0, 6, 0 )
272                ENDIF
273!
274!--             Calculate error
275                err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4 
276                errmax  = 0.0
277                errmax  = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
278!
279!--             Leave loop if accuracy is sufficient, otherwise try again
280!--             with a reduced stepsize
281                IF ( errmax <= 1.0 )  THEN
282                   EXIT
283                ELSE
284                   dt_ros = SIGN( MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), &
285                                       shrnk * ABS( dt_ros ) ), dt_ros )
286                ENDIF
287
288             ENDDO  ! loop for stepsize adjustment
289
290!
291!--          Calculate next internal timestep
292             IF ( errmax > errcon )  THEN
293                dt_ros_next = 0.9 * dt_ros * errmax**pgrow
294             ELSE
295                dt_ros_next = grow * dt_ros
296             ENDIF
297
298!
299!--          Estimated timestep is reduced if the PALM time step is exceeded
300             dt_ros_last = dt_ros
301             IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d )  THEN
302                dt_ros = dt_3d - dt_ros_sum
303             ELSE
304                dt_ros = dt_ros_next
305             ENDIF
306
307             drdt_m = drdt
308
309          ENDDO
310!
311!--       Store derivative and internal timestep values for next PALM step
312          particles(n)%rvar1 = dt_ros_last
313          particles(n)%rvar2 = dt_ros_next
314          particles(n)%rvar3 = drdt_m
315
316          new_r = r_ros
317!
318!--       Radius should not fall below 1E-8 because Rosenbrock method may
319!--       lead to errors otherwise
320          new_r = MAX( new_r, 1.0E-8 )
321
322       ENDIF
323
324       delta_r = new_r - particles(n)%radius
325
326!
327!--    Sum up the change in volume of liquid water for the respective grid
328!--    volume (this is needed later in lpm_calc_liquid_water_content for
329!--    calculating the release of latent heat)
330       i = ( particles(n)%x + 0.5 * dx ) * ddx
331       j = ( particles(n)%y + 0.5 * dy ) * ddy
332       k = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
333           ! only exact if equidistant
334
335       ql_c(k,j,i) = ql_c(k,j,i) + particles(n)%weight_factor *            &
336                                   rho_l * 1.33333333 * pi *               &
337                                   ( new_r**3 - particles(n)%radius**3 ) / &
338                                   ( rho_surface * dx * dy * dz )
339       IF ( ql_c(k,j,i) > 100.0 )  THEN
340          WRITE( message_string, * ) 'k=',k,' j=',j,' i=',i,      &
341                       ' ql_c=',ql_c(k,j,i), ' &part(',n,')%wf=', &
342                       particles(n)%weight_factor,' delta_r=',delta_r
343          CALL message( 'lpm_droplet_condensation', 'PA0143', 2, 2, -1, 6, 1 )
344       ENDIF
345
346!
347!--    Change the droplet radius
348       IF ( ( new_r - particles(n)%radius ) < 0.0  .AND.  new_r < 0.0 ) &
349       THEN
350          WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i,   &
351                       ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int,  &
352                       ' &delta_r=',delta_r,                      &
353                       ' particle_radius=',particles(n)%radius
354          CALL message( 'lpm_droplet_condensation', 'PA0144', 2, 2, -1, 6, 1 )
355       ENDIF
356
357!
358!--    Sum up the total volume of liquid water (needed below for
359!--    re-calculating the weighting factors)
360       ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * new_r**3
361
362       particles(n)%radius = new_r
363
364!
365!--    Determine radius class of the particle needed for collision
366       IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  use_kernel_tables ) &
367       THEN
368          particles(n)%class = ( LOG( new_r ) - rclass_lbound ) /  &
369                               ( rclass_ubound - rclass_lbound ) * &
370                               radius_classes
371          particles(n)%class = MIN( particles(n)%class, radius_classes )
372          particles(n)%class = MAX( particles(n)%class, 1 )
373       ENDIF
374
375    ENDDO
376
377    CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'stop' )
378
379
380 END SUBROUTINE lpm_droplet_condensation
Note: See TracBrowser for help on using the repository browser.