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

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

Addition of the microphysics model in moments.
File size: 9.2 KB
Line 
1MODULE MP2M_METHODS
2    !============================================================================
3    !
4    !     Purpose
5    !     -------
6    !     Model miscellaneous methods module.
7    !
8    !     The module contains miscellaneous methods used in the haze of the model.
9    !     The module contains nine methods:
10    !        - mm_lambda_air
11    !        - mm_eta_air
12    !        - mm_ps2s
13    !        - mm_qmean
14    !        - mm_get_btk
15    !        - mm_get_kco
16    !        - mm_get_kfm
17    !
18    !     Authors
19    !     -------
20    !     B. de Batz de Trenquelléon, J. Burgalat (11/2024)
21    !
22    !============================================================================
23
24    USE MP2M_MPREC
25    USE MP2M_GLOBALS
26    USE LINT_DSET
27    USE LINT_LOCATORS
28    IMPLICIT NONE
29
30    PRIVATE
31
32    PUBLIC  :: mm_lambda_air, mm_eta_air, mm_ps2s, mm_qmean, mm_get_btk, mm_get_kfm, mm_get_kco
33
34
35    CONTAINS
36
37    !============================================================================
38    !                          GENERAL METHODS
39    !============================================================================
40
41    ELEMENTAL FUNCTION mm_lambda_air(T,P) RESULT(res)
42      !! Get the air mean free path at given temperature and pressure.
43      !!
44      REAL(kind=mm_wp), INTENT(in) :: T ! Temperature (K).
45      REAL(kind=mm_wp), INTENT(in) :: P ! Pressure level (Pa).
46      REAL(kind=mm_wp) :: res           ! Air mean free path (m).
47
48      res = (mm_kboltz*T) / (dsqrt(2._mm_wp)*mm_pi*(2._mm_wp*mm_air_rad)**2*P)
49     
50      RETURN
51    END FUNCTION mm_lambda_air
52
53   
54    ELEMENTAL FUNCTION mm_eta_air(T) RESULT (res)
55      !! Get the air dynamical viscosity at a given temperature using Sutherland method.
56      !!
57      REAL(kind=mm_wp), INTENT(in) :: T ! Temperature (K).
58      REAL(kind=mm_wp) :: res           ! Air viscosity at given temperature (Pa.s-1).
59
60      REAL(kind=mm_wp), PARAMETER :: eta0 = 1.74e-5_mm_wp
61      REAL(kind=mm_wp), PARAMETER :: Tsut = 109._mm_wp
62      REAL(kind=mm_wp), PARAMETER :: Tref = 293._mm_wp
63     
64      res = eta0 * dsqrt(T/Tref) * ((1._mm_wp + Tsut/Tref) / (1._mm_wp + Tsut/T))
65     
66      RETURN
67    END FUNCTION mm_eta_air
68
69
70    !============================================================================
71    !                      AEROSOL COAGULATION METHODS
72    !============================================================================
73   
74    FUNCTION mm_ps2s(rcs,k,flow) RESULT(res)
75      !! Get the proportion of aerosols that remains in the spherical mode during SS coagulation.
76      !!
77      !! From __k__ and __flow__ values, the method selects one of the four probability datasets
78      !! in mm_globals(module) module (for instance mm_pco0p) and interpolates linearly probability
79      !! for the given value of __rcs__, __T__ and __P__.
80      !!
81      !! @Warning
82      !! Here, the method assumes the datasets define the probability for __spherical__ particles to
83      !! be transferred in the __fractal__ mode, but returns the proportion of particles that remains
84      !! in the mode (which is expected by MP2M model).
85      !!
86      ! Characteristic radius of the spherical size-distribution (m).
87      REAL(kind=mm_wp), INTENT(in) :: rcs
88      ! Order of the moment (0 or 3 expected).
89      INTEGER, INTENT(in)          :: k
90      ! Flow regime indicator (0: Continuous - Kn << 1, 1: Free-Molecular - Kn >> 1).
91      INTEGER, INTENT(in)          :: flow
92     
93      ! Proportion of spherical particles that stay in the spherical mode during SS coagulation.
94      REAL(kind=mm_wp) :: res
95     
96      ! Local variable.
97      TYPE(dset1d), POINTER :: pp
98     
99      res = 1._mm_wp
100      IF (rcs <= 0.0_mm_wp .OR. .NOT.mm_w_ps2s) RETURN
101     
102      SELECT CASE(k+flow)
103        CASE(0)      ; pp => mm_pco0p ! 0 = 0 + 0 -> M0 / CO
104        CASE(1)      ; pp => mm_pfm0p ! 1 = 0 + 1 -> M0 / FM
105        CASE(3)      ; pp => mm_pco3p ! 3 = 3 + 0 -> M3 / CO
106        CASE(4)      ; pp => mm_pfm3p ! 4 = 3 + 1 -> M3 / FM
107        CASE DEFAULT ; RETURN
108      END SELECT
109
110      IF (.NOT.hdcd_lint_dset(rcs,pp,locate_reg_ext,res)) THEN
111        WRITE(*,'(a)') "mp2m_methods:ps2s_sc: Cannot interpolate transfert probability"
112        call EXIT(10)
113      ELSE
114        ! Sanity check: bound probability value between 0 and 1.
115        res = MAX(0.0_mm_wp,MIN(res,1.0_mm_wp))
116        ! We have interpolated f = 1 - p and we need p !
117        res = 1._mm_wp - res
118      ENDIF
119    END FUNCTION mm_ps2s
120
121
122    FUNCTION mm_qmean(rc1,rc2,order,modes,T) RESULT(res)
123      !! Get the electric correction for coagulation kernel.
124      !!
125      !! The method computes the eletric charging correction to apply to the coagulation
126      !! kernel as a function of the temperature and the characteristic radius of the
127      !! mode involved in the coagulation.
128      !! Here the electric charging correction is computed using linear interpolation from
129      !! pre-tabulated values.
130      !!
131      !! @Warning:
132      !! Modes are referred by a two letters uppercase string with the combination of:
133      !!    - S : spherical mode
134      !!    - F : fractal mode
135      !!
136      REAL(kind=mm_wp), INTENT(in) :: rc1   ! Characteristic radius of the first mode (m).
137      REAL(kind=mm_wp), INTENT(in) :: rc2   ! Characteristic radius of the the second mode (m).
138      INTEGER, INTENT(in)          :: order ! Moment's order (0 or 3 expected).
139      CHARACTER(len=2), INTENT(in) :: modes ! Interaction mode (combination of [S,F]).
140      REAL(kind=mm_wp), INTENT(in) :: T     ! Temperature (K).
141     
142      ! Electric charging correction.
143      REAL(kind=mm_wp) :: res
144
145      ! Local variable.
146      INTEGER          :: chx
147      REAL(kind=mm_wp) :: r_tmp, t_tmp
148     
149      chx = 0
150      IF (.NOT.mm_w_qe) THEN
151        res = 1._mm_wp
152        RETURN
153      ENDIF
154   
155      IF (SCAN(modes(1:1),"sS") /= 0) chx = chx + 1
156      IF (SCAN(modes(2:2),"sS") /= 0) chx = chx + 1
157      IF (SCAN(modes(1:1),"fF") /= 0) chx = chx + 3
158      IF (SCAN(modes(2:2),"fF") /= 0) chx = chx + 3
159     
160      chx = chx + order
161     
162      SELECT CASE(chx)
163        ! M0/SS:
164        CASE(2)
165          res = 1._mm_wp
166        ! M0/SF:
167        CASE(4)     
168          ! Fix max values of input parameters
169          r_tmp = MAX(MIN(log(rc1),mm_qbsf0_e(2,2)),mm_qbsf0_e(2,1))
170          t_tmp = MAX(MIN(T,mm_qbsf0_e(1,2)),mm_qbsf0_e(1,1))
171          ! Interpolates values
172          IF (.NOT.hdcd_lint_dset(t_tmp,r_tmp,mm_qbsf0,locate_reg,res)) THEN
173            WRITE(*,'(a)') "mp2m_methods:mm_qmean: Cannot interpolate mean Qelec"
174            call EXIT(10)
175          ENDIF
176        ! M3/SS:
177        CASE(5)
178          res = 1._mm_wp
179        ! M0/FF:
180        CASE(6)
181          ! Fix max values of input parameters
182          r_tmp = MAX(MIN(log(rc1),mm_qbff0_e(2,2)),mm_qbff0_e(2,1))
183          t_tmp = MAX(MIN(T,mm_qbff0_e(1,2)),mm_qbff0_e(1,1))
184          ! Interpolates values
185          IF (.NOT.hdcd_lint_dset(t_tmp,r_tmp,mm_qbff0,locate_reg,res)) THEN
186            WRITE(*,'(a)') "mp2m_methods:mm_qmean: Cannot interpolate mean Qelec"
187            call EXIT(10)
188          ENDIF
189        ! M3/SF:
190        CASE(7)
191          ! Fix max values of input parameters
192          r_tmp = MAX(MIN(log(rc1),mm_qbsf3_e(2,2)),mm_qbsf3_e(2,1))
193          t_tmp = MAX(MIN(T,mm_qbsf3_e(1,2)),mm_qbsf3_e(1,1))
194          ! Interpolates values
195          IF (.NOT.hdcd_lint_dset(t_tmp,r_tmp,mm_qbsf3,locate_reg,res)) THEN
196            WRITE(*,'(a)') "mp2m_methods:mm_qmean: Cannot interpolate mean Qelec"
197            call EXIT(10)
198          ENDIF
199        ! Anything else:
200        CASE DEFAULT
201          res = 1._mm_wp
202      END SELECT
203     
204      RETURN
205    END FUNCTION mm_qmean
206
207
208    PURE FUNCTION mm_get_btk(t,k) RESULT(res)
209      !! Get the value of the Free-Molecular regime coagulation pre-factor b^t_k.
210      !!
211      !! @Note
212      !! __k__ can only be one of the following value: 0 or 3. __t__ ranges only from 1 to 5.
213      !!
214      INTEGER, INTENT(in) :: t ! Index of the b^t_k coefficient.
215      INTEGER, INTENT(in) :: k ! Moment order of the b^t_k coefficient.
216     
217      ! b^t_k coefficient.
218      REAL(kind=mm_wp) :: res
219
220      ! Sanity check:
221      IF (.NOT.(k == 3 .OR. k == 0)) THEN
222        res = 0._mm_wp
223      ENDIF
224     
225      ! Sanity check:
226      IF (t > 5 .OR. t < 1) THEN
227        res = 0._mm_wp
228      ENDIF
229     
230      IF (k == 0) THEN
231        res = mm_bt0(t)
232
233      ELSE IF (k == 3) THEN
234        res = mm_bt3(t)
235      ENDIF
236     
237      RETURN
238    END FUNCTION mm_get_btk
239
240
241    ELEMENTAL FUNCTION mm_get_kco(T) RESULT(res)
242      !! Get the Continuous regime (Kn << 1) thermodynamics pre-factor of the coagulation kernel.
243      !!
244      REAL(kind=mm_wp), INTENT(in) :: T ! Temperature (K).
245      REAL(kind=mm_wp) :: res           ! Continuous regime thermodynamics pre-factor (m3.s-1).
246
247      res = (2._mm_wp*mm_kboltz*T) / (3._mm_wp*mm_eta_air(T))
248     
249      RETURN
250    END FUNCTION mm_get_kco
251
252
253    ELEMENTAL FUNCTION mm_get_kfm(T) RESULT(res)
254      !! Get the Free-Molecular regime (Kn >> 1) thermodynamics pre-factor of the coagulation kernel.
255      !!
256      REAL(kind=mm_wp), INTENT(in) :: T ! Temperature (K).
257      REAL(kind=mm_wp) :: res           ! Free-Molecular regime thermodynamics pre-factor (m^(5/2).s-1).
258
259      res = (6._mm_wp*mm_kboltz*T / mm_rhoaer)**(0.5_mm_wp)
260     
261      RETURN
262    END FUNCTION mm_get_kfm
263
264END MODULE MP2M_METHODS
Note: See TracBrowser for help on using the repository browser.