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

Last change on this file since 2263 was 2263, checked in by schwenkel, 4 years ago

Implemented splitting and merging algorithm

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