Context Navigation

calc_ysat.F90 @ 3497

Last change on this file since 3497 was 3318, checked in by slebonnois, 7 months ago
Titan PCM update : optics + microphysics
File size: 8.8 KB

Line
1	SUBROUTINE 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.15e0temp(:,:) - 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.15e0temp(:,:) - 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
215	END SUBROUTINE calc_ysat

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format