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

Last change on this file since 3036 was 2932, checked in by maronga, 7 years ago

renamed all Fortran NAMELISTS

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