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

Last change on this file since 3026 was 2109, checked in by jvatant, 6 years ago

Fix some problems for the microphysics :
+ Altitude of the last level at 1e7m from physics was certainly source of divergence
+ Sanity check for negative is moved from within mm_microphysic to the end of calmufi avoiding rounding pbs
--JVO
File size: 38.1 KB
Line 
1! Copyright 2013-2015,2017 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
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
76  CONTAINS
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)                            :: m,n,p
548    REAL(kind=mm_wp), PARAMETER                 :: fac = 4._mm_wp/3._mm_wp * mm_pi
549
550    ALLOCATE(ft(mm_nle),wth(mm_nle),fdcor(mm_nle))
551
552    !mm_aer_s_flux(:) = 0._mm_wp ; mm_aer_f_flux(:) = 0._mm_wp
553    IF (mm_wsed_m0) THEN
554      ! Spherical particles
555      ! M0
556      call get_weff(mm_m0aer_s,0._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth,fdcor)
557      ft(:) = wth(:) * fdcor(:)  ; mm_m0as_vsed(:) = ft(1:mm_nla)
558      call let_me_fall_in_peace(mm_m0aer_s,-ft,mm_dt,dm0s)
559      ! M3
560      mm_m3as_vsed(:) = ft(1:mm_nla)
561      call let_me_fall_in_peace(mm_m3aer_s,-ft,mm_dt,dm3s)
562      ! get mass flux
563      mm_aer_s_flux(:) = fac * mm_rhoaer * ft(2:) * mm_m3aer_s
564      ! Fractal particles
565      ! M0
566      call get_weff(mm_m0aer_f,0._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth,fdcor)
567      ft(:) = wth(:) * fdcor(:)  ; mm_m0af_vsed(:) = ft(1:mm_nla)
568      call let_me_fall_in_peace(mm_m0aer_f,-ft,mm_dt,dm0f)
569      ! M3
570      mm_m3af_vsed(:) = ft(1:mm_nla)
571      call let_me_fall_in_peace(mm_m3aer_f,-ft,mm_dt,dm3f)
572      ! get mass flux
573      mm_aer_f_flux(:) = fac * mm_rhoaer * ft(2:) * mm_m3aer_f
574    ELSEIF (mm_wsed_m3) THEN
575      ! Spherical particles
576      ! M3
577      call get_weff(mm_m3aer_s,3._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth,fdcor)
578      ft(:) = wth(:) * fdcor(:)  ; mm_m3as_vsed(:) = ft(1:mm_nla)
579      call let_me_fall_in_peace(mm_m3aer_s,-ft,mm_dt,dm3s)
580      ! get mass flux
581      mm_aer_s_flux(:) = fac * mm_rhoaer * ft(2:) * mm_m3aer_s
582      ! M0
583      mm_m0as_vsed(:) = ft(1:mm_nla)
584      call let_me_fall_in_peace(mm_m0aer_s,-ft,mm_dt,dm0s)
585      ! Fractal particles
586      ! M3
587      call get_weff(mm_m3aer_f,3._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth,fdcor)
588      ft(:) = wth(:) * fdcor(:)  ; mm_m3af_vsed(:) = ft(1:mm_nla)
589      call let_me_fall_in_peace(mm_m3aer_f,-ft,mm_dt,dm3f)
590      ! get mass flux
591      mm_aer_f_flux(:) = fac * mm_rhoaer * ft(2:) * mm_m3aer_f
592      ! M0
593      mm_m0af_vsed(:) = ft(1:mm_nla)
594      call let_me_fall_in_peace(mm_m0aer_f,-ft,mm_dt,dm0f)
595    ELSE
596      ! Spherical particles
597      ! M0
598      call get_weff(mm_m0aer_s,0._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth,fdcor)
599      ft(:) = wth(:) * fdcor(:)  ; mm_m0as_vsed(:) = ft(1:mm_nla)
600      call let_me_fall_in_peace(mm_m0aer_s,-ft,mm_dt,dm0s)
601      ! M3
602      call get_weff(mm_m3aer_s,3._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth,fdcor)
603      ft(:) = wth(:) * fdcor(:)  ; mm_m3as_vsed(:) = ft(1:mm_nla)
604      call let_me_fall_in_peace(mm_m3aer_s,-ft,mm_dt,dm3s)
605      ! get mass flux
606      mm_aer_s_flux(:) = fac * mm_rhoaer * ft(2:) * mm_m3aer_s
607      ! Fractal particles
608      ! M0
609      call get_weff(mm_m0aer_f,0._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth,fdcor)
610      ft(:) = wth(:) * fdcor(:)  ; mm_m0af_vsed(:) = ft(1:mm_nla)
611      call let_me_fall_in_peace(mm_m0aer_f,-ft,mm_dt,dm0f)
612      ! M3
613      call get_weff(mm_m3aer_f,3._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth,fdcor)
614      ft(:) = wth(:) * fdcor(:)  ; mm_m3af_vsed(:) = ft(1:mm_nla)
615      call let_me_fall_in_peace(mm_m3aer_f,-ft,mm_dt,dm3f)
616      ! get mass flux
617      mm_aer_f_flux(:) = fac * mm_rhoaer * ft(2:) * mm_m3aer_f
618    ENDIF
619    DEALLOCATE(ft,wth,fdcor)
620    RETURN
621  END SUBROUTINE mm_haze_sedimentation
622
623  SUBROUTINE let_me_fall_in_peace(mk,ft,dt,dmk)
624    !! Compute the tendency of the moment through sedimentation process.
625    !!
626    !!
627    !! The method computes the time evolution of the \(k^{th}\) order moment through sedimentation:
628    !!
629    !! $$ \dfrac{dM_{k}}{dt} = \dfrac{\Phi_{k}}{dz} $$
630    !!
631    !! The equation is resolved using a [Crank-Nicolson algorithm](http://en.wikipedia.org/wiki/Crank-Nicolson_method).
632    !!
633    !! Sedimentation algorithm is quite messy. It appeals to the dark side of the Force and uses evil black magic spells
634    !! from ancient times. It is based on \cite{toon1988b,fiadeiro1977,turco1979} and is an update of the algorithm
635    !! originally implemented in the LMDZ-Titan 2D GCM.
636    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: mk  !! \(k^{th}\) order moment to sediment (in \(m^{k}\)).
637    REAL(kind=mm_wp), INTENT(in)                :: dt  !! Time step (s).
638    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: ft  !! Downward sedimentation flux  (effective velocity of the moment).
639    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dmk !! Tendency of \(k^{th}\) order moment (in \(m^{k}.m^{-3}\)).
640    INTEGER                                 :: i
641    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: as,bs,cs,mko
642    ALLOCATE(as(mm_nla), bs(mm_nla), cs(mm_nla), mko(mm_nla))
643    mko(1:mm_nla) = 0._mm_wp
644    cs(1:mm_nla) = ft(2:mm_nla+1) - mm_dzlay(1:mm_nla)/dt
645    IF (ANY(cs > 0._mm_wp)) THEN
646      ! implicit case
647      as(1:mm_nla) = ft(1:mm_nla)
648      bs(1:mm_nla) = -(ft(2:mm_nle)+mm_dzlay(1:mm_nla)/dt)
649      cs(1:mm_nla) = -mm_dzlay(1:mm_nla)/dt*mk(1:mm_nla)
650      ! (Tri)diagonal matrix inversion
651      mko(1) = cs(1)/bs(1)
652      DO i=2,mm_nla ; mko(i) = (cs(i)-mko(i-1)*as(i))/bs(i) ; ENDDO
653    ELSE
654      ! explicit case
655      as(1:mm_nla)=-mm_dzlay(1:mm_nla)/dt
656      bs(1:mm_nla)=-ft(1:mm_nla)
657      ! boundaries
658      mko(1)=cs(1)*mk(1)/as(1)
659      mko(mm_nla)=(bs(mm_nla)*mk(mm_nla-1)+cs(mm_nla)*mk(mm_nla))/as(mm_nla)
660      ! interior
661      mko(2:mm_nla-1)=(bs(2:mm_nla-1)*mk(1:mm_nla-2) + &
662                       cs(2:mm_nla-1)*mk(2:mm_nla-1)   &
663                      )/as(2:mm_nla-1)
664    ENDIF
665    dmk = mko - mk
666    DEALLOCATE(as,bs,cs,mko)
667    RETURN
668  END SUBROUTINE let_me_fall_in_peace
669
670  SUBROUTINE get_weff(mk,k,df,rc,dt,afun,wth,corf)
671    !! Get the effective settling velocity for aerosols moments.
672    !!
673    !! This method computes the effective settling velocity of the \(k^{th}\) order moment of aerosol
674    !! tracers. The basic settling velocity (\(v^{eff}_{M_{k}}\)) is computed using the following
675    !! equation:
676    !!
677    !! $$
678    !! \begin{eqnarray*}
679    !! \Phi^{sed}_{M_{k}} &=& \int_{0}^{\infty} n(r) r^{k} \times w(r) dr
680    !!                     == M_{k}  \times v^{eff}_{M_{k}} \\
681    !!     v^{eff}_{M_{k} &=& \dfrac{2 \rho g r_{c}^{\dfrac{3D_{f}-3}{D_{f}}}}
682    !!     {r_{m}^{D_{f}-3}/D_{f}} \times \alpha(k)} \times \left( \alpha \right(
683    !!     \frac{D_{f}(k+3)-3}{D_{f}}\left) + \dfrac{A_{kn}\lambda_{g}}{r_{c}^{
684    !!     3/D_{f}}} \alpha \right( \frac{D_{f}(k+3)-6}{D_{f}}\left)\left)
685    !! \end{eqnarray*}
686    !! $$
687    !!
688    !! \(v^{eff}_{M_{k}\) is then corrected to reduce numerical diffusion of the sedimentation algorithm
689    !! as defined in \cite{toon1988b}.
690    !!
691    !! @warning
692    !! Both __df__, __rc__ and __afun__ must be consistent with each other otherwise wrong values will
693    !! be computed.
694    REAL(kind=mm_wp), INTENT(in), DIMENSION(mm_nla)            :: mk
695      !! Moment of order __k__ (\(m^{k}.m^{-3}\)) at each layer.
696    REAL(kind=mm_wp), INTENT(in), DIMENSION(mm_nla)            :: rc
697      !! Characteristic radius associated to the moment at each layer.
698    REAL(kind=mm_wp), INTENT(in)                               :: k
699      !! The order of the moment.
700    REAL(kind=mm_wp), INTENT(in)                               :: df
701      !! Fractal dimension of the aersols.
702    REAL(kind=mm_wp), INTENT(in)                               :: dt
703      !! Time step (s).
704    REAL(kind=mm_wp), INTENT(out), DIMENSION(mm_nle)           :: wth
705      !! Theoretical Settling velocity at each vertical __levels__ (\( wth \times corf = weff\)).
706    REAL(kind=mm_wp), INTENT(out), DIMENSION(mm_nle), OPTIONAL :: corf
707      !! _Fiadero_ correction factor applied to the theoretical settling velocity at each vertical __levels__.
708    INTERFACE
709      FUNCTION afun(order)
710        !! Inter-moment relation function (see [[mm_interfaces(module):mm_alpha_s(interface)]]).
711        IMPORT mm_wp
712        REAL(kind=mm_wp), INTENT(in) :: order !! Order of the moment.
713        REAL(kind=mm_wp) :: afun              !! Alpha value.
714      END FUNCTION afun
715    END INTERFACE
716    INTEGER                             :: i
717    REAL(kind=mm_wp)                    :: af1,af2,ar1,ar2
718    REAL(kind=mm_wp)                    :: csto,cslf,ratio,wdt,dzb
719    REAL(kind=mm_wp)                    :: rb2ra
720    REAL(kind=mm_wp), DIMENSION(mm_nle) :: zcorf
721    ! ------------------
722   
723    wth(:) = 0._mm_wp ; zcorf(:) = 1._mm_wp
724
725    ar1 = (3._mm_wp*df -3._mm_wp)/df    ; ar2 = (3._mm_wp*df -6._mm_wp)/df
726    af1 = (df*(k+3._mm_wp)-3._mm_wp)/df ; af2 = (df*(k+3._mm_wp)-6._mm_wp)/df
727    rb2ra = mm_rm**((df-3._mm_wp)/df)
728    DO i=2,mm_nla
729      IF (rc(i-1) <= 0._mm_wp) CYCLE
730      dzb = (mm_dzlay(i)+mm_dzlay(i-1))/2._mm_wp
731      csto = 2._mm_wp*mm_rhoaer*mm_effg(mm_zlev(i))/(9._mm_wp*mm_eta_g(mm_btemp(i)))
732      cslf = mm_akn * mm_lambda_g(mm_btemp(i),mm_plev(i))
733      wth(i) = - csto/(rb2ra*afun(k)) * (rc(i-1)**ar1 * afun(af1) + cslf/rb2ra * rc(i-1)**ar2 * afun(af2))
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
806  END SUBROUTINE mm_haze_production
807
808END MODULE MM_HAZE
Note: See TracBrowser for help on using the repository browser.