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

Last change on this file since 4080 was 4032, checked in by debatzbr, 6 weeks ago

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