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

Last change on this file since 849 was 849, checked in by raasch, 10 years ago

Changed:


Original routine advec_particles split into several new subroutines and renamed
lpm.
init_particles renamed lpm_init
user_advec_particles renamed user_lpm_advec,
particle_boundary_conds renamed lpm_boundary_conds,
set_particle_attributes renamed lpm_set_attributes,
user_init_particles renamed user_lpm_init,
user_particle_attributes renamed user_lpm_set_attributes
(Makefile, lpm_droplet_collision, lpm_droplet_condensation, init_3d_model, modules, palm, read_var_list, time_integration, write_var_list, deleted: advec_particles, init_particles, particle_boundary_conds, set_particle_attributes, user_advec_particles, user_init_particles, user_particle_attributes, new: lpm, lpm_advec, lpm_boundary_conds, lpm_calc_liquid_water_content, lpm_data_output_particles, lpm_droplet_collision, lpm_drollet_condensation, lpm_exchange_horiz, lpm_extend_particle_array, lpm_extend_tails, lpm_extend_tail_array, lpm_init, lpm_init_sgs_tke, lpm_pack_arrays, lpm_read_restart_file, lpm_release_set, lpm_set_attributes, lpm_sort_arrays, lpm_write_exchange_statistics, lpm_write_restart_file, user_lpm_advec, user_lpm_init, user_lpm_set_attributes

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