Context Navigation

getk_V09.F @ 426

Last change on this file since 426 was 414, checked in by aslmd, 14 years ago

LMDZ.MARS : NEW NLTE MODEL FROM GRANADA AMIGOS

23/11/11 == FGG + MALV

New parameterization of the NLTE 15 micron cooling. The old parameterization is kept as an option, including or not variable atomic oxygen concentration. A new flag is introduced in callphys.def, nltemodel, to select which parameterization wants to be used (new one, old one with variable [O], or old one with fixed [O], see below). Includes many new subroutines and commons in phymars. Some existing routines are also modified:

-physiq.F. Call to the new subroutine NLTE_leedat in first call. Call to nltecool modified to allow for variable atomic oxygen. Depending on the value of nltemodel, the new subroutine NLTEdlvr09_TCOOL is called instead of nltecool.

-inifis.F. Reading of nltemodel is added.

-callkeys.h Declaration of nltemodel is added.

The following lines should be added to callphys.def (ideally after setting callnlte):

# NLTE 15um scheme to use.
# 0-> Old scheme, static oxygen
# 1-> Old scheme, dynamic oxygen
# 2-> New scheme
nltemodel = 2

A new directory, NLTEDAT, has to be included in datagcm.

Improvements into NLTE NIR heating parameterization to take into account variability with O/CO2 ratio and SZA. A new subroutine, NIR_leedat.F, and a new common, NIRdata.h, are included. A new flag, nircorr, is added in callphys.def, to include or not these corrections. The following files are modified:

-nirco2abs.F: nq and pq are added in the arguments. The corrections factors are interpolated to the GCM grid and included in the clculation. A new subroutine, interpnir, is added at the end of the file.

-physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs

-inifis: Reading new flag nircorr.

-callkeys.h: Declaration of nircorr.

The following lines have to be added to callphys.def (ideally after callnirco2):

# NIR NLTE correction for variable SZA and O/CO2?
# matters only if callnirco2=T
# 0-> no correction
# 1-> include correction
nircorr=1

A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm.

Small changes to the molecular diffusion scheme to fix the number of species considered, to avoid problems when compiling with more than 15 tracers (for example, when CH4 is included). Modified subroutines: aeronomars/moldiff.F and aeronomars/moldiffcoeff.F

File size: 40.8 KB

Line
1	c***********************************************************************
2
3	subroutine getk (tt)
4
5	c jan 98 malv version for mz1d. copied from solar10/getk.f
6	c jul 2011 malv+fgg adapted to LMD-MGCM
7	c***********************************************************************
8
9	implicit none
10
11	include 'nltedefs.h'
12	include 'tcr_15um.h'
13	include 'nlte_data.h'
14	include 'nlte_rates.h'
15
16	c arguments
17	real tt ! i. temperature
18
19	!! local variables:
20	real*8 k1x, k7x,k7xa,k7xb
21	real*8 k3x, k3xaa,k3xac,k3xab,k3xbb,k3xba,k3xbc
22	real*8 k3xco2,k3xn2,k3xco, k6x,k6x1,k6x2
23	real*8 k20x,k20xa,k20xb,k20xc
24	real*8 k19xca,k19xcb,k19xcc
25	real*8 k19xba,k19xbb,k19xbc
26	real*8 k19xaa,k19xab,k19xac
27	real*8 k21x,k21xa,k21xb,k21xc
28	real*8 anu, factor , k21xa_626
29	real tt1,tt2, de
30	integer i
31
32	c***********************************************************************
33
34	co2i(001) + n2 ---> co2i + n2(1) k1 (not considered in the model
35	c k1(i) i = 1 --- 626
36	! 2 --- 636
37	! 3 --- 628
38	! 4 --- 627
39
40	! k1x = 5.d-13 * sqrt(300.d0/tt)
41	! do i=1,nisot
42	! k1(i) = k1x * rf1
43	! k1p(i) = k1(i) * exp( -ee/tt 1.d0 (-nun2+nu(i,4)) )
44	! end do
45
46	co2(001) + co2i ---> co2 + co2i(001) k2 (vv1 in table 4, paper i)
47	c k2i i = x --- 628
48	! y --- 636
49	! z --- 627
50
51	k2a = 6.8d-12 * sqrt(tt) ! delta(e)< 42 cm-1
52	k2b = 3.6d-11 * sqrt(tt) * exp(-65.6557/26.3) ! > 42 cm-1 see table 3
53	k2a = k2a * rf2desac
54	k2b = k2b * rf2desac
55
56	k2x = 6.8d-12 * sqrt(tt)
57	k2y = 3.6d-11 * sqrt(tt) * exp(-65.6557/26.3)
58	k2z = 6.8d-12 * sqrt(tt)
59	k2x = k2x * rf2iso
60	k2y = k2y * rf2iso
61	k2z = k2z * rf2iso
62	k2xp = k2x * exp( dble( -ee/tt * (nu(1,4)-nu(2,4)) ) )
63	k2yp = k2y * exp( dble( -ee/tt * (nu(1,4)-nu(3,4)) ) )
64	k2zp = k2z * exp( dble( -ee/tt * (nu(1,4)-nu(4,4)) ) )
65
66	! these are vt1 in table 3, paper i
67	co2i(001) + m ---> co2i(0v0) + m k3
68	co2i(001) + o3p ---> co2i(0v0) + o3p k7
69	c k3vm(i) m = a --- co2 v = a --- 3 i = 1,2,3,4
70	! b --- n2,o2 b --- 2
71	! c --- co
72	c k7v(i) v = a --- 3 i = 1,2,3,4
73	! b --- 2
74
75	k7x = 2.0d-13*sqrt(tt/300.d0)
76	k7x = k7x * rf7
77
78	if (iopt3.eq.0) then
79
80	k3x = 2.2d-15 + 1.14d-10 * exp( dble(-76.75/tt**(1.d0/3.d0)) )
81	k3xaa = 3.0d-15 + 1.72d-10 * exp( dble(-76.75/tt**(1.d0/3.d0)) )
82	k3xac = 2.2d-15 + 9.66d-11 * exp( dble(-76.75/tt**(1.d0/3.d0)) )
83	k3x = k3x * rf3
84	k3xaa = k3xaa * rf3
85	k3xac = k3xac * rf3
86
87	tt1=220.0 !500.0 !220.0
88	tt2=250.0 !250.0
89	if(tt.le.tt1)then
90	k3xab = k3x
91	k3xbb = 0.d0
92	k7xa=k7x
93	k7xb=0.d0
94	else if(tt.gt.tt2)then
95	k3xab = 0.d0
96	k3xbb = k3x
97	k7xa=0.d0
98	k7xb=k7x
99	else
100	k3xab = k3x/30.d0*(tt2-tt)
101	k3xbb = k3x/30.d0*(tt-tt1)
102	k7xa=k7x/30.d0*(tt2-tt)
103	k7xb=k7x/30.d0*(tt-tt1)
104	end if
105	k3xba = k3xbb
106	k3xbc = k3xbb
107
108	elseif (iopt3.gt.0) then ! bauer et. al., 1987
109
110	if (tt.ge.190. .and. tt.le.250.) then
111	factor = 0.9d0 + dble( (0.1-0.9)/(250.-190.) * (tt-190.) )
112	elseif (tt.lt.190.) then
113	factor = 0.9d0
114	elseif (tt.gt.250.) then
115	factor = 0.1d0
116	end if
117
118	k3xn2 = 2.2d-15 + 1.14d-10 * exp(dble( -76.75/tt**(1.d0/3.d0)) )
119	k3xn2 = k3xn2 * rf3
120	k3xab = k3xn2 * factor
121	k3xbb = k3xn2 * (1.-factor)
122
123	k7xa = k7x * factor
124	k7xb = k7x * (1.-factor)
125
126	if (iopt3.eq.1) then
127
128	if (tt.le.148.8) then
129	k3xco2 = 13.8 - 0.1807 * (tt-140.)
130	elseif (tt.ge.148.8 .and. tt.le.159.6) then
131	k3xco2 = 12.21 - 0.1787 * (tt-148.8)
132	elseif (tt.ge.159.6 .and. tt.le.171.0) then
133	k3xco2 = 10.28 - 0.04035 * (tt-159.6)
134	elseif (tt.ge.171.0 .and. tt.le.186.4) then
135	k3xco2 = 9.82 - 0.027273 * (tt-171.0)
136	elseif (tt.ge.186.4 .and. tt.le.244.1) then
137	k3xco2 = 9.4 + 0.002253 * (tt-186.4)
138	elseif (tt.ge.244.1 .and. tt.le.300) then
139	k3xco2 = 9.53 + 0.02129 * (tt-244.1)
140	elseif (tt.ge.300 .and. tt.le.336.1) then
141	k3xco2 = 10.72 + 0.052632 * (tt-300)
142	elseif (tt.ge.336.1 .and. tt.le.397.0) then
143	k3xco2 = 12.62 + 0.0844 * (tt-336.1)
144	elseif (tt.ge.397.0 .and. tt.le.523.4) then
145	k3xco2 = 17.76 + 0.2615 * (tt-397.)
146	end if
147	k3xco2 = k3xco2 * 1.d-15 * rf3
148	k3xaa = 0.82d0 * k3xco2
149	k3xba = 0.18d0 * k3xco2
150
151	if (tt.le.163.) then
152	k3xco = 10.58 - 0.093 * (tt-140)
153	elseif (tt.ge.163. .and. tt.le.180.) then
154	k3xco = 8.44 - 0.05353 * (tt-163.)
155	elseif (tt.ge.180. .and. tt.le.196.) then
156	k3xco = 7.53 - 0.034375 * (tt-180.)
157	elseif (tt.ge.196. .and. tt.le.248.) then
158	k3xco = 6.98 - 0.0108 * (tt-196.)
159	elseif (tt.ge.248. .and. tt.le.301.) then
160	k3xco = 6.42 + 0.01415 * (tt-248.)
161	elseif (tt.ge.301. .and. tt.le.353.) then
162	k3xco = 7.17 + 0.02865 * (tt-301.)
163	end if
164	k3xac = k3xco * 1.d-15 * rf3
165	k3xbc = 0.d0
166
167	elseif (iopt3.eq.2) then ! revision for the papers (feb 93)
168
169	k3xco2 = 7.3d-14 * exp( -850.3/tt + 86523/tt**2.d0 )
170	k3xco2 = k3xco2 * rf3
171	k3xaa = 0.82d0 * k3xco2
172	k3xba = 0.18d0 * k3xco2
173
174	k3xco = 1.7d-14 * exp( -448.3/tt + 53636/tt**2.d0 )
175	k3xac = k3xco * rf3
176	k3xbc = 0.d0
177
178	end if
179
180	end if
181
182	do i=1,nisot
183	k3aa(i) = k3xaa
184	k3ba(i) = k3xba
185	k3ab(i) = k3xab
186	k3bb(i) = k3xbb
187	k3ac(i) = k3xac
188	k3bc(i) = k3xbc
189	anu = nu(i,4)-nu(i,3)
190	! anu = nu(i,4)-nu(i,3)+70
191	k3aap(i) = k3aa(i) * exp( -ee/tt * anu )/6.d0
192	k3abp(i) = k3ab(i) * exp( -ee/tt * anu )/6.d0
193	k3acp(i) = k3ac(i) * exp( -ee/tt * anu )/6.d0
194	anu = nu(i,4)-nu(i,2)
195	! anu = nu(i,4)-nu(i,2)+40
196	k3bap(i) = k3ba(i) * exp( -ee/tt * anu )/4.d0
197	k3bbp(i) = k3bb(i) * exp( -ee/tt * anu )/4.d0
198	k3bcp(i) = k3bc(i) * exp( -ee/tt * anu )/4.d0
199
200	k7a(i) = k7xa
201	k7b(i) = k7xb
202	k7ap(i) = k7a(i) * exp(dble( -ee/tt*(nu(i,4)-nu(i,3)) ))/6.d0
203	k7bp(i) = k7b(i) * exp(dble( -ee/tt*(nu(i,4)-nu(i,2)) ))/4.d0
204	end do
205
206
207	! the next ones correspond to vv2 in table 4, paper i
208	co2i(001) + co2 ---> co2i(020) + co2(010) k6
209	! k6(i) i = 1,2,3,4
210	! we need a new index for the inverse rates due to both fractions :
211	c k6a(i) i=2,3,4 co2i(001) + co2 ---> co2i(020) + co2(010)
212	c k6b(i) " co2i(001) + co2 ---> co2i(010) + co2(020)
213
214	if (iopt6.eq.1) then
215
216	if(tt.le.175.d0)then
217	k6x=8.6d-15
218	elseif(tt.gt.175.0.and.tt.le.200.d0)then
219	k6x=8.6d-15+9.d-16*(175.d0-tt)/25.d0
220	elseif(tt.gt.200.0.and.tt.le.225.d0)then
221	k6x=7.7d-15+5.d-16*(200.d0-tt)/25.d0
222	elseif(tt.gt.225.0.and.tt.le.250.d0)then
223	k6x=7.20d-15+6.d-16*(tt-225.d0)/25.d0
224	elseif(tt.gt.250.0.and.tt.le.275.d0)then
225	k6x=7.80d-15+1.d-15*(tt-250.d0)/25.d0
226	elseif(tt.gt.275.0.and.tt.le.300.d0)then
227	k6x=8.80d-15+1.3d-15*(tt-275.d0)/25.d0
228	elseif(tt.gt.300.0.and.tt.le.325.d0)then
229	k6x=10.1d-15+1.54d-15*(tt-300.d0)/25.d0
230	elseif(tt.gt.325.0)then
231	k6x=11.6d-15
232	end if
233
234	elseif (iopt6.eq.2) then
235
236	k6x = 3.6d-13 * exp( -1660/tt + 176948/tt**2.d0 )
237	if (tt.lt.175) k6x = 8.8d-15
238
239	end if
240
241	k6x1 = k6x * rf6 * frac6
242	k6x2 = k6x * rf6 * (1.-frac6)
243
244	k6 = k6x * rf6
245	k6p = k6 / 8.d0 * exp(dble( -ee/tt * (nu(1,4)-nu(1,2)-nu(1,1)) ))
246	do i=2,nisot
247	k6a(i) = k6x1
248	k6b(i) = k6x2
249	anu = nu(i,4)-nu(i,2)-nu(1,1)
250	k6ap(i) = k6a(i) / 8.d0 * exp(dble( -ee*anu/tt ))
251	anu = nu(i,4)-nu(i,1)-nu(1,2)
252	k6bp(i) = k6b(i) / 8.d0 * exp(dble( -ee*anu/tt ))
253	end do
254
255
256	co2i(0v0) + co2i ---> co2i(0v-10) + co2i(010)
257	c k5 esta reaccion es desdenable frente a co2 como colisionante.
258	! k5=3.0d-15+6.0d-17*(tt-210.d0)
259	! k5=k5*rf5
260	! k5p=k5/2.d0exp(-ee125.77/tt)
261
262
263	co2i(0v0) + m ---> co2i(0v-10) + m k19 (vt2,k5,k6 in table 3, paper
264	co2i(0v0) + o3p ---> co2i(0v-10) + o3p k20 (vt2,k7 in table 3, paper
265	c k19vm(i) m = a --- co2 v = a --- 3 i=1,2,3,4
266	! b --- n2 b --- 2
267	! c --- co c --- 1
268	c k20v(i) v = a --- 3 i = 1,2,3,4
269	! b --- 2
270	! c --- 1
271	!
272	! k20x=1.9d-8exp(-76.75/(tt*(1.d0/3.d0))) ! taylor,74 reajusted
273	! k20x=2.32d-9exp(-76.75/(tt(1./3.)))+1.0d-14sqrt(tt) ! k&j, 83
274	! k20x = 1.43d-12*(tt/300.d0) ! shved et al, 90
275
276	if (iopt20.eq.1) then ! first version of pap1
277	k20x=2.32d-9exp(-76.75/(tt(1./3.)))+3.5d-13sqrt(tt) ! s&w, 91
278	k20xb = k20x / 2.d0 * rf20
279	k20xc = k20xb
280	k20xa = 3.d0/2.d0 * k20xb
281	elseif(iopt20.eq.2) then ! revision for the papers in feb 93
282	k20x=3.d-12 ! minimum value of lopez-puertas et al., 92
283	k20xc = k20x * rf20
284	k20xb = 2.d0 * k20xc
285	k20xa = 3.d0/2.d0 * k20xb
286	elseif(iopt20.eq.3) then ! values from boug & roble '91
287	k20x=1.d-12/sqrt(300.) * sqrt(tt)
288	k20xc = k20x * rf20
289	k20xb = 2.d0 * k20xc
290	k20xa = 3.d0/2.d0 * k20xb
291	elseif(iopt20.eq.4) then ! values from boug & dick '88 case b
292	k20x=7.d-13
293	k20xc = k20x * rf20
294	k20xb = 2.d0 * k20xc
295	k20xa = 3.d0/2.d0 * k20xb
296	elseif(iopt20.eq.5) then ! values from s.bougher (oct-98)
297	k20x = 1.732d-13 * sqrt(tt) ! 1/sqrt(300) * sqrt(t)
298	k20xc = k20x * rf20
299	k20xb = 2.d0 * k20xc
300	k20xa = 3.d0/2.d0 * k20xb
301	end if
302
303	if (iopt19.eq.0) then
304
305	k19xca = 4.64d-10 * exp(dble( - 74.75 / tt**(1.d0/3.d0) ))
306	k19xcb = 6.69d-10 * exp(dble( - 84.07 / tt**(1.d0/3.d0) ))
307	k19xcc = k19xcb
308
309	if ( tt.le.250 ) then
310	k19xba = 181.25d0
311	elseif ( tt.ge.310 ) then
312	k19xba = 200.d0 + 0.9d0 * ( tt - 310.d0 )
313	else
314	k19xba = 181.25d0 + 0.3125d0 * ( tt - 250.d0 )
315	end if
316	k19xba = k19xba * 1.03558d-19 * tt ! cm-1 s-1
317	k19xbb = 1.24d-14 * ( tt / 273.3d0 )**2.d0 ! taine & lepoutre 1979
318	k19xbc = k19xbb
319
320	k19xaa = 3.d0/2.d0 * k19xba
321	k19xab = 3.d0/2.d0 * k19xbb
322	k19xac = 3.d0/2.d0 * k19xbc
323
324	k19xaa = k19xaa * rf19
325	k19xab = k19xab * rf19
326	k19xac = k19xac * rf19
327	k19xba = k19xba * rf19
328	k19xbb = k19xbb * rf19
329	k19xbc = k19xbc * rf19
330	k19xca = k19xca * rf19
331	k19xcb = k19xcb * rf19
332	k19xcc = k19xcc * rf19
333
334	elseif (iopt19.ge.1) then
335
336	if (iopt19.eq.1) then ! lunt et. al., 1985 (thesis values)
337
338	if (tt.le.175.) then
339	k19xca = 4.d0 - 0.02d0 * (tt-140.d0)
340	k19xcb = 0.494d0 + 0.0076 * (tt-140.d0)
341	elseif (tt.ge.175. .and. tt.le.200.) then
342	k19xca = 3.3d0 - 0.02d0 * (tt-175.)
343	k19xcb = 0.76d0 + 0.0076d0 * (tt-175.d0)
344	elseif (tt.ge.200. .and. tt.le.225.) then
345	k19xca = 2.8d0 + 0.004d0 * (tt-200.d0)
346	k19xcb = 0.95d0 + 0.014d0 * (tt-200.d0)
347	elseif (tt.ge.225. .and. tt.le.250.) then
348	k19xca = 2.9d0 + 0.024d0 * (tt-225.d0)
349	k19xcb = 1.3d0 + 0.016d0 * (tt-225.d0)
350	elseif (tt.ge.250. .and. tt.le.275.) then
351	k19xca = 3.5d0 + 0.04d0 * (tt-250.d0)
352	k19xcb = 1.7d0 + 0.032d0 * (tt-250.d0)
353	elseif (tt.ge.275. .and. tt.le.295.) then
354	k19xca = 4.5d0 + 0.055d0 * (tt-275.d0)
355	k19xcb = 2.5d0 + 0.045d0 * (tt-275.d0)
356	elseif (tt.ge.295. .and. tt.le.320.) then
357	k19xca = 5.6d0 + 0.54d0 * (tt-295.d0)
358	k19xcb = 3.4d0 + 0.045d0 * (tt-295.d0)
359	end if
360	k19xca = k19xca * 1.d-15 * rf19
361	k19xcb = k19xcb * 1.d-15 * rf19
362	k19xcc = k19xcb
363
364	elseif (iopt19.eq.2) then ! revision for the papers, feb 1993
365
366	! k19xca = 7.3d-14 * exp( -850.3d0/tt + 86523.d0/tt**2.d0 )
367	k19xca = 4.2d-12 * exp( -2988.d0/tt + 303930.d0/tt**2.d0 )
368	if (tt.le.175.) k19xca = 3.3d-15
369	k19xcb = 2.1d-12 * exp( -2659.d0/tt + 223052.d0/tt**2.d0 )
370	if (tt.le.175.) k19xcb = 7.6d-16
371	k19xca = k19xca * rf19
372	k19xcb = k19xcb * rf19
373	k19xcc = k19xcb
374
375	elseif (iopt19.eq.3) then ! values from dick'72 for k19xc
376	! k19xcb is not modified
377	if (tt.le.158.) then
378	k19xca = 0.724d-15
379	elseif (tt.le.190.) then
380	k19xca = 0.724d-15 +
381	@ (1.1d-15-0.724d-15) * (tt-158.) / (190.-158.)
382	elseif (tt.le.250.) then
383	k19xca = 1.1d-15 +
384	@ (3.45d-15-1.1d-15) * (tt-190.) / (250.-190.)
385	elseif (tt.gt.250.) then
386	k19xca = 3.45d-15
387	end if
388	k19xcb = 2.1d-12 * exp( -2659.d0/tt + 223052.d0/tt**2.d0 )
389	if (tt.le.175.) k19xcb = 7.6d-16
390	k19xca = k19xca * rf19
391	k19xcb = k19xcb * rf19
392	k19xcc = k19xcb
393
394	elseif (iopt19.eq.5) then
395
396	k19xca = 5.2d-15 ! s.bougher, nov-98
397	k19xcb = 7.6d-16 ! nuestro, de feb-93
398	k19xcc = k19xcb
399
400	k19xca = k19xca * rf19
401	k19xcb = k19xcb * rf19
402
403	end if
404
405	factor = 2.5d0
406	k19xba = factor * k19xca
407	k19xbb = factor * k19xcb
408	k19xbc = factor * k19xcc
409	factor = 3.d0/2.d0
410	k19xaa = factor * k19xba
411	k19xab = factor * k19xbb
412	k19xac = factor * k19xbc
413
414	end if
415
416	do i = 1, nisot
417
418	k19aa(i) = k19xaa
419	k19ba(i) = k19xba
420	k19ca(i) = k19xca
421	k19ab(i) = k19xab
422	k19bb(i) = k19xbb
423	k19cb(i) = k19xcb
424	k19ac(i) = k19xac
425	k19bc(i) = k19xbc
426	k19cc(i) = k19xcc
427	anu = nu(i,3)-nu(i,2)
428	k19aap(i) = k19aa(i) * 6.d0/4.d0 * exp(dble( -ee*anu/tt))
429	k19abp(i) = k19ab(i) * 6.d0/4.d0 * exp(dble( -ee*anu/tt))
430	k19acp(i) = k19ac(i) * 6.d0/4.d0 * exp(dble( -ee*anu/tt))
431	anu = nu(i,2)-nu(i,1)
432	k19bap(i) = k19ba(i) * 2.d0 * exp(dble( -ee*anu/tt))
433	k19bbp(i) = k19bb(i) * 2.d0 * exp(dble( -ee*anu/tt))
434	k19bcp(i) = k19bc(i) * 2.d0 * exp(dble( -ee*anu/tt))
435	anu = nu(i,1)
436	k19cap(i) = k19ca(i) * 2.d0 * exp(dble( -ee*anu/tt))
437	k19cbp(i) = k19cb(i) * 2.d0 * exp(dble( -ee*anu/tt))
438	k19ccp(i) = k19cc(i) * 2.d0 * exp(dble( -ee*anu/tt))
439
440	k20a(i) = k20xa
441	k20b(i) = k20xb
442	k20c(i) = k20xc
443	k20ap(i) = k20a(i)6.d0/4.d0
444	@ exp(dble( -ee/tt * (nu(i,3)-nu(i,2)) ))
445	k20bp(i) = k20b(i)4.d0/2.d0
446	@ exp(dble( -ee/tt * (nu(i,2)-nu(i,1)) ))
447	k20cp(i) = k20c(i)2.d0/1.d0
448	@ exp(dble( -ee/tt * nu(i,1) ))
449	end do
450
451	! write(1,*) tt,k19cap(1),k19ac(1)
452
453	! the next ones correspond to vv3 in table 4 (paper i)
454	co2i(0v0) + co2 ---> co2i(0v-10) + co2(010) k21 also see k33
455	c k21v(i) v = a --- 3 i = 1,2,3,4
456	! b --- 2
457	! c --- 1
458	! we need a new index for the 030i rates due to both fractions :
459	c k21a1 co2i(030) + co2 ---> co2i(020) + co2(010)
460	c k21a2 co2i(030) + co2 ---> co2i(010) + co2(020)
461	co2i(010) + co2j(000) ---> co2i(000) + co2j(010) kijk21c see pag.22-s
462	! k23k21c i=628,j=636
463	! k24k21c i=628,j=627
464	! k34k21c i=636,j=627
465
466	if (iopt21.eq.0) then
467	k21x = 1.2d-11
468	k21xb = k21x
469	k21xa = 3.d0/2.d0 * k21x
470	k21xc = k21x ! esta ultima no se usa con 626
471	elseif (iopt21.eq.1) then
472	k21x = 2.49d-11 ! orr & smith, 1987
473	k21xb = k21x
474	k21xa = 3.d0/2.d0 * k21xb ! oscilador armonico
475	k21xc = k21xb / 2.d0 ! novedad mia
476	elseif (iopt21.eq.2) then
477	k21x = 100.d0*k19xca ! dickinson'76 (icarus)
478	k21xb = k21x
479	k21xa = 3.d0/2.d0 * k21xb ! oscilador armonico
480	k21xc = k21xb / 2.d0 ! novedad mia
481	end if
482	k21xa_626 = k21xa * rf21a !* 0.01d0 !* 10.d0
483	k21xa = k21xa * rf21a !* 0.01d0
484	k21xb = k21xb * rf21b
485	k21xc = k21xc * rf21c
486
487	k21a = k21xa_626
488	k21ap = k21a * 6.d0/8.d0 *
489	@ exp( dble( -ee/tt * (nu(1,3)-nu(1,2)-nu(1,1)) ) )
490	do i = 2, nisot
491	k21a1(i) = k21xa * frac21
492	k21a2(i) = k21xa * (1.d0-frac21)
493	k21a1p(i) = k21a1(i) * 6.d0/8.d0 *
494	@ exp(dble( -ee/tt* (nu(i,3)-nu(i,2)-nu(1,1)) ))
495	k21a2p(i) = k21a2(i) * 6.d0/8.d0 *
496	@ exp(dble( -ee/tt* (nu(i,3)-nu(i,1)-nu(1,2)) ))
497	end do
498
499
500	do i = 1, nisot
501	k21b(i) = k21xb
502	k21c(i) = k21xc
503	k21bp(i) = k21b(i) *
504	@ exp(dble( -ee/tt* (nu(i,2)-nu(i,1)-nu(1,1)) ))
505	k21cp(i) = k21c(i) *
506	@ exp(dble( -ee/tt * (nu(i,1)-nu(1,1)) ))
507	end do
508
509	k23k21c = k21xc
510	k24k21c = k21xc
511	k34k21c = k21xc
512	k23k21cp = k23k21c2.d0/2.d0
513	@ exp(dble( -ee/tt* (nu(2,1)-nu(3,1)) ))
514	k24k21cp = k24k21c2.d0/2.d0
515	@ exp(dble( -ee/tt* (nu(2,1)-nu(4,1)) ))
516	k34k21cp = k34k21c2.d0/2.d0
517	@ exp(dble( -ee/tt* (nu(3,1)-nu(4,1)) ))
518
519	!these are also vv3 in table 4, paper i
520	c k31 & k32
521
522	k31 = k21x * rf31 ! we're suposing thar the rate for the deactivation
523	! v-v from high combinational levels is the same
524	k32 = k21x * rf32 ! that the one for : (020) --> (010) + (010)
525
526	co2(*) + co2i ---> co2() + co2i(**) k33
527	c k33a1 : co2(001) + co2i ---> co2(020) + co2i(010) (vv2, table 4,
528	! a2 : co2(001) + co2i ---> co2(010) + co2i(020) "
529	! b1 : co2(030) + co2i ---> co2(020) + co2i(010) (vv3, table 4,
530	! b2 : co2(030) + co2i ---> co2(010) + co2i(020) "
531	! c : co2(020) + co2i ---> co2(010) + co2i(010) "
532	! we have to add an index to the inverse rates, depending on the isotope
533
534	k33c = k21x * rf33bc
535	k33b1 = 3.d0/3.d0 * k33c * frac33
536	k33b2 = 3.d0/3.d0 * k33c * (1.d0-frac33)
537	k33a1 = k6x * rf33a * frac33
538	k33a2 = k6x * rf33a * (1.d0-frac33)
539
540	do i=2,nisot
541	k33a1p(i)=k33a1*
542	@ 1.d0/8.d0exp(dble( -ee/tt (nu(1,4)-nu(1,2)-nu(i,1)) ))
543	k33a2p(i)=k33a2*
544	@ 1.d0/8.d0exp(dble( -ee/tt (nu(1,4)-nu(1,1)-nu(i,2)) ))
545	k33b1p(i)=k33b1*
546	@ 6.d0/8.d0exp(dble( -ee/tt (nu(1,3)-nu(1,2)-nu(i,1)) ))
547	k33b2p(i)=k33b2*
548	@ 6.d0/8.d0exp(dble( -ee/tt (nu(1,4)-nu(1,1)-nu(i,2)) ))
549	k33cp(i) =k33c * exp(dble( -ee/tt * (nu(1,2)-nu(1,1)-nu(i,1)) ))
550	end do
551
552	! here they are the vt3 in table 3, paper i
553	co2(2.7um) + m ---> co2(2.7um) + m k27
554	c k27a : n8 + m ---> n6 + m
555	! b : n7 + m ---> n6 + m
556	! c : n8 + m ---> n7 + m
557
558	if (iopt27.eq.0) then
559	k27a = 3.d-11 !between fermi levels
560	k27b = 3.d-13 !between side levels
561	k27c = 2.d0 * k27b !between side levels
562	elseif (iopt27.ge.1) then ! orr & smith, 1987
563	k27a = 1.55d-12
564	k27c = 4.97d-12
565	k27b = k27c
566	end if
567	k27a = k27a * rf27f
568	k27b = k27b * rf27s
569	k27c = k27c * rf27s
570
571	k27ap = k27a * exp(dble( -ee/tt * (nu(1,8)-nu(1,6)) ))
572	k27bp = k27b * exp(dble( -ee/tt * (nu(1,7)-nu(1,6)) ))
573	k27cp = k27c * exp(dble( -ee/tt * (nu(1,8)-nu(1,7)) ))
574
575
576	! the next two are not used in the model:
577
578	c k28 : n* + n2 ---> n*low + n2(1)
579	! k28v v = a --- n8
580	! b --- n7
581	! c --- n6
582	! k28a = 5.d-13 * sqrt(300.d0/tt) * rf28 ! = k1
583	! k28b = k28a
584	! k28c = k28a
585	! k28ap = k28a * exp( -ee/tt * (nu(1,8)-1388.1847-nun2) )
586	! k28bp = k28b * exp( -ee/tt * (nu(1,7)-1335.1317-nun2) )
587	! k28cp = k28c * exp( -ee/tt * (nu(1,6)-1285.4087-nun2) )
588
589	c k29 : n* + co ---> n*low + co(1)
590	! k29v v = a --- n8
591	! b --- n7
592	! c --- n6
593
594	c k29a = ?????????? * rf29
595	c k29b = k29a
596	c k29c = k29a
597
598	c k29ap = k29a * exp( -ee/tt * (nu(1,8)-1388.1847-nuco) )
599	c k29bp = k29b * exp( -ee/tt * (nu(1,7)-1335.1317-nuco) )
600	c k29cp = k29c * exp( -ee/tt * (nu(1,6)-1285.4087-nuco) )
601
602
603	! these are also vv1 processes in table 4, paper i
604	c k26 :
605	! 1. deactivation of the 626 isotope:
606	! reaction: n* + co2i ---> n*low + co2i(001) ; i=1-4
607	! nomenclature: k26v ; v=a,b,c,d for n8,n7,n6,n5 respectively
608	! inverse rates: k26vp(i) ; i=1-4
609	! 2. deactivation of the minor isotopes:
610	! reaction: n_i + co2j ---> n_i_low + (001)_j ; i=2-4 ; j=1-4
611	! nomenclature: k26vij ; v=a,c,d for n8,n6,n5 respectively
612	! inverse rates: k26vijp
613	! 3. notes:
614	! a * it is clear that : k26vij=/=k26vjip,
615	! b * at the moment we do not include inverse rates for the case 2., o
616	! the deactivation of the main isotope (pg. 32b, sn1).
617	! c * not 0221 level for minor isotopes is considered.
618	! d * only a value is known for these rates, so all of the deactivatio
619	! the same, but not the inverse rates.
620	! e * although all the direct deactivation constants have the same val
621	! it is useful to distinguish between them with the present nam
622
623	k26a = 6.8d-12 * sqrt(tt) * rf26 ! = k2
624	k26b = k26a
625	k26c = k26a
626	if (iopt26.eq.0 .or. iopt26.eq.2) then
627	k26d = k26a
628	elseif( iopt26.eq.1) then
629	k26d = 1.15d-10 * rf26
630	end if
631
632	do i=1,4
633	k26ap(i) = k26a *
634	@ exp(dble( -ee/tt * (nu(1,8)-nu12_1000-nu(i,4)) ))
635	k26bp(i) = k26b *
636	@ exp(dble( -ee/tt * (nu(1,7)- nu(1,2) -nu(i,4)) ))
637	k26cp(i) = k26c *
638	@ exp(dble( -ee/tt * (nu(1,6)-nu12_0200-nu(i,4)) ))
639	k26dp(i) = k26d *
640	@ exp(dble( -ee/tt * (nu(1,5)- nu(1,1) -nu(i,4)) ))
641	end do
642
643	k26a21 = k26a
644	k26c21 = k26c
645	k26d21 = k26d
646	k26a22 = k26a
647	k26c22 = k26c
648	k26d22 = k26d
649	k26a23 = k26a
650	k26c23 = k26c
651	k26d23 = k26d
652	k26a24 = k26a
653	k26c24 = k26c
654	k26d24 = k26d
655
656	k26a31 = k26a
657	k26c31 = k26c
658	k26d31 = k26d
659	k26a32 = k26a
660	k26c32 = k26c
661	k26d32 = k26d
662	k26a33 = k26a
663	k26c33 = k26c
664	k26d33 = k26d
665	k26a34 = k26a
666	k26c34 = k26c
667	k26d34 = k26d
668
669	k26a41 = k26a
670	k26c41 = k26c
671	k26d41 = k26d
672	k26a42 = k26a
673	k26c42 = k26c
674	k26d42 = k26d
675	k26a43 = k26a
676	k26c43 = k26c
677	k26d43 = k26d
678	k26a44 = k26a
679	k26c44 = k26c
680	k26d44 = k26d
681
682	!! some examples of inverse rates, although not used at the moment
683	! k26a32p = k26a32 * exp( -ee* (nu(3,8)-nu32_1000-nu(2,4)) / tt )*1./1.
684	! k26c32p = k26c32 * exp( -ee* (nu(3,6)-nu32_0200-nu(2,4)) / tt )*1./1.
685	! k26d32p = k26d32 * exp( -ee* (nu(3,5)- nu(3,1) -nu(2,4)) / tt )*2./2.
686	!
687	! k26a43 = k26a34 * exp( -ee* (nu(3,8)-nu32_1000-nu(4,4)) / tt )*1./1.
688	! k26c43 = k26c34 * exp( -ee* (nu(3,6)-nu32_0200-nu(4,4)) / tt )*1./1.
689	! k26d43 = k26d34 * exp( -ee* (nu(3,5)- nu(3,1) -nu(4,4)) / tt )*2./2.
690	!
691	! k26a24p = k26a24 * exp( -ee* (nu(2,8)-nu22_1000-nu(4,4)) / tt )*1./1.
692	! k26c24p = k26c24 * exp( -ee* (nu(2,6)-nu22_0200-nu(4,4)) / tt )*1./1.
693	! k26d24p = k26d24 * exp( -ee* (nu(2,5)- nu(2,1) -nu(4,4)) / tt )*2./2.
694
695
696	! this is taken as vv4 in table 4, paper i (in the inverse direction)
697	c k41 : co(v) + co2 ---> co(v-1) + co2(001) + de
698	! k41_v v=1,2,3,4
699	!
700	! de = -205.9 cm-1 when v=1
701	! de = -232.9 cm-1 when v=2
702	! de = -258.6 cm-1 when v=3
703	! de = -285.0 cm-1 when v=4
704
705	k41p_taylor = 1.56d-11 * exp( -30.12/tt**0.333333 )! [ s-1 cm+3 ]
706	k41p_shved = 7.5d7/sqrt(tt) ! [ s-1 atm-1 ]
707	k41p_shved = k41p_shved * 1.38d-16/1013250. * tt! [ s-1 cm+3 ]
708	k41p_starr_hannock = 6.27d3 ! [ s-1 torr-1 ]
709
710	if (iopt41.eq.1) then
711	k41p_1 = k41p_starr_hannock *
712	@ 760.1.38d-16/1013250. tt ! [ s-1 cm+3 ]
713	elseif (iopt41.eq.2) then
714	k41p_1 = 1.6d-12 * exp( -1169/tt + 77601/tt**2.d0 )
715	end if
716	k41p_1 = k41p_1 * rf41
717	k41_1 = k41p_1 * exp(dble( -ee * 205.9/tt ))
718	k41_2 = k41_1
719	k41p_2 = k41_2 * exp(dble( -ee * (-232.9)/tt ))
720	k41_3 = k41_1
721	k41p_3 = k41_3 * exp(dble( -ee * (-258.6)/tt ))
722	k41_4 = k41_1
723	k41p_4 = k41_4 * exp(dble( -ee * (-285.0)/tt ))
724
725	!k41p_1 = k41p_1 * 1.d-6
726	!k41p_2 = k41p_2 * 1.d-6
727
728	c k41iso : 63(v) + co2 ---> 63(v-1) + co2(001) + de
729	! k41iso_v v=1,2,3
730	! de = -253 cm-1 when v=1
731	! de = -278 cm-1 when v=2
732	! de = -303 cm-1 when v=3
733
734	k41iso_1 = k41_1
735	k41iso_1p = k41iso_1 * exp(dble( -ee * (-253.)/tt ))
736	k41iso_2 = k41iso_1
737	k41iso_2p = k41iso_2 * exp(dble( -ee * (-278.)/tt ))
738	k41iso_3 = k41iso_1
739	k41iso_3p = k41iso_3 * exp(dble( -ee * (-303.)/tt ))
740
741
742
743	c k42 : co(v) + co ---> co(v-1) + co(1) + de=-26.481 si v=2 K42
744	! -52.8940 si v=3 k42b
745	! -79.2402 si v=4 k42c
746	! tomado de stepanova & shved (ellos de powell, 1975), ver pg .. l5
747	! solo para v=2 :
748
749	k42 = 2.89d-10 * (1./sqrt(tt) + 67.4/tt*1.5)
750	@ exp(dble(24.78/tt))
751	k42 = k42 * rf42
752	k42b = k42
753	k42c = k42
754	k42p = k42 * exp(dble( -ee * (-26.481)/tt ))
755	k42bp = k42b * exp(dble( -ee * (-52.894)/tt ))
756	k42cp = k42c * exp(dble( -ee * (-79.24)/tt ))
757
758	c k42iso : 63(v) + 63 ---> 63(v-1) + 63(1) + de=-25.31 si v=2 K42iso
759	! -50.57 si v=3 k42isob
760	! tomado de stepanova & shved (ellos de powell, 1975), ver pg .. l5
761	! solo para v=2 :
762
763	k42iso = k42
764	k42isop = k42iso * exp(dble( -ee * (-25.31)/tt ))
765	k42isob = k42
766	k42isobp = k42isob * exp(dble( -ee * (-50.57)/tt ))
767
768
769	c k43 : co(v) + o3p ---> co(v-1) + o3p + de=2143
770	! tomado de lewittes et. al, 1978 para v=1
771
772	if (iopt43.eq.1) then
773	tt1 = tt - 300.
774	k43 = 2.85d-14 * exp( dble( 9.5e-3tt1 + 1.11e-4tt1**2. ) )
775	elseif (iopt43.eq.2) then
776	k43 = 1.4d-5 * exp( -10957.d0 / tt + 1.486d6 / tt**2.d0 )
777	if ( tt.lt.265.0 ) k43 = 2.3d-14
778	end if
779	k43 = k43 * rf43
780	k43p = k43 * exp( -ee * dble(2143.3 / tt) )
781
782	c k43iso : co63(v) + o3p ---> co63(v-1) + o3p + de=2096
783	! Por similitud con el anterior
784
785	k43iso = k43
786	k43isop = k43iso * exp( -ee * dble(2096. / tt) )
787
788
789	c k44 : co62(v) + co63 ---> co62(v-1) + co63(1) + de
790	! basado en Lopez-Valverde et al para el caso v=1, solo usamos este
791	! k44x x = a --- v=1 de= 147.33
792	! b --- v=2 de= 20.7241
793	! c --- v=3 de= -5.7
794	! d --- v=4 de= -32.0361
795
796	k44a = 2.0d-12 * rf44 ! Solo vamos a usar este, no los b,c,d
797	k44b = k44a
798	k44c = k44a
799	k44d = k44a
800
801	de = 147.33
802	k44ap = k44a * exp(dble( -ee * de/tt ))
803	de = 20.7241
804	k44bp = k44b * exp(dble( -ee * de/tt ))
805	de = -5.7
806	k44cp = k44c * exp(dble( -ee * de/tt ))
807	de = -32.0361
808	k44dp = k44d * exp(dble( -ee * de/tt ))
809
810
811	co2(hcl) + co2 --> co2 + co2 + de(hcl)
812	! este rate tambien lo usamos para los high combination levels (para tra
813	! al lte. cualquier valor vale, supongo. es k_vthcl
814
815	k_vthcl = 3.3d-15 ! similar al valor pequenho del vt2
816	k_vthcl = k_vthcl * rf_hcl
817
818	return
819	end

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

Download in other formats:

Original Format