PLUME-MoM  1.0
Integralvolcanicplumemodel
 All Classes Files Functions Variables Pages
particles.f90
Go to the documentation of this file.
1 !********************************************************************************
3 !
10 !********************************************************************************
11 MODULE particles_module
12  !
13  USE moments_module, ONLY : n_mom , n_nodes
14 
15  IMPLICIT NONE
16 
17  INTEGER :: n_part
18 
20  REAL*8, ALLOCATABLE, DIMENSION(:) :: solid_partial_mass_fraction
21 
23  REAL*8, ALLOCATABLE, DIMENSION(:) :: solid_partial_volume_fraction
24 
26  REAL*8, ALLOCATABLE, DIMENSION(:) :: solid_mass_fraction
27 
29  REAL*8, ALLOCATABLE, DIMENSION(:) :: solid_volume_fraction
30 
32  REAL*8, ALLOCATABLE, DIMENSION(:,:) :: mom
33 
35  REAL*8, ALLOCATABLE, DIMENSION(:,:) :: set_mom
36 
38  REAL*8, ALLOCATABLE, DIMENSION(:,:) :: rhop_mom
39 
41  REAL*8, ALLOCATABLE, DIMENSION(:,:) :: set_rhop_mom
42 
44  REAL*8, ALLOCATABLE, DIMENSION(:,:) :: cp_rhop_mom
45 
47  REAL*8, ALLOCATABLE, DIMENSION(:,:) :: set_cp_rhop_mom
48 
50  REAL*8, ALLOCATABLE, DIMENSION(:,:) :: set_cp_mom
51 
53  REAL*8, ALLOCATABLE, DIMENSION(:,:) :: birth_term
54 
56  REAL*8, ALLOCATABLE, DIMENSION(:,:) :: death_term
57 
59  REAL*8 :: shape_factor
60 
62  REAL*8, ALLOCATABLE :: diam1(:)
63 
65  REAL*8, ALLOCATABLE :: rho1(:)
66 
68  REAL*8, ALLOCATABLE :: diam2(:)
69 
71  REAL*8, ALLOCATABLE :: rho2(:)
72 
73  REAL*8, ALLOCATABLE :: cp_part(:)
74 
76  REAL*8, DIMENSION(1:50,0:100) :: mom0
77 
82  CHARACTER*10 :: settling_model
83 
89  CHARACTER(LEN=20) :: distribution
90 
91  CHARACTER(LEN=20) :: distribution_variable
92 
96  LOGICAL :: aggregation
97 
98  SAVE
99 
100 CONTAINS
101 
102  !******************************************************************************
104  !
111  !******************************************************************************
112 
113  SUBROUTINE allocate_particles
114 
115  IMPLICIT NONE
116 
117  ALLOCATE ( solid_partial_mass_fraction(1:n_part) )
118  ALLOCATE ( solid_partial_volume_fraction(1:n_part) )
119  ALLOCATE ( solid_mass_fraction(1:n_part) )
120  ALLOCATE ( solid_volume_fraction(1:n_part) )
121 
122  ! Allocation of the arrays for the moments
123  ALLOCATE ( mom(1:n_part,0:n_mom-1) )
124  ALLOCATE ( set_mom(1:n_part,0:n_mom-1) )
125  ALLOCATE ( rhop_mom(1:n_part,0:n_mom-1) )
126  ALLOCATE ( set_rhop_mom(1:n_part,0:n_mom-1) )
127  ALLOCATE ( set_cp_rhop_mom(1:n_part,0:n_mom-1) )
128  ALLOCATE ( set_cp_mom(1:n_part,0:n_mom-1) )
129  ALLOCATE ( cp_rhop_mom(1:n_part,0:n_mom-1) )
130  ALLOCATE ( birth_term(1:n_part,0:n_mom-1) )
131  ALLOCATE ( death_term(1:n_part,0:n_mom-1) )
132 
133  ! Allocation of the parameters for the variable density
134  ALLOCATE ( diam1(n_part) )
135  ALLOCATE ( rho1(n_part) )
136  ALLOCATE ( diam2(n_part) )
137  ALLOCATE ( rho2(n_part) )
138 
139  ALLOCATE ( cp_part(n_part) )
140 
141  END SUBROUTINE allocate_particles
142 
143  !******************************************************************************
145  !
152  !******************************************************************************
153 
154  SUBROUTINE initialize_particles
155 
156  USE moments_module, ONLY : wheeler_algorithm , moments_correction
157  USE variables, ONLY : verbose_level
158 
159  IMPLICIT NONE
160 
161  REAL*8, DIMENSION(n_part,n_nodes) :: xi
162  REAL*8, DIMENSION(n_part,n_nodes) :: w
163 
164  INTEGER :: iter
165 
166  INTEGER :: i_part
167 
168  DO i_part=1,n_part
169 
170  mom(i_part,:) = mom0(i_part,0:n_mom-1)
171 
172  ! CALL moments_correction(mom(i_part,:),iter)
173 
174  IF ( distribution .EQ. 'constant' ) THEN
175 
176  CALL wheeler_algorithm( mom(i_part,0:1) , xi(i_part,:) , &
177  w(i_part,:) )
178 
179  ELSE
180 
181  CALL wheeler_algorithm( mom(i_part,:) , xi(i_part,:) , w(i_part,:) )
182 
183  END IF
184 
185  IF ( verbose_level .GE. 1 ) THEN
186 
187  WRITE(*,*) 'part ',i_part
188  WRITE(*,*) 'mom',mom(i_part,:)
189  WRITE(*,*) 'abscissas',xi(i_part,:)
190  WRITE(*,*) 'weights',w(i_part,:)
191  READ(*,*)
192 
193  END IF
194 
195 
196  END DO
197 
198  CALL eval_particles_moments( xi(1:n_part,:) , w(1:n_part,:) )
199 
200 
201  END SUBROUTINE initialize_particles
202 
203 
204  !******************************************************************************
206  !
216  !******************************************************************************
217 
218  FUNCTION particles_settling_velocity(i_part,diam_in)
219  !
220  USE meteo_module, ONLY : rho_atm , rho_atm0 , visc_atm
221 
222  USE variables, ONLY : gi , pi_g , verbose_level
223 
224  IMPLICIT NONE
225 
226  REAL*8 :: particles_settling_velocity
227 
228  INTEGER, INTENT(IN) :: i_part
229  REAL*8, INTENT(IN) :: diam_in
230 
231  REAL*8 :: diam
232 
233  REAL*8 :: rhop
234 
235  REAL*8 :: k1 , k2 , k3
236  REAL*8 :: cd
237 
239  REAL*8 :: a_cs
240 
242  REAL*8 :: cd_100 , cd_1000
243 
245  REAL*8 :: cd_interp
246 
248  REAL*8 :: us_100 , us_1000
249 
251  REAL*8 :: mass
252 
254  REAL*8 :: us , us_1 ,us_2
255 
257  REAL*8 :: rey1 , rey2
258 
260  REAL*8 :: c0 , c1 , c2
261 
263  REAL*8 :: sqrt_delta
264 
265  IF ( distribution_variable .EQ. 'mass_fraction' ) THEN
266 
267  diam = 1.d-3 * 2.d0 ** ( - diam_in )
268 
269  ELSE
270 
271  diam = diam_in
272 
273  END IF
274 
275  rhop = particles_density(i_part,diam_in)
276 
277  IF ( settling_model .EQ. 'textor' ) THEN
278 
279  ! Textor et al. 2006
280 
281  IF ( diam .LE. 1.d-4 ) THEN
282 
283  k1 = 1.19d5 ! (m^2 kg^-1 s^-1 )
284 
285  particles_settling_velocity = k1 * rhop * sqrt( rho_atm0 / rho_atm ) * &
286  ( 0.5d0 * diam )**2
287 
288  ELSEIF ( diam .LE. 1.d-3 ) THEN
289 
290  k2 = 8.d0 ! (m^3 kg^-1 s^-1 )
291 
292  particles_settling_velocity = k2 * rhop * sqrt( rho_atm0 / rho_atm ) * &
293  ( 0.5d0 * diam )
294 
295  ELSE
296 
297  k3 = 4.833d0 ! (m^2 kg^-0.5 s^-1 )
298  cd = 0.75d0
299 
300  particles_settling_velocity = k3 * sqrt( rhop / cd ) * sqrt( rho_atm0 &
301  / rho_atm ) * sqrt( 0.5d0 * diam )
302 
303  END IF
304 
305  ELSEIF ( settling_model .EQ. 'pfeiffer' ) THEN
306 
307  k1 = shape_factor**(-0.828)
308  k2 = 2.d0 * dsqrt( 1.07 - shape_factor )
309 
310  mass = rhop * 4.d0/3.d0 * pi_g * ( 0.5*diam )**3
311 
312  a_cs = pi_g * ( 0.5*diam )**2
313 
314  c0 = -2.d0 * diam * mass * gi
315  c1 = 24.d0 * visc_atm * k1 * a_cs
316  c2 = rho_atm * diam * k2 * a_cs
317 
318  sqrt_delta = sqrt( c1**2 - 4 * c0*c2 )
319 
320  us_1 = ( - c1 + sqrt_delta ) / ( 2 * c2 )
321  us_2 = ( - c1 - sqrt_delta ) / ( 2 * c2 )
322 
323 
324  cd_100 = 24.d0/100.d0 * k1 + k2
325  us_100 = sqrt( 2 * mass * gi / ( cd_100*rho_atm * a_cs ) )
326 
327  cd_1000 = 1.d0
328  us_1000 = sqrt( 2 * mass * gi / ( cd_1000*rho_atm * a_cs ) )
329 
330  rey1 = rho_atm * diam * us_1 / visc_atm
331  rey2 = rho_atm * diam * us_2 / visc_atm
332 
333  IF ( verbose_level .GE. 4 ) THEN
334 
335  WRITE(*,*) 'rho_atm , diam , Us_1 , visc_atm',rho_atm , diam , us_1 , &
336  visc_atm
337  WRITE(*,*) 'Rey1,Rey2',rey1,rey2
338  READ(*,*)
339 
340  END IF
341 
342  IF ( ( rey1 .GT. 0.d0 ) .AND. ( rey1 .LE. 100.d0 ) ) THEN
343 
344  ! For small Reynolds numbers the drag coefficient is given by Eq.8
345  ! of Pfeiffer et al. 2005 and the settling velocity is Us_1
346 
347  us = us_1
348 
349  ELSEIF ( ( rey1 .GT. 100.d0 ) .AND. ( rey1 .LE. 1000.d0 ) ) THEN
350 
351  ! For intermediate Reyonlds numbers, 100<Re<1000, the drag coefficient
352  ! is linearly interpolated between Cd_100 and Cd_1000
353 
354  cd_interp = cd_100 + ( rey1 - 100 ) / ( 1000 - 100 ) * &
355  ( cd_1000 - cd_100)
356  us = sqrt( 2 * mass * gi / ( cd_interp * rho_atm * a_cs ) )
357 
358  ELSEIF ( rey1 .GT. 1000.d0 ) THEN
359 
360  ! For large Reynolds numbers the drag coefficient is taken as Cd=1,
361  ! as in Pfeiffer et al. 2005 with the settling velocity is Us_1000
362 
363  us = us_1000
364 
365  END IF
366 
367  IF ( ( rey2 .GT. 0.d0 ) .AND. ( rey2 .LE. 100.d0 ) ) THEN
368 
369  us = us_2
370 
371  ELSEIF ( ( rey2 .GT. 100.d0 ) .AND. ( rey2 .LE. 1000.d0 ) ) THEN
372 
373  cd_interp = cd_100 + ( rey2 - 100 ) / ( 1000 - 100 ) * ( cd_1000 - cd_100)
374  us = sqrt( 2 * mass * gi / ( cd_interp * rho_atm * a_cs ) )
375 
376  ELSEIF ( rey2 .GT. 1000.d0 ) THEN
377 
378  us = us_1000
379 
380  END IF
381 
382  particles_settling_velocity = us
383 
384  ELSE
385 
386  WRITE(*,*) 'wrong settling model'
387  stop
388 
389  END IF
390 
391  END FUNCTION particles_settling_velocity
392 
393  !******************************************************************************
395  !
403  !******************************************************************************
404 
405  FUNCTION particles_heat_capacity(i_part,diam_in)
406  !
407  IMPLICIT NONE
408 
409  REAL*8 :: particles_heat_capacity
410  INTEGER, INTENT(IN) :: i_part
411  REAL*8, INTENT(IN) :: diam_in
412  REAL*8 :: diam
413 
414  IF ( distribution_variable .EQ. 'mass_fraction' ) THEN
415 
416  diam = 1.d-3 * 2.d0 ** ( - diam_in )
417 
418  ELSE
419 
420  diam = diam_in
421 
422  END IF
423 
424  particles_heat_capacity = cp_part(i_part)
425 
426  END FUNCTION particles_heat_capacity
427 
428  !******************************************************************************
430  !
438  !******************************************************************************
439 
440  FUNCTION particles_density(i_part,diam_in)
441  !
442  IMPLICIT NONE
443 
444  REAL*8 :: particles_density
445 
446  INTEGER, INTENT(IN) :: i_part
447  REAL*8, INTENT(IN) :: diam_in
448 
449  REAL*8 :: diam
450 
451  REAL*8 :: diam_phi , diam1_phi , diam2_phi
452 
453 
454  IF ( distribution_variable .EQ. 'mass_fraction' ) THEN
455 
456  diam = 1.d-3 * 2.d0 ** ( - diam_in )
457 
458  ELSE
459 
460  diam = diam_in
461 
462  END IF
463 
464  IF ( diam .LE. diam1(i_part) ) THEN
465 
466  particles_density = rho1(i_part)
467 
468  ELSEIF ( diam .LE. diam2(i_part) ) THEN
469 
470  diam_phi = -log(diam*1000)/log(2.d0)
471  diam1_phi = -log(diam1(i_part)*1000)/log(2.d0)
472  diam2_phi = -log(diam2(i_part)*1000)/log(2.d0)
473 
474  particles_density = rho1(i_part) + ( diam_phi - diam1_phi ) / &
475  ( diam2_phi - diam1_phi ) * ( rho2(i_part) - rho1(i_part) )
476 
477  ELSE
478 
479  particles_density = rho2(i_part)
480 
481  END IF
482 
483  RETURN
484 
485  END FUNCTION particles_density
486 
487  !******************************************************************************
489  !
498  !******************************************************************************
499 
500  FUNCTION particles_beta(i_part,diam_i,diam_j)
501  !
502  IMPLICIT NONE
503 
504  REAL*8 :: particles_beta
505 
506  INTEGER, INTENT(IN) :: i_part
507  REAL*8, INTENT(IN) :: diam_i
508  REAL*8, INTENT(IN) :: diam_j
509 
510  ! Diameters in meters
511  REAL*8 :: diam_im , diam_jm
512 
513 
514  IF ( distribution_variable .EQ. 'mass_fraction' ) THEN
515 
516  diam_im = 1.d-3 * 2.d0 ** ( - diam_i )
517  diam_jm = 1.d-3 * 2.d0 ** ( - diam_j )
518 
519  ELSE
520 
521  diam_im = diam_i
522  diam_jm = diam_j
523 
524  END IF
525 
526  particles_beta = ( diam_im + diam_jm ) ** 2 / ( diam_im + diam_jm )
527 
528  RETURN
529 
530  END FUNCTION particles_beta
531 
532 
533  !******************************************************************************
535  !
544  !******************************************************************************
545 
546  FUNCTION aggregation_kernel(i_part,diam1,diam2)
547 
548  IMPLICIT NONE
549 
550  REAL*8 :: aggregation_kernel
551 
552  INTEGER, INTENT(IN) :: i_part
553  REAL*8, INTENT(IN) :: diam1
554  REAL*8, INTENT(IN) :: diam2
555 
556  REAL*8 :: beta
557  REAL*8 :: alfa
558 
559  beta = collision_kernel(i_part,diam1,diam2)
560 
561  alfa = coalescence_efficiency(i_part,diam1,diam2)
562 
563  aggregation_kernel = beta * alfa
564 
565  END FUNCTION aggregation_kernel
566 
567  !******************************************************************************
569  !
576  !******************************************************************************
577 
578  FUNCTION collision_kernel(i_part,diam1,diam2)
579 
580  USE meteo_module, ONLY : visc_atm
581 
582  USE variables, ONLY: pi_g
583 
584  IMPLICIT NONE
585 
586  REAL*8 :: collision_kernel
587 
588  INTEGER, INTENT(IN) :: i_part
589  REAL*8, INTENT(IN) :: diam1
590  REAL*8, INTENT(IN) :: diam2
591 
593  REAL*8 :: beta_b
594 
596  REAL*8 :: beta_s
597 
599  REAL*8 :: beta_ds
600 
602  REAL*8 :: k_b
603 
605  REAL*8 :: vs_1 , vs_2
606 
608  REAL*8 :: e_coll
609 
611  REAL*8 :: epsilon
612 
614  REAL*8 :: gamma_s
615 
617  REAL*8 :: air_kin_viscosity
618 
619  REAL*8 :: tp
620 
621  k_b =1.3806488d-23
622 
623  beta_b = 2.d0 / 3.d0 * k_b * tp / visc_atm * ( diam1 + diam2 )**2 / ( diam1*diam2 )
624 
625  gamma_s = dsqrt( 1.3d0 * epsilon * air_kin_viscosity )
626 
627  e_coll = collision_efficiency(i_part,diam1,diam2)
628 
629  beta_s = 1.d0 / 6.d0 * gamma_s * ( diam1 + diam2 )**3
630 
631  vs_1 = particles_settling_velocity(i_part,diam1)
632 
633  vs_2 = particles_settling_velocity(i_part,diam2)
634 
635  beta_ds = pi_g / 4.d0 * ( diam1 + diam2 )**2 * abs( vs_2 - vs_1 )
636 
637  collision_kernel = beta_b + beta_s + beta_ds
638 
639  END FUNCTION collision_kernel
640 
641  !******************************************************************************
643  !
650  !******************************************************************************
651 
652  FUNCTION collision_efficiency(i_part,diam1,diam2)
653 
654  USE variables, ONLY : gi
655 
656  IMPLICIT NONE
657 
658  REAL*8 :: collision_efficiency
659 
660  INTEGER, INTENT(IN) :: i_part
661  REAL*8, INTENT(IN) :: diam1
662  REAL*8, INTENT(IN) :: diam2
663 
664  REAL*8 :: e_v , e_a
665 
666  REAL*8 :: re
667 
668  REAL*8 :: stokes
669 
670  REAL*8 :: kin_visc_air
671 
673  REAL*8 :: vs_1 , vs_2
674 
675  vs_1 = particles_settling_velocity(i_part,diam1)
676 
677  vs_2 = particles_settling_velocity(i_part,diam2)
678 
679  IF ( diam1 .GT. diam2 ) THEN
680 
681  re = diam1 * vs_1 / kin_visc_air
682 
683  stokes = 2.d0 * vs_2 * abs( vs_1 - vs_2 ) / diam1 * gi
684 
685  ELSE
686 
687  re = diam2 * vs_2 / kin_visc_air
688 
689  stokes = 2.d0 * vs_1 * abs( vs_2 - vs_1 ) / diam2 * gi
690 
691  END IF
692 
693 
694 
695  IF ( stokes > 1.214 ) THEN
696 
697  e_v = ( 1.d0 + ( 0.75 * log( 2.d0 * stokes ) / ( stokes - 1.214 ) ) )** &
698  ( -2.d0 )
699 
700  ELSE
701 
702  e_v = 0.d0
703 
704  END IF
705 
706 
707 
708  collision_efficiency = ( 60.d0 * e_v + e_a * re ) / ( 60.d0 * re )
709 
710  END FUNCTION collision_efficiency
711 
712 
713  !******************************************************************************
715  !
722  !******************************************************************************
723 
724  FUNCTION coalescence_efficiency(i_part,diam1,diam2)
725 
726  USE variables, ONLY: gi
727 
728  IMPLICIT NONE
729 
730  REAL*8 :: coalescence_efficiency
731 
732  INTEGER :: i_part
733  REAL*8, INTENT(IN) :: diam1
734  REAL*8, INTENT(IN) :: diam2
735 
737  REAL*8 :: stokes
738 
740  REAL*8 :: stokes_cr
741 
743  REAL*8 :: q
744 
746  REAL*8 :: vs_1 , vs_2
747 
748  REAL*8 :: tp
749 
750  IF ( tp .LE. 273 ) THEN
751 
752  coalescence_efficiency = 0.09d0
753 
754  ELSE
755 
756  vs_1 = particles_settling_velocity(i_part,diam1)
757 
758  vs_2 = particles_settling_velocity(i_part,diam2)
759 
760  IF ( diam1 .GT. diam2 ) THEN
761 
762  stokes = 2.d0 * vs_2 * abs( vs_1 - vs_2 ) / diam1 * gi
763 
764  ELSE
765 
766  stokes = 2.d0 * vs_1 * abs( vs_2 - vs_1 ) / diam2 * gi
767 
768  END IF
769 
770  stokes_cr = 1.3d0
771 
772  q = 0.8d0
773 
774  coalescence_efficiency = 1.d0 / ( 1.d0 + ( stokes / stokes_cr ) ) ** q
775 
776  END IF
777 
778  END FUNCTION coalescence_efficiency
779 
780 
781 
782  !******************************************************************************
784  !
792  !******************************************************************************
793 
794  SUBROUTINE eval_particles_moments( xi , w )
795 
796  ! external variables
797  USE variables, ONLY : verbose_level
798 
799  ! external procedures
800  USE meteo_module, ONLY : zmet
801 
802  IMPLICIT NONE
803 
804  REAL*8, DIMENSION(n_part,n_nodes), INTENT(IN) :: xi
805  REAL*8, DIMENSION(n_part,n_nodes), INTENT(IN) :: w
806 
807  REAL*8, DIMENSION(n_part,n_nodes) :: part_dens_array
808  REAL*8, DIMENSION(n_part,n_nodes) :: part_set_vel_array
809  REAL*8, DIMENSION(n_part,n_nodes) :: part_cp_array
810  REAL*8, DIMENSION(n_part,n_nodes,n_nodes) :: part_beta_array
811 
812  INTEGER :: i , j , j1 , j2
813  INTEGER :: i_part
814 
815  CALL zmet
816 
817  DO i_part=1,n_part
818 
819  DO j=1,n_nodes
820 
821  part_dens_array(i_part,j) = particles_density( i_part , xi(i_part,j) )
822 
823  part_set_vel_array(i_part,j) = particles_settling_velocity( i_part , &
824  xi(i_part,j) )
825 
826  part_cp_array(i_part,j) = particles_heat_capacity( i_part,xi(i_part,j))
827 
828  IF ( aggregation) THEN
829 
830  DO j2=1,n_nodes
831 
832  part_beta_array(i_part,j,j2) = particles_beta( i_part , &
833  xi(i_part,j) , xi(i_part,j2) )
834 
835  END DO
836 
837  END IF
838 
839  END DO
840 
841  IF ( verbose_level .GE. 2 ) THEN
842 
843  WRITE(*,*) 'i_part',i_part
844  WRITE(*,*) 'abscissas', xi(i_part,1:n_nodes)
845  WRITE(*,*) 'weights', w(i_part,1:n_nodes)
846  WRITE(*,*) 'part_dens_array',part_dens_array(i_part,:)
847  WRITE(*,*) 'part_set_vel_array',part_set_vel_array(i_part,:)
848  WRITE(*,*) 'part_cp_array',part_cp_array(i_part,:)
849 
850  END IF
851 
852  END DO
853 
854 
855  DO i_part=1,n_part
856 
857  DO i=0,n_mom-1
858 
859  set_mom(i_part,i) = sum( part_set_vel_array(i_part,:) * w(i_part,:) &
860  * xi(i_part,:)**i ) / mom(i_part,i)
861 
862  rhop_mom(i_part,i) = sum( part_dens_array(i_part,:) * w(i_part,:) &
863  * xi(i_part,:)**i ) / mom(i_part,i)
864 
865  cp_rhop_mom(i_part,i) = sum( part_cp_array(i_part,:) &
866  * part_dens_array(i_part,:) * w(i_part,:) * xi(i_part,:)**i ) &
867  / mom(i_part,i)
868 
869  set_rhop_mom(i_part,i) = sum( part_set_vel_array(i_part,:) &
870  * part_dens_array(i_part,:) * w(i_part,:) * xi(i_part,:)**i ) &
871  / mom(i_part,i)
872 
873  set_cp_rhop_mom(i_part,i) = sum( part_set_vel_array(i_part,:) &
874  * part_cp_array(i_part,:) * part_dens_array(i_part,:) * &
875  w(i_part,:) * xi(i_part,:)**i ) / mom(i_part,i)
876 
877  set_cp_mom(i_part,i) = sum( part_set_vel_array(i_part,:) &
878  * part_cp_array(i_part,:) * w(i_part,:) * xi(i_part,:)**i ) &
879  / mom(i_part,i)
880 
881  IF ( aggregation ) THEN
882 
883  birth_term(i_part,i) = 0.d0
884  death_term(i_part,i) = 0.d0
885 
886  DO j1=1,n_nodes
887 
888  DO j2=1,n_nodes
889 
890  birth_term(i_part,i) = birth_term(i_part,i) + w(i_part,j1) &
891  * w(i_part,j2) * part_beta_array(i_part,j1,j2) &
892  *( xi(i_part,j1)**3 + xi(i_part,j2)**3 ) ** ( i / 3.d0 )
893 
894  death_term(i_part,i) = death_term(i_part,i) - w(i_part,j1) &
895  * xi(i_part,j1) * part_beta_array(i_part,j1,j2) &
896  * w(i_part,j1) * w(i_part,j2)
897 
898  END DO
899 
900  END DO
901 
902  birth_term(i_part,i) = 0.5d0 * birth_term(i_part,i)
903 
904  END IF
905 
906  END DO
907 
908  END DO
909 
910  RETURN
911 
912  END SUBROUTINE eval_particles_moments
913 
914 END MODULE particles_module
915 
particles_module::particles_density
real *8 function particles_density(i_part, diam_in)
Particle density.
Definition: particles.f90:440
meteo_module
Meteo module.
Definition: meteo.f90:11
particles_module::particles_beta
real *8 function particles_beta(i_part, diam_i, diam_j)
Brownian aggregation.
Definition: particles.f90:500
meteo_module::zmet
subroutine zmet
Meteo parameters.
Definition: meteo.f90:131
variables
Global variables.
Definition: variables.f90:10
particles_module::initialize_particles
subroutine initialize_particles
Particles variables inizialization.
Definition: particles.f90:154
particles_module::particles_heat_capacity
real *8 function particles_heat_capacity(i_part, diam_in)
Heat capacity.
Definition: particles.f90:405
particles_module::collision_kernel
real *8 function collision_kernel(i_part, diam1, diam2)
Collision kernel.
Definition: particles.f90:578
particles_module::coalescence_efficiency
real *8 function coalescence_efficiency(i_part, diam1, diam2)
Coalescence efficiency.
Definition: particles.f90:724
moments_module::moments_correction
subroutine moments_correction(m, iter)
Moments correction.
Definition: moments.f90:171
particles_module::aggregation_kernel
real *8 function aggregation_kernel(i_part, diam1, diam2)
Aggregation kernel.
Definition: particles.f90:546
particles_module::allocate_particles
subroutine allocate_particles
Particles variables inizialization.
Definition: particles.f90:113
particles_module::particles_settling_velocity
real *8 function particles_settling_velocity(i_part, diam_in)
Settling velocity.
Definition: particles.f90:218
moments_module
Method of Moments module.
Definition: moments.f90:11
particles_module::collision_efficiency
real *8 function collision_efficiency(i_part, diam1, diam2)
Collision efficiency.
Definition: particles.f90:652
particles_module::eval_particles_moments
subroutine eval_particles_moments(xi, w)
Particles moments computation.
Definition: particles.f90:794
particles_module
Particles module.
Definition: particles.f90:11
moments_module::wheeler_algorithm
subroutine wheeler_algorithm(m, xi, w)
Wheeler algorithm.
Definition: moments.f90:41
Generated on Wed Jun 1 2016 10:07:08 for PLUME-MoM by   doxygen 1.8.7