source: trunk/LMDZ.TITAN/libf/phytitan/calc_ysat.F90

Last change on this file was 3318, checked in by slebonnois, 7 months ago

Titan PCM update : optics + microphysics
File size: 8.8 KB
Line 
1SUBROUTINE calc_ysat(nlon,nlay,press,temp,ysat)
2
3  !     ==============================================================================
4  !     Purpose
5  !     -------
6  !     Compute saturation profiles (mol/mol) for chemical tracers
7  !     
8  !     Authors
9  !     -------
10  !     S. Lebonnois
11  !     -> inicondens.F in old Titan, with T,P in planetary average
12  !     J. Vatant d'Ollone
13  !       -> 2017 : Adapt to new physics, move to F90 and compute every grid point
14  !       -> 2018 : Update saturation profiles from recent litterature (Vuitton 2019)
15  !     Modified
16  !     --------
17  !     B. de Batz de Trenquelléon
18  !       -> 2022 : Update saturation profiles and correct the pressure unit
19  !     ==============================================================================
20
21
22  !-----------------------------------------------------------------------
23  !     Declarations:
24  !     -------------
25
26  USE comchem_h, only: nkim, cnames
27
28  IMPLICIT NONE
29
30  !     Arguments :
31  !     -----------
32  INTEGER, INTENT(IN)                           :: nlon  ! # of grid points in the chunk
33  INTEGER, INTENT(IN)                           :: nlay  ! # of vertical layes
34  REAL, DIMENSION(nlon,nlay),      INTENT(IN)   :: press ! Mid-layers pressure (!!! mbar !!!)
35  REAL, DIMENSION(nlon,nlay),      INTENT(IN)   :: temp  ! Mid-layers temperature (K)
36  REAL, DIMENSION(nlon,nlay,nkim), INTENT(OUT)  :: ysat  ! Saturation profiles (mol/mol)
37
38  !     Local variables :
39  !     -----------------
40  INTEGER                           ::  ic
41  REAL,DIMENSION(:,:), ALLOCATABLE  ::  x
42  !     -------------------------------------------------------------------------
43
44  ALLOCATE(x(nlon,nlay))
45
46  ! By default, ysat = 1, meaning we do not condense
47  ysat(:,:,:) = 1.0
48
49  ! Main loop
50  do ic = 1, nkim
51     
52      ! CH4 :
53      !------
54      if(trim(cnames(ic)).eq."CH4") then
55         !where (temp(:,:).lt.90.65)                                     
56         !   ysat(:,:,ic) = 10.0**(4.42507e0 - ( ( ( 1165560.7e0 / temp(:,:)                 &
57         !        -  115352.19e0 ) / temp(:,:) + 4055.6016e0 ) / temp(:,:)                   &
58         !        + 453.92414e0 ) / temp(:,:) ) / press(:,:) * 1013.25e0
59         !elsewhere
60         !   ysat(:,:,ic) = 10.0**(3.901408e0 - ( (154567.02e0 / temp(:,:) - 1598.8512e0)    &
61         !        / temp(:,:) + 437.54809e0 ) / temp(:,:)) / press(:,:) * 1013.25e0
62         !endwhere
63
64         ! Fray and Schmidt (2009)
65         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp(1.051e1 - 1.110e3/temp(:,:)                &
66                        - 4.341e3/temp(:,:)**2 + 1.035e5/temp(:,:)**3 - 7.910e5/temp(:,:)**4)
67
68         ! Forcing CH4 to 1.4% minimum               
69         where (ysat(:,:,ic).lt.0.014) ysat(:,:,ic) = 0.014
70
71      ! C2H2 :
72      !-------
73      else if(trim(cnames(ic)).eq."C2H2") then
74         !ysat(:,:,ic) = 10.0**(6.09748e0-1644.1e0/temp(:,:)+7.42346e0                       &
75         !     * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
76
77         ! Fray and Schmidt (2009)
78         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp(1.340e1 - 2.536e3/temp(:,:))
79
80      ! C2H4 :
81      !-------
82      else if(trim(cnames(ic)).eq."C2H4") then
83         ! not far from Fray and Schmidt (2009)
84         !where (temp(:,:).lt.89.0)
85         !   ysat(:,:,ic) = 10.0**(1.5477e0 + (1.0e0/temp(:,:) - 0.011e0)                     &
86         !         * (16537.0e0*(1.0e0/temp(:,:) - 0.011e0) - 1038.1e0))                      &
87         !         / press(:,:) * 1.01325e0 / 760.0
88         !elsewhere (temp(:,:).lt.104.0)
89         !   ysat(:,:,ic) = 10.0**(8.724e0 - 901.6e0/(temp(:,:) - 2.555e0) )                  &
90         !         / press(:,:) * 1013.25e0 / 760.0
91         !elsewhere (temp(:,:).lt.120.0)
92         !   ysat(:,:,ic) = 10.0**(50.79e0 - 1703.0e0/temp(:,:) - 17.141e0                    &
93         !         * alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0
94         !elsewhere (temp(:,:).lt.155.0)
95         !   ysat(:,:,ic) = 10.0**(6.74756e0 - 585.0e0/(temp(:,:) - 18.16e0) )                &
96         !         / press(:,:) * 1013.25e0 / 760.0
97         !endwhere
98
99         ! Fray and Schmidt (2009)
100         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp(1.540e1 - 2.206e3/temp(:,:)                &
101                         - 1.216e4/temp(:,:)**2 + 2.843e5/temp(:,:)**3 - 2.203e6/temp(:,:)**4)
102     
103      ! C2H6 :
104      !-------
105      else if(trim(cnames(ic)).eq."C2H6") then
106         where (temp(:,:).lt.90.)
107            !ysat(:,:,ic) = 10.0**(10.01e0-1085.0e0/(temp(:,:)-0.561e0) )                    &
108            !      / press(:,:) * 1013.25e0 / 760.0e0
109            ! Fray and Schmidt (2009)
110            ysat(:,:,ic) = (1.e3 / press(:,:)) * exp ( 15.11 - 2207./temp(:,:)               &
111                           - 24110./(temp(:,:)**2) + 7.744E5/(temp(:,:)**3)                  &
112                           - 1.161E7/(temp(:,:)**4) + 6.763E7/(temp(:,:)**5))
113         elsewhere
114            ysat(:,:,ic) = 10.0**(5.9366e0 - 1086.17e0/temp(:,:) + 3.83464e0                 &
115                  * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
116         endwhere
117     
118      ! CH3CCH :
119      !---------
120      else if((trim(cnames(ic)).eq."CH3CCH") .or. (trim(cnames(ic)).eq."CH2CCH2")) then
121         ysat(:,:,ic) = 10.0**(2.8808e0 - 4.5e0*(249.9e0 - temp(:,:))                        &
122             / (1.15e0*temp(:,:) - 37.485e0) ) / press(:,:) * 1013.25e0 / 760.0e0
123
124      ! C3H6 :
125      !-------
126      else if(trim(cnames(ic)).eq."C3H6")  then
127         ysat(:,:,ic) = 10.0**(7.4463e0 - 1028.5654e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
128
129      ! C3H8 :
130      !-------
131     else if(trim(cnames(ic)).eq."C3H8")  then
132        ysat(:,:,ic) = 10.0**(7.217e0 - 994.30251e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
133
134      ! C4H2 :
135      !-------
136      else if((trim(cnames(ic)).eq."C4H2").or.(trim(cnames(ic)).eq."C4H2s")) then
137         ysat(:,:,ic) = 10.0**(96.26781e0 - 4651.872e0/temp(:,:) - 31.68595e0                &
138             * alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
139     
140      ! C4H4 :
141      !-------
142      else if(trim(cnames(ic)).eq."C4H4")  then
143         ysat(:,:,ic)= 1.0e3 * exp(9.3898e0 - 2203.57/(temp(:,:)-43.15e0) ) / press(:,:)
144
145      ! C4H6 :
146      !-------
147      else if(trim(cnames(ic)).eq."C4H6")  then
148         ysat(:,:,ic)= 10.0**(2.8808e0 - 4.6e0*(262.3e0 - temp(:,:))                         &
149             /(1.15e0*temp(:,:) - 39.345e0) ) / press(:,:) * 1013.25e0 / 760.0e0
150
151      ! C4H10 :
152      !--------
153      else if(trim(cnames(ic)).eq."C4H10")  then
154         ysat(:,:,ic)= 10.0**(8.446e0 - 1461.2e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
155
156      ! C6H2 :
157      !-------
158      else if(trim(cnames(ic)).eq."C6H2")  then
159         ysat(:,:,ic)= 10.0**(4.666e0 - 4956e0/temp(:,:) + 25.845e0 *                        &
160             alog10(1.0e3/temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
161
162      ! C8H2 :
163      !-------
164      else if(trim(cnames(ic)).eq."C8H2")  then
165         ysat(:,:,ic)= 10.0**(3.95e0 - 6613e0/temp(:,:) + 35.055e0 *                         &
166             alog10(1.0e3/temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
167
168      ! C6H6 :
169      !-------
170      else if(trim(cnames(ic)).eq."AC6H6")  then
171         !x = 1.0e0 - temp(:,:) / 562.2e0
172         !ysat(:,:,ic)= 48.9e3 * exp( ( 1.33213 * x**1.5 - 6.98273 * x                        &
173         !     - x**3 * (2.62863 + 3.33399 * x**3) ) * 562.2e0/temp(:,:) ) / press(:,:)
174
175         ! Fray and Schmidt (2009)
176         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp (17.35-5663./temp(:,:))
177
178      ! HCN :
179      !-------
180      else if(trim(cnames(ic)).eq."HCN")  then
181         !ysat(:,:,ic)= 10.0**(8.6165e0 - 1516.5e0/(temp(:,:) - 26.2e0) ) / press(:,:) * 1013.25e0 / 760.0e0
182         
183         ! Fray and Schmidt (2009)
184         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp (1.393e1 - 3.624e3/temp(:,:)               &
185                        - 1.325e5/temp(:,:)**2 + 6.314e6/temp(:,:)**3 - 1.128e8/temp(:,:)**4)
186
187      ! CH3CN :
188      !--------
189      else if(trim(cnames(ic)).eq."CH3CN")  then
190         ysat(:,:,ic)= 10.0**(8.458e0 - 1911.7e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
191   
192      ! C2H3CN :
193      !-------
194      else if(trim(cnames(ic)).eq."C2H3CN")  then
195         ysat(:,:,ic)= 10.0**(9.3051e0 - 2782.21/(temp(:,:) - 51.15e0) ) / press(:,:) * 1013.25e0 / 760.0e0
196
197      ! NCCN :
198      !-------
199      else if(trim(cnames(ic)).eq."NCCN")  then
200         ysat(:,:,ic)=  10.0**(7.454e0 - 1832e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
201   
202      ! HC3N :
203      !-------
204      else if(trim(cnames(ic)).eq."HC3N")  then
205         ysat(:,:,ic)= 10.0**(7.7446e0 - 1453.5609e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
206   
207      ! C4N2 :
208      !-------
209      else if(trim(cnames(ic)).eq."C4N2")  then
210         ysat(:,:,ic)= 10.0**(8.269e0 - 2155.0e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
211
212      endif
213  enddo
214
215END SUBROUTINE calc_ysat
Note: See TracBrowser for help on using the repository browser.