Context Navigation

calc_ysat.F90 @ 3094

Last change on this file since 3094 was 3090, checked in by slebonnois, 15 months ago
BdeBatz? : Cleans microphysics and makes few corrections for physics
File size: 8.5 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	! C2H6 :
100	!-------
101	else if(trim(cnames(ic)).eq."C2H6") then
102	where (temp(:,:).lt.90.)
103	!ysat(:,:,ic) = 10.0**(10.01e0-1085.0e0/(temp(:,:)-0.561e0) ) &
104	! / press(:,:) * 1013.25e0 / 760.0e0
105	! Fray and Schmidt (2009)
106	ysat(:,:,ic) = (1.e3 / press(:,:)) * exp ( 15.11 - 2207./temp(:,:) &
107	- 24110./(temp(:,:)2) + 7.744E5/(temp(:,:)3) &
108	- 1.161E7/(temp(:,:)4) + 6.763E7/(temp(:,:)5))
109	elsewhere
110	ysat(:,:,ic) = 10.0**(5.9366e0 - 1086.17e0/temp(:,:) + 3.83464e0 &
111	* alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
112	endwhere
113
114	! CH3CCH :
115	!---------
116	else if((trim(cnames(ic)).eq."CH3CCH") .or. (trim(cnames(ic)).eq."CH2CCH2")) then
117	ysat(:,:,ic) = 10.0*(2.8808e0 - 4.5e0(249.9e0 - temp(:,:)) &
118	/ (1.15e0temp(:,:) - 37.485e0) ) / press(:,:) 1013.25e0 / 760.0e0
119
120	! C3H6 :
121	!-------
122	else if(trim(cnames(ic)).eq."C3H6") then
123	ysat(:,:,ic) = 10.0*(7.4463e0 - 1028.5654e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
124
125	! C3H8 :
126	!-------
127	else if(trim(cnames(ic)).eq."C3H8") then
128	ysat(:,:,ic) = 10.0*(7.217e0 - 994.30251e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
129
130	! C4H2 :
131	!-------
132	else if((trim(cnames(ic)).eq."C4H2").or.(trim(cnames(ic)).eq."C4H2s")) then
133	ysat(:,:,ic) = 10.0**(96.26781e0 - 4651.872e0/temp(:,:) - 31.68595e0 &
134	* alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
135
136	! C4H4 :
137	!-------
138	else if(trim(cnames(ic)).eq."C4H4") then
139	ysat(:,:,ic)= 1.0e3 * exp(9.3898e0 - 2203.57/(temp(:,:)-43.15e0) ) / press(:,:)
140
141	! C4H6 :
142	!-------
143	else if(trim(cnames(ic)).eq."C4H6") then
144	ysat(:,:,ic)= 10.0*(2.8808e0 - 4.6e0(262.3e0 - temp(:,:)) &
145	/(1.15e0temp(:,:) - 39.345e0) ) / press(:,:) 1013.25e0 / 760.0e0
146
147	! C4H10 :
148	!--------
149	else if(trim(cnames(ic)).eq."C4H10") then
150	ysat(:,:,ic)= 10.0*(8.446e0 - 1461.2e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
151
152	! C6H2 :
153	!-------
154	else if(trim(cnames(ic)).eq."C6H2") then
155	ysat(:,:,ic)= 10.0*(4.666e0 - 4956e0/temp(:,:) + 25.845e0 &
156	alog10(1.0e3/temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
157
158	! C8H2 :
159	!-------
160	else if(trim(cnames(ic)).eq."C8H2") then
161	ysat(:,:,ic)= 10.0*(3.95e0 - 6613e0/temp(:,:) + 35.055e0 &
162	alog10(1.0e3/temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
163
164	! C6H6 :
165	!-------
166	else if(trim(cnames(ic)).eq."AC6H6") then
167	!x = 1.0e0 - temp(:,:) / 562.2e0
168	!ysat(:,:,ic)= 48.9e3 * exp( ( 1.33213 * x*1.5 - 6.98273 x &
169	! - x*3 (2.62863 + 3.33399 * x*3) ) 562.2e0/temp(:,:) ) / press(:,:)
170
171	! Fray and Schmidt (2009)
172	ysat(:,:,ic) = (1.e3 / press(:,:)) * exp (17.35-5663./temp(:,:))
173
174	! HCN :
175	!-------
176	else if(trim(cnames(ic)).eq."HCN") then
177	!ysat(:,:,ic)= 10.0*(8.6165e0 - 1516.5e0/(temp(:,:) - 26.2e0) ) / press(:,:) 1013.25e0 / 760.0e0
178
179	! Fray and Schmidt (2009)
180	ysat(:,:,ic) = (1.e3 / press(:,:)) * exp (1.393e1 - 3.624e3/temp(:,:) &
181	- 1.325e5/temp(:,:)2 + 6.314e6/temp(:,:)3 - 1.128e8/temp(:,:)**4)
182
183	! CH3CN :
184	!--------
185	else if(trim(cnames(ic)).eq."CH3CN") then
186	ysat(:,:,ic)= 10.0*(8.458e0 - 1911.7e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
187
188	! C2H3CN :
189	!-------
190	else if(trim(cnames(ic)).eq."C2H3CN") then
191	ysat(:,:,ic)= 10.0*(9.3051e0 - 2782.21/(temp(:,:) - 51.15e0) ) / press(:,:) 1013.25e0 / 760.0e0
192
193	! NCCN :
194	!-------
195	else if(trim(cnames(ic)).eq."NCCN") then
196	ysat(:,:,ic)= 10.0*(7.454e0 - 1832e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
197
198	! HC3N :
199	!-------
200	else if(trim(cnames(ic)).eq."HC3N") then
201	ysat(:,:,ic)= 10.0*(7.7446e0 - 1453.5609e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
202
203	! C4N2 :
204	!-------
205	else if(trim(cnames(ic)).eq."C4N2") then
206	ysat(:,:,ic)= 10.0*(8.269e0 - 2155.0e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
207
208	endif
209	enddo
210
211	END SUBROUTINE calc_ysat

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

Download in other formats:

Original Format