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

Last change on this file since 2707 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

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