source: trunk/LMDZ.TITAN/libf/muphytitan/mm_haze.f90 @ 3094

Last change on this file since 3094 was 3090, checked in by slebonnois, 15 months ago

BdeBatz? : Cleans microphysics and makes few corrections for physics
File size: 37.8 KB
Line 
1! Copyright 2013-2015,2017,2022 Université de Reims Champagne-Ardenne
2! Contributor: J. Burgalat (GSMA, URCA)
3! email of the author : jeremie.burgalat@univ-reims.fr
4!
5! This software is a computer program whose purpose is to compute
6! microphysics processes using a two-moments scheme.
7!
8! This library is governed by the CeCILL-B license under French law and
9! abiding by the rules of distribution of free software.  You can  use,
10! modify and/ or redistribute the software under the terms of the CeCILL-B
11! license as circulated by CEA, CNRS and INRIA at the following URL
12! "http://www.cecill.info".
13!
14! As a counterpart to the access to the source code and  rights to copy,
15! modify and redistribute granted by the license, users are provided only
16! with a limited warranty  and the software's author,  the holder of the
17! economic rights,  and the successive licensors  have only  limited
18! liability.
19!
20! In this respect, the user's attention is drawn to the risks associated
21! with loading,  using,  modifying and/or developing or reproducing the
22! software by the user in light of its specific status of free software,
23! that may mean  that it is complicated to manipulate,  and  that  also
24! therefore means  that it is reserved for developers  and  experienced
25! professionals having in-depth computer knowledge. Users are therefore
26! encouraged to load and test the software's suitability as regards their
27! requirements in conditions enabling the security of their systems and/or
28! data to be ensured and,  more generally, to use and operate it in the
29! same conditions as regards security.
30!
31! The fact that you are presently reading this means that you have had
32! knowledge of the CeCILL-B license and that you accept its terms.
33
34!! file: mm_haze.f90
35!! summary: Haze microphysics module.
36!! author: J. Burgalat
37!! date: 2013-2015,2017,2022
38MODULE MM_HAZE
39  !! Haze microphysics module.
40  !!
41  !! This module contains all definitions of the microphysics processes related to aerosols:
42  !!
43  !! - [coagulation](page/haze.html#coagulation)
44  !! - [sedimentation](page/haze.html#sedimentation)
45  !! - [production](page/haze.html#production)
46  !!
47  !! @note
48  !! The production function is specific to Titan, where aerosols are created above the detached
49  !! haze layer. No other source is taken into account.  This process is controled by two parameters,
50  !! the pressure level of production and the production rate. Then both M0 and M3 of the aerosols
51  !! distribution are updated in the production zone by addition of the production rate along a
52  !! gaussian shape.
53  !!
54  !! @note
55  !! The interface methods always uses the global variables defined in [[mm_globals(module)]] when
56  !! values (any kind, temperature, pressure, moments...) over the vertical grid are required.
57  !!
58  !! @warning
59  !! The tendencies returned by the method are always defined over the vertical grid from __TOP__
60  !! to __GROUND__.
61  !!
62  !! @todo
63  !! Modify tests on tendencies vectors to get sure that allocation is done:
64  !! Currently, we assume the compiler handles automatic allocation of arrays.
65  USE MM_MPREC
66  USE MM_GLOBALS
67  USE MM_INTERFACES
68  USE MM_METHODS
69  IMPLICIT NONE
70
71  PRIVATE
72
73  PUBLIC :: mm_haze_microphysics, mm_haze_coagulation, mm_haze_sedimentation, &
74    mm_haze_production
75
76CONTAINS
77
78  !============================================================================
79  ! HAZE MICROPHYSICS INTERFACE SUBROUTINE
80  !============================================================================
81
82  SUBROUTINE mm_haze_microphysics(dm0a_s,dm3a_s,dm0a_f,dm3a_f)
83    !! Get the evolution of moments tracers through haze microphysics processes.
84    !!
85    !! The subroutine is a wrapper to the haze microphysics methods. It computes the tendencies
86    !! of moments tracers for coagulation, sedimentation and production processes for the
87    !! atmospheric column.
88    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0a_s
89    !! Tendency of the 0th order moment of the spherical mode distribution (\(m^{-3}\)).
90    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3a_s
91    !! Tendency of the 3rd order moment of the spherical mode distribution (\(m^{3}.m^{-3}\)).
92    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0a_f
93    !! Tendency of the 0th order moment of the fractal mode distribution (\(m^{-3}\)).
94    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3a_f
95    !! Tendency of the 3rd order moment of the fractal mode distribution (\(m^{3}.m^{-3}\)).
96    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm0as
97    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm3as
98    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm0af
99    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm3af
100
101    dm0a_s = 0._mm_wp ; dm3a_s = 0._mm_wp ; dm0a_f = 0._mm_wp ; dm3a_f = 0._mm_wp
102
103    ALLOCATE(zdm0as(mm_nla),zdm3as(mm_nla),zdm0af(mm_nla),zdm3af(mm_nla))
104    zdm0as(1:mm_nla) = 0._mm_wp
105    zdm3as(1:mm_nla) = 0._mm_wp
106    zdm0af(1:mm_nla) = 0._mm_wp
107    zdm3af(1:mm_nla) = 0._mm_wp
108
109    IF (mm_w_haze_coag) THEN
110      ! Calls coagulation
111      call mm_haze_coagulation(dm0a_s,dm3a_s,dm0a_f,dm3a_f)
112    ENDIF
113
114    IF (mm_w_haze_sed) THEN
115      ! Calls sedimentation
116      call mm_haze_sedimentation(zdm0as,zdm3as,zdm0af,zdm3af)
117
118      ! Computes precipitations
119      mm_aer_prec = SUM(zdm3as*mm_dzlev) + SUM(zdm3af*mm_dzlev)
120
121      ! Updates tendencies
122      dm0a_s=dm0a_s+zdm0as ; dm3a_s=dm3a_s+zdm3as
123      dm0a_f=dm0a_f+zdm0af ; dm3a_f=dm3a_f+zdm3af
124    ENDIF
125
126    IF (mm_w_haze_prod) THEN
127      call mm_haze_production(zdm0as,zdm3as)
128      ! We only produce spherical aerosols
129      dm0a_s=dm0a_s+zdm0as ; dm3a_s=dm3a_s+zdm3as
130    ENDIF
131
132    RETURN
133  END SUBROUTINE mm_haze_microphysics
134
135
136  !============================================================================
137  ! COAGULATION PROCESS RELATED METHODS
138  !============================================================================
139
140  SUBROUTINE mm_haze_coagulation(dM0s,dM3s,dM0f,dM3f)
141    !! Get the evolution of the aerosols moments vertical column due to coagulation process.
142    !!
143    !! This is main method of the coagulation process:
144    !!
145    !! 1. Computes gamma pre-factor for each parts of the coagulation equation(s)
146    !! 2. Applies the electic correction on the gamma pre-factor
147    !! 3. Computes the specific flow regime "kernels"
148    !! 4. Computes the harmonic mean of the kernels
149    !! 5. Finally computes the tendencies of the moments.
150    !!
151    !! All arguments are assumed vectors of __N__ elements where __N__ is the total number of
152    !! vertical __layers__.
153    !!
154    !! @note
155    !! The method uses directly the global variables related to the vertical atmospheric structure
156    !! stored in [[mm_globals(module)]]. Consequently they must be updated before calling the subroutine.
157    !!
158    !! @bug
159    !! If the transfert probabilities are set to 1 for the two flow regimes (pco and pfm),
160    !! a floating point exception occured (i.e. a NaN) as we perform a division by zero
161    !!
162    !! @todo
163    !! Get rid of the fu\*\*\*\* STOP statement...
164    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM0s
165    !! Tendency of the 0th order moment of the spherical size-distribution over a time step (\(m^{-3}\)).
166    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM3s
167    !! Tendency of the 3rd order moment of the spherical size-distribution (\(m^{3}.m^{-3}\)).
168    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM0f
169    !! Tendency of the 0th order moment of the fractal size-distribution over a time step (\(m^{-3}\)).
170    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM3f
171    !! Tendency of the 3rd order moment of the fractal size-distribution over a time step (\(m^{3}.m^{-3}\)).
172    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: c_kco,c_kfm,c_slf,tmp, &
173      kco,kfm,pco,pfm,mq
174    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: a_ss,a_sf,b_ss,b_ff,c_ss,c_sf
175    INTEGER                                     :: i
176
177    IF (mm_coag_choice < 0 .OR. mm_coag_choice > mm_coag_ss+mm_coag_sf+mm_coag_ff) &
178      STOP "invalid choice for coagulation mode interaction activation"
179
180    ! Alloctes local arrays
181    ALLOCATE(kco(mm_nla),kfm(mm_nla),c_slf(mm_nla), &
182      c_kco(mm_nla),c_kfm(mm_nla),mq(mm_nla), &
183      pco(mm_nla),pfm(mm_nla))
184    ALLOCATE(a_ss(mm_nla),a_sf(mm_nla), &
185      b_ss(mm_nla),b_ff(mm_nla), &
186      c_ss(mm_nla),c_sf(mm_nla))
187
188    a_ss(:) = 0._mm_wp ; a_sf(:) = 0._mm_wp
189    b_ss(:) = 0._mm_wp ; b_ff(:) = 0._mm_wp
190    c_ss(:) = 0._mm_wp ; c_sf(:) = 0._mm_wp
191
192    ! gets kco, kfm pre-factors
193    c_kco(:) = mm_get_kco(mm_temp) ; c_kfm(:) = mm_get_kfm(mm_temp)
194    ! get slf (slip-flow factor)
195    c_slf(:) = mm_akn * mm_lambda_g(mm_temp,mm_play)
196
197    DO i=1,mm_nla
198      ! SS interactions
199      IF (mm_rcs(i) > mm_rc_min .AND. IAND(mm_coag_choice,mm_coag_ss) /= 0) THEN
200        ! compute probability for M0/CO and M0/FM (resp.)
201        pco(i) = mm_ps2s(mm_rcs(i),0,0,mm_temp(i),mm_play(i))
202        pfm(i) = mm_ps2s(mm_rcs(i),0,1,mm_temp(i),mm_play(i))
203        ! (A_SS_CO x A_SS_FM) / (A_SS_CO + A_SS_FM)
204        kco(i) = g0ssco(mm_rcs(i),c_slf(i),c_kco(i))
205        kfm(i) = g0ssfm(mm_rcs(i),c_kfm(i))
206        IF (kco(i)*(pco(i)-2._mm_wp)+kfm(i)*(pfm(i)-2._mm_wp) /=0) THEN
207          a_ss(i) = (kco(i)*(pco(i)-2._mm_wp)*kfm(i)*(pfm(i)-2._mm_wp))/(kco(i)*(pco(i)-2._mm_wp)+kfm(i)*(pfm(i)-2._mm_wp))
208        ENDIF
209        ! (B_SS_CO x B_SS_FM) / (B_SS_CO + B_SS_FM)
210        kco(i) = kco(i) * (1._mm_wp-pco(i)) ; kfm(i) = kfm(i) * (1._mm_wp-pfm(i))
211        IF (kco(i) + kfm(i) /= 0._mm_wp) THEN
212          b_ss(i) = kco(i)*kfm(i)/(kco(i)+kfm(i))
213        ENDIF
214        ! compute and apply eletric charge correction for M0/SS interactions
215        mq(i) = mm_qmean(mm_rcs(i),mm_rcs(i),0,'SS',mm_temp(i),mm_play(i))
216
217        a_ss(i) = a_ss(i) * mq(i)
218        b_ss(i) = b_ss(i) * mq(i)
219        kco(i) = 0._mm_wp ; kfm(i) = 0._mm_wp ; mq(i) = 1._mm_wp
220        ! compute probability for M3/CO and M3/FM (resp.)
221        pco(i) = mm_ps2s(mm_rcs(i),3,0,mm_temp(i),mm_play(i))
222        pfm(i) = mm_ps2s(mm_rcs(i),3,1,mm_temp(i),mm_play(i))
223        ! (C_SS_CO x C_SS_FM) / (C_SS_CO + C_SS_FM)
224        kco(i) = g3ssco(mm_rcs(i),c_slf(i),c_kco(i))*(pco(i)-1._mm_wp)
225        kfm(i) = g3ssfm(mm_rcs(i),c_kfm(i))*(pfm(i)-1._mm_wp)
226        IF (kco(i) + kfm(i) /= 0._mm_wp) THEN
227          c_ss(i) = (kco(i)*kfm(i))/(kco(i)+kfm(i))
228        ENDIF
229        IF (b_ss(i) <= 0._mm_wp) c_ss(i) = 0._mm_wp
230        ! compute and apply eletric charge correction for M3/SS interactions
231        mq(i) = mm_qmean(mm_rcs(i),mm_rcs(i),3,'SS',mm_temp(i),mm_play(i))
232        c_ss(i) = c_ss(i) * mq(i)
233      ENDIF
234      kco(i) = 0._mm_wp ; kfm(i) = 0._mm_wp ; mq(i) = 1._mm_wp
235
236      ! SF interactions
237      IF (mm_rcs(i) > mm_rc_min .AND. mm_rcf(i) > mm_rc_min .AND. IAND(mm_coag_choice,mm_coag_sf) /= 0) THEN
238        ! (A_SF_CO x A_SF_FM) / (A_SF_CO + A_SF_FM)
239        kco(i) = g0sfco(mm_rcs(i),mm_rcf(i),c_slf(i),c_kco(i))
240        kfm(i) = g0sffm(mm_rcs(i),mm_rcf(i),c_kfm(i))
241        IF(kco(i)+kfm(i) /= 0._mm_wp) THEN
242          a_sf(i) = (kco(i)*kfm(i))/(kco(i)+kfm(i))
243        ENDIF
244        ! compute and apply eletric charge correction for M0/SF interactions
245        mq(i) = mm_qmean(mm_rcs(i),mm_rcf(i),0,'SF',mm_temp(i),mm_play(i))
246        a_sf(i) = a_sf(i) * mq(i)
247        ! (C_SF_CO x C_SF_FM) / (C_SF_CO + C_SF_FM)
248        kco(i) = g3sfco(mm_rcs(i),mm_rcf(i),c_slf(i),c_kco(i))
249        kfm(i) = g3sffm(mm_rcs(i),mm_rcf(i),c_kfm(i))
250        IF (kco(i)+kfm(i) /= 0._mm_wp) THEN
251          c_sf(i) = (kco(i)*kfm(i))/(kco(i)+kfm(i))
252        ENDIF
253        ! compute and apply eletric charge correction for M3/SF interactions
254        mq(i) = mm_qmean(mm_rcs(i),mm_rcf(i),3,'SF',mm_temp(i),mm_play(i))
255        c_sf(i) = c_sf(i) * mq(i)
256      ENDIF
257      kco(i) = 0._mm_wp ; kfm(i) = 0._mm_wp ; mq(i) = 1._mm_wp
258      ! FF interactions
259      IF(mm_rcf(i) > mm_rc_min .AND. IAND(mm_coag_choice,mm_coag_sf) /= 0) THEN
260        ! (B_FF_CO x B_FF_FM) / (B_FF_CO + B_FF_FM)
261        kco(i) = g0ffco(mm_rcf(i),c_slf(i),c_kco(i))
262        kfm(i) = g0fffm(mm_rcf(i),c_kfm(i))
263        b_ff(i) = (kco(i)*kfm(i))/(kco(i)+kfm(i))
264        ! compute and apply eletric charge correction for M0/FF interactions
265        mq(i) = mm_qmean(mm_rcf(i),mm_rcf(i),0,'FF',mm_temp(i),mm_play(i))
266        b_ff(i) = b_ff(i) * mq(i)
267      ENDIF
268    ENDDO
269
270    DEALLOCATE(kco,kfm,c_kco,c_kfm,pco,pfm,c_slf)
271
272    ! Now we will use the kharm two by two to compute :
273    ! dm_0_S/mm_dt = kharm(1) * m_0_S^2 - kharm(2) * m_0_S * m_0_F
274    ! dm_0_F/mm_dt = kharm(3) * m_0_S^2 - kharm(4) * m_0_F^2
275    ! dm_3_S/mm_dt = kharm(5) * m_3_S^2 - kharm(6) * m_3_S * m_3_F
276    ! ... and finally :
277    ! dm_3_F/mm_dt = - dm_3_S/mm_dt
278    !
279    ! We use a (semi) implicit scheme : when X appears as square we set one X
280    ! at t+1, the other a t
281    ALLOCATE(tmp(mm_nla))
282    ! --- dm0s
283    tmp(:) = mm_dt*(a_ss*mm_m0aer_s - a_sf*mm_m0aer_f)
284    dm0s(:) =  mm_m0aer_s * (tmp/(1._mm_wp - tmp))
285    ! --- dm0f
286    tmp(:) = b_ff*mm_dt*mm_m0aer_f
287    dm0f(:) = (b_ss*mm_dt*mm_m0aer_s**2 - tmp*mm_m0aer_f)/(1._mm_wp + tmp)
288    ! --- dm3s
289    tmp(:) = mm_dt*(c_ss*mm_m3aer_s - c_sf*mm_m3aer_f)
290    dm3s(:) =  mm_m3aer_s * (tmp/(1._mm_wp - tmp))
291    ! --- dmm3f
292    dm3f(:) = -dm3s
293
294    ! Deallocates memory explicitly ... another obsolete statement :)
295    DEALLOCATE(a_ss,a_sf,b_ss,b_ff,c_ss,c_sf,tmp)
296
297    ! Time to do something else !
298    RETURN
299  END SUBROUTINE mm_haze_coagulation
300
301  ELEMENTAL FUNCTION g0ssco(rcs,c_slf,c_kco) RESULT(res)
302    !! Get &gamma; pre-factor for the 0th order moment with SS interactions in the continuous flow regime.
303    !!
304    !! @note
305    !! If __rcs__ is 0, the function returns 0.
306    REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
307    REAL(kind=mm_wp), INTENT(in) :: c_slf !! Slip-Flow correction pre-factor.
308    REAL(kind=mm_wp), INTENT(in) :: c_kco !! Thermodynamic continuous flow regime pre-factor.
309    REAL(kind=mm_wp) :: res               !! &gamma; coagulation kernel pre-factor.
310    REAL(kind=mm_wp) :: a1, a2, a3
311    res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN
312    ! computes mm_alpha coefficients
313    a1=mm_alpha_s(1._mm_wp) ; a2=mm_alpha_s(-1._mm_wp) ; a3=mm_alpha_s(-2._mm_wp)
314    ! Computes gamma pre-factor
315    res = (1._mm_wp + a1*a2 + c_slf/rcs *(a2+a1*a3))*c_kco
316    RETURN
317  END FUNCTION g0ssco
318
319  ELEMENTAL FUNCTION g0sfco(rcs,rcf,c_slf,c_kco) RESULT(res)
320    !! Get &gamma; pre-factor for the 0th order moment with SF interactions in the continuous flow regime.
321    !!
322    !! @note
323    !! If __rcs__ or __rcf__ is 0, the function returns 0.
324    REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
325    REAL(kind=mm_wp), INTENT(in) :: rcf   !! Characteristic radius of the fractal size-distribution.
326    REAL(kind=mm_wp), INTENT(in) :: c_slf !! Slip-Flow correction pre-factor.
327    REAL(kind=mm_wp), INTENT(in) :: c_kco !! Thermodynamic continuous flow regime pre-factor.
328    REAL(kind=mm_wp) :: res               !! &gamma; coagulation kernel pre-factor.
329    REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, e, rcff
330    res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN
331    e=3._mm_wp/mm_df ; rcff = rcf**e*mm_rb2ra
332    ! computes mm_alpha coefficients
333    a1=mm_alpha_s(1._mm_wp)   ; a2=mm_alpha_f(-e)   ; a3=mm_alpha_f(e)
334    a4=mm_alpha_s(-1._mm_wp) ; a5=mm_alpha_s(-2._mm_wp) ; a6=mm_alpha_f(-2._mm_wp*e)
335    ! Computes gamma pre-factor
336    res = c_kco*( 2._mm_wp + a1*a2*rcs/rcff + a4*a3*rcff/rcs + c_slf*( a4/rcs + &
337      a2/rcff + a5*a3*rcff/rcs**2 + a1*a6*rcs/rcff**2))
338    RETURN
339  END FUNCTION g0sfco
340
341  ELEMENTAL FUNCTION g0ffco(rcf,c_slf,c_kco) RESULT(res)
342    !! Get &gamma; pre-factor for the 0th order moment with FF interactions in the continuous flow regime.
343    !!
344    !! @note
345    !! If __rcf__ is 0, the function returns 0.
346    REAL(kind=mm_wp), INTENT(in) :: rcf   !! Characteristic radius of the fractal size-distribution.
347    REAL(kind=mm_wp), INTENT(in) :: c_slf !! Slip-Flow correction pre-factor.
348    REAL(kind=mm_wp), INTENT(in) :: c_kco !! Thermodynamic continuous flow regime pre-factor.
349    REAL(kind=mm_wp) :: res               !! &gamma; coagulation kernel pre-factor.
350    REAL(kind=mm_wp) :: a1, a2, a3, e, rcff
351    res = 0._mm_wp ; IF (rcf <= 0._mm_wp) RETURN
352    ! computes mm_alpha coefficients
353    e = 3._mm_wp/mm_df ; rcff = rcf**e*mm_rb2ra
354    a1=mm_alpha_f(e) ; a2=mm_alpha_f(-e) ; a3=mm_alpha_s(-2._mm_wp*e)
355    ! Computes gamma pre-factor
356    res = (1._mm_wp + a1*a2 + c_slf/rcff *(a2+a1*a3))*c_kco
357    RETURN
358  END FUNCTION g0ffco
359
360  ELEMENTAL FUNCTION g3ssco(rcs, c_slf, c_kco) RESULT(res)
361    !! Get &gamma; pre-factor for the 3rd order moment with SS interactions in the continuous flow regime.
362    !!
363    !! @note
364    !! If __rcs__ is 0, the function returns 0.
365    REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
366    REAL(kind=mm_wp), INTENT(in) :: c_slf !! Slip-Flow correction pre-factor.
367    REAL(kind=mm_wp), INTENT(in) :: c_kco !! Thermodynamic continuous flow regime pre-factor.
368    REAL(kind=mm_wp) :: res               !! &gamma; coagulation kernel pre-factor.
369    REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6
370    res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN
371    ! computes mm_alpha coefficients
372    a1=mm_alpha_s(3._mm_wp) ; a2=mm_alpha_s(2._mm_wp)  ; a3=mm_alpha_s(1._mm_wp)
373    a4=mm_alpha_s(4._mm_wp) ; a5=mm_alpha_s(-1._mm_wp) ; a6=mm_alpha_s(-2._mm_wp)
374
375    ! Computes gamma pre-factor
376    res = (2._mm_wp*a1 + a2*a3 + a4*a5 + c_slf/rcs*(a3**2 + a4*a6 + a1*a5 + a2))* &
377      c_kco/(a1**2*rcs**3)
378    RETURN
379  END FUNCTION g3ssco
380
381  ELEMENTAL FUNCTION g3sfco(rcs, rcf, c_slf, c_kco) RESULT(res)
382    !! Get &gamma; pre-factor for the 3rd order moment with SF interactions in the continuous flow regime.
383    !!
384    !! @note
385    !! If __rcs__ or __rcf__ is 0, the function returns 0.
386    REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
387    REAL(kind=mm_wp), INTENT(in) :: rcf   !! Characteristic radius of the fractal size-distribution.
388    REAL(kind=mm_wp), INTENT(in) :: c_slf !! Slip-Flow correction pre-factor.
389    REAL(kind=mm_wp), INTENT(in) :: c_kco !! Thermodynamic continuous flow regime pre-factor.
390    REAL(kind=mm_wp) :: res               !! &gamma; coagulation kernel pre-factor.
391    REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, a7, a8, e, rcff
392    res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN
393    ! computes mm_alpha coefficients
394    e=3._mm_wp/mm_df ; rcff = rcf**e*mm_rb2ra
395    a1=mm_alpha_s(3._mm_wp)    ; a2=mm_alpha_s(4._mm_wp) ; a3=mm_alpha_f(-e)
396    a4=mm_alpha_s(2._mm_wp)    ; a5=mm_alpha_f(e)    ; a6=mm_alpha_s(1._mm_wp)
397    a7=mm_alpha_f(-2._mm_wp*e) ; a8=mm_alpha_f(3._mm_wp)
398    ! Computes gamma pre-factor
399    res = (2._mm_wp*a1*rcs**3 + a2*rcs**4*a3/rcff + a4*rcs**2*a5*rcff + &
400      c_slf *( a4*rcs**2 + a1*rcs**3*a3/rcff + a6*rcs*a5*rcff + &
401      a2*rcs**4*a7/rcff**2))* c_kco/(a1*a8*(rcs*rcf)**3)
402    RETURN
403  END FUNCTION g3sfco
404
405  ELEMENTAL FUNCTION g0ssfm(rcs, c_kfm) RESULT(res)
406    !! Get &gamma; pre-factor for the 0th order moment with SS interactions in the Free Molecular flow regime.
407    !!
408    !! @note
409    !! If __rcs__ is 0, the function returns 0.
410    REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
411    REAL(kind=mm_wp), INTENT(in) :: c_kfm !! Thermodynamic free molecular flow regime pre-factor.
412    REAL(kind=mm_wp) :: res               !! &gamma; coagulation kernel pre-factor.
413    REAL(kind=mm_wp) :: a1, a2, a3, a4, a5
414    res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN
415    ! computes mm_alpha coefficients
416    a1=mm_alpha_s(0.5_mm_wp) ; a2=mm_alpha_s(1._mm_wp) ; a3=mm_alpha_s(-0.5_mm_wp)
417    a4=mm_alpha_s(2._mm_wp)   ; a5=mm_alpha_s(-1.5_mm_wp)
418    ! Computes gamma pre-factor
419    res = (a1 + 2._mm_wp*a2*a3 + a4*a5)*rcs**0.5_mm_wp*mm_get_btk(1,0)*c_kfm
420    RETURN
421  END FUNCTION g0ssfm
422
423  ELEMENTAL FUNCTION g0sffm(rcs, rcf, c_kfm) RESULT(res)
424    !> Get &gamma; pre-factor for the 0th order moment with SF interactions in the Free Molecular flow regime.
425    !!
426    !! @note
427    !! If __rcs__ or __rcf__ is 0, the function returns 0.
428    REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
429    REAL(kind=mm_wp), INTENT(in) :: rcf   !! Characteristic radius of the fractal size-distribution.
430    REAL(kind=mm_wp), INTENT(in) :: c_kfm !! Thermodynamic free molecular flow regime pre-factor.
431    REAL(kind=mm_wp) :: res               !! &gamma; coagulation kernel pre-factor.
432    REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, a7, a8, a9, a10
433    REAL(kind=mm_wp) :: e1, e2, e3, e4
434    REAL(kind=mm_wp) :: rcff1, rcff2, rcff3, rcff4
435    res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN
436    ! computes mm_alpha coefficients
437    e1 = 3._mm_wp/mm_df
438    e2 = 6._mm_wp/mm_df
439    e3 = (6._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
440    e4 = (12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
441
442    rcff1 = mm_rb2ra * rcf**e1
443    rcff2 = rcff1**2
444    rcff3 = mm_rb2ra * rcf**e3
445    rcff4 = mm_rb2ra**2 * rcf**e4
446
447    a1=mm_alpha_s(0.5_mm_wp)  ; a2=mm_alpha_s(-0.5_mm_wp) ; a3=mm_alpha_f(e1)
448    a4=mm_alpha_s(-1.5_mm_wp) ; a5=mm_alpha_f(e2)      ; a6=mm_alpha_s(2._mm_wp)
449    a7=mm_alpha_f(-1.5_mm_wp) ; a8=mm_alpha_s(1._mm_wp)   ; a9=mm_alpha_f(e3)
450    a10=mm_alpha_f(e4)
451
452    ! Computes gamma pre-factor
453    res = (a1*rcs**0.5_mm_wp + 2._mm_wp*rcff1*a2*a3/rcs**0.5_mm_wp + a4*a5*rcff2/rcs**1.5_mm_wp + &
454      a6*a7*rcs**2/rcf**1.5_mm_wp + 2._mm_wp*a8*a9*rcs*rcff3 + a10*rcff4 &
455      )*mm_get_btk(4,0)*c_kfm
456    RETURN
457  END FUNCTION g0sffm
458
459  ELEMENTAL FUNCTION g0fffm(rcf, c_kfm) RESULT(res)
460    !! Get &gamma; pre-factor for the 0th order moment with FF interactions in the Free Molecular flow regime.
461    !!
462    !! @note
463    !! If __rcf__ is 0, the function returns 0.
464    REAL(kind=mm_wp), INTENT(in) :: rcf   !! Characteristic radius of the fractal size-distribution.
465    REAL(kind=mm_wp), INTENT(in) :: c_kfm !! Thermodynamic free molecular flow regime pre-factor.
466    REAL(kind=mm_wp) :: res               !! &gamma; coagulation kernel pre-factor.
467    REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, e1, e2, e3, rcff
468    res = 0._mm_wp ; IF (rcf <= 0._mm_wp) RETURN
469    ! computes mm_alpha coefficients
470    e1=3._mm_wp/mm_df ; e2=(6_mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
471    e3=(12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
472    rcff=mm_rb2ra**2*rcf**e3
473    a1=mm_alpha_f(e3)      ; a2=mm_alpha_f(e1) ; a3=mm_alpha_f(e2)
474    a4=mm_alpha_f(-1.5_mm_wp) ; a5=mm_alpha_f(2._mm_wp*e1)
475    ! Computes gamma pre-factor
476    res = (a1 + 2._mm_wp*a2*a3 + a4*a5)*rcff*mm_get_btk(3,0)*c_kfm
477    RETURN
478  END FUNCTION g0fffm
479
480  ELEMENTAL FUNCTION g3ssfm(rcs, c_kfm) RESULT(res)
481    !! Get &gamma; pre-factor for the 3rd order moment with SS interactions in the Free Molecular flow regime.
482    !!
483    !! @note
484    !! If __rcs__ is 0, the function returns 0.
485    REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
486    REAL(kind=mm_wp), INTENT(in) :: c_kfm !! Thermodynamic free molecular flow regime pre-factor.
487    REAL(kind=mm_wp) :: res               !! &gamma; coagulation kernel pre-factor.
488    REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11
489    res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN
490    ! computes mm_alpha coefficients
491    a1=mm_alpha_s(3.5_mm_wp)  ; a2=mm_alpha_s(1._mm_wp)  ; a3=mm_alpha_s(2.5_mm_wp)
492    a4=mm_alpha_s(2._mm_wp)   ; a5=mm_alpha_s(1.5_mm_wp) ; a6=mm_alpha_s(5._mm_wp)
493    a7=mm_alpha_s(-1.5_mm_wp) ; a8=mm_alpha_s(4._mm_wp)  ; a9=mm_alpha_s(-0.5_mm_wp)
494    a10=mm_alpha_s(3._mm_wp) ; a11=mm_alpha_s(0.5_mm_wp)
495    ! Computes gamma pre-factor
496    res = (a1 + 2._mm_wp*a2*a3 + a4*a5 + a6*a7 + 2._mm_wp*a8*a9 + a10*a11) &
497      *mm_get_btk(1,3)*c_kfm/(a10**2*rcs**2.5_mm_wp)
498    RETURN
499  END FUNCTION g3ssfm
500
501  ELEMENTAL FUNCTION g3sffm(rcs, rcf, c_kfm) RESULT(res)
502    !! Get &gamma; pre-factor for the 3rd order moment with SF interactions in the Free Molecular flow regime.
503    !!
504    !! @note
505    !! If __rcs__ or __rcf__ is 0, the function returns 0.
506    REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
507    REAL(kind=mm_wp), INTENT(in) :: rcf   !! Characteristic radius of the fractal size-distribution.
508    REAL(kind=mm_wp), INTENT(in) :: c_kfm !! Thermodynamic free molecular flow regime pre-factor.
509    REAL(kind=mm_wp) :: res               !! &gamma; coagulation kernel pre-factor.
510    REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12
511    REAL(kind=mm_wp) :: e1, e2, e3, rcff1, rcff2, rcff3
512    res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN
513    ! computes mm_alpha coefficients
514    e1=3._mm_wp/mm_df
515    e2=(6._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
516    e3=(12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
517    rcff1=mm_rb2ra*rcf**e1 ; rcff2=mm_rb2ra*rcf**e2 ; rcff3=mm_rb2ra**2*rcf**e3
518    a1=mm_alpha_s(3.5_mm_wp)  ; a2=mm_alpha_s(2.5_mm_wp)   ; a3=mm_alpha_f(e1)
519    a4=mm_alpha_s(1.5_mm_wp)  ; a5=mm_alpha_f(2._mm_wp*e1) ; a6=mm_alpha_s(5._mm_wp)
520    a7=mm_alpha_f(-1.5_mm_wp) ; a8=mm_alpha_s(4._mm_wp)    ; a9=mm_alpha_f(e2)
521    a10=mm_alpha_s(3._mm_wp)  ; a11=mm_alpha_f(e3)         ; a12=mm_alpha_f(3._mm_wp)
522    ! Computes gamma pre-factor
523    res = (a1*rcs**3.5_mm_wp + 2._mm_wp*a2*a3*rcs**2.5_mm_wp*rcff1 + a4*a5*rcs**1.5_mm_wp*rcff1**2 + &
524      a6*a7*rcs**5/rcf**1.5_mm_wp + 2._mm_wp*a8*a9*rcs**4*rcff2 + &
525      a10*a11*rcs**3*rcff3)*mm_get_btk(4,3)*c_kfm/(a10*a12*(rcs*rcf)**3)
526    RETURN
527  END FUNCTION g3sffm
528
529  !============================================================================
530  ! SEDIMENTATION PROCESS RELATED METHODS
531  !============================================================================
532
533  SUBROUTINE mm_haze_sedimentation(dm0s,dm3s,dm0f,dm3f)
534    !! Interface to sedimentation algorithm.
535    !!
536    !! The subroutine computes the evolution of each moment of the aerosols tracers
537    !! through sedimentation process and returns their tendencies for a timestep.
538    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0s
539    !! Tendency of the 0th order moment of the spherical mode (\(m^{-3}\)).
540    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3s
541    !! Tendency of the 3rd order moment of the spherical mode (\(m^{3}.m^{-3}\)).
542    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0f
543    !! Tendency of the 0th order moment of the fractal mode (\(m^{-3}\)).
544    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3f
545    !! Tendency of the 3rd order moment of the fractal mode (\(m^{3}.m^{-3}\)).
546    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: ft,fdcor,wth
547    REAL(kind=mm_wp), PARAMETER                 :: fac = 4._mm_wp/3._mm_wp * mm_pi
548
549    ALLOCATE(ft(mm_nle),wth(mm_nle),fdcor(mm_nle))
550
551    !mm_aer_s_flux(:) = 0._mm_wp ; mm_aer_f_flux(:) = 0._mm_wp
552    IF (mm_wsed_m0) THEN
553      ! Spherical particles
554      ! M0
555      call get_weff(mm_m0aer_s,0._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth,fdcor)
556      ft(:) = wth(:) * fdcor(:)  ; mm_m0as_vsed(:) = ft(1:mm_nla)
557      call let_me_fall_in_peace(mm_m0aer_s,-ft,mm_dt,dm0s)
558      ! M3
559      mm_m3as_vsed(:) = ft(1:mm_nla)
560      call let_me_fall_in_peace(mm_m3aer_s,-ft,mm_dt,dm3s)
561      ! get mass flux
562      mm_aer_s_flux(:) = fac * mm_rhoaer * ft(2:) * mm_m3aer_s
563      ! Fractal particles
564      ! M0
565      call get_weff(mm_m0aer_f,0._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth,fdcor)
566      ft(:) = wth(:) * fdcor(:)  ; mm_m0af_vsed(:) = ft(1:mm_nla)
567      call let_me_fall_in_peace(mm_m0aer_f,-ft,mm_dt,dm0f)
568      ! M3
569      mm_m3af_vsed(:) = ft(1:mm_nla)
570      call let_me_fall_in_peace(mm_m3aer_f,-ft,mm_dt,dm3f)
571      ! get mass flux
572      mm_aer_f_flux(:) = fac * mm_rhoaer * ft(2:) * mm_m3aer_f
573    ELSEIF (mm_wsed_m3) THEN
574      ! Spherical particles
575      ! M3
576      call get_weff(mm_m3aer_s,3._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth,fdcor)
577      ft(:) = wth(:) * fdcor(:)  ; mm_m3as_vsed(:) = ft(1:mm_nla)
578      call let_me_fall_in_peace(mm_m3aer_s,-ft,mm_dt,dm3s)
579      ! get mass flux
580      mm_aer_s_flux(:) = fac * mm_rhoaer * ft(2:) * mm_m3aer_s
581      ! M0
582      mm_m0as_vsed(:) = ft(1:mm_nla)
583      call let_me_fall_in_peace(mm_m0aer_s,-ft,mm_dt,dm0s)
584      ! Fractal particles
585      ! M3
586      call get_weff(mm_m3aer_f,3._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth,fdcor)
587      ft(:) = wth(:) * fdcor(:)  ; mm_m3af_vsed(:) = ft(1:mm_nla)
588      call let_me_fall_in_peace(mm_m3aer_f,-ft,mm_dt,dm3f)
589      ! get mass flux
590      mm_aer_f_flux(:) = fac * mm_rhoaer * ft(2:) * mm_m3aer_f
591      ! M0
592      mm_m0af_vsed(:) = ft(1:mm_nla)
593      call let_me_fall_in_peace(mm_m0aer_f,-ft,mm_dt,dm0f)
594    ELSE
595      ! Spherical particles
596      ! M0
597      call get_weff(mm_m0aer_s,0._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth,fdcor)
598      ft(:) = wth(:) * fdcor(:)  ; mm_m0as_vsed(:) = ft(1:mm_nla)
599      call let_me_fall_in_peace(mm_m0aer_s,-ft,mm_dt,dm0s)
600      ! M3
601      call get_weff(mm_m3aer_s,3._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth,fdcor)
602      ft(:) = wth(:) * fdcor(:)  ; mm_m3as_vsed(:) = ft(1:mm_nla)
603      call let_me_fall_in_peace(mm_m3aer_s,-ft,mm_dt,dm3s)
604      ! get mass flux
605      mm_aer_s_flux(:) = fac * mm_rhoaer * ft(2:) * mm_m3aer_s
606      ! Fractal particles
607      ! M0
608      call get_weff(mm_m0aer_f,0._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth,fdcor)
609      ft(:) = wth(:) * fdcor(:)  ; mm_m0af_vsed(:) = ft(1:mm_nla)
610      call let_me_fall_in_peace(mm_m0aer_f,-ft,mm_dt,dm0f)
611      ! M3
612      call get_weff(mm_m3aer_f,3._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth,fdcor)
613      ft(:) = wth(:) * fdcor(:)  ; mm_m3af_vsed(:) = ft(1:mm_nla)
614      call let_me_fall_in_peace(mm_m3aer_f,-ft,mm_dt,dm3f)
615      ! get mass flux
616      mm_aer_f_flux(:) = fac * mm_rhoaer * ft(2:) * mm_m3aer_f
617    ENDIF
618    DEALLOCATE(ft,wth,fdcor)
619    RETURN
620  END SUBROUTINE mm_haze_sedimentation
621
622  SUBROUTINE let_me_fall_in_peace(mk,ft,dt,dmk)
623    !! Compute the tendency of the moment through sedimentation process.
624    !!
625    !!
626    !! The method computes the time evolution of the \(k^{th}\) order moment through sedimentation:
627    !!
628    !! $$ \dfrac{dM_{k}}{dt} = \dfrac{\Phi_{k}}{dz} $$
629    !!
630    !! The equation is resolved using a [Crank-Nicolson algorithm](http://en.wikipedia.org/wiki/Crank-Nicolson_method).
631    !!
632    !! Sedimentation algorithm is quite messy. It appeals to the dark side of the Force and uses evil black magic spells
633    !! from ancient times. It is based on \cite{toon1988b,fiadeiro1977,turco1979} and is an update of the algorithm
634    !! originally implemented in the LMDZ-Titan 2D GCM.
635    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: mk  !! \(k^{th}\) order moment to sediment (in \(m^{k}\)).
636    REAL(kind=mm_wp), INTENT(in)                :: dt  !! Time step (s).
637    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: ft  !! Downward sedimentation flux  (effective velocity of the moment).
638    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dmk !! Tendency of \(k^{th}\) order moment (in \(m^{k}.m^{-3}\)).
639    INTEGER                                 :: i
640    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: as,bs,cs,mko
641    ALLOCATE(as(mm_nla), bs(mm_nla), cs(mm_nla), mko(mm_nla))
642    mko(1:mm_nla) = 0._mm_wp
643    cs(1:mm_nla) = ft(2:mm_nla+1) - mm_dzlay(1:mm_nla)/dt
644    IF (ANY(cs > 0._mm_wp)) THEN
645      ! implicit case
646      as(1:mm_nla) = ft(1:mm_nla)
647      bs(1:mm_nla) = -(ft(2:mm_nle)+mm_dzlay(1:mm_nla)/dt)
648      cs(1:mm_nla) = -mm_dzlay(1:mm_nla)/dt*mk(1:mm_nla)
649      ! (Tri)diagonal matrix inversion
650      mko(1) = cs(1)/bs(1)
651      DO i=2,mm_nla ; mko(i) = (cs(i)-mko(i-1)*as(i))/bs(i) ; ENDDO
652    ELSE
653      ! explicit case
654      as(1:mm_nla)=-mm_dzlay(1:mm_nla)/dt
655      bs(1:mm_nla)=-ft(1:mm_nla)
656      ! boundaries
657      mko(1)=cs(1)*mk(1)/as(1)
658      mko(mm_nla)=(bs(mm_nla)*mk(mm_nla-1)+cs(mm_nla)*mk(mm_nla))/as(mm_nla)
659      ! interior
660      mko(2:mm_nla-1)=(bs(2:mm_nla-1)*mk(1:mm_nla-2) + &
661        cs(2:mm_nla-1)*mk(2:mm_nla-1)   &
662        )/as(2:mm_nla-1)
663    ENDIF
664    dmk = mko - mk
665    DEALLOCATE(as,bs,cs,mko)
666    RETURN
667  END SUBROUTINE let_me_fall_in_peace
668
669  SUBROUTINE get_weff(mk,k,df,rc,dt,afun,wth,corf)
670    !! Get the effective settling velocity for aerosols moments.
671    !!
672    !! This method computes the effective settling velocity of the \(k^{th}\) order moment of aerosol
673    !! tracers. The basic settling velocity (\(v^{eff}_{M_{k}}\)) is computed using the following
674    !! equation:
675    !!
676    !! $$
677    !! \begin{eqnarray*}
678    !! \Phi^{sed}_{M_{k}} &=& \int_{0}^{\infty} n(r) r^{k} \times w(r) dr
679    !!                     == M_{k}  \times v^{eff}_{M_{k}} \\
680    !!     v^{eff}_{M_{k} &=& \dfrac{2 \rho g r_{c}^{\dfrac{3D_{f}-3}{D_{f}}}}
681    !!     {r_{m}^{D_{f}-3}/D_{f}} \times \alpha(k)} \times \left( \alpha \right(
682    !!     \frac{D_{f}(k+3)-3}{D_{f}}\left) + \dfrac{A_{kn}\lambda_{g}}{r_{c}^{
683    !!     3/D_{f}}} \alpha \right( \frac{D_{f}(k+3)-6}{D_{f}}\left)\left)
684    !! \end{eqnarray*}
685    !! $$
686    !!
687    !! \(v^{eff}_{M_{k}\) is then corrected to reduce numerical diffusion of the sedimentation algorithm
688    !! as defined in \cite{toon1988b}.
689    !!
690    !! @warning
691    !! Both __df__, __rc__ and __afun__ must be consistent with each other otherwise wrong values will
692    !! be computed.
693    REAL(kind=mm_wp), INTENT(in), DIMENSION(mm_nla)            :: mk
694    !! Moment of order __k__ (\(m^{k}.m^{-3}\)) at each layer.
695    REAL(kind=mm_wp), INTENT(in), DIMENSION(mm_nla)            :: rc
696    !! Characteristic radius associated to the moment at each layer.
697    REAL(kind=mm_wp), INTENT(in)                               :: k
698    !! The order of the moment.
699    REAL(kind=mm_wp), INTENT(in)                               :: df
700    !! Fractal dimension of the aersols.
701    REAL(kind=mm_wp), INTENT(in)                               :: dt
702    !! Time step (s).
703    REAL(kind=mm_wp), INTENT(out), DIMENSION(mm_nle)           :: wth
704    !! Theoretical Settling velocity at each vertical __levels__ (\( wth \times corf = weff\)).
705    REAL(kind=mm_wp), INTENT(out), DIMENSION(mm_nle), OPTIONAL :: corf
706    !! _Fiadero_ correction factor applied to the theoretical settling velocity at each vertical __levels__.
707    INTERFACE
708      FUNCTION afun(order)
709        !! Inter-moment relation function (see [[mm_interfaces(module):mm_alpha_s(interface)]]).
710        IMPORT mm_wp
711        REAL(kind=mm_wp), INTENT(in) :: order !! Order of the moment.
712        REAL(kind=mm_wp) :: afun              !! Alpha value.
713      END FUNCTION afun
714    END INTERFACE
715    INTEGER                             :: i
716    REAL(kind=mm_wp)                    :: af1,af2,ar1,ar2
717    REAL(kind=mm_wp)                    :: csto,cslf,ratio,wdt,dzb
718    REAL(kind=mm_wp)                    :: rb2ra
719    REAL(kind=mm_wp), DIMENSION(mm_nle) :: zcorf
720    ! ------------------
721
722    wth(:) = 0._mm_wp ; zcorf(:) = 1._mm_wp
723
724    ar1 = (3._mm_wp*df -3._mm_wp)/df    ; ar2 = (3._mm_wp*df -6._mm_wp)/df
725    af1 = (df*(k+3._mm_wp)-3._mm_wp)/df ; af2 = (df*(k+3._mm_wp)-6._mm_wp)/df
726    rb2ra = mm_rm**((df-3._mm_wp)/df)
727    DO i=2,mm_nla
728      IF (rc(i-1) <= 0._mm_wp) CYCLE
729      dzb = (mm_dzlay(i)+mm_dzlay(i-1))/2._mm_wp
730      csto = 2._mm_wp*mm_rhoaer*mm_effg(mm_zlev(i))/(9._mm_wp*mm_eta_g(mm_btemp(i)))
731      cslf = mm_akn * mm_lambda_g(mm_btemp(i),mm_plev(i))
732      wth(i) = - csto/(rb2ra*afun(k)) * (rc(i-1)**ar1 * afun(af1) + cslf/rb2ra * rc(i-1)**ar2 * afun(af2))
733     
734      ! now correct velocity to reduce numerical diffusion
735      IF (.NOT.mm_no_fiadero_w) THEN
736        IF (mk(i) <= 0._mm_wp) THEN
737          ratio = mm_fiadero_max
738        ELSE
739          ratio = MAX(MIN(mk(i-1)/mk(i),mm_fiadero_max),mm_fiadero_min)
740        ENDIF
741        ! apply correction
742        IF ((ratio <= 0.9_mm_wp .OR. ratio >= 1.1_mm_wp) .AND. wth(i) /= 0._mm_wp) THEN
743          wdt = wth(i)*dt
744          ! bugfix: max exponential arg to 30)
745          zcorf(i) = dzb/wdt * (exp(MIN(30._mm_wp,-wdt*log(ratio)/dzb))-1._mm_wp) / (1._mm_wp-ratio)
746          !zcorf(i) = dzb/wdt * (exp(-wdt*log(ratio)/dzb)-1._mm_wp) / (1._mm_wp-ratio)
747        ENDIF
748      ENDIF
749    ENDDO
750    ! last value (ground) set to first layer value: arbitrary :)
751    wth(i) = wth(i-1)
752    IF (PRESENT(corf)) corf(:) = zcorf(:)
753    RETURN
754  END SUBROUTINE get_weff
755
756  !============================================================================
757  ! PRODUCTION PROCESS RELATED METHOD
758  !============================================================================
759
760  SUBROUTINE mm_haze_production(dm0s,dm3s)
761    !! Compute the production of aerosols moments.
762    !!
763    !! The method computes the tendencies of M0 and M3 for the spherical mode through production process.
764    !! Production values are distributed along a normal law in altitude, centered  at
765    !! [[mm_globals(module):mm_p_prod(variable)]] pressure level with a fixed sigma of 20km.
766    !!
767    !! First M3 tendency is computed and M0 is retrieved using the inter-moments relation a spherical
768    !! characteristic radius set to [[mm_globals(module):mm_rc_prod(variable)]].
769    !!
770    !! If [[mm_globals(module):mm_var_prod(variable)]] is set to .true., the method computes time-dependent
771    !! tendencies using a sine wave of anuglar frequency [[mm_globals(module):mm_w_prod(variable)]]
772    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0s !! Tendency of M0 (\(m^{-3}\)).
773    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3s !! Tendency of M3 (\(m^{3}.m^{-3}\)).
774    INTEGER                     :: i
775    REAL(kind=mm_wp)            :: zprod,cprod,timefact
776    REAL(kind=mm_wp), PARAMETER :: sigz  = 20e3_mm_wp, &
777                                   fnorm = 1._mm_wp/(dsqrt(2._mm_wp*mm_pi)*sigz), &
778                                   znorm = dsqrt(2._mm_wp)*sigz
779    REAL(kind=mm_wp), SAVE      :: ctime = 0._mm_wp
780    !$OMP THREADPRIVATE (ctime)
781    zprod = -1._mm_wp
782    ! locate production altitude
783    DO i=1, mm_nla-1
784      IF (mm_plev(i) < mm_p_prod.AND.mm_plev(i+1) >= mm_p_prod) THEN
785        zprod = mm_zlay(i) ; EXIT
786      ENDIF
787    ENDDO
788    IF (zprod < 0._mm_wp) THEN
789      WRITE(*,'(a)') "cannot find aerosols production altitude"
790      call EXIT(11)
791    ENDIF
792
793    dm3s(:)= mm_tx_prod *0.75_mm_wp/mm_pi *mm_dt / mm_rhoaer / 2._mm_wp / mm_dzlev(1:mm_nla) * &
794             (erf((mm_zlev(1:mm_nla)-zprod)/znorm) - &
795             erf((mm_zlev(2:mm_nla+1)-zprod)/znorm))
796    dm0s(:) = dm3s(:)/(mm_rc_prod**3*mm_alpha_s(3._mm_wp))
797
798    IF (mm_var_prod) THEN
799      timefact = mm_d_prod*sin(mm_w_prod*ctime)+1._mm_wp
800      dm3s(:) = timefact*dm3s(:)
801      dm0s(:) = timefact*dm0s(:)
802      ctime = ctime + mm_dt
803    ENDIF
804
805  END SUBROUTINE mm_haze_production
806
807END MODULE MM_HAZE
Note: See TracBrowser for help on using the repository browser.