Context Navigation

calc_ysat.F90 @ 3026

Last change on this file since 3026 was 2326, checked in by jvatant, 5 years ago
Update Titan reference photochemistry (reaction constants,branching ratios, condensation rates) according to Vuitton et al 2019. --JVO
File size: 7.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	! ==============================================================================
16
17
18	!-----------------------------------------------------------------------
19	! Declarations:
20	! -------------
21
22	USE comchem_h, only: nkim, cnames
23
24	IMPLICIT NONE
25
26	! Arguments :
27	! -----------
28	INTEGER, INTENT(IN) :: nlon ! # of grid points in the chunk
29	INTEGER, INTENT(IN) :: nlay ! # of vertical layes
30	REAL, DIMENSION(nlon,nlay), INTENT(IN) :: press ! Mid-layers pressure (Pa)
31	REAL, DIMENSION(nlon,nlay), INTENT(IN) :: temp ! Mid-layers temperature (K)
32	REAL, DIMENSION(nlon,nlay,nkim), INTENT(OUT) :: ysat ! Saturation profiles (mol/mol)
33
34	! Local variables :
35	! -----------------
36	INTEGER :: ic
37	REAL,DIMENSION(:,:), ALLOCATABLE :: x
38	! -------------------------------------------------------------------------
39
40	ALLOCATE(x(nlon,nlay))
41
42	! By default, ysat=1, meaning we do not condense
43	ysat(:,:,:) = 1.0
44
45	! Main loop
46
47	do ic=1,nkim
48
49	if(trim(cnames(ic)).eq."CH4") then
50
51	where (temp(:,:).lt.90.65)
52	ysat(:,:,ic) = 10.0**(4.42507e0 - ( ( ( 1165560.7e0 / temp(:,:) &
53	- 115352.19e0 ) / temp(:,:) + 4055.6016e0 ) / temp(:,:) &
54	+ 453.92414e0 ) / temp(:,:) ) / press(:,:) * 1013.25e0
55	elsewhere
56	ysat(:,:,ic) = 10.0**(3.901408e0 - ( ( 154567.02e0 / temp(:,:) - 1598.8512e0 ) &
57	/ temp(:,:) + 437.54809e0 ) / temp(:,:)) / press(:,:) * 1013.25e0
58	endwhere
59
60	! Forcing CH4 to 1.4% minimum
61	where (ysat(:,:,ic).lt.0.014) ysat(:,:,ic) = 0.014
62
63	else if(trim(cnames(ic)).eq."C2H2") then
64
65	! ysat(:,:,ic) = 10.0**(6.09748e0-1644.1e0/temp(:,:)+7.42346e0 &
66	! * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
67
68	ysat(:,:,ic) = 1.E5 * exp(13.4-2536./temp(:,:)) / press(:,:) ! Fray and Schmidt (2009)
69
70	else if(trim(cnames(ic)).eq."C2H4") then
71
72	where (temp(:,:).lt.89.0) ! not far from Fray and Schmidt, 2009
73	ysat(:,:,ic) = 10.0**(1.5477e0 + (1.0e0/temp(:,:) - 0.011e0) &
74	* (16537.0e0*(1.0e0/temp(:,:) - 0.011e0) - 1038.1e0)) &
75	/ press(:,:) * 1.01325e0 / 760.0
76	elsewhere (temp(:,:).lt.104.0) ! not far from Fray and Schmidt, 2009
77	ysat(:,:,ic) = 10.0**(8.724e0 - 901.6e0/(temp(:,:) - 2.555e0) ) &
78	/ press(:,:) * 1013.25e0 / 760.0
79	elsewhere (temp(:,:).lt.120.0)
80	ysat(:,:,ic) = 10.0**(50.79e0 - 1703.0e0/temp(:,:) - 17.141e0 &
81	* alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0
82	elsewhere (temp(:,:).lt.155.0)
83	ysat(:,:,ic) = 10.0**(6.74756e0 - 585.0e0/(temp(:,:) - 18.16e0) ) &
84	/ press(:,:) * 1013.25e0 / 760.0
85	endwhere
86
87	else if(trim(cnames(ic)).eq."C2H6") then
88
89	where (temp(:,:).lt.90.)
90	! ysat(:,:,ic) = 10.0**(10.01e0-1085.0e0/(temp(:,:)-0.561e0) ) &
91	! / press(:,:) * 1013.25e0 / 760.0e0
92	ysat(:,:,ic) = 1.E5 * exp ( 15.11 - 2207./temp(:,:) - 24110./(temp(:,:)**2) &
93	+ 7.744E5/(temp(:,:)3) - 1.161E7/(temp(:,:)4) &
94	+ 6.763E7/(temp(:,:)**5) ) / press(:,:) ! Fray and Schmidt (2009)
95	elsewhere
96	ysat(:,:,ic) = 10.0**(5.9366e0 - 1086.17e0/temp(:,:) + 3.83464e0 &
97	* alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
98	endwhere
99
100	else if((trim(cnames(ic)).eq."CH3CCH") .or. (trim(cnames(ic)).eq."CH2CCH2")) then
101
102	ysat(:,:,ic) = 10.0*(2.8808e0 - 4.5e0(249.9e0 - temp(:,:)) &
103	/ (1.15e0temp(:,:) - 37.485e0) ) / press(:,:) 1013.25e0 / 760.0e0
104
105	else if(trim(cnames(ic)).eq."C3H6") then
106
107	ysat(:,:,ic) = 10.0*(7.4463e0 - 1028.5654e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
108
109	else if(trim(cnames(ic)).eq."C3H8") then
110
111	ysat(:,:,ic) = 10.0*(7.217e0 - 994.30251e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
112
113	else if((trim(cnames(ic)).eq."C4H2").or.(trim(cnames(ic)).eq."C4H2s")) then
114
115	ysat(:,:,ic) = 10.0**(96.26781e0 - 4651.872e0/temp(:,:) - 31.68595e0 &
116	* alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
117
118	else if(trim(cnames(ic)).eq."C4H4") then
119
120	ysat(:,:,ic)= 1.0e3 * exp(9.3898e0 - 2203.57/(temp(:,:)-43.15e0) ) / press(:,:)
121
122	else if(trim(cnames(ic)).eq."C4H6") then
123
124	ysat(:,:,ic)= 10.0*(2.8808e0 - 4.6e0(262.3e0 - temp(:,:)) &
125	/(1.15e0temp(:,:) - 39.345e0) ) / press(:,:) 1013.25e0 / 760.0e0
126
127	else if(trim(cnames(ic)).eq."C4H10") then
128
129	ysat(:,:,ic)= 10.0*(8.446e0 - 1461.2e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
130
131	else if(trim(cnames(ic)).eq."C6H2") then
132
133	ysat(:,:,ic)= 10.0*(4.666e0 - 4956e0/temp(:,:) + 25.845e0 &
134	alog10(1.0e3/temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
135
136	else if(trim(cnames(ic)).eq."C8H2") then
137
138	ysat(:,:,ic)= 10.0*(3.95e0 - 6613e0/temp(:,:) + 35.055e0 &
139	alog10(1.0e3/temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
140
141	else if(trim(cnames(ic)).eq."AC6H6") then
142
143	!x = 1.0e0 - temp(:,:) / 562.2e0
144	!ysat(:,:,ic)= 48.9e3 * exp( ( 1.33213 * x*1.5 - 6.98273 x &
145	! - x*3 (2.62863 + 3.33399 * x*3) ) 562.2e0/temp(:,:) ) / press(:,:)
146
147	ysat(:,:,ic) = 1.E5 * exp (17.35-5663./temp(:,:)) / press(:,:) ! Fray and Schmidt (2009)
148
149	else if(trim(cnames(ic)).eq."HCN") then
150
151	!ysat(:,:,ic)= 10.0*(8.6165e0 - 1516.5e0/(temp(:,:) - 26.2e0) ) / press(:,:) 1013.25e0 / 760.0e0
152
153	ysat(:,:,ic) = 1.E5 * exp ( 13.93 - 3624./temp(:,:) - 1.325E5/(temp(:,:)**2) &
154	+ 6.314E6/(temp(:,:)3) - 1.128E8/(temp(:,:)4) ) / press(:,:) ! Fray and Schmidt (2009)
155
156	else if(trim(cnames(ic)).eq."CH3CN") then
157
158	ysat(:,:,ic)= 10.0*(8.458e0 - 1911.7e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
159
160	else if(trim(cnames(ic)).eq."C2H3CN") then
161
162	ysat(:,:,ic)= 10.0*(9.3051e0 - 2782.21/(temp(:,:) - 51.15e0) ) / press(:,:) 1013.25e0 / 760.0e0
163
164	else if(trim(cnames(ic)).eq."NCCN") then
165
166	ysat(:,:,ic)= 10.0*(7.454e0 - 1832e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
167
168	else if(trim(cnames(ic)).eq."HC3N") then
169
170	ysat(:,:,ic)= 10.0*(7.7446e0 - 1453.5609e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
171
172	else if(trim(cnames(ic)).eq."C4N2") then
173
174	ysat(:,:,ic)= 10.0*(8.269e0 - 2155.0e0/temp(:,:) ) / press(:,:) 1013.25e0 / 760.0e0
175
176	endif
177	enddo
178
179	END SUBROUTINE calc_ysat

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

Download in other formats:

Original Format