source: trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_haze.F90 @ 3590

Last change on this file since 3590 was 3560, checked in by debatzbr, 5 weeks ago

Addition of the microphysics model in moments.
File size: 38.9 KB
Line 
1MODULE MP2M_HAZE
2    !============================================================================
3    !
4    !     Purpose
5    !     -------
6    !     Haze microphysics module.
7    !
8    !     This module contains all definitions of the microphysics processes related to
9    !     aerosols (production, coagulation, sedimentation).
10    !
11    !     @Warning
12    !     The tendencies returned by the method are always defined over the vertical grid
13    !     from __TOP__ to __GROUND__.
14    !
15    !     The module also contains sixteen method:
16    !        - mm_haze_microphysics
17    !        - mm_haze_production
18    !        - mm_haze_sedimentation
19    !           * get_weff
20    !           * let_me_fall_in_peace
21    !        - mm_haze_coagulation
22    !           * g0ssco
23    !           * g0sfco
24    !           * g0ffco
25    !           * g3ssco
26    !           * g3sfco
27    !           * g0ssfm
28    !           * g0sffm
29    !           * g0fffm
30    !           * g3ssfm
31    !           * g3sffm
32    !
33    !     Authors
34    !     -------
35    !     B. de Batz de Trenquelléon, J. Burgalat (11/2024)
36    !
37    !============================================================================
38
39    USE MP2M_MPREC
40    USE MP2M_GLOBALS
41    USE MP2M_METHODS
42    IMPLICIT NONE
43
44    PRIVATE
45
46    PUBLIC :: mm_haze_microphysics,  &
47              mm_haze_production,    &
48              mm_haze_sedimentation, &
49              mm_haze_coagulation
50
51    CONTAINS
52 
53    !============================================================================
54    ! HAZE MICROPHYSICS INTERFACE SUBROUTINE
55    !============================================================================
56 
57    SUBROUTINE mm_haze_microphysics(m3a_s_prod,dm0a_s,dm3a_s,dm0a_f,dm3a_f)
58      !! Get the evolution of moments tracers through haze microphysics processes.
59      !!
60      !! The subroutine is a wrapper to the haze microphysics methods. It computes the tendencies
61      !! of moments tracers for production, sedimentation, and coagulation processes for the
62      !! atmospheric column.
63      !!
64
65      ! Production of the 3rd order moment of the spherical mode distribution from CH4 photolysis (m3.m-3).
66      REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: m3a_s_prod
67      ! Tendency of the 0th order moment of the spherical mode distribution (m-3).
68      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0a_s
69      ! Tendency of the 3rd order moment of the spherical mode distribution (m3.m-3).
70      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3a_s
71      ! Tendency of the 0th order moment of the fractal mode distribution (m-3).
72      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0a_f
73      ! Tendency of the 3rd order moment of the fractal mode distribution (m3.m-3).
74      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3a_f
75     
76      ! Local variables:
77      !~~~~~~~~~~~~~~~~~
78      REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm0as
79      REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm3as
80      REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm0af
81      REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm3af
82 
83      ! Initialization:
84      !~~~~~~~~~~~~~~~~
85      dm0a_s = 0._mm_wp
86      dm3a_s = 0._mm_wp
87      dm0a_f = 0._mm_wp
88      dm3a_f = 0._mm_wp
89 
90      ALLOCATE(zdm0as(mm_nla),zdm3as(mm_nla),zdm0af(mm_nla),zdm3af(mm_nla))
91      zdm0as(1:mm_nla) = 0._mm_wp
92      zdm3as(1:mm_nla) = 0._mm_wp
93      zdm0af(1:mm_nla) = 0._mm_wp
94      zdm3af(1:mm_nla) = 0._mm_wp
95 
96
97      ! Microphysical processes:
98      !~~~~~~~~~~~~~~~~~~~~~~~~~
99     
100      ! Coagulation
101      IF (mm_w_haze_coag) THEN
102        call mm_haze_coagulation(dm0a_s,dm3a_s,dm0a_f,dm3a_f)
103      ENDIF
104 
105      ! Sedimentation
106      IF (mm_w_haze_sed) THEN
107        call mm_haze_sedimentation(zdm0as,zdm3as,zdm0af,zdm3af)
108 
109        ! Computes precipitations
110        mm_aers_prec = SUM(zdm3as*mm_dzlev*mm_rhoaer)
111        mm_aerf_prec = SUM(zdm3af*mm_dzlev*mm_rhoaer)
112 
113        ! Updates tendencies
114        dm0a_s = dm0a_s + zdm0as
115        dm3a_s = dm3a_s + zdm3as
116        dm0a_f = dm0a_f + zdm0af
117        dm3a_f = dm3a_f + zdm3af
118      ENDIF
119 
120      ! Production
121      IF (mm_w_haze_prod) THEN
122        ! Production by photolysis of CH4
123        IF (mm_haze_prod_pCH4) THEN
124          dm0a_s = dm0a_s + (m3a_s_prod / (mm_rc_prod**3 * mm_alpha_s(3._mm_wp)))
125          dm3a_s = dm3a_s + m3a_s_prod
126        ELSE
127          call mm_haze_production(zdm0as,zdm3as)
128          dm0a_s = dm0a_s + zdm0as
129          dm3a_s = dm3a_s + zdm3as
130        ENDIF
131      ENDIF
132 
133      RETURN
134    END SUBROUTINE mm_haze_microphysics
135
136
137    !============================================================================
138    ! PRODUCTION PROCESS RELATED METHOD
139    !============================================================================
140   
141    SUBROUTINE mm_haze_production(dm0s,dm3s)
142      !! Compute the production of aerosols.
143      !!
144      !! @warning
145      !! Spherical aerosols are created at one pressure level. No other source is taken into account.
146      !! This process is controled by two parameters, the pressure level of production and the production rate.
147      !! Then both M0 and M3 of the spherical aerosols distribution are updated in the production zone by addition
148      !! of the production rate along a gaussian shape.
149      !!
150      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0s !! Tendency of M0 (m-3).
151      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3s !! Tendency of M3 (m3.m-3).
152     
153      ! Local variables:
154      !~~~~~~~~~~~~~~~~~
155      INTEGER                     :: i
156      REAL(kind=mm_wp)            :: zprod
157      REAL(kind=mm_wp), PARAMETER :: sigz  = 20e3_mm_wp, &
158                                     znorm = dsqrt(2._mm_wp)*sigz
159     
160      ! Locates production altitude
161      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
162      zprod = -1._mm_wp
163      DO i=1,mm_nla-1
164        IF (mm_plev(i) < mm_p_prod.AND.mm_plev(i+1) >= mm_p_prod) THEN
165          zprod = mm_zlay(i)
166          EXIT
167        ENDIF
168      ENDDO
169     
170      ! Sanity check
171      !~~~~~~~~~~~~~
172      IF (zprod < 0._mm_wp) THEN
173        WRITE(*,'(a)') "cannot find aerosols production altitude"
174        call EXIT(11)
175      ENDIF
176 
177      ! Computes production rate
178      !~~~~~~~~~~~~~~~~~~~~~~~~~
179      dm3s(:) = mm_tx_prod  * (0.75_mm_wp * mm_dt) / (mm_pi * mm_rhoaer) / (2._mm_wp * dsqrt(2._mm_wp) * mm_dzlev(1:mm_nla)) * &
180               (erf((mm_zlev(1:mm_nla)-zprod)/znorm) - erf((mm_zlev(2:mm_nla+1)-zprod)/znorm))
181      dm0s(:) = dm3s(:) / (mm_rc_prod**3*mm_alpha_s(3._mm_wp))
182    END SUBROUTINE mm_haze_production
183
184
185    !============================================================================
186    ! SEDIMENTATION PROCESS RELATED METHODS
187    !============================================================================
188 
189    SUBROUTINE mm_haze_sedimentation(dm0s,dm3s,dm0f,dm3f)
190      !! Interface to sedimentation algorithm.
191      !!
192      !! The subroutine computes the evolution of each moment of the aerosols tracers
193      !! through sedimentation process and returns their tendencies for a timestep.
194      !!
195      ! Tendency of the 0th order moment of the spherical mode (m-3).
196      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0s
197      ! Tendency of the 3rd order moment of the spherical mode (m3.m-3).
198      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3s
199      ! Tendency of the 0th order moment of the fractal mode (m-3).
200      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0f
201      ! Tendency of the 3rd order moment of the fractal mode (m3.m-3).
202      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3f
203
204      ! Local variables:
205      !~~~~~~~~~~~~~~~~~
206      REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: wth
207 
208      ALLOCATE(wth(mm_nle))
209 
210      ! Force settling velocity to M0
211      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
212      IF (mm_wsed_m0) THEN
213        ! Spherical particles
214        call get_weff(0._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth)
215       
216        call let_me_fall_in_peace(mm_m0aer_s,-wth,mm_dt,dm0s)
217        call let_me_fall_in_peace(mm_m3aer_s,-wth,mm_dt,dm3s)
218
219        ! Get settling velocity and mass flux
220        mm_m0as_vsed(:)  = wth(1:mm_nla)
221        mm_m3as_vsed(:)  = wth(1:mm_nla)
222        mm_aer_s_flux(:) = (mm_m3aer_s * mm_rhoaer) * wth(2:)
223       
224        ! Fractal particles
225        call get_weff(0._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth)
226       
227        call let_me_fall_in_peace(mm_m0aer_f,-wth,mm_dt,dm0f)
228        call let_me_fall_in_peace(mm_m3aer_f,-wth,mm_dt,dm3f)
229       
230        ! Get settling velocity and mass flux
231        mm_m0af_vsed(:)  = wth(1:mm_nla)
232        mm_m3af_vsed(:)  = wth(1:mm_nla)
233        mm_aer_f_flux(:) = (mm_m3aer_f * mm_rhoaer) * wth(2:)
234     
235      ! Force settling velocity to M3
236      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
237      ELSEIF (mm_wsed_m3) THEN
238        ! Spherical particles
239        call get_weff(3._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth)
240       
241        call let_me_fall_in_peace(mm_m3aer_s,-wth,mm_dt,dm3s)
242        call let_me_fall_in_peace(mm_m0aer_s,-wth,mm_dt,dm0s)
243       
244        ! Get settling velocity and mass flux
245        mm_m3as_vsed(:)  = wth(1:mm_nla)
246        mm_m0as_vsed(:)  = wth(1:mm_nla)
247        mm_aer_s_flux(:) = (mm_m3aer_s * mm_rhoaer) * wth(2:)
248       
249        ! Fractal particles
250        call get_weff(3._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth)
251       
252        call let_me_fall_in_peace(mm_m3aer_f,-wth,mm_dt,dm3f)
253        call let_me_fall_in_peace(mm_m0aer_f,-wth,mm_dt,dm0f)
254
255        ! Get settling velocity and mass flux
256        mm_m3af_vsed(:) = wth(1:mm_nla)
257        mm_m0af_vsed(:) = wth(1:mm_nla)
258        mm_aer_f_flux(:) = (mm_m3aer_f * mm_rhoaer) * wth(2:)
259     
260      ! No forcing of settling velocity (might be a problem...)
261      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
262      ELSE
263        ! Spherical particles
264        call get_weff(0._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth)
265        call let_me_fall_in_peace(mm_m0aer_s,-wth,mm_dt,dm0s)
266       
267        ! Get settling velocity
268        mm_m0as_vsed(:) = wth(1:mm_nla)
269       
270        call get_weff(3._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth)
271        call let_me_fall_in_peace(mm_m3aer_s,-wth,mm_dt,dm3s)
272       
273        ! Get settling velocity
274        mm_m3as_vsed(:) = wth(1:mm_nla)
275       
276        ! Get mass flux
277        mm_aer_s_flux(:) = (mm_m3aer_s * mm_rhoaer) * wth(2:)
278       
279        ! Fractal aerosols
280        call get_weff(0._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth)
281        call let_me_fall_in_peace(mm_m0aer_f,-wth,mm_dt,dm0f)
282       
283        ! Get settling velocity
284        mm_m0af_vsed(:) = wth(1:mm_nla)
285       
286        call get_weff(3._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth)
287        call let_me_fall_in_peace(mm_m3aer_f,-wth,mm_dt,dm3f)
288       
289        ! Get settling velocity
290        mm_m3af_vsed(:) = wth(1:mm_nla)
291       
292        ! Get mass flux
293        mm_aer_f_flux(:) = (mm_m3aer_f * mm_rhoaer) * wth(2:)
294      ENDIF
295
296      DEALLOCATE(wth)
297     
298      RETURN
299    END SUBROUTINE mm_haze_sedimentation
300
301
302    SUBROUTINE get_weff(k,df,rc,dt,afun,wth)
303      !! Get the effective settling velocity for aerosols moments.
304      !!
305      !! This method computes the effective settling velocity of the kth order moment of aerosol tracers.
306      !! The basic settling velocity (weff(Mk)) is computed using the following equation:
307      !!     Phi_sed (Mk) = \int_0^infty (n(r) . r**k . w(r) . dr)
308      !!                  = Mk . weff (Mk)
309      !!
310      ! The order of the moment.
311      REAL(kind=mm_wp), INTENT(in)                     :: k
312      ! Fractal dimension of the aersols.
313      REAL(kind=mm_wp), INTENT(in)                     :: df
314      ! Characteristic radius associated to the moment at each layer (m).
315      REAL(kind=mm_wp), INTENT(in), DIMENSION(mm_nla)  :: rc
316      ! Time step (s).
317      REAL(kind=mm_wp), INTENT(in)                     :: dt
318      ! Theoretical effective settling velocity at each vertical __levels__ (m.s-1).
319      REAL(kind=mm_wp), INTENT(out), DIMENSION(mm_nle) :: wth
320     
321      ! Inter-moment relation function.
322      INTERFACE
323        FUNCTION afun(order)
324          IMPORT mm_wp
325          REAL(kind=mm_wp), INTENT(in) :: order ! Order of the moment.
326          REAL(kind=mm_wp)             :: afun  ! Alpha value.
327        END FUNCTION afun
328      END INTERFACE
329     
330      ! Local variables:
331      !~~~~~~~~~~~~~~~~~
332      INTEGER          :: i
333      REAL(kind=mm_wp) :: af1,af2,ar1,ar2
334      REAL(kind=mm_wp) :: mrcoef, brhoair, thermal_w
335      REAL(kind=mm_wp) :: csto,cslf
336      REAL(kind=mm_wp) :: rb2ra
337     
338      wth(:) = 0._mm_wp
339
340      ! Molecular's case
341      !~~~~~~~~~~~~~~~~~
342      mrcoef = 0.74_mm_wp
343
344      ar1 = (3._mm_wp*df - 6._mm_wp) / df
345      af1 = (df*(k+3._mm_wp) - 6._mm_wp) / df
346     
347      rb2ra = mm_rm**((6._mm_wp-2._mm_wp*df)/df)
348
349      DO i=1,mm_nle
350        thermal_w = sqrt((8._mm_wp * mm_kboltz * mm_btemp(i)) / (mm_pi * mm_air_mmas))
351     
352        IF (i == 1) THEN
353          IF (rc(i) <= 0._mm_wp) CYCLE
354          brhoair = mm_rhoair(i)
355          wth(i) = - (mrcoef * mm_effg(mm_zlev(i)) * (mm_rhoaer/brhoair) / thermal_w) * rb2ra * &
356                   (afun(af1)/afun(k)) * rc(i)**ar1
357        ELSEIF (i == mm_nle) THEN
358          IF (rc(i-1) <= 0._mm_wp) CYCLE
359          brhoair = mm_rhoair(mm_nla) + (mm_rhoair(mm_nla) - mm_rhoair(mm_nla-1)) / 2._mm_wp
360          wth(i) = - (mrcoef * mm_effg(mm_zlev(i)) * (mm_rhoaer/brhoair) / thermal_w) * rb2ra * &
361                   (afun(af1)/afun(k)) * rc(i-1)**ar1
362        ELSE
363          IF (rc(i-1) <= 0._mm_wp) CYCLE
364          brhoair = (mm_rhoair(i-1) + mm_rhoair(i)) / 2._mm_wp
365          wth(i) = - (mrcoef * mm_effg(mm_zlev(i)) * (mm_rhoaer/brhoair) / thermal_w) * rb2ra * &
366                   (afun(af1)/afun(k)) * rc(i-1)**ar1
367        ENDIF
368      ENDDO
369
370      ! Hydrodynamical's case
371      ! In fact: transitory case which is the general case
372      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
373      !ar1 = (3._mm_wp*df - 3._mm_wp) / df
374      !ar2 = (3._mm_wp*df - 6._mm_wp) / df
375      !
376      !af1 = (df*(k+3._mm_wp) - 3._mm_wp) / df
377      !af2 = (df*(k+3._mm_wp) - 6._mm_wp) / df
378      !
379      !rb2ra = mm_rm**((df-3._mm_wp)/df)
380      !
381      !DO i=1,mm_nle
382      !  csto = (2._mm_wp*mm_rhoaer*mm_effg(mm_zlev(i))) / (9._mm_wp*mm_eta_air(mm_btemp(i)))
383      !  cslf = mm_akn * mm_lambda_air(mm_btemp(i),mm_plev(i)) / rb2ra
384      ! 
385      !  IF (i == 1) THEN
386      !    IF (rc(i) <= 0._mm_wp) CYCLE
387      !    wth(i) = - csto/(rb2ra*afun(k)) * &
388      !             ((afun(af1)*rc(i)**ar1) + (cslf*afun(af2)*rc(i)**ar2))
389      !  ELSE
390      !    IF (rc(i-1) <= 0._mm_wp) CYCLE
391      !    wth(i) = - csto/(rb2ra*afun(k)) * &
392      !             ((afun(af1)*rc(i-1)**ar1) + (cslf*afun(af2)*rc(i-1)**ar2))
393      !  ENDIF
394      !ENDDO
395
396      RETURN
397    END SUBROUTINE get_weff
398 
399   
400    SUBROUTINE let_me_fall_in_peace(mk,ft,dt,dmk)
401      !! Compute the tendency of the moment through sedimentation process.
402      !!
403      !! The method computes the time evolution of the kth order moment through sedimentation:
404      !! dM(k) / dt = Phi(k) / dz
405      !!
406      !! The equation is resolved using a [Crank-Nicolson algorithm](http://en.wikipedia.org/wiki/Crank-Nicolson_method).
407      !!
408      !! Sedimentation algorithm is quite messy. It appeals to the dark side of the Force and uses evil black magic spells
409      !! from ancient times. It is based on \cite{toon1988b,fiadeiro1977,turco1979} and is an update of the algorithm
410      !! originally implemented in the LMDZ-Titan 2D GCM.
411      !!
412      REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: mk  !! kth order moment to sediment (m^k.m^-3).
413      REAL(kind=mm_wp), INTENT(in)                :: dt  !! Time step (s).
414      REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: ft  !! Downward sedimentation flux (effective velocity of the moment).
415      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dmk !! Tendency of kth order moment (m^k.m^-3).
416     
417      ! Local variables:
418      !~~~~~~~~~~~~~~~~~
419      INTEGER :: i
420      REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: as,bs,cs,mko
421     
422      ALLOCATE(as(mm_nla), bs(mm_nla), cs(mm_nla), mko(mm_nla))
423
424      mko(:) = 0._mm_wp
425      cs(1:mm_nla) = ft(2:mm_nla+1) - mm_dzlay(1:mm_nla)/dt
426     
427      IF (ANY(cs > 0._mm_wp)) THEN
428        ! Implicit case
429        as(1:mm_nla) = ft(1:mm_nla)
430        bs(1:mm_nla) = -(ft(2:mm_nle) + mm_dzlay(1:mm_nla)/dt)
431        cs(1:mm_nla) = -mm_dzlay(1:mm_nla) / dt * mk(1:mm_nla)
432        ! (Tri)diagonal matrix inversion
433        mko(1) = cs(1) / bs(1)
434        DO i=2,mm_nla ; mko(i) = (cs(i)-mko(i-1)*as(i))/bs(i) ; ENDDO
435     
436      ELSE
437        ! Explicit case
438        as(1:mm_nla) = -mm_dzlay(1:mm_nla) / dt
439        bs(1:mm_nla) = -ft(1:mm_nla)
440        ! Boundaries
441        mko(1) = cs(1) * mk(1) / as(1)
442        mko(mm_nla) = (bs(mm_nla)*mk(mm_nla-1) + cs(mm_nla)*mk(mm_nla)) / as(mm_nla)
443        ! Interior
444        mko(2:mm_nla-1) = (bs(2:mm_nla-1)*mk(1:mm_nla-2) + cs(2:mm_nla-1)*mk(2:mm_nla-1)) &
445                          /as(2:mm_nla-1)
446      ENDIF
447     
448      dmk = mko - mk
449      DEALLOCATE(as,bs,cs,mko)
450     
451      RETURN
452    END SUBROUTINE let_me_fall_in_peace
453 
454   
455    !============================================================================
456    ! COAGULATION PROCESS RELATED METHODS
457    !============================================================================
458 
459    SUBROUTINE mm_haze_coagulation(dM0s,dM3s,dM0f,dM3f)
460      !! Get the evolution of the aerosols moments vertical column due to coagulation process.
461      !!
462      !! This is the main method of the coagulation process:
463      !!    1. Computes gamma pre-factor for each parts of the coagulation equation(s).
464      !!    2. Applies the electic correction on the gamma pre-factor.
465      !!    3. Computes the specific flow regime kernels (molecular or hydrodynamic).
466      !!    4. Computes the harmonic mean of the kernels (transitory regime).
467      !!    5. Finally computes the tendencies of the moments.
468      !!
469      !! @Warning
470      !! If the transfert probabilities are set to 1 for the two flow regimes (pco and pfm),
471      !! a floating point exception occured (i.e. a NaN) as we perform a division by zero
472      !!
473      ! Tendency of the 0th order moment of the spherical size-distribution over a time step (m-3).
474      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM0s
475      ! Tendency of the 3rd order moment of the spherical size-distribution (m3.m-3).
476      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM3s
477      ! Tendency of the 0th order moment of the fractal size-distribution over a time step (m-3).
478      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM0f
479      ! Tendency of the 3rd order moment of the fractal size-distribution over a time step (m3.m-3).
480      REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM3f
481
482      ! Local variables:
483      !~~~~~~~~~~~~~~~~~
484      INTEGER :: i
485
486      REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: kco,kfm,slf
487      REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: g_co,g_fm,pco,pfm,mq,tmp
488      REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: a_ss,a_sf,b_ss,b_ff,c_ss,c_sf
489
490      ! Sanity check:
491      IF (mm_coag_choice < 0 .OR. mm_coag_choice > mm_coag_ss+mm_coag_sf+mm_coag_ff) &
492          STOP "invalid choice for coagulation mode interaction activation"
493 
494      ! Allocates local arrays:
495      ALLOCATE(kco(mm_nla),kfm(mm_nla),slf(mm_nla))
496      ALLOCATE(g_co(mm_nla),g_fm(mm_nla), &
497               pco(mm_nla),pfm(mm_nla),   &
498               mq(mm_nla), tmp(mm_nla))
499      ALLOCATE(a_ss(mm_nla),a_sf(mm_nla), &
500               b_ss(mm_nla),b_ff(mm_nla), &
501               c_ss(mm_nla),c_sf(mm_nla))
502 
503      a_ss(:) = 0._mm_wp ; a_sf(:) = 0._mm_wp
504      b_ss(:) = 0._mm_wp ; b_ff(:) = 0._mm_wp
505      c_ss(:) = 0._mm_wp ; c_sf(:) = 0._mm_wp
506 
507
508      ! 1. Gamma pre-factor:
509      !~~~~~~~~~~~~~~~~~~~~~
510      ! gets kco, kfm pre-factors
511      kco(:) = mm_get_kco(mm_temp) ; kfm(:) = mm_get_kfm(mm_temp)
512      ! get slf (slip-flow factor)
513      slf(:) = mm_akn * mm_lambda_air(mm_temp,mm_play)
514 
515
516      ! 2,3,4. Electic correction, Flow regime kernels, Harmonic mean:
517      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
518      DO i=1,mm_nla
519       
520        ! SS interactions
521        !~~~~~~~~~~~~~~~~
522        IF (mm_rcs(i) > mm_rcs_min .AND. IAND(mm_coag_choice,mm_coag_ss) /= 0) THEN
523          ! Probability S --> S for M0/CO and M0/FM
524          pco(i) = mm_ps2s(mm_rcs(i),0,0)
525          pfm(i) = mm_ps2s(mm_rcs(i),0,1)
526         
527          ! Gamma (kernel dependent)
528          g_co(i) = g0ssco(mm_rcs(i),slf(i),kco(i))
529          g_fm(i) = g0ssfm(mm_rcs(i),kfm(i))
530         
531          ! Molecular's case
532          !~~~~~~~~~~~~~~~~~
533          a_ss(i) = (pfm(i) - 2._mm_wp) * g_fm(i)
534          b_ss(i) = (1._mm_wp - pfm(i)) * g_fm(i)
535
536          ! Hydrodynamical's case
537          ! In fact: transitory case which is the general case
538          !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
539          ! (A_SS_CO x A_SS_FM) / (A_SS_CO + A_SS_FM)
540          !IF (g_co(i)*(pco(i)-2._mm_wp)+g_fm(i)*(pfm(i)-2._mm_wp) /= 0._mm_wp) THEN
541          !  a_ss(i) = (g_co(i)*(pco(i)-2._mm_wp)*g_fm(i)*(pfm(i)-2._mm_wp))/(g_co(i)*(pco(i)-2._mm_wp)+g_fm(i)*(pfm(i)-2._mm_wp))
542          !ENDIF
543          ! (B_SS_CO x B_SS_FM) / (B_SS_CO + B_SS_FM)
544          !IF (g_co(i)*(1._mm_wp-pco(i))+g_fm(i)*(1._mm_wp-pfm(i)) /= 0._mm_wp) THEN
545          !  b_ss(i) = (g_co(i)*(1._mm_wp-pco(i))*g_fm(i)*(1._mm_wp-pfm(i)))/(g_co(i)*(1._mm_wp-pco(i))+g_fm(i)*(1._mm_wp-pfm(i)))
546          !ENDIF
547
548          ! Eletric charge correction for M0/SS interactions
549          mq(i) = mm_qmean(mm_rcs(i),mm_rcs(i),0,'SS',mm_temp(i))
550          a_ss(i) = a_ss(i) * mq(i)
551          b_ss(i) = b_ss(i) * mq(i)
552         
553          ! Probability S --> S for M3/CO and M3/FM
554          pco(i) = mm_ps2s(mm_rcs(i),3,0)
555          pfm(i) = mm_ps2s(mm_rcs(i),3,1)
556
557          ! Gamma (kernel dependent)
558          g_co(i) = g3ssco(mm_rcs(i),slf(i),kco(i)) * (pco(i)-1._mm_wp)
559          g_fm(i) = g3ssfm(mm_rcs(i),kfm(i)) * (pfm(i)-1._mm_wp)
560         
561          ! Molecular's case
562          !~~~~~~~~~~~~~~~~~
563          c_ss(i) = g_fm(i)
564
565          ! Hydrodynamical's case
566          ! In fact: transitory case which is the general case
567          !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
568          ! (C_SS_CO x C_SS_FM) / (C_SS_CO + C_SS_FM)
569          !IF (g_co(i) + g_fm(i) /= 0._mm_wp) THEN
570          !  c_ss(i) = (g_co(i) * g_fm(i)) / (g_co(i) + g_fm(i))
571          !ENDIF
572
573          IF (b_ss(i) <= 0._mm_wp) c_ss(i) = 0._mm_wp
574         
575          ! Eletric charge correction for M3/SS interactions
576          mq(i) = mm_qmean(mm_rcs(i),mm_rcs(i),3,'SS',mm_temp(i))
577          c_ss(i) = c_ss(i) * mq(i)
578        ENDIF ! end SS interactions
579 
580        ! SF interactions
581        !~~~~~~~~~~~~~~~~
582        IF (mm_rcs(i) > mm_rcs_min .AND. mm_rcf(i) > mm_rcf_min .AND. IAND(mm_coag_choice,mm_coag_sf) /= 0) THEN
583          ! Gamma (kernel dependent)
584          g_co(i) = g0sfco(mm_rcs(i),mm_rcf(i),slf(i),kco(i))
585          g_fm(i) = g0sffm(mm_rcs(i),mm_rcf(i),kfm(i))
586         
587          ! Molecular's case
588          !~~~~~~~~~~~~~~~~~
589          a_sf(i) = g_fm(i)
590
591          ! Hydrodynamical's case
592          ! In fact: transitory case which is the general case
593          !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
594          ! (A_SF_CO x A_SF_FM) / (A_SF_CO + A_SF_FM)
595          !IF(g_co(i) + g_fm(i) /= 0._mm_wp) THEN
596          !  a_sf(i) = (g_co(i) * g_fm(i)) / (g_co(i) + g_fm(i))
597          !ENDIF
598
599          ! Eletric charge correction for M0/SF interactions
600          mq(i) = mm_qmean(mm_rcs(i),mm_rcf(i),0,'SF',mm_temp(i))
601          a_sf(i) = a_sf(i) * mq(i)
602         
603          ! Gamma (kernel dependent)
604          g_co(i) = g3sfco(mm_rcs(i),mm_rcf(i),slf(i),kco(i))
605          g_fm(i) = g3sffm(mm_rcs(i),mm_rcf(i),kfm(i))
606         
607          ! Molecular's case
608          !~~~~~~~~~~~~~~~~~
609          c_sf(i) = g_fm(i)
610
611          ! Hydrodynamical's case
612          ! In fact: transitory case which is the general case
613          !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
614          ! (C_SF_CO x C_SF_FM) / (C_SF_CO + C_SF_FM)
615          !IF (g_co(i) + g_fm(i) /= 0._mm_wp) THEN
616          !  c_sf(i) = (g_co(i) * g_fm(i)) / (g_co(i) + g_fm(i))
617          !ENDIF
618         
619          ! Eletric charge correction for M3/SF interactions
620          mq(i) = mm_qmean(mm_rcs(i),mm_rcf(i),3,'SF',mm_temp(i))
621          c_sf(i) = c_sf(i) * mq(i)
622        ENDIF ! end SF interactions
623       
624        ! FF interactions
625        !~~~~~~~~~~~~~~~~
626        IF(mm_rcf(i) > mm_rcf_min .AND. IAND(mm_coag_choice,mm_coag_sf) /= 0) THEN
627          ! Gamma (kernel dependent)
628          g_co(i) = g0ffco(mm_rcf(i),slf(i),kco(i))
629          g_fm(i) = g0fffm(mm_rcf(i),kfm(i))
630
631          ! Molecular's case
632          !~~~~~~~~~~~~~~~~~
633          b_ff(i) = g_fm(i)
634
635          ! Hydrodynamical's case
636          ! In fact: transitory case which is the general case
637          !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
638          ! (B_FF_CO x B_FF_FM) / (B_FF_CO + B_FF_FM)
639          !IF (g_co(i) + g_fm(i) /= 0._mm_wp) THEN
640          !  b_ff(i) = (g_co(i) * g_fm(i)) / (g_co(i) + g_fm(i))
641          !ENDIF
642         
643          ! Eletric charge correction for M0/FF interactions
644          mq(i) = mm_qmean(mm_rcf(i),mm_rcf(i),0,'FF',mm_temp(i))
645          b_ff(i) = b_ff(i) * mq(i)
646        ENDIF ! end FF interactions
647     
648      ENDDO ! end n_lay
649 
650      DEALLOCATE(g_co,g_fm,kco,kfm,pco,pfm,slf)
651 
652      ! 5. Tendencies of the moments:
653      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
654      ! Now we will use the kharm two by two to compute :
655      ! dM_0_S/mm_dt = kharm(1) * M_0_S^2 - kharm(2) * M_0_S * m_0_F
656      ! dM_0_F/mm_dt = kharm(3) * M_0_S^2 - kharm(4) * M_0_F^2
657      ! dM_3_S/mm_dt = kharm(5) * M_3_S^2 - kharm(6) * M_3_S * m_3_F
658      ! dM_3_F/mm_dt = - dM_3_S/mm_dt
659      !
660      ! We use a (semi) implicit scheme:
661      ! When X appears as square we set one X at t+1, the other a t.
662      ! --- dM0s
663      tmp(:) = mm_dt * (a_ss*mm_m0aer_s - a_sf*mm_m0aer_f)
664      dm0s(:) =  mm_m0aer_s * (tmp / (1._mm_wp - tmp))
665      ! --- dM0f
666      tmp(:) = mm_dt * b_ff * mm_m0aer_f
667      dm0f(:) = (b_ss*mm_dt*mm_m0aer_s**2 - tmp*mm_m0aer_f) / (1._mm_wp + tmp)
668      ! --- dM3s
669      tmp(:) = mm_dt * (c_ss*mm_m3aer_s - c_sf*mm_m3aer_f)
670      dm3s(:) =  mm_m3aer_s * (tmp / (1._mm_wp - tmp))
671      ! --- dM3f
672      dm3f(:) = -dm3s
673 
674      DEALLOCATE(a_ss,a_sf,b_ss,b_ff,c_ss,c_sf,tmp)
675 
676      RETURN
677    END SUBROUTINE mm_haze_coagulation
678 
679 
680    ! =====================
681    ! Hydrodynamical's case
682    ! =====================
683    ELEMENTAL FUNCTION g0ssco(rcs,slf,kco) RESULT(res)
684      !! Get gamma pre-factor for the 0th order moment with SS interactions in the Continuous flow regime.
685      !!
686      REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution.
687      REAL(kind=mm_wp), INTENT(in) :: slf !! Slip-Flow correction pre-factor.
688      REAL(kind=mm_wp), INTENT(in) :: kco !! Thermodynamic continuous flow regime pre-factor.
689      REAL(kind=mm_wp) :: res             !! Gamma coagulation kernel pre-factor.
690      REAL(kind=mm_wp) :: a1, a2, a3
691      res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN
692      ! Computes mm_alpha coefficients
693      a1 = mm_alpha_s(1._mm_wp)
694      a2 = mm_alpha_s(-1._mm_wp)
695      a3 = mm_alpha_s(-2._mm_wp)
696      ! Computes gamma pre-factor
697      res = kco * (1._mm_wp + a1*a2 + slf/rcs*(a2+a1*a3))
698      RETURN
699    END FUNCTION g0ssco
700 
701
702    ELEMENTAL FUNCTION g0sfco(rcs,rcf,slf,kco) RESULT(res)
703      !! Get gamma pre-factor for the 0th order moment with SF interactions in the Continuous flow regime.
704      !!
705      REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution.
706      REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution.
707      REAL(kind=mm_wp), INTENT(in) :: slf !! Slip-Flow correction pre-factor.
708      REAL(kind=mm_wp), INTENT(in) :: kco !! Thermodynamic continuous flow regime pre-factor.
709      REAL(kind=mm_wp) :: res             !! Gamma coagulation kernel pre-factor.
710      REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, e, rcff
711      res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN
712      ! Computes mm_alpha and other coefficients
713      e = 3._mm_wp/mm_df
714      rcff = rcf**e * mm_rb2ra
715      a1 = mm_alpha_s(1._mm_wp)  ; a2 = mm_alpha_f(-e)        ; a3 = mm_alpha_f(e)
716      a4 = mm_alpha_s(-1._mm_wp) ; a5 = mm_alpha_s(-2._mm_wp) ; a6 = mm_alpha_f(-2._mm_wp*e)
717      ! Computes gamma pre-factor
718      res = kco * (2._mm_wp + a1*a2*rcs/rcff + &
719                   a4*a3*rcff/rcs +            &
720                   slf * (a4/rcs + a2/rcff +   &
721                          a5*a3*rcff/rcs**2 +  &
722                          a1*a6*rcs/rcff**2))
723      RETURN
724    END FUNCTION g0sfco
725 
726
727    ELEMENTAL FUNCTION g0ffco(rcf,slf,kco) RESULT(res)
728      !! Get gamma pre-factor for the 0th order moment with FF interactions in the Continuous flow regime.
729      !!
730      REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution.
731      REAL(kind=mm_wp), INTENT(in) :: slf !! Slip-Flow correction pre-factor.
732      REAL(kind=mm_wp), INTENT(in) :: kco !! Thermodynamic continuous flow regime pre-factor.
733      REAL(kind=mm_wp) :: res             !! Gamma coagulation kernel pre-factor.
734      REAL(kind=mm_wp) :: a1, a2, a3, e, rcff
735      res = 0._mm_wp ; IF (rcf <= 0._mm_wp) RETURN
736      ! Computes mm_alpha and other coefficients
737      e = 3._mm_wp/mm_df
738      rcff = rcf**e * mm_rb2ra
739      a1 = mm_alpha_f(e) ; a2 = mm_alpha_f(-e) ; a3 = mm_alpha_f(-2._mm_wp*e)
740      ! Computes gamma pre-factor
741      res = kco * (1._mm_wp + a1*a2 + slf/rcff *(a2+a1*a3))
742      RETURN
743    END FUNCTION g0ffco
744
745
746    ELEMENTAL FUNCTION g3ssco(rcs, slf, kco) RESULT(res)
747      !! Get gamma pre-factor for the 3rd order moment with SS interactions in the Continuous flow regime.
748      !!
749      REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution.
750      REAL(kind=mm_wp), INTENT(in) :: slf !! Slip-Flow correction pre-factor.
751      REAL(kind=mm_wp), INTENT(in) :: kco !! Thermodynamic continuous flow regime pre-factor.
752      REAL(kind=mm_wp) :: res             !! Gamma coagulation kernel pre-factor.
753      REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6
754      res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN
755      ! Computes mm_alpha coefficients
756      a1 = mm_alpha_s(3._mm_wp) ; a2 = mm_alpha_s(2._mm_wp)  ; a3 = mm_alpha_s(1._mm_wp)
757      a4 = mm_alpha_s(4._mm_wp) ; a5 = mm_alpha_s(-1._mm_wp) ; a6 = mm_alpha_s(-2._mm_wp)
758      ! Computes gamma pre-factor
759      res = kco/(a1**2*rcs**3) * (2._mm_wp*a1 + a2*a3 + a4*a5 + &
760                                  slf/rcs*(a3**2 + a4*a6 + a1*a5 + a2))
761      RETURN
762    END FUNCTION g3ssco
763
764
765    ELEMENTAL FUNCTION g3sfco(rcs, rcf, slf, kco) RESULT(res)
766      !! Get gamma pre-factor for the 3rd order moment with SF interactions in the Continuous flow regime.
767      !!
768      REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution.
769      REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution.
770      REAL(kind=mm_wp), INTENT(in) :: slf !! Slip-Flow correction pre-factor.
771      REAL(kind=mm_wp), INTENT(in) :: kco !! Thermodynamic continuous flow regime pre-factor.
772      REAL(kind=mm_wp) :: res             !! Gamma coagulation kernel pre-factor.
773      REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, a7, a8, e, rcff
774      res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN
775      ! Computes mm_alpha and other coefficients
776      e = 3._mm_wp/mm_df
777      rcff = rcf**e * mm_rb2ra
778      a1 = mm_alpha_s(3._mm_wp)    ; a2 = mm_alpha_s(4._mm_wp) ; a3 = mm_alpha_f(-e)
779      a4 = mm_alpha_s(2._mm_wp)    ; a5 = mm_alpha_f(e)        ; a6 = mm_alpha_s(1._mm_wp)
780      a7 = mm_alpha_f(-2._mm_wp*e) ; a8 = mm_alpha_f(3._mm_wp)
781      ! Computes gamma pre-factor
782      res = kco/(a1*a8*(rcs*rcf)**3) * (2._mm_wp*a1*rcs**3 + a2*rcs**4*a3/rcff + a4*rcs**2*a5*rcff + &
783                                        slf * (a4*rcs**2 + a1*rcs**3*a3/rcff +                       &
784                                               a6*rcs*a5*rcff +                                      &
785                                               a2*rcs**4*a7/rcff**2))
786      RETURN
787    END FUNCTION g3sfco
788
789
790    ! ================
791    ! Molecular's case
792    ! ================
793    ELEMENTAL FUNCTION g0ssfm(rcs, kfm) RESULT(res)
794      !! Get gamma pre-factor for the 0th order moment with SS interactions in the Free Molecular flow regime.
795      !!
796      REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution.
797      REAL(kind=mm_wp), INTENT(in) :: kfm !! Thermodynamic free molecular flow regime pre-factor.
798      REAL(kind=mm_wp) :: res             !! Gamma coagulation kernel pre-factor.
799      REAL(kind=mm_wp) :: a1, a2, a3, a4, a5
800      res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN
801      ! Computes mm_alpha coefficients
802      a1 = mm_alpha_s(0.5_mm_wp) ; a2 = mm_alpha_s(1._mm_wp)   ; a3 = mm_alpha_s(-0.5_mm_wp)
803      a4 = mm_alpha_s(2._mm_wp)  ; a5 = mm_alpha_s(-1.5_mm_wp)
804      ! Computes gamma pre-factor
805      res = kfm * (rcs**0.5_mm_wp*mm_get_btk(1,0)) * (a1 + 2._mm_wp*a2*a3 + a4*a5)
806      RETURN
807    END FUNCTION g0ssfm
808 
809
810    ELEMENTAL FUNCTION g0sffm(rcs, rcf, kfm) RESULT(res)
811      !! Get gamma pre-factor for the 0th order moment with SF interactions in the Free Molecular flow regime.
812      !!
813      REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution.
814      REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution.
815      REAL(kind=mm_wp), INTENT(in) :: kfm !! Thermodynamic free molecular flow regime pre-factor.
816      REAL(kind=mm_wp) :: res             !! Gamma coagulation kernel pre-factor.
817      REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, a7, a8, a9, a10
818      REAL(kind=mm_wp) :: e1, e2, e3, e4
819      REAL(kind=mm_wp) :: rcff1, rcff2, rcff3, rcff4
820      res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN
821      ! Computes mm_alpha and other coefficients
822      e1 = 3._mm_wp/mm_df
823      e2 = 6._mm_wp/mm_df
824      e3 = (6._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
825      e4 = (12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
826 
827      rcff1 = rcf**e1 * mm_rb2ra
828      rcff2 = rcff1**2
829      rcff3 = rcf**e3 * mm_rb2ra
830      rcff4 = rcf**e4 * mm_rb2ra**2
831 
832      a1 = mm_alpha_s(0.5_mm_wp)  ; a2 = mm_alpha_s(-0.5_mm_wp) ; a3 = mm_alpha_f(e1)
833      a4 = mm_alpha_s(-1.5_mm_wp) ; a5 = mm_alpha_f(e2)         ; a6 = mm_alpha_s(2._mm_wp)
834      a7 = mm_alpha_f(-1.5_mm_wp) ; a8 = mm_alpha_s(1._mm_wp)   ; a9 = mm_alpha_f(e3)
835      a10 = mm_alpha_f(e4)
836 
837      ! Computes gamma pre-factor
838      res = kfm * mm_get_btk(4,0) * (a1*rcs**0.5_mm_wp + 2._mm_wp*a2*a3*rcff1/rcs**0.5_mm_wp + &
839                                     a4*a5*rcff2/rcs**1.5_mm_wp +                              &
840                                     a6*a7*rcs**2/rcf**1.5_mm_wp +                             &
841                                     2._mm_wp*a8*a9*rcs*rcff3 +                                &
842                                     a10*rcff4)
843      RETURN
844    END FUNCTION g0sffm
845
846
847    ELEMENTAL FUNCTION g0fffm(rcf, kfm) RESULT(res)
848      !! Get gamma pre-factor for the 0th order moment with FF interactions in the Free Molecular flow regime.
849      !!
850      REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution.
851      REAL(kind=mm_wp), INTENT(in) :: kfm !! Thermodynamic free molecular flow regime pre-factor.
852      REAL(kind=mm_wp) :: res             !! Gamma coagulation kernel pre-factor.
853      REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, e1, e2, e3, rcff
854      res = 0._mm_wp ; IF (rcf <= 0._mm_wp) RETURN
855      ! Computes mm_alpha and other coefficients
856      e1 = 3._mm_wp/mm_df
857      e2 = (6_mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
858      e3 = (12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
859      rcff = rcf**e3 * mm_rb2ra**2
860      a1 = mm_alpha_f(e3)         ; a2 = mm_alpha_f(e1)          ; a3 = mm_alpha_f(e2)
861      a4 = mm_alpha_f(-1.5_mm_wp) ; a5 = mm_alpha_f(2._mm_wp*e1)
862      ! Computes gamma pre-factor
863      res = kfm * (rcff*mm_get_btk(3,0)) * (a1 + 2._mm_wp*a2*a3 + a4*a5)
864      RETURN
865    END FUNCTION g0fffm
866
867
868    ELEMENTAL FUNCTION g3ssfm(rcs, kfm) RESULT(res)
869      !! Get gamma pre-factor for the 3rd order moment with SS interactions in the Free Molecular flow regime.
870      !!
871      REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution.
872      REAL(kind=mm_wp), INTENT(in) :: kfm !! Thermodynamic free molecular flow regime pre-factor.
873      REAL(kind=mm_wp) :: res             !! Gamma coagulation kernel pre-factor.
874      REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11
875      res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN
876      ! Computes mm_alpha coefficients
877      a1 = mm_alpha_s(3.5_mm_wp)  ; a2 = mm_alpha_s(1._mm_wp)  ; a3 = mm_alpha_s(2.5_mm_wp)
878      a4 = mm_alpha_s(2._mm_wp)   ; a5 = mm_alpha_s(1.5_mm_wp) ; a6 = mm_alpha_s(5._mm_wp)
879      a7 = mm_alpha_s(-1.5_mm_wp) ; a8 = mm_alpha_s(4._mm_wp)  ; a9 = mm_alpha_s(-0.5_mm_wp)
880      a10 = mm_alpha_s(3._mm_wp)  ; a11 = mm_alpha_s(0.5_mm_wp)
881      ! Computes gamma pre-factor
882      res = kfm * (mm_get_btk(1,3)/(a10**2*rcs**2.5_mm_wp)) * (a1 + 2._mm_wp*a2*a3 + &
883                                                               a4*a5 + a6*a7 +       &
884                                                               2._mm_wp*a8*a9 + a10*a11)
885      RETURN
886    END FUNCTION g3ssfm
887
888
889    ELEMENTAL FUNCTION g3sffm(rcs, rcf, kfm) RESULT(res)
890      !! Get gamma pre-factor for the 3rd order moment with SF interactions in the Free Molecular flow regime.
891      !!
892      REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution.
893      REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution.
894      REAL(kind=mm_wp), INTENT(in) :: kfm !! Thermodynamic free molecular flow regime pre-factor.
895      REAL(kind=mm_wp) :: res             !! Gamma coagulation kernel pre-factor.
896      REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12
897      REAL(kind=mm_wp) :: e1, e2, e3, rcff1, rcff2, rcff3
898      res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN
899      ! Computes mm_alpha and other coefficients
900      e1 = 3._mm_wp/mm_df
901      e2 = (6._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
902      e3 = (12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
903
904      rcff1 = rcf**e1 * mm_rb2ra
905      rcff2 = rcf**e2 * mm_rb2ra
906      rcff3 = rcf**e3 * mm_rb2ra**2
907     
908      a1 = mm_alpha_s(3.5_mm_wp)  ; a2 = mm_alpha_s(2.5_mm_wp)   ; a3 = mm_alpha_f(e1)
909      a4 = mm_alpha_s(1.5_mm_wp)  ; a5 = mm_alpha_f(2._mm_wp*e1) ; a6 = mm_alpha_s(5._mm_wp)
910      a7 = mm_alpha_f(-1.5_mm_wp) ; a8 = mm_alpha_s(4._mm_wp)    ; a9 = mm_alpha_f(e2)
911      a10 = mm_alpha_s(3._mm_wp)  ; a11 = mm_alpha_f(e3)         ; a12 = mm_alpha_f(3._mm_wp)
912     
913      ! Computes gamma pre-factor
914      res = kfm * (mm_get_btk(4,3)/(a10*a12*(rcs*rcf)**3)) * (a1*rcs**3.5_mm_wp + &
915                                                              2._mm_wp*a2*a3*rcs**2.5_mm_wp*rcff1 + &
916                                                              a4*a5*rcs**1.5_mm_wp*rcff1**2 +       &
917                                                              a6*a7*rcs**5/rcf**1.5_mm_wp +         &
918                                                              2._mm_wp*a8*a9*rcs**4*rcff2 +         &
919                                                              a10*a11*rcs**3*rcff3)
920      RETURN
921    END FUNCTION g3sffm
922
923END MODULE MP2M_HAZE
Note: See TracBrowser for help on using the repository browser.