source: trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_clouds_methods.F90 @ 3985

Last change on this file since 3985 was 3957, checked in by debatzbr, 6 weeks ago

Pluto PCM: Implementation of a microphysical model for clouds in moments.
BBT
File size: 12.9 KB
Line 
1MODULE MP2M_CLOUDS_METHODS
2    !============================================================================
3    !
4    !     Purpose
5    !     -------
6    !     Cloud model miscellaneous methods module.
7    !
8    !     The module contains miscellaneous methods used in the clouds of the model.
9    !     The module contains four interfaces (8 methods):
10    !        - mm_sigX   | Compute surface tension
11    !        - mm_fshape | Compute shape factor
12    !        - mm_qsatX  | Compute saturation molar mixing ratio
13    !        - mm_LheatX | Compute latent heat released
14    !
15    !     Authors
16    !     -------
17    !     B. de Batz de Trenquelléon (10/2025)
18    !
19    !============================================================================
20
21    USE MP2M_MPREC
22    USE MP2M_GLOBALS
23    USE LINT_DSET
24    USE LINT_LOCATORS
25    IMPLICIT NONE
26
27    PRIVATE
28
29    PUBLIC :: mm_sigX, mm_fshape, mm_qsatX, mm_LheatX
30
31    !============================================================================
32    !                             INTERFACES
33    !============================================================================
34
35    !! Interface to surface tension computation functions.
36    !! The method computes the surface tension of a given species at given temperature(s).
37    INTERFACE mm_sigX
38      MODULE PROCEDURE sigX_sc,sigX_ve
39    END INTERFACE mm_sigX
40
41    !! Interface to shape factor computation functions.
42    !! The method computes the shape factor for the heterogeneous nucleation.
43    INTERFACE mm_fshape
44      MODULE PROCEDURE fshape_sc,fshape_ve
45    END INTERFACE mm_fshape
46
47    !! Interface to saturation molar mixing ratio computation functions.
48    !! The method computes the molar mixing ratio at saturation of a given species at given temperature(s)
49    !! and pressure level(s).
50    INTERFACE mm_qsatX
51      MODULE PROCEDURE qsatX_sc,qsatX_ve
52    END INTERFACE mm_qsatX
53   
54    !! Interface to latent heat computation functions.
55    !! The method computes the latent heat released of a given species at given temperature(s).
56    INTERFACE mm_LheatX
57      MODULE PROCEDURE LheatX_sc,LheatX_ve
58    END INTERFACE mm_LheatX
59
60    CONTAINS
61
62    !============================================================================
63    !                   CLOUD CONDENSATION NUCLEI METHODS
64    !============================================================================
65
66    ! FUNCTION mm_sigX(temp,xESP):
67    ! xESP must always be given as a scalar. If temp is given as a vector, then the method
68    ! computes the result for all the temperatures and returns a vector of same size than temp.
69    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
70    FUNCTION sigX_sc(temp,xESP) RESULT(res)
71        !! Get the surface tension between a given species and the air (scalar).
72        !! The method computes the surface tension equation as given in Reid et al. (1986) p. 637 (eq. 12-3.6).
73        !!
74        REAL(kind=mm_wp), INTENT(in) :: temp ! Temperature (K).
75        TYPE(mm_esp), INTENT(in)     :: xESP ! Specie properties.
76        REAL(kind=mm_wp)             :: res  ! Surface tension (N.m-1).
77
78        ! Local variables:
79        REAL(kind=mm_wp) :: Tr, Tbr, sig0, sig
80
81        Tr = MIN(temp/xESP%Tc,0.99_mm_wp)
82        Tbr = xESP%Tb/xESP%Tc
83
84        sig0 = 0.1196_mm_wp*(1._mm_wp+(Tbr*dlog(xESP%pc/1.01325_mm_wp))/(1._mm_wp-Tbr))-0.279_mm_wp
85
86        sig = xESP%pc**(2._mm_wp/3._mm_wp) * xESP%Tc**(1._mm_wp/3._mm_wp) * sig0 * (1._mm_wp-Tr)**(11._mm_wp/9._mm_wp)
87
88        ! Convertion (dyn/cm) --> (N/m):
89        res = sig*1e-3_mm_wp
90        RETURN
91    END FUNCTION sigX_sc
92
93    FUNCTION sigX_ve(temp,xESP) RESULT(res)
94        !! Get the surface tension between a given species and the air (vector).
95        !! The method computes the surface tension equation as given in Reid et al. (1986) p. 637 (eq. 12-3.6).
96        !!
97        REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp ! Temperatures (K).
98        TYPE(mm_esp), INTENT(in)                   :: xESP ! Specie properties.
99        REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: res  ! Surface tensions (N.m-1).
100
101        ! Local variables:
102        INTEGER          :: i
103        REAL(kind=mm_wp) :: Tr, Tbr, sig0, sig
104
105        Tbr = xESP%Tb/xESP%Tc
106
107        sig0 = 0.1196_mm_wp*(1._mm_wp+(Tbr*dlog(xESP%pc/1.01325_mm_wp))/(1._mm_wp-Tbr))-0.279_mm_wp
108
109        DO i = 1, SIZE(temp)
110          Tr  = MIN(temp(i)/xESP%Tc,0.99_mm_wp)
111          sig = xESP%pc**(2._mm_wp/3._mm_wp) * xESP%Tc**(1._mm_wp/3._mm_wp) * sig0 * (1._mm_wp-Tr)**(11._mm_wp/9._mm_wp)
112         
113          ! Convertion (dyn/cm) --> (N/m):
114          res(i) = sig*1e-3_mm_wp
115        ENDDO
116        RETURN
117    END FUNCTION sigX_ve
118
119
120    ! FUNCTION mm_fshape(m,x):
121    ! Where m is cosine of the contact angle and x the curvature radius. m must always be
122    ! given as a scalar. If x is given as a vector, then the method compute the result for each
123    ! value of x and and returns a vector of same size than x.
124    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
125    FUNCTION fshape_sc(m,x) RESULT(res)
126        !! Get the shape factor of a ccn (scalar).
127        !! The method computes the shape factor for the heterogeneous nucleation on a fractal particle.
128        !! Details about the shape factor can be found in Fletcher et al. (1958).
129        !!
130        REAL(kind=mm_wp), INTENT(in) :: m   ! Cosine of the contact angle.
131        REAL(kind=mm_wp), INTENT(in) :: x   ! Curvature radius (r_particle/r*).
132        REAL(kind=mm_wp)             :: res ! Shape factor value.
133
134        ! Local variables:
135        REAL(kind=mm_wp) :: phi, a, b, c
136
137        IF (x > 3000._mm_wp) THEN
138          res = ((2._mm_wp+m) * (1._mm_wp-m)**2) / 4._mm_wp
139        ELSE
140          phi = dsqrt(1._mm_wp-2._mm_wp*m*x+x**2)
141          a = 1._mm_wp + ((1._mm_wp-m*x)/phi)**3
142          b = (x**3) * (2._mm_wp - 3._mm_wp*(x-m)/phi + ((x-m)/phi)**3)
143          c = 3._mm_wp*m*(x**2) * ((x-m)/phi-1._mm_wp)
144          res = 0.5_mm_wp*(a + b + c)
145        ENDIF
146        RETURN
147    END FUNCTION fshape_sc
148
149    FUNCTION fshape_ve(m,x) RESULT(res)
150        !! Get the shape factor of a ccn (vector).
151        !! The method computes the shape factor for the heterogeneous nucleation on a fractal particle.
152        !! Details about the shape factor can be found in Fletcher et al. (1958).
153        !!
154        REAL(kind=mm_wp), INTENT(in)               :: m   ! Cosine of the contact angle.
155        REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: x   ! Curvature radii (r_particle/r*).
156        REAL(kind=mm_wp), DIMENSION(SIZE(x))       :: res ! Shape factor value.
157
158        ! Local variables:
159        REAL(kind=mm_wp), DIMENSION(SIZE(x)) :: phi, a, b, c
160
161        WHERE(x > 3000._mm_wp)
162          res = ((2._mm_wp+m) * (1._mm_wp-m)**2) / 4._mm_wp
163        ELSEWHERE
164          phi = dsqrt(1._mm_wp-2._mm_wp*m*x+x**2)
165          a = 1._mm_wp + ((1._mm_wp-m*x)/phi)**3
166          b = (x**3) * (2._mm_wp - 3._mm_wp*(x-m)/phi + ((x-m)/phi)**3)
167          c = 3._mm_wp*m*(x**2) * ((x-m)/phi-1._mm_wp)
168          res = 0.5_mm_wp*(a + b + c)
169        ENDWHERE
170        RETURN
171    END FUNCTION fshape_ve
172
173    !============================================================================
174    !                       CONDENSATION METHODS
175    !============================================================================
176
177    ! FUNCTION mm_qsatX(temp,pres,xESP)
178    ! xESP must always be given as a scalar. If temp and pres  are given as a vector,
179    ! then the method computes the result for each couple of (temperature, pressure) and returns
180    ! a vector of same size than temp.
181    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
182    FUNCTION qsatX_sc(temp,pres,xESP) RESULT(res)
183        !! Get the molar mixing ratio of a given species at saturation (scalar).
184        !! Compute saturation molar mixing ratio for condensable tracers.
185        !!
186        !! @warning:
187        !! The formula depends on the species (and the reference)!
188        !!
189        REAL(kind=mm_wp), INTENT(in) :: temp ! Temperature (K).
190        REAL(kind=mm_wp), INTENT(in) :: pres ! Pressure level (Pa).
191        TYPE(mm_esp), INTENT(in)     :: xESP ! Specie properties.
192        REAL(kind=mm_wp)             :: res  ! Saturation molar mixing ratio (mol/mol).
193
194        ! Local variables:
195        REAL(kind=mm_wp) :: fp, fsat
196
197        fp = (1.0e5 / pres)
198
199        ! C2H2, C6H6, HCN: Fray & Schmitt (2009)
200        !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
201        if ((xESP%name == 'C2H2').OR.(xESP%name == 'C6H6').OR.(xESP%name == 'HCN')) then
202            fsat = xESP%a0_sat + xESP%a1_sat/temp + xESP%a2_sat/temp**2 + xESP%a3_sat/temp**3 + &
203                   xESP%a4_sat/temp**4 + xESP%a5_sat/temp**5 + xESP%a6_sat/temp**6
204
205            fsat = exp(fsat)
206
207        ! C2H6: Dykyj et al. (1999)
208        ! C4H2: Orton et al. (2014)
209        !~~~~~~~~~~~~~~~~~~~~~~~~~~
210        else if ((xESP%name == 'C2H6').OR.(xESP%name == 'C4H2')) then
211            fsat = xESP%a0_sat + xESP%a1_sat/(temp+xESP%a2_sat)
212
213            fsat = 10.**(fsat)
214       
215        ! Otherwise: error
216        !~~~~~~~~~~~~~~~~~
217        else
218            write(*,'(a)') "[qsatX_sc] This is a fatal error..."
219            write(*,'(a)') "Species no found for condensation!"
220            call exit(111)
221        endif
222
223        res = fp * fsat
224        RETURN
225    END FUNCTION qsatX_sc
226
227    FUNCTION qsatX_ve(temp,pres,xESP) RESULT(res)
228        !! Get the molar mixing ratio of a given species at saturation (vector).
229        !! Compute saturation molar mixing ratio for condensable tracers.
230        !!
231        !! @warning:
232        !! The formula depends on the species (and the reference)!
233        !!
234        REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp ! Temperatures (K).
235        REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pres ! Pressure levels (Pa).
236        TYPE(mm_esp), INTENT(in)                   :: xESP ! Specie properties.
237        REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: res  ! Saturation molar mixing ratio (mol/mol).
238
239        REAL(kind=mm_wp), DIMENSION(SIZE(temp)) :: fp, fsat
240
241        fp = (1.0e5 / pres)
242
243        ! C2H2, C6H6, HCN: Fray & Schmitt (2009)
244        !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
245        if ((xESP%name == 'C2H2').OR.(xESP%name == 'C6H6').OR.(xESP%name == 'HCN')) then
246            fsat = xESP%a0_sat + xESP%a1_sat/temp + xESP%a2_sat/temp**2 + xESP%a3_sat/temp**3 + &
247                   xESP%a4_sat/temp**4 + xESP%a5_sat/temp**5 + xESP%a6_sat/temp**6
248
249            fsat = exp(fsat)
250       
251        ! C2H6: Dykyj et al. (1999)
252        ! C4H2: Orton et al. (2014)
253        !~~~~~~~~~~~~~~~~~~~~~~~~~~
254        else if ((xESP%name == 'C2H6').OR.(xESP%name == 'C4H2')) then
255            fsat = xESP%a0_sat + xESP%a1_sat/(temp+xESP%a2_sat)
256
257            fsat = 10.**(fsat)
258       
259        ! Otherwise: error
260        !~~~~~~~~~~~~~~~~~
261        else
262            write(*,'(a)') "[qsatX_ve] This is a fatal error..."
263            write(*,'(a)') "Species no found for condensation!"
264            call exit(111)
265        endif
266
267        res = fp * fsat
268        RETURN
269    END FUNCTION qsatX_ve
270
271
272    ! FUNCTION mm_LheatX(temp,xESP):
273    ! xESP must always be given as a scalar. If temp is given as a vector, then the method
274    ! computes the result for all the temperatures and returns a vector of same size than temp.
275    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
276    FUNCTION LheatX_sc(temp,xESP) RESULT(res)
277        !! Compute latent heat of a given species at given temperature (scalar).
278        !! The method computes the latent heat equation as given in Reid et al. (1986) p. 220 (eq. 7-9.4).
279        !!
280        REAL(kind=mm_wp), INTENT(in) :: temp ! Temperature (K).
281        TYPE(mm_esp), INTENT(in)     :: xESP ! Specie properties.
282        REAL(kind=mm_wp)             :: res  ! Latent heat of given species at given temperature (J.kg-1).
283   
284        ! Local variables:
285        REAL(kind=mm_wp) :: Tr, ftm
286
287        Tr = temp / xESP%Tc
288        ftm = MAX(1._mm_wp - Tr,1.e-3_mm_wp)
289
290        res = (7.08_mm_wp*ftm**0.354_mm_wp + 10.95_mm_wp*xESP%w*ftm**0.456_mm_wp) * mm_rgas * xESP%Tc / xESP%masmol
291        RETURN
292    END FUNCTION LheatX_sc
293
294    FUNCTION LheatX_ve(temp,xESP) RESULT(res)
295        !! Compute latent heat of a given species at given temperature (vector).
296        !! The method computes the latent heat equation as given in Reid et al. (1986) p. 220 (eq. 7-9.4).
297        !!
298        REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp ! Temperatures (K).
299        TYPE(mm_esp), INTENT(in)                   :: xESP ! Specie properties.
300        REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: res  ! Latent heat of given species at given temperatures (J.kg-1).
301
302        ! Local variables:
303        INTEGER :: i
304        REAL(kind=mm_wp) :: Tr, ftm
305
306        DO i=1,SIZE(temp)
307          Tr = temp(i) / xESP%Tc
308          ftm = MAX(1._mm_wp - Tr,1.e-3_mm_wp)
309         
310          res(i) = (7.08_mm_wp*ftm**0.354_mm_wp + 10.95_mm_wp*xESP%w*ftm**0.456_mm_wp) * &
311                    mm_rgas * xESP%Tc / xESP%masmol
312        ENDDO
313        RETURN
314    END FUNCTION LheatX_ve
315
316END MODULE MP2M_CLOUDS_METHODS
Note: See TracBrowser for help on using the repository browser.