source: palm/trunk/SOURCE/lpm_splitting.f90 @ 3655

Last change on this file since 3655 was 3655, checked in by knoop, 5 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

  • Property svn:keywords set to Id
File size: 30.8 KB
Line 
1!> @file lpm_splitting.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
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.
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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_splitting.f90 3655 2019-01-07 16:51:22Z knoop $
27! Modularization of all bulk cloud physics code components
28!
29! 3241 2018-09-12 15:02:00Z raasch
30! unused variables removed
31!
32! 2932 2018-03-26 09:39:22Z maronga
33! renamed particles_par to particle_parameters
34!
35! 2718 2018-01-02 08:49:38Z maronga
36! Corrected "Former revisions" section
37!
38!
39! Change in file header (GPL part)
40!
41! Added comments
42!
43!
44! 2263 2017-06-08 14:59:01Z schwenkel
45! Initial revision
46!
47!
48!
49! Description:
50! ------------
51! This routine is a part of the Lagrangian particle model. Super droplets which
52! fulfill certain criterion's (e.g. a big weighting factor and a large radius)
53! can be split into several super droplets with a reduced number of
54! represented particles of every super droplet. This mechanism ensures an
55! improved representation of the right tail of the drop size distribution with
56! a feasible amount of computational costs. The limits of particle creation
57! should be chosen carefully! The idea of this algorithm is based on
58! Unterstrasser and Soelch, 2014.
59!------------------------------------------------------------------------------!
60 SUBROUTINE lpm_splitting
61
62
63    USE arrays_3d,                                                             &
64        ONLY:  ql
65
66    USE basic_constants_and_equations_mod,                                     &
67        ONLY:  pi, rho_l
68
69    USE cpulog,                                                                &
70        ONLY:  cpu_log, log_point_s
71
72    USE indices,                                                               &
73        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
74
75    USE kinds
76
77    USE lpm_exchange_horiz_mod,                                                &
78        ONLY:  realloc_particles_array
79
80    USE particle_attributes,                                                   &
81        ONLY:  grid_particles, initial_weighting_factor, isf, i_splitting_mode,&
82               max_number_particles_per_gridbox, new_particles, n_max,         &
83               number_of_particles, particles, particle_type, prt_count,       &
84               radius_split, splitting_factor, splitting_factor_max,           &
85               sum_new_particles, weight_factor_split                       
86
87    USE pegrid
88
89    IMPLICIT NONE
90
91    INTEGER(iwp) ::  i                !<
92    INTEGER(iwp) ::  j                !<
93    INTEGER(iwp) ::  jpp              !<
94    INTEGER(iwp) ::  k                !<
95    INTEGER(iwp) ::  n                !<
96    INTEGER(iwp) ::  new_particles_gb !< counter of created particles within one grid box
97    INTEGER(iwp) ::  new_size         !< new particle array size
98    INTEGER(iwp) ::  np               !<
99    INTEGER(iwp) ::  old_size         !< old particle array size
100   
101    LOGICAL ::  first_loop_stride = .TRUE. !< flag to calculate constants only once
102
103    REAL(wp) ::  diameter                 !< diameter of droplet
104    REAL(wp) ::  dlog                     !< factor for DSD calculation
105    REAL(wp) ::  factor_volume_to_mass    !< pre calculate factor volume to mass
106    REAL(wp) ::  lambda                   !< slope parameter of gamma-distribution
107    REAL(wp) ::  lwc                      !< liquid water content of grid box
108    REAL(wp) ::  lwc_total                !< average liquid water content of cloud
109    REAL(wp) ::  m1                       !< first moment of DSD
110    REAL(wp) ::  m1_total                 !< average over all PEs of first moment of DSD
111    REAL(wp) ::  m2                       !< second moment of DSD
112    REAL(wp) ::  m2_total                 !< average average over all PEs second moment of DSD
113    REAL(wp) ::  m3                       !< third moment of DSD
114    REAL(wp) ::  m3_total                 !< average average over all PEs third moment of DSD
115    REAL(wp) ::  mu                       !< spectral shape parameter of gamma distribution
116    REAL(wp) ::  nrclgb                   !< number of cloudy grid boxes (ql >= 1.0E-5 kg/kg)
117    REAL(wp) ::  nrclgb_total             !< average over all PEs of number of cloudy grid boxes
118    REAL(wp) ::  nr                       !< number concentration of cloud droplets
119    REAL(wp) ::  nr_total                 !< average over all PEs of number of cloudy grid boxes
120    REAL(wp) ::  nr0                      !< intercept parameter of gamma distribution
121    REAL(wp) ::  pirho_l                  !< pi * rho_l / 6.0
122    REAL(wp) ::  ql_crit = 1.0E-5_wp      !< threshold lwc for cloudy grid cells
123                                          !< (Siebesma et al 2003, JAS, 60)
124    REAL(wp) ::  rm                       !< volume averaged mean radius
125    REAL(wp) ::  rm_total                 !< average over all PEs of volume averaged mean radius
126    REAL(wp) ::  r_min = 1.0E-6_wp        !< minimum radius of approximated spectra
127    REAL(wp) ::  r_max = 1.0E-3_wp        !< maximum radius of approximated spectra
128    REAL(wp) ::  sigma_log = 1.5_wp       !< standard deviation of the LOG-distribution
129    REAL(wp) ::  zeta                     !< Parameter for DSD calculation of Seifert
130
131    REAL(wp), DIMENSION(0:n_max-1) ::  an_spl     !< size dependent critical weight factor
132    REAL(wp), DIMENSION(0:n_max-1) ::  r_bin_mid  !< mass weighted mean radius of a bin
133    REAL(wp), DIMENSION(0:n_max)   ::  r_bin      !< boundaries of a radius bin
134   
135    TYPE(particle_type) ::  tmp_particle   !< temporary particle TYPE
136
137    CALL cpu_log( log_point_s(80), 'lpm_splitting', 'start' )
138
139    IF ( first_loop_stride )  THEN
140       IF ( i_splitting_mode == 2  .OR.  i_splitting_mode == 3 )  THEN
141          dlog   = ( LOG10(r_max) - LOG10(r_min) ) / ( n_max - 1 )
142          DO  i = 0, n_max-1
143             r_bin(i) = 10.0_wp**( LOG10(r_min) + i * dlog - 0.5_wp * dlog )
144             r_bin_mid(i) = 10.0_wp**( LOG10(r_min) + i * dlog )
145          ENDDO
146          r_bin(n_max) = 10.0_wp**( LOG10(r_min) + n_max * dlog - 0.5_wp * dlog )
147       ENDIF   
148       factor_volume_to_mass =  4.0_wp / 3.0_wp * pi * rho_l
149       pirho_l  = pi * rho_l / 6.0_wp
150       IF ( weight_factor_split == -1.0_wp )  THEN
151          weight_factor_split = 0.1_wp * initial_weighting_factor 
152       ENDIF
153    ENDIF
154
155    new_particles  = 0
156
157    IF ( i_splitting_mode == 1 )  THEN
158
159       DO  i = nxl, nxr
160          DO  j = nys, nyn
161             DO  k = nzb+1, nzt
162
163                new_particles_gb = 0
164                number_of_particles = prt_count(k,j,i)
165                IF ( number_of_particles <= 0  .OR.                            & 
166                     ql(k,j,i) < ql_crit )  CYCLE
167                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
168!               
169!--             Start splitting operations. Each particle is checked if it
170!--             fulfilled the splitting criterion's. In splitting mode 'const'   
171!--             a critical radius  (radius_split) a critical weighting factor
172!--             (weight_factor_split) and a splitting factor (splitting_factor)
173!--             must  be prescribed (see particle_parameters). Super droplets
174!--             which have a larger radius and larger weighting factor are split
175!--             into 'splitting_factor' super droplets. Therefore, the weighting
176!--             factor of  the super droplet and all created clones is reduced
177!--             by the factor of 'splitting_factor'.
178                DO  n = 1, number_of_particles
179                   IF ( particles(n)%particle_mask  .AND.                      &
180                        particles(n)%radius >= radius_split  .AND.             & 
181                        particles(n)%weight_factor >= weight_factor_split )    &
182                   THEN         
183!
184!--                   Calculate the new number of particles.
185                      new_size = prt_count(k,j,i) + splitting_factor - 1                                   
186!
187!--                   Cycle if maximum number of particles per grid box
188!--                   is greater than the allowed maximum number.
189                      IF ( new_size >= max_number_particles_per_gridbox )  CYCLE                     
190!
191!--                   Reallocate particle array if necessary.
192                      IF ( new_size > SIZE(particles) )  THEN
193                         CALL realloc_particles_array(i,j,k,new_size)
194                      ENDIF
195                      old_size = prt_count(k,j,i)
196!
197!--                   Calculate new weighting factor.
198                      particles(n)%weight_factor =  & 
199                         particles(n)%weight_factor / splitting_factor
200                      tmp_particle = particles(n)
201!
202!--                   Create splitting_factor-1 new particles.
203                      DO  jpp = 1, splitting_factor-1
204                         grid_particles(k,j,i)%particles(jpp+old_size) =       & 
205                            tmp_particle         
206                      ENDDO 
207                      new_particles_gb = new_particles_gb + splitting_factor - 1
208!   
209!--                   Save the new number of super droplets for every grid box.
210                      prt_count(k,j,i) = prt_count(k,j,i) +                    &
211                                         splitting_factor - 1         
212                   ENDIF
213                ENDDO
214               
215                new_particles       = new_particles     + new_particles_gb
216                sum_new_particles   = sum_new_particles + new_particles_gb 
217             ENDDO
218          ENDDO
219       ENDDO
220
221    ELSEIF ( i_splitting_mode == 2 )  THEN 
222!
223!--    Initialize summing variables.
224       lwc          = 0.0_wp
225       lwc_total    = 0.0_wp 
226       m1           = 0.0_wp
227       m1_total     = 0.0_wp
228       m2           = 0.0_wp
229       m2_total     = 0.0_wp
230       m3           = 0.0_wp
231       m3_total     = 0.0_wp
232       nr           = 0.0_wp   
233       nrclgb       = 0.0_wp
234       nrclgb_total = 0.0_wp 
235       nr_total     = 0.0_wp
236       rm           = 0.0_wp
237       rm_total     = 0.0_wp
238       
239       DO  i = nxl, nxr
240          DO  j = nys, nyn
241             DO  k = nzb+1, nzt
242                number_of_particles = prt_count(k,j,i)
243                IF ( number_of_particles <= 0  .OR.                            & 
244                     ql(k,j,i) < ql_crit )  CYCLE
245                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
246                nrclgb = nrclgb + 1.0_wp
247!               
248!--             Calculate moments of DSD.               
249                DO  n = 1, number_of_particles
250                   IF ( particles(n)%particle_mask  .AND.                      &
251                        particles(n)%radius >= r_min )                         &
252                   THEN
253                      nr  = nr  + particles(n)%weight_factor
254                      rm  = rm  + factor_volume_to_mass  *                     &
255                                 particles(n)%radius**3  *                     &
256                                 particles(n)%weight_factor
257                      IF ( isf == 1 )  THEN           
258                         diameter   = particles(n)%radius * 2.0_wp           
259                         lwc = lwc + factor_volume_to_mass *                   &
260                                     particles(n)%radius**3 *                  & 
261                                     particles(n)%weight_factor 
262                         m1  = m1  + particles(n)%weight_factor * diameter                                               
263                         m2  = m2  + particles(n)%weight_factor * diameter**2           
264                         m3  = m3  + particles(n)%weight_factor * diameter**3
265                      ENDIF   
266                   ENDIF                       
267                ENDDO 
268             ENDDO
269          ENDDO
270       ENDDO
271
272#if defined( __parallel )
273       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
274       CALL MPI_ALLREDUCE( nr, nr_total, 1 , &
275       MPI_REAL, MPI_SUM, comm2d, ierr )
276       CALL MPI_ALLREDUCE( rm, rm_total, 1 , &
277       MPI_REAL, MPI_SUM, comm2d, ierr )
278       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
279       CALL MPI_ALLREDUCE( nrclgb, nrclgb_total, 1 , &
280       MPI_REAL, MPI_SUM, comm2d, ierr )
281       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
282       CALL MPI_ALLREDUCE( lwc, lwc_total, 1 , &
283       MPI_REAL, MPI_SUM, comm2d, ierr )
284       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
285       CALL MPI_ALLREDUCE( m1, m1_total, 1 , &
286       MPI_REAL, MPI_SUM, comm2d, ierr )
287       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
288       CALL MPI_ALLREDUCE( m2, m2_total, 1 , &
289       MPI_REAL, MPI_SUM, comm2d, ierr )
290       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
291       CALL MPI_ALLREDUCE( m3, m3_total, 1 , &
292       MPI_REAL, MPI_SUM, comm2d, ierr )
293#endif 
294
295!
296!--    Calculate number concentration and mean volume averaged radius.
297       nr_total = MERGE( nr_total / nrclgb_total,                              &
298                         0.0_wp, nrclgb_total > 0.0_wp                         &
299                       )
300       rm_total = MERGE( ( rm_total /                                          &
301                            ( nr_total * factor_volume_to_mass )               &
302                          )**0.3333333_wp, 0.0_wp, nrclgb_total > 0.0_wp       &
303                       )                         
304!
305!--    Check which function should be used to approximate the DSD.
306       IF ( isf == 1 )  THEN
307          lwc_total = MERGE( lwc_total / nrclgb_total,                         &
308                             0.0_wp, nrclgb_total > 0.0_wp                     &
309                           )
310          m1_total  = MERGE( m1_total / nrclgb_total,                          &
311                             0.0_wp, nrclgb_total > 0.0_wp                     &
312                           )
313          m2_total  = MERGE( m2_total / nrclgb_total,                          &
314                             0.0_wp, nrclgb_total > 0.0_wp                     &
315                           )
316          m3_total  = MERGE( m3_total / nrclgb_total,                          &
317                             0.0_wp, nrclgb_total > 0.0_wp                     &
318                           )
319          zeta = m1_total * m3_total / m2_total**2                             
320          mu   = MAX( ( ( 1.0_wp - zeta ) * 2.0_wp + 1.0_wp ) /                &
321                        ( zeta - 1.0_wp ), 0.0_wp                              &
322                    )
323
324          lambda = ( pirho_l * nr_total / lwc_total *                          &
325                     ( mu + 3.0_wp ) * ( mu + 2.0_wp ) * ( mu + 1.0_wp )       &                                         
326                   )**0.3333333_wp
327          nr0 = nr_total / gamma( mu + 1.0_wp ) * lambda**( mu + 1.0_wp ) 
328         
329          DO  n = 0, n_max-1
330             diameter  = r_bin_mid(n) * 2.0_wp           
331             an_spl(n) = nr0 * diameter**mu * EXP( -lambda * diameter ) *      & 
332                         ( r_bin(n+1) - r_bin(n) ) * 2.0_wp 
333          ENDDO
334       ELSEIF ( isf == 2 )  THEN
335          DO  n = 0, n_max-1
336             an_spl(n) = nr_total / ( SQRT( 2.0_wp * pi ) *                    &
337                                     LOG(sigma_log) * r_bin_mid(n)             &
338                                     ) *                                       &
339                         EXP( -( LOG( r_bin_mid(n) / rm_total )**2 ) /         &
340                               ( 2.0_wp * LOG(sigma_log)**2 )                  & 
341                             ) *                                               & 
342                         ( r_bin(n+1) - r_bin(n) )
343          ENDDO
344       ELSEIF( isf == 3 )  THEN
345          DO  n = 0, n_max-1 
346             an_spl(n) = 3.0_wp * nr_total * r_bin_mid(n)**2 / rm_total**3  *  &
347                         EXP( - ( r_bin_mid(n)**3 / rm_total**3 ) )         *  &
348                         ( r_bin(n+1) - r_bin(n) )
349          ENDDO
350       ENDIF
351!
352!--    Criterion to avoid super droplets with a weighting factor < 1.0.
353       an_spl = MAX(an_spl, 1.0_wp)
354 
355       DO  i = nxl, nxr
356          DO  j = nys, nyn
357             DO  k = nzb+1, nzt
358                number_of_particles = prt_count(k,j,i)
359                IF ( number_of_particles <= 0  .OR.                            & 
360                     ql(k,j,i) < ql_crit )  CYCLE
361                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
362                new_particles_gb = 0         
363!               
364!--             Start splitting operations. Each particle is checked if it
365!--             fulfilled the splitting criterion's. In splitting mode 'cl_av'   
366!--             a critical radius (radius_split) and a splitting function must
367!--             be prescribed (see particles_par). The critical weighting factor
368!--             is calculated while approximating a 'gamma', 'log' or 'exp'-
369!--             drop size distribution. In this mode the DSD is calculated as
370!--             an average over all cloudy grid boxes. Super droplets which
371!--             have a larger radius and larger weighting factor are split into
372!--             'splitting_factor' super droplets. In this case the splitting 
373!--             factor is calculated of weighting factor of the super droplet 
374!--             and the approximated number concentration for droplet of such
375!--             a size. Due to the splitting, the weighting factor of the 
376!--             super droplet and all created clones is reduced by the factor 
377!--             of 'splitting_facor'.
378                DO  n = 1, number_of_particles
379                   DO  np = 0, n_max-1
380                      IF ( r_bin(np) >= radius_split  .AND.                    &
381                           particles(n)%particle_mask  .AND.                   &
382                           particles(n)%radius >= r_bin(np)  .AND.             &
383                           particles(n)%radius < r_bin(np+1)  .AND.            &
384                           particles(n)%weight_factor >= an_spl(np)  )         &
385                      THEN
386!
387!--                      Calculate splitting factor
388                         splitting_factor =                                    & 
389                             MIN( INT( particles(n)%weight_factor /            &
390                                        an_spl(np)                             &
391                                     ), splitting_factor_max                   &
392                                )
393                         IF ( splitting_factor < 2 )  CYCLE
394!
395!--                      Calculate the new number of particles.                                                           
396                         new_size = prt_count(k,j,i) + splitting_factor - 1
397!
398!--                      Cycle if maximum number of particles per grid box
399!--                      is greater than the allowed maximum number.                         
400                         IF ( new_size >= max_number_particles_per_gridbox )   & 
401                         CYCLE
402!
403!--                      Reallocate particle array if necessary.
404                         IF ( new_size > SIZE(particles) )  THEN
405                            CALL realloc_particles_array(i,j,k,new_size)
406                         ENDIF
407                         old_size  = prt_count(k,j,i)                             
408                         new_particles_gb = new_particles_gb +                 &
409                                            splitting_factor - 1
410!
411!--                      Calculate new weighting factor.
412                         particles(n)%weight_factor =                          & 
413                            particles(n)%weight_factor / splitting_factor
414                         tmp_particle = particles(n)
415!
416!--                      Create splitting_factor-1 new particles.                                                 
417                         DO  jpp = 1, splitting_factor-1
418                            grid_particles(k,j,i)%particles(jpp+old_size) =    &
419                                                                    tmp_particle         
420                         ENDDO
421!   
422!--                      Save the new number of super droplets.
423                         prt_count(k,j,i) = prt_count(k,j,i) +                 &
424                                            splitting_factor - 1
425                      ENDIF
426                   ENDDO
427                ENDDO 
428               
429                new_particles       = new_particles     + new_particles_gb
430                sum_new_particles   = sum_new_particles + new_particles_gb                   
431             ENDDO
432          ENDDO
433       ENDDO
434
435    ELSEIF ( i_splitting_mode == 3 )  THEN
436
437       DO  i = nxl, nxr
438          DO  j = nys, nyn
439             DO  k = nzb+1, nzt
440             
441!
442!--             Initialize summing variables.             
443                lwc = 0.0_wp
444                m1  = 0.0_wp
445                m2  = 0.0_wp
446                m3  = 0.0_wp
447                nr  = 0.0_wp
448                rm  = 0.0_wp 
449               
450                new_particles_gb = 0
451                number_of_particles = prt_count(k,j,i)
452                IF ( number_of_particles <= 0  .OR.                            & 
453                     ql(k,j,i) < ql_crit )  CYCLE
454                particles => grid_particles(k,j,i)%particles
455!               
456!--             Calculate moments of DSD.               
457                DO  n = 1, number_of_particles
458                   IF ( particles(n)%particle_mask  .AND.                      &
459                        particles(n)%radius >= r_min )                         &
460                   THEN
461                      nr  = nr + particles(n)%weight_factor
462                      rm  = rm + factor_volume_to_mass  *                      &
463                                 particles(n)%radius**3  *                     &
464                                 particles(n)%weight_factor
465                      IF ( isf == 1 )  THEN           
466                         diameter   = particles(n)%radius * 2.0_wp           
467                         lwc = lwc + factor_volume_to_mass *                   &
468                                     particles(n)%radius**3 *                  & 
469                                     particles(n)%weight_factor 
470                         m1  = m1 + particles(n)%weight_factor * diameter
471                         m2  = m2 + particles(n)%weight_factor * diameter**2
472                         m3  = m3 + particles(n)%weight_factor * diameter**3
473                      ENDIF     
474                   ENDIF                                           
475                ENDDO 
476
477                IF ( nr <= 0.0  .OR.  rm <= 0.0_wp )  CYCLE 
478!
479!--             Calculate mean volume averaged radius.               
480                rm = ( rm / ( nr * factor_volume_to_mass ) )**0.3333333_wp
481!
482!--             Check which function should be used to approximate the DSD.             
483                IF ( isf == 1 )  THEN
484!
485!--                Gamma size distribution to calculate 
486!--                critical weight_factor (e.g. Marshall + Palmer, 1948).
487                   zeta = m1 * m3 / m2**2
488                   mu   = MAX( ( ( 1.0_wp - zeta ) * 2.0_wp + 1.0_wp ) /       &
489                                ( zeta - 1.0_wp ), 0.0_wp                      &
490                             )   
491                   lambda = ( pirho_l * nr / lwc *                             &
492                              ( mu + 3.0_wp ) * ( mu + 2.0_wp ) *              &
493                              ( mu + 1.0_wp )                                  &                                 
494                            )**0.3333333_wp
495                   nr0 =  ( nr / (gamma( mu + 1.0_wp ) ) ) *                   &
496                          lambda**( mu + 1.0_wp ) 
497
498                   DO  n = 0, n_max-1
499                      diameter         = r_bin_mid(n) * 2.0_wp           
500                      an_spl(n) = nr0 * diameter**mu *                         &
501                                  EXP( -lambda * diameter ) *                  & 
502                                  ( r_bin(n+1) - r_bin(n) ) * 2.0_wp 
503                   ENDDO
504                ELSEIF ( isf == 2 )  THEN
505!
506!--                Lognormal size distribution to calculate critical
507!--                weight_factor (e.g. Levin, 1971, Bradley + Stow, 1974).   
508                   DO  n = 0, n_max-1
509                      an_spl(n) = nr / ( SQRT( 2.0_wp * pi ) *                 &
510                                              LOG(sigma_log) * r_bin_mid(n)    &
511                                        ) *                                    &
512                                  EXP( -( LOG( r_bin_mid(n) / rm )**2 ) /      &
513                                        ( 2.0_wp * LOG(sigma_log)**2 )         & 
514                                      ) *                                      & 
515                                  ( r_bin(n+1) - r_bin(n) )
516                   ENDDO
517                ELSEIF ( isf == 3 )  THEN
518!
519!--                Exponential size distribution to calculate critical
520!--                weight_factor (e.g. Berry + Reinhardt, 1974). 
521                   DO  n = 0, n_max-1
522                      an_spl(n) = 3.0_wp * nr * r_bin_mid(n)**2 / rm**3 *     &
523                                  EXP( - ( r_bin_mid(n)**3 / rm**3 ) ) *      &
524                                  ( r_bin(n+1) - r_bin(n) )
525                   ENDDO
526                ENDIF
527               
528!
529!--             Criterion to avoid super droplets with a weighting factor < 1.0.                                   
530                an_spl = MAX(an_spl, 1.0_wp)
531!               
532!--             Start splitting operations. Each particle is checked if it
533!--             fulfilled the splitting criterion's. In splitting mode 'gb_av'   
534!--             a critical radius (radius_split) and a splitting function must
535!--             be prescribed (see particles_par). The critical weighting factor
536!--             is calculated while appoximating a 'gamma', 'log' or 'exp'-
537!--             drop size distribution. In this mode a DSD is calculated for
538!--             every cloudy grid box. Super droplets which have a larger
539!--             radius and larger weighting factor are split into
540!--             'splitting_factor' super droplets. In this case the splitting 
541!--             factor is calculated of weighting factor of the super droplet 
542!--             and theapproximated number concentration for droplet of such
543!--             a size. Due to the splitting, the weighting factor of the 
544!--             super droplet and all created clones is reduced by the factor 
545!--             of 'splitting_facor'.
546                DO  n = 1, number_of_particles
547                   DO  np = 0, n_max-1
548                      IF ( r_bin(np) >= radius_split  .AND.                    &
549                           particles(n)%particle_mask  .AND.                   &
550                           particles(n)%radius >= r_bin(np)    .AND.           &
551                           particles(n)%radius < r_bin(np+1)   .AND.           &
552                           particles(n)%weight_factor >= an_spl(np) )          &
553                      THEN
554!
555!--                      Calculate splitting factor.
556                         splitting_factor =                                    & 
557                             MIN( INT( particles(n)%weight_factor /            &
558                                        an_spl(np)                             &
559                                     ), splitting_factor_max                   &
560                                )
561                         IF ( splitting_factor < 2 )  CYCLE
562
563!
564!--                      Calculate the new number of particles.                                                                                         
565                         new_size = prt_count(k,j,i) + splitting_factor - 1
566!
567!--                      Cycle if maximum number of particles per grid box
568!--                      is greater than the allowed maximum number.                                                 
569                         IF ( new_size >= max_number_particles_per_gridbox )   &
570                         CYCLE
571!
572!--                      Reallocate particle array if necessary.                         
573                         IF ( new_size > SIZE(particles) )  THEN
574                            CALL realloc_particles_array(i,j,k,new_size)
575                         ENDIF
576!
577!--                      Calculate new weighting factor.
578                         particles(n)%weight_factor = & 
579                            particles(n)%weight_factor / splitting_factor
580                         tmp_particle               = particles(n)
581                         old_size                   = prt_count(k,j,i)
582!
583!--                      Create splitting_factor-1 new particles.
584                         DO jpp = 1, splitting_factor-1
585                            grid_particles(k,j,i)%particles(jpp+old_size) =    &
586                               tmp_particle                 
587                         ENDDO
588!
589!--                      Save the new number of droplets for every grid box.
590                         prt_count(k,j,i)    = prt_count(k,j,i) +              &
591                                               splitting_factor - 1
592                         new_particles_gb    = new_particles_gb +              &
593                                               splitting_factor - 1                                       
594                      ENDIF
595                   ENDDO 
596                ENDDO
597
598                new_particles       = new_particles + new_particles_gb
599                sum_new_particles   = sum_new_particles + new_particles_gb                                                 
600             ENDDO
601          ENDDO
602       ENDDO
603    ENDIF
604       
605    CALL cpu_log( log_point_s(80), 'lpm_splitting', 'stop' )
606
607 END SUBROUTINE lpm_splitting
608 
Note: See TracBrowser for help on using the repository browser.