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

Last change on this file since 1036 was 1036, checked in by raasch, 11 years ago

code has been put under the GNU General Public License (v3)

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