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

Last change on this file since 927 was 850, checked in by raasch, 12 years ago

last commit documented

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