Context Navigation

nlte_leedat.F @ 652

Last change on this file since 652 was 498, checked in by emillour, 14 years ago
Mars GCM: Cleanup of the NLTE routines which have been packed together to limit the number of files. Also enforced that file names are non-capitalized (needed by the create_make_gcm scripts to better evaluate dependencies when building the makefile). FGG+EM
File size: 54.0 KB

Rev	Line
[498]	1	c***********************************************************************
	2	subroutine NLTE_leedat
	3
	4	c reads planetary and molecular parameters
	5
	6	c jan 98 malv first version
	7	c jul 2011 malv+fgg adapted to LMD-MGCM
	8	c***********************************************************************
	9
	10	implicit none
	11
	12	include 'datafile.h'
	13	include 'nlte_paramdef.h'
	14	include 'nlte_commons.h'
	15
	16
	17	c local variables
	18	integer i,j, k,lun1, lun2
	19	integer :: ib = 0
	20	integer isot
	21	character isotcode*2
	22	character ibcode1*1
	23
	24	c formats
	25	132 format (i2)
	26	101 format(i1)
	27	c***********************************************************************
	28
	29	lun1 = 1
	30	lun2 = 2
	31
	32	do k=1,nisot
	33	! encode (2,132,isotcode) indexisot(k)
	34	write (isotcode,132) indexisot(k)
	35	open (lun1,
	36	@ file=trim(datafile)//'/NLTEDAT/enelow'//isotcode//'.dat',
	37	@ status='old')
	38	open (lun2,
	39	@ file=trim(datafile)//'/NLTEDAT/deltanu'//isotcode//'.dat',
	40	@ status='old')
	41	read (lun1,*)
	42	read (lun2,*)
	43	read (lun1,*) (elow(k,i), i=1,nb)
	44	read (lun2,*) (deltanu(k,i), i=1,nb)
	45	close (lun1)
	46	close (lun2)
	47	end do
	48
	49
	50	c Call to rhist
	51
	52	hfile1='hid'
	53	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
	54	if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
	55
	56	ib=1
	57	do isot=1,4
	58	c Cases 1,2,3,4
	59	if (isot.eq.1) hisfile = hfile1//'26-1.dat'
	60	if (isot.eq.2) hisfile = hfile1//'28-1.dat'
	61	if (isot.eq.3) hisfile = hfile1//'36-1.dat'
	62	if (isot.eq.4) hisfile = hfile1//'27-1.dat'
	63	call rhist(1.0)
	64	if (isot.eq.1) then !Case 1
	65	mm_c1=mm
	66	nbox_c1=nbox
	67	tmin_c1=tmin
	68	tmax_c1=tmax
	69	do i=1,nbox_max
	70	no_c1(i)=no(i)
	71	dist_c1(i)=dist(i)
	72	do j=1,nhist
	73	sk1_c1(j,i)=sk1(j,i)
	74	xls1_c1(j,i)=xls1(j,i)
	75	xln1_c1(j,i)=xln1(j,i)
	76	xld1_c1(j,i)=xld1(j,i)
	77	enddo
	78	enddo
	79	do j=1,nhist
	80	thist_c1(j)=thist(j)
	81	enddo
	82	else if(isot.eq.2) then !Case 2
	83	mm_c2=mm
	84	nbox_c2=nbox
	85	tmin_c2=tmin
	86	tmax_c2=tmax
	87	do i=1,nbox_max
	88	no_c2(i)=no(i)
	89	dist_c2(i)=dist(i)
	90	do j=1,nhist
	91	sk1_c2(j,i)=sk1(j,i)
	92	xls1_c2(j,i)=xls1(j,i)
	93	xln1_c2(j,i)=xln1(j,i)
	94	xld1_c2(j,i)=xld1(j,i)
	95	enddo
	96	enddo
	97	do j=1,nhist
	98	thist_c2(j)=thist(j)
	99	enddo
	100	else if (isot.eq.3) then ! Case 3
	101	mm_c3=mm
	102	nbox_c3=nbox
	103	tmin_c3=tmin
	104	tmax_c3=tmax
	105	do i=1,nbox_max
	106	no_c3(i)=no(i)
	107	dist_c3(i)=dist(i)
	108	do j=1,nhist
	109	sk1_c3(j,i)=sk1(j,i)
	110	xls1_c3(j,i)=xls1(j,i)
	111	xln1_c3(j,i)=xln1(j,i)
	112	xld1_c3(j,i)=xld1(j,i)
	113	enddo
	114	enddo
	115	do j=1,nhist
	116	thist_c3(j)=thist(j)
	117	enddo
	118	else if (isot.eq.4) then ! Case 4
	119	mm_c4=mm
	120	nbox_c4=nbox
	121	tmin_c4=tmin
	122	tmax_c4=tmax
	123	do i=1,nbox_max
	124	no_c4(i)=no(i)
	125	dist_c4(i)=dist(i)
	126	do j=1,nhist
	127	sk1_c4(j,i)=sk1(j,i)
	128	xls1_c4(j,i)=xls1(j,i)
	129	xln1_c4(j,i)=xln1(j,i)
	130	xld1_c4(j,i)=xld1(j,i)
	131	enddo
	132	enddo
	133	do j=1,nhist
	134	thist_c4(j)=thist(j)
	135	enddo
	136	endif
	137	enddo
	138	do ib=2,4
	139	isot=1
	140	write (ibcode1,101) ib
	141	hisfile = hfile1//'26-'//ibcode1//'.dat'
	142	call rhist (1.0)
	143	if (ib.eq.2) then !Case 5
	144	mm_c5=mm
	145	nbox_c5=nbox
	146	tmin_c5=tmin
	147	tmax_c5=tmax
	148	do i=1,nbox_max
	149	no_c5(i)=no(i)
	150	dist_c5(i)=dist(i)
	151	do j=1,nhist
	152	sk1_c5(j,i)=sk1(j,i)
	153	xls1_c5(j,i)=xls1(j,i)
	154	xln1_c5(j,i)=xln1(j,i)
	155	xld1_c5(j,i)=xld1(j,i)
	156	enddo
	157	enddo
	158	do j=1,nhist
	159	thist_c5(j)=thist(j)
	160	enddo
	161	else if (ib.eq.3) then !Case 6
	162	mm_c6=mm
	163	nbox_c6=nbox
	164	tmin_c6=tmin
	165	tmax_c6=tmax
	166	do i=1,nbox_max
	167	no_c6(i)=no(i)
	168	dist_c6(i)=dist(i)
	169	do j=1,nhist
	170	sk1_c6(j,i)=sk1(j,i)
	171	xls1_c6(j,i)=xls1(j,i)
	172	xln1_c6(j,i)=xln1(j,i)
	173	xld1_c6(j,i)=xld1(j,i)
	174	enddo
	175	enddo
	176	do j=1,nhist
	177	thist_c6(j)=thist(j)
	178	enddo
	179	else if (ib.eq.4) then !Case 7
	180	mm_c7=mm
	181	nbox_c7=nbox
	182	tmin_c7=tmin
	183	tmax_c7=tmax
	184	do i=1,nbox_max
	185	no_c7(i)=no(i)
	186	dist_c7(i)=dist(i)
	187	do j=1,nhist
	188	sk1_c7(j,i)=sk1(j,i)
	189	xls1_c7(j,i)=xls1(j,i)
	190	xln1_c7(j,i)=xln1(j,i)
	191	xld1_c7(j,i)=xld1(j,i)
	192	enddo
	193	enddo
	194	do j=1,nhist
	195	thist_c7(j)=thist(j)
	196	enddo
	197	endif
	198	enddo
	199
	200	c end
	201	return
	202
	203	end
	204
	205
	206
	207	c *******************************************************************
	208
	209	subroutine rhist (factor_comp)
	210
	211	c reads histogram data arrays created by ~/spectral/his.for
	212	c malv nov-98 add average distance between lines for overlapp
	213
	214	c *******************************************************************
	215
	216
	217	implicit none
	218
	219	include 'nlte_paramdef.h'
	220	include 'nlte_commons.h'
	221	include 'datafile.h'
	222
	223	c arguments
	224	real factor_comp
	225
	226	c local variables
	227	integer j, r
	228	real*8 sk1_aux, xls1_aux, xln1_aux, xld1_aux,weight,nu0
	229	character tonto*50
	230
	231	c formats
	232	! 100 format(80a1) ! Esto es si fuese byte dummy(80)
	233	100 format(a80) ! Esto es si fuese character dummy*80
	234	150 format(a50) ! Esto es si fuese character dummy*80
	235
	236	c ***************
	237
	238	open(unit=3,
	239	$ file=trim(datafile)//'/NLTEDAT/'
	240	$ //hisfile(1:len_trim(hisfile)),status='old')
	241	!read(3,100) dummy
	242	read(3,150) tonto
	243	read(3,*) weight
	244	read(3,*) mm
	245	read(3,*) nu0
	246	read(3,*) nbox
	247	!read(3,'(a)') dumm
	248	read(3,'(a)') tonto
	249
	250	if ( nbox .gt. nbox_max ) then
	251	write (,) ' nbox too large in input file ', hisfile
	252	stop ' Check maximum number nbox_max in mz1d.par '
	253	endif
	254	do 1 j=1,mm ! for each temperature
	255	read(3,*) thist(j)
	256	do r=1,nbox ! for each box
	257	read(3,*) no(r), sk1(j,r), xls1(j,r),xln1(j,r),xld1(j,r),
	258	@ dist(r)
	259	c xld1(j,r)=xld1(j,r)*0.83255 !0.83255=sqrt(log(2))
	260	enddo
	261	1 continue
	262	tmax=thist(1)
	263	tmin=thist(mm)
	264
	265	!close(unit=3,dispose='save')
	266	close(unit=3)
	267
	268
	269	do 2 j=1,mm
	270	do r=1,nbox
	271	sk1(j,r) = sk1(j,r) * factor_comp
	272	enddo
	273	2 continue
	274
	275
	276	return
	277	end
	278
	279
	280
	281
	282	c***********************************************************************
	283	subroutine leetvt
	284
	285	c reads input vibr. temps. from external files or sets lte values
	286	c according to the driver table
	287
	288	c jul 2011 malv+fgg adapted to LMD-MGCM
	289	c malv Jan 07 Add new vertical fine-grid for NLTE
	290	c jan 98 malv based on solar10sub
	291	c***********************************************************************
	292
	293	implicit none
	294
	295	include 'nlte_paramdef.h'
	296	include 'nlte_commons.h'
	297
	298	c local variables
	299	integer i
	300	real*8 zld(nl), zyd(nzy)
	301	real*8 xvt11(nzy), xvt21(nzy), xvt31(nzy), xvt41(nzy)
	302
	303	c***********************************************************************
	304
	305	do i=1,nzy
	306	zyd(i) = dble(zy(i))
	307	xvt11(i)= dble( ty(i) )
	308	xvt21(i)= dble( ty(i) )
	309	xvt31(i)= dble( ty(i) )
	310	xvt41(i)= dble( ty(i) )
	311	end do
	312
	313
	314	c interpolate to the nlte subgrid
	315
	316	do i=1,nl
	317	zld(i) = dble( zl(i) )
	318	enddo
	319	call interdp ( v626t1,zld,nl, xvt11,zyd,nzy, 1)
	320	call interdp ( v628t1,zld,nl, xvt21,zyd,nzy, 1)
	321	call interdp ( v636t1,zld,nl, xvt31,zyd,nzy, 1)
	322	call interdp ( v627t1,zld,nl, xvt41,zyd,nzy, 1)
	323
	324
	325	c end
	326	return
	327	end
	328
	329
	330
	331	c***********************************************************************
	332
	333	subroutine getk (tt)
	334
	335	c jan 98 malv version for mz1d. copied from solar10/getk.f
	336	c jul 2011 malv+fgg adapted to LMD-MGCM
	337	c***********************************************************************
	338
	339	implicit none
	340
	341	include 'nlte_paramdef.h'
	342	include 'nlte_commons.h'
	343
	344	c arguments
	345	real tt ! i. temperature
	346
	347	!! local variables:
	348	real*8 k1x, k7x,k7xa,k7xb
	349	real*8 k3x, k3xaa,k3xac,k3xab,k3xbb,k3xba,k3xbc
	350	real*8 k3xco2,k3xn2,k3xco, k6x,k6x1,k6x2
	351	real*8 k20x,k20xa,k20xb,k20xc
	352	real*8 k19xca,k19xcb,k19xcc
	353	real*8 k19xba,k19xbb,k19xbc
	354	real*8 k19xaa,k19xab,k19xac
	355	real*8 k21x,k21xa,k21xb,k21xc
	356	real*8 anu, factor , k21xa_626
	357	real tt1,tt2, de
	358	integer i
	359
	360	c***********************************************************************
	361
	362	co2i(001) + n2 ---> co2i + n2(1) k1 (not considered in the model
	363	c k1(i) i = 1 --- 626
	364	! 2 --- 636
	365	! 3 --- 628
	366	! 4 --- 627
	367
	368	! k1x = 5.d-13 * sqrt(300.d0/tt)
	369	! do i=1,nisot
	370	! k1(i) = k1x * rf1
	371	! k1p(i) = k1(i) * exp( -ee/tt 1.d0 (-nun2+nu(i,4)) )
	372	! end do
	373
	374	co2(001) + co2i ---> co2 + co2i(001) k2 (vv1 in table 4, paper i)
	375	c k2i i = x --- 628
	376	! y --- 636
	377	! z --- 627
	378
	379	k2a = 6.8d-12 * sqrt(tt) ! delta(e)< 42 cm-1
	380	k2b = 3.6d-11 * sqrt(tt) * exp(-65.6557/26.3) ! > 42 cm-1 see table 3
	381	k2a = k2a * rf2desac
	382	k2b = k2b * rf2desac
	383
	384	k2x = 6.8d-12 * sqrt(tt)
	385	k2y = 3.6d-11 * sqrt(tt) * exp(-65.6557/26.3)
	386	k2z = 6.8d-12 * sqrt(tt)
	387	k2x = k2x * rf2iso
	388	k2y = k2y * rf2iso
	389	k2z = k2z * rf2iso
	390	k2xp = k2x * exp( dble( -ee/tt * (nu(1,4)-nu(2,4)) ) )
	391	k2yp = k2y * exp( dble( -ee/tt * (nu(1,4)-nu(3,4)) ) )
	392	k2zp = k2z * exp( dble( -ee/tt * (nu(1,4)-nu(4,4)) ) )
	393
	394	! these are vt1 in table 3, paper i
	395	co2i(001) + m ---> co2i(0v0) + m k3
	396	co2i(001) + o3p ---> co2i(0v0) + o3p k7
	397	c k3vm(i) m = a --- co2 v = a --- 3 i = 1,2,3,4
	398	! b --- n2,o2 b --- 2
	399	! c --- co
	400	c k7v(i) v = a --- 3 i = 1,2,3,4
	401	! b --- 2
	402
	403	k7x = 2.0d-13*sqrt(tt/300.d0)
	404	k7x = k7x * rf7
	405
	406	if (iopt3.eq.0) then
	407
	408	k3x = 2.2d-15 + 1.14d-10 * exp( dble(-76.75/tt**(1.d0/3.d0)) )
	409	k3xaa = 3.0d-15 + 1.72d-10 * exp( dble(-76.75/tt**(1.d0/3.d0)))
	410	k3xac = 2.2d-15 + 9.66d-11 * exp( dble(-76.75/tt**(1.d0/3.d0)))
	411	k3x = k3x * rf3
	412	k3xaa = k3xaa * rf3
	413	k3xac = k3xac * rf3
	414
	415	tt1=220.0 !500.0 !220.0
	416	tt2=250.0 !250.0
	417	if(tt.le.tt1)then
	418	k3xab = k3x
	419	k3xbb = 0.d0
	420	k7xa=k7x
	421	k7xb=0.d0
	422	else if(tt.gt.tt2)then
	423	k3xab = 0.d0
	424	k3xbb = k3x
	425	k7xa=0.d0
	426	k7xb=k7x
	427	else
	428	k3xab = k3x/30.d0*(tt2-tt)
	429	k3xbb = k3x/30.d0*(tt-tt1)
	430	k7xa=k7x/30.d0*(tt2-tt)
	431	k7xb=k7x/30.d0*(tt-tt1)
	432	end if
	433	k3xba = k3xbb
	434	k3xbc = k3xbb
	435
	436	elseif (iopt3.gt.0) then ! bauer et. al., 1987
	437
	438	if (tt.ge.190. .and. tt.le.250.) then
	439	factor = 0.9d0 + dble( (0.1-0.9)/(250.-190.) * (tt-190.) )
	440	elseif (tt.lt.190.) then
	441	factor = 0.9d0
	442	elseif (tt.gt.250.) then
	443	factor = 0.1d0
	444	end if
	445
	446	k3xn2 = 2.2d-15 + 1.14d-10 * exp(dble( -76.75/tt**(1.d0/3.d0)))
	447	k3xn2 = k3xn2 * rf3
	448	k3xab = k3xn2 * factor
	449	k3xbb = k3xn2 * (1.-factor)
	450
	451	k7xa = k7x * factor
	452	k7xb = k7x * (1.-factor)
	453
	454	if (iopt3.eq.1) then
	455
	456	if (tt.le.148.8) then
	457	k3xco2 = 13.8 - 0.1807 * (tt-140.)
	458	elseif (tt.ge.148.8 .and. tt.le.159.6) then
	459	k3xco2 = 12.21 - 0.1787 * (tt-148.8)
	460	elseif (tt.ge.159.6 .and. tt.le.171.0) then
	461	k3xco2 = 10.28 - 0.04035 * (tt-159.6)
	462	elseif (tt.ge.171.0 .and. tt.le.186.4) then
	463	k3xco2 = 9.82 - 0.027273 * (tt-171.0)
	464	elseif (tt.ge.186.4 .and. tt.le.244.1) then
	465	k3xco2 = 9.4 + 0.002253 * (tt-186.4)
	466	elseif (tt.ge.244.1 .and. tt.le.300) then
	467	k3xco2 = 9.53 + 0.02129 * (tt-244.1)
	468	elseif (tt.ge.300 .and. tt.le.336.1) then
	469	k3xco2 = 10.72 + 0.052632 * (tt-300)
	470	elseif (tt.ge.336.1 .and. tt.le.397.0) then
	471	k3xco2 = 12.62 + 0.0844 * (tt-336.1)
	472	elseif (tt.ge.397.0 .and. tt.le.523.4) then
	473	k3xco2 = 17.76 + 0.2615 * (tt-397.)
	474	end if
	475	k3xco2 = k3xco2 * 1.d-15 * rf3
	476	k3xaa = 0.82d0 * k3xco2
	477	k3xba = 0.18d0 * k3xco2
	478
	479	if (tt.le.163.) then
	480	k3xco = 10.58 - 0.093 * (tt-140)
	481	elseif (tt.ge.163. .and. tt.le.180.) then
	482	k3xco = 8.44 - 0.05353 * (tt-163.)
	483	elseif (tt.ge.180. .and. tt.le.196.) then
	484	k3xco = 7.53 - 0.034375 * (tt-180.)
	485	elseif (tt.ge.196. .and. tt.le.248.) then
	486	k3xco = 6.98 - 0.0108 * (tt-196.)
	487	elseif (tt.ge.248. .and. tt.le.301.) then
	488	k3xco = 6.42 + 0.01415 * (tt-248.)
	489	elseif (tt.ge.301. .and. tt.le.353.) then
	490	k3xco = 7.17 + 0.02865 * (tt-301.)
	491	end if
	492	k3xac = k3xco * 1.d-15 * rf3
	493	k3xbc = 0.d0
	494
	495	elseif (iopt3.eq.2) then ! revision for the papers (feb 93)
	496
	497	k3xco2 = 7.3d-14 * exp( -850.3/tt + 86523/tt**2.d0 )
	498	k3xco2 = k3xco2 * rf3
	499	k3xaa = 0.82d0 * k3xco2
	500	k3xba = 0.18d0 * k3xco2
	501
	502	k3xco = 1.7d-14 * exp( -448.3/tt + 53636/tt**2.d0 )
	503	k3xac = k3xco * rf3
	504	k3xbc = 0.d0
	505
	506	end if
	507
	508	end if
	509
	510	do i=1,nisot
	511	k3aa(i) = k3xaa
	512	k3ba(i) = k3xba
	513	k3ab(i) = k3xab
	514	k3bb(i) = k3xbb
	515	k3ac(i) = k3xac
	516	k3bc(i) = k3xbc
	517	anu = nu(i,4)-nu(i,3)
	518	! anu = nu(i,4)-nu(i,3)+70
	519	k3aap(i) = k3aa(i) * exp( -ee/tt * anu )/6.d0
	520	k3abp(i) = k3ab(i) * exp( -ee/tt * anu )/6.d0
	521	k3acp(i) = k3ac(i) * exp( -ee/tt * anu )/6.d0
	522	anu = nu(i,4)-nu(i,2)
	523	! anu = nu(i,4)-nu(i,2)+40
	524	k3bap(i) = k3ba(i) * exp( -ee/tt * anu )/4.d0
	525	k3bbp(i) = k3bb(i) * exp( -ee/tt * anu )/4.d0
	526	k3bcp(i) = k3bc(i) * exp( -ee/tt * anu )/4.d0
	527
	528	k7a(i) = k7xa
	529	k7b(i) = k7xb
	530	k7ap(i) = k7a(i) * exp(dble( -ee/tt*(nu(i,4)-nu(i,3)) ))/6.d0
	531	k7bp(i) = k7b(i) * exp(dble( -ee/tt*(nu(i,4)-nu(i,2)) ))/4.d0
	532	end do
	533
	534
	535	! the next ones correspond to vv2 in table 4, paper i
	536	co2i(001) + co2 ---> co2i(020) + co2(010) k6
	537	! k6(i) i = 1,2,3,4
	538	! we need a new index for the inverse rates due to both fractions :
	539	c k6a(i) i=2,3,4 co2i(001) + co2 ---> co2i(020) + co2(010)
	540	c k6b(i) " co2i(001) + co2 ---> co2i(010) + co2(020)
	541
	542	if (iopt6.eq.1) then
	543
	544	if(tt.le.175.d0)then
	545	k6x=8.6d-15
	546	elseif(tt.gt.175.0.and.tt.le.200.d0)then
	547	k6x=8.6d-15+9.d-16*(175.d0-tt)/25.d0
	548	elseif(tt.gt.200.0.and.tt.le.225.d0)then
	549	k6x=7.7d-15+5.d-16*(200.d0-tt)/25.d0
	550	elseif(tt.gt.225.0.and.tt.le.250.d0)then
	551	k6x=7.20d-15+6.d-16*(tt-225.d0)/25.d0
	552	elseif(tt.gt.250.0.and.tt.le.275.d0)then
	553	k6x=7.80d-15+1.d-15*(tt-250.d0)/25.d0
	554	elseif(tt.gt.275.0.and.tt.le.300.d0)then
	555	k6x=8.80d-15+1.3d-15*(tt-275.d0)/25.d0
	556	elseif(tt.gt.300.0.and.tt.le.325.d0)then
	557	k6x=10.1d-15+1.54d-15*(tt-300.d0)/25.d0
	558	elseif(tt.gt.325.0)then
	559	k6x=11.6d-15
	560	end if
	561
	562	elseif (iopt6.eq.2) then
	563
	564	k6x = 3.6d-13 * exp( -1660/tt + 176948/tt**2.d0 )
	565	if (tt.lt.175) k6x = 8.8d-15
	566
	567	end if
	568
	569	k6x1 = k6x * rf6 * frac6
	570	k6x2 = k6x * rf6 * (1.-frac6)
	571
	572	k6 = k6x * rf6
	573	k6p = k6 / 8.d0 * exp(dble( -ee/tt * (nu(1,4)-nu(1,2)-nu(1,1)) ))
	574	do i=2,nisot
	575	k6a(i) = k6x1
	576	k6b(i) = k6x2
	577	anu = nu(i,4)-nu(i,2)-nu(1,1)
	578	k6ap(i) = k6a(i) / 8.d0 * exp(dble( -ee*anu/tt ))
	579	anu = nu(i,4)-nu(i,1)-nu(1,2)
	580	k6bp(i) = k6b(i) / 8.d0 * exp(dble( -ee*anu/tt ))
	581	end do
	582
	583
	584	co2i(0v0) + co2i ---> co2i(0v-10) + co2i(010)
	585	c k5 esta reaccion es desdenable frente a co2 como colisionante.
	586	! k5=3.0d-15+6.0d-17*(tt-210.d0)
	587	! k5=k5*rf5
	588	! k5p=k5/2.d0exp(-ee125.77/tt)
	589
	590
	591	co2i(0v0) + m ---> co2i(0v-10) + m k19 (vt2,k5,k6 in table 3, paper
	592	co2i(0v0) + o3p ---> co2i(0v-10) + o3p k20 (vt2,k7 in table 3, paper
	593	c k19vm(i) m = a --- co2 v = a --- 3 i=1,2,3,4
	594	! b --- n2 b --- 2
	595	! c --- co c --- 1
	596	c k20v(i) v = a --- 3 i = 1,2,3,4
	597	! b --- 2
	598	! c --- 1
	599	!
	600	! k20x=1.9d-8exp(-76.75/(tt*(1.d0/3.d0))) ! taylor,74 reajusted
	601	! k20x=2.32d-9exp(-76.75/(tt(1./3.)))+1.0d-14sqrt(tt) ! k&j, 83
	602	! k20x = 1.43d-12*(tt/300.d0) ! shved et al, 90
	603
	604	if (iopt20.eq.1) then ! first version of pap1
	605	k20x=2.32d-9exp(-76.75/(tt(1./3.)))+3.5d-13sqrt(tt) ! s&w, 91
	606	k20xb = k20x / 2.d0 * rf20
	607	k20xc = k20xb
	608	k20xa = 3.d0/2.d0 * k20xb
	609	elseif(iopt20.eq.2) then ! revision for the papers in feb 93
	610	k20x=3.d-12 ! minimum value of lopez-puertas et al., 92
	611	k20xc = k20x * rf20
	612	k20xb = 2.d0 * k20xc
	613	k20xa = 3.d0/2.d0 * k20xb
	614	elseif(iopt20.eq.3) then ! values from boug & roble '91
	615	k20x=1.d-12/sqrt(300.) * sqrt(tt)
	616	k20xc = k20x * rf20
	617	k20xb = 2.d0 * k20xc
	618	k20xa = 3.d0/2.d0 * k20xb
	619	elseif(iopt20.eq.4) then ! values from boug & dick '88 case b
	620	k20x=7.d-13
	621	k20xc = k20x * rf20
	622	k20xb = 2.d0 * k20xc
	623	k20xa = 3.d0/2.d0 * k20xb
	624	elseif(iopt20.eq.5) then ! values from s.bougher (oct-98)
	625	k20x = 1.732d-13 * sqrt(tt) ! 1/sqrt(300) * sqrt(t)
	626	k20xc = k20x * rf20
	627	k20xb = 2.d0 * k20xc
	628	k20xa = 3.d0/2.d0 * k20xb
	629	end if
	630
	631	if (iopt19.eq.0) then
	632
	633	k19xca = 4.64d-10 * exp(dble( - 74.75 / tt**(1.d0/3.d0) ))
	634	k19xcb = 6.69d-10 * exp(dble( - 84.07 / tt**(1.d0/3.d0) ))
	635	k19xcc = k19xcb
	636
	637	if ( tt.le.250 ) then
	638	k19xba = 181.25d0
	639	elseif ( tt.ge.310 ) then
	640	k19xba = 200.d0 + 0.9d0 * ( tt - 310.d0 )
	641	else
	642	k19xba = 181.25d0 + 0.3125d0 * ( tt - 250.d0 )
	643	end if
	644	k19xba = k19xba * 1.03558d-19 * tt ! cm-1 s-1
	645	k19xbb = 1.24d-14 * ( tt / 273.3d0 )**2.d0 ! taine & lepoutre 1979
	646	k19xbc = k19xbb
	647
	648	k19xaa = 3.d0/2.d0 * k19xba
	649	k19xab = 3.d0/2.d0 * k19xbb
	650	k19xac = 3.d0/2.d0 * k19xbc
	651
	652	k19xaa = k19xaa * rf19
	653	k19xab = k19xab * rf19
	654	k19xac = k19xac * rf19
	655	k19xba = k19xba * rf19
	656	k19xbb = k19xbb * rf19
	657	k19xbc = k19xbc * rf19
	658	k19xca = k19xca * rf19
	659	k19xcb = k19xcb * rf19
	660	k19xcc = k19xcc * rf19
	661
	662	elseif (iopt19.ge.1) then
	663
	664	if (iopt19.eq.1) then ! lunt et. al., 1985 (thesis values)
	665
	666	if (tt.le.175.) then
	667	k19xca = 4.d0 - 0.02d0 * (tt-140.d0)
	668	k19xcb = 0.494d0 + 0.0076 * (tt-140.d0)
	669	elseif (tt.ge.175. .and. tt.le.200.) then
	670	k19xca = 3.3d0 - 0.02d0 * (tt-175.)
	671	k19xcb = 0.76d0 + 0.0076d0 * (tt-175.d0)
	672	elseif (tt.ge.200. .and. tt.le.225.) then
	673	k19xca = 2.8d0 + 0.004d0 * (tt-200.d0)
	674	k19xcb = 0.95d0 + 0.014d0 * (tt-200.d0)
	675	elseif (tt.ge.225. .and. tt.le.250.) then
	676	k19xca = 2.9d0 + 0.024d0 * (tt-225.d0)
	677	k19xcb = 1.3d0 + 0.016d0 * (tt-225.d0)
	678	elseif (tt.ge.250. .and. tt.le.275.) then
	679	k19xca = 3.5d0 + 0.04d0 * (tt-250.d0)
	680	k19xcb = 1.7d0 + 0.032d0 * (tt-250.d0)
	681	elseif (tt.ge.275. .and. tt.le.295.) then
	682	k19xca = 4.5d0 + 0.055d0 * (tt-275.d0)
	683	k19xcb = 2.5d0 + 0.045d0 * (tt-275.d0)
	684	elseif (tt.ge.295. .and. tt.le.320.) then
	685	k19xca = 5.6d0 + 0.54d0 * (tt-295.d0)
	686	k19xcb = 3.4d0 + 0.045d0 * (tt-295.d0)
	687	end if
	688	k19xca = k19xca * 1.d-15 * rf19
	689	k19xcb = k19xcb * 1.d-15 * rf19
	690	k19xcc = k19xcb
	691
	692	elseif (iopt19.eq.2) then ! revision for the papers, feb 1993
	693
	694	! k19xca = 7.3d-14 * exp( -850.3d0/tt + 86523.d0/tt**2.d0 )
	695	k19xca = 4.2d-12 * exp( -2988.d0/tt + 303930.d0/tt**2.d0 )
	696	if (tt.le.175.) k19xca = 3.3d-15
	697	k19xcb = 2.1d-12 * exp( -2659.d0/tt + 223052.d0/tt**2.d0 )
	698	if (tt.le.175.) k19xcb = 7.6d-16
	699	k19xca = k19xca * rf19
	700	k19xcb = k19xcb * rf19
	701	k19xcc = k19xcb
	702
	703	elseif (iopt19.eq.3) then ! values from dick'72 for k19xc
	704	! k19xcb is not modified
	705	if (tt.le.158.) then
	706	k19xca = 0.724d-15
	707	elseif (tt.le.190.) then
	708	k19xca = 0.724d-15 +
	709	@ (1.1d-15-0.724d-15) * (tt-158.) / (190.-158.)
	710	elseif (tt.le.250.) then
	711	k19xca = 1.1d-15 +
	712	@ (3.45d-15-1.1d-15) * (tt-190.) / (250.-190.)
	713	elseif (tt.gt.250.) then
	714	k19xca = 3.45d-15
	715	end if
	716	k19xcb = 2.1d-12 * exp( -2659.d0/tt + 223052.d0/tt**2.d0 )
	717	if (tt.le.175.) k19xcb = 7.6d-16
	718	k19xca = k19xca * rf19
	719	k19xcb = k19xcb * rf19
	720	k19xcc = k19xcb
	721
	722	elseif (iopt19.eq.5) then
	723
	724	k19xca = 5.2d-15 ! s.bougher, nov-98
	725	k19xcb = 7.6d-16 ! nuestro, de feb-93
	726	k19xcc = k19xcb
	727
	728	k19xca = k19xca * rf19
	729	k19xcb = k19xcb * rf19
	730
	731	end if
	732
	733	factor = 2.5d0
	734	k19xba = factor * k19xca
	735	k19xbb = factor * k19xcb
	736	k19xbc = factor * k19xcc
	737	factor = 3.d0/2.d0
	738	k19xaa = factor * k19xba
	739	k19xab = factor * k19xbb
	740	k19xac = factor * k19xbc
	741
	742	end if
	743
	744	do i = 1, nisot
	745
	746	k19aa(i) = k19xaa
	747	k19ba(i) = k19xba
	748	k19ca(i) = k19xca
	749	k19ab(i) = k19xab
	750	k19bb(i) = k19xbb
	751	k19cb(i) = k19xcb
	752	k19ac(i) = k19xac
	753	k19bc(i) = k19xbc
	754	k19cc(i) = k19xcc
	755	anu = nu(i,3)-nu(i,2)
	756	k19aap(i) = k19aa(i) * 6.d0/4.d0 * exp(dble( -ee*anu/tt))
	757	k19abp(i) = k19ab(i) * 6.d0/4.d0 * exp(dble( -ee*anu/tt))
	758	k19acp(i) = k19ac(i) * 6.d0/4.d0 * exp(dble( -ee*anu/tt))
	759	anu = nu(i,2)-nu(i,1)
	760	k19bap(i) = k19ba(i) * 2.d0 * exp(dble( -ee*anu/tt))
	761	k19bbp(i) = k19bb(i) * 2.d0 * exp(dble( -ee*anu/tt))
	762	k19bcp(i) = k19bc(i) * 2.d0 * exp(dble( -ee*anu/tt))
	763	anu = nu(i,1)
	764	k19cap(i) = k19ca(i) * 2.d0 * exp(dble( -ee*anu/tt))
	765	k19cbp(i) = k19cb(i) * 2.d0 * exp(dble( -ee*anu/tt))
	766	k19ccp(i) = k19cc(i) * 2.d0 * exp(dble( -ee*anu/tt))
	767
	768	k20a(i) = k20xa
	769	k20b(i) = k20xb
	770	k20c(i) = k20xc
	771	k20ap(i) = k20a(i)6.d0/4.d0
	772	@ exp(dble( -ee/tt * (nu(i,3)-nu(i,2)) ))
	773	k20bp(i) = k20b(i)4.d0/2.d0
	774	@ exp(dble( -ee/tt * (nu(i,2)-nu(i,1)) ))
	775	k20cp(i) = k20c(i)2.d0/1.d0
	776	@ exp(dble( -ee/tt * nu(i,1) ))
	777	end do
	778
	779	! write(1,*) tt,k19cap(1),k19ac(1)
	780
	781	! the next ones correspond to vv3 in table 4 (paper i)
	782	co2i(0v0) + co2 ---> co2i(0v-10) + co2(010) k21 also see k33
	783	c k21v(i) v = a --- 3 i = 1,2,3,4
	784	! b --- 2
	785	! c --- 1
	786	! we need a new index for the 030i rates due to both fractions :
	787	c k21a1 co2i(030) + co2 ---> co2i(020) + co2(010)
	788	c k21a2 co2i(030) + co2 ---> co2i(010) + co2(020)
	789	co2i(010) + co2j(000) ---> co2i(000) + co2j(010) kijk21c see pag.22-s
	790	! k23k21c i=628,j=636
	791	! k24k21c i=628,j=627
	792	! k34k21c i=636,j=627
	793
	794	if (iopt21.eq.0) then
	795	k21x = 1.2d-11
	796	k21xb = k21x
	797	k21xa = 3.d0/2.d0 * k21x
	798	k21xc = k21x ! esta ultima no se usa con 626
	799	elseif (iopt21.eq.1) then
	800	k21x = 2.49d-11 ! orr & smith, 1987
	801	k21xb = k21x
	802	k21xa = 3.d0/2.d0 * k21xb ! oscilador armonico
	803	k21xc = k21xb / 2.d0 ! novedad mia
	804	elseif (iopt21.eq.2) then
	805	k21x = 100.d0*k19xca ! dickinson'76 (icarus)
	806	k21xb = k21x
	807	k21xa = 3.d0/2.d0 * k21xb ! oscilador armonico
	808	k21xc = k21xb / 2.d0 ! novedad mia
	809	end if
	810	k21xa_626 = k21xa * rf21a !* 0.01d0 !* 10.d0
	811	k21xa = k21xa * rf21a !* 0.01d0
	812	k21xb = k21xb * rf21b
	813	k21xc = k21xc * rf21c
	814
	815	k21a = k21xa_626
	816	k21ap = k21a * 6.d0/8.d0 *
	817	@ exp( dble( -ee/tt * (nu(1,3)-nu(1,2)-nu(1,1)) ) )
	818	do i = 2, nisot
	819	k21a1(i) = k21xa * frac21
	820	k21a2(i) = k21xa * (1.d0-frac21)
	821	k21a1p(i) = k21a1(i) * 6.d0/8.d0 *
	822	@ exp(dble( -ee/tt* (nu(i,3)-nu(i,2)-nu(1,1)) ))
	823	k21a2p(i) = k21a2(i) * 6.d0/8.d0 *
	824	@ exp(dble( -ee/tt* (nu(i,3)-nu(i,1)-nu(1,2)) ))
	825	end do
	826
	827
	828	do i = 1, nisot
	829	k21b(i) = k21xb
	830	k21c(i) = k21xc
	831	k21bp(i) = k21b(i) *
	832	@ exp(dble( -ee/tt* (nu(i,2)-nu(i,1)-nu(1,1)) ))
	833	k21cp(i) = k21c(i) *
	834	@ exp(dble( -ee/tt * (nu(i,1)-nu(1,1)) ))
	835	end do
	836
	837	k23k21c = k21xc
	838	k24k21c = k21xc
	839	k34k21c = k21xc
	840	k23k21cp = k23k21c2.d0/2.d0
	841	@ exp(dble( -ee/tt* (nu(2,1)-nu(3,1)) ))
	842	k24k21cp = k24k21c2.d0/2.d0
	843	@ exp(dble( -ee/tt* (nu(2,1)-nu(4,1)) ))
	844	k34k21cp = k34k21c2.d0/2.d0
	845	@ exp(dble( -ee/tt* (nu(3,1)-nu(4,1)) ))
	846
	847	! these are also vv3 in table 4, paper i
	848	c k31 & k32
	849
	850	k31 = k21x * rf31 ! we're suposing thar the rate for the deactivation
	851	! v-v from high combinational levels is the same
	852	k32 = k21x * rf32 ! that the one for : (020) --> (010) + (010)
	853
	854	c o2(*) + co2i ---> co2() + co2i(**) k33
	855	c k33a1 : co2(001) + co2i ---> co2(020) + co2i(010) (vv2, table 4,
	856	! a2 : co2(001) + co2i ---> co2(010) + co2i(020) "
	857	! b1 : co2(030) + co2i ---> co2(020) + co2i(010) (vv3, table 4,
	858	! b2 : co2(030) + co2i ---> co2(010) + co2i(020) "
	859	! c : co2(020) + co2i ---> co2(010) + co2i(010) "
	860	! we have to add an index to the inverse rates, depending on the isotope
	861
	862	k33c = k21x * rf33bc
	863	k33b1 = 3.d0/3.d0 * k33c * frac33
	864	k33b2 = 3.d0/3.d0 * k33c * (1.d0-frac33)
	865	k33a1 = k6x * rf33a * frac33
	866	k33a2 = k6x * rf33a * (1.d0-frac33)
	867
	868	do i=2,nisot
	869	k33a1p(i)=k33a1*
	870	@ 1.d0/8.d0exp(dble( -ee/tt (nu(1,4)-nu(1,2)-nu(i,1)) ))
	871	k33a2p(i)=k33a2*
	872	@ 1.d0/8.d0exp(dble( -ee/tt (nu(1,4)-nu(1,1)-nu(i,2)) ))
	873	k33b1p(i)=k33b1*
	874	@ 6.d0/8.d0exp(dble( -ee/tt (nu(1,3)-nu(1,2)-nu(i,1)) ))
	875	k33b2p(i)=k33b2*
	876	@ 6.d0/8.d0exp(dble( -ee/tt (nu(1,4)-nu(1,1)-nu(i,2)) ))
	877	k33cp(i) =k33c * exp(dble( -ee/tt * (nu(1,2)-nu(1,1)-nu(i,1)) ))
	878	end do
	879
	880	! here they are the vt3 in table 3, paper i
	881	co2(2.7um) + m ---> co2(2.7um) + m k27
	882	c k27a : n8 + m ---> n6 + m
	883	! b : n7 + m ---> n6 + m
	884	! c : n8 + m ---> n7 + m
	885
	886	if (iopt27.eq.0) then
	887	k27a = 3.d-11 !between fermi levels
	888	k27b = 3.d-13 !between side levels
	889	k27c = 2.d0 * k27b !between side levels
	890	elseif (iopt27.ge.1) then ! orr & smith, 1987
	891	k27a = 1.55d-12
	892	k27c = 4.97d-12
	893	k27b = k27c
	894	end if
	895	k27a = k27a * rf27f
	896	k27b = k27b * rf27s
	897	k27c = k27c * rf27s
	898
	899	k27ap = k27a * exp(dble( -ee/tt * (nu(1,8)-nu(1,6)) ))
	900	k27bp = k27b * exp(dble( -ee/tt * (nu(1,7)-nu(1,6)) ))
	901	k27cp = k27c * exp(dble( -ee/tt * (nu(1,8)-nu(1,7)) ))
	902
	903
	904	! the next two are not used in the model:
	905
	906	c k28 : n* + n2 ---> n*low + n2(1)
	907	! k28v v = a --- n8
	908	! b --- n7
	909	! c --- n6
	910	! k28a = 5.d-13 * sqrt(300.d0/tt) * rf28 ! = k1
	911	! k28b = k28a
	912	! k28c = k28a
	913	! k28ap = k28a * exp( -ee/tt * (nu(1,8)-1388.1847-nun2) )
	914	! k28bp = k28b * exp( -ee/tt * (nu(1,7)-1335.1317-nun2) )
	915	! k28cp = k28c * exp( -ee/tt * (nu(1,6)-1285.4087-nun2) )
	916
	917	c k29 : n* + co ---> n*low + co(1)
	918	! k29v v = a --- n8
	919	! b --- n7
	920	! c --- n6
	921
	922	c k29a = ?????????? * rf29
	923	c k29b = k29a
	924	c k29c = k29a
	925
	926	c k29ap = k29a * exp( -ee/tt * (nu(1,8)-1388.1847-nuco) )
	927	c k29bp = k29b * exp( -ee/tt * (nu(1,7)-1335.1317-nuco) )
	928	c k29cp = k29c * exp( -ee/tt * (nu(1,6)-1285.4087-nuco) )
	929
	930
	931	! these are also vv1 processes in table 4, paper i
	932	c k26 :
	933	! 1. deactivation of the 626 isotope:
	934	! reaction: n* + co2i ---> n*low + co2i(001) ; i=1-4
	935	! nomenclature: k26v ; v=a,b,c,d for n8,n7,n6,n5 respectively
	936	! inverse rates: k26vp(i) ; i=1-4
	937	! 2. deactivation of the minor isotopes:
	938	! reaction: n_i + co2j ---> n_i_low + (001)_j ; i=2-4 ; j=1-4
	939	! nomenclature: k26vij ; v=a,c,d for n8,n6,n5 respectively
	940	! inverse rates: k26vijp
	941	! 3. notes:
	942	! a * it is clear that : k26vij=/=k26vjip,
	943	! b * at the moment we do not include inverse rates for the case 2., o
	944	! the deactivation of the main isotope (pg. 32b, sn1).
	945	! c * not 0221 level for minor isotopes is considered.
	946	! d * only a value is known for these rates, so all of the deactivatio
	947	! the same, but not the inverse rates.
	948	! e * although all the direct deactivation constants have the same val
	949	! it is useful to distinguish between them with the present nam
	950
	951	k26a = 6.8d-12 * sqrt(tt) * rf26 ! = k2
	952	k26b = k26a
	953	k26c = k26a
	954	if (iopt26.eq.0 .or. iopt26.eq.2) then
	955	k26d = k26a
	956	elseif( iopt26.eq.1) then
	957	k26d = 1.15d-10 * rf26
	958	end if
	959
	960	do i=1,4
	961	k26ap(i) = k26a *
	962	@ exp(dble( -ee/tt * (nu(1,8)-nu12_1000-nu(i,4)) ))
	963	k26bp(i) = k26b *
	964	@ exp(dble( -ee/tt * (nu(1,7)- nu(1,2) -nu(i,4)) ))
	965	k26cp(i) = k26c *
	966	@ exp(dble( -ee/tt * (nu(1,6)-nu12_0200-nu(i,4)) ))
	967	k26dp(i) = k26d *
	968	@ exp(dble( -ee/tt * (nu(1,5)- nu(1,1) -nu(i,4)) ))
	969	end do
	970
	971	k26a21 = k26a
	972	k26c21 = k26c
	973	k26d21 = k26d
	974	k26a22 = k26a
	975	k26c22 = k26c
	976	k26d22 = k26d
	977	k26a23 = k26a
	978	k26c23 = k26c
	979	k26d23 = k26d
	980	k26a24 = k26a
	981	k26c24 = k26c
	982	k26d24 = k26d
	983
	984	k26a31 = k26a
	985	k26c31 = k26c
	986	k26d31 = k26d
	987	k26a32 = k26a
	988	k26c32 = k26c
	989	k26d32 = k26d
	990	k26a33 = k26a
	991	k26c33 = k26c
	992	k26d33 = k26d
	993	k26a34 = k26a
	994	k26c34 = k26c
	995	k26d34 = k26d
	996
	997	k26a41 = k26a
	998	k26c41 = k26c
	999	k26d41 = k26d
	1000	k26a42 = k26a
	1001	k26c42 = k26c
	1002	k26d42 = k26d
	1003	k26a43 = k26a
	1004	k26c43 = k26c
	1005	k26d43 = k26d
	1006	k26a44 = k26a
	1007	k26c44 = k26c
	1008	k26d44 = k26d
	1009
	1010	!! some examples of inverse rates, although not used at the moment
	1011	! k26a32p = k26a32 * exp( -ee* (nu(3,8)-nu32_1000-nu(2,4)) / tt )*1./1.
	1012	! k26c32p = k26c32 * exp( -ee* (nu(3,6)-nu32_0200-nu(2,4)) / tt )*1./1.
	1013	! k26d32p = k26d32 * exp( -ee* (nu(3,5)- nu(3,1) -nu(2,4)) / tt )*2./2.
	1014	!
	1015	! k26a43 = k26a34 * exp( -ee* (nu(3,8)-nu32_1000-nu(4,4)) / tt )*1./1.
	1016	! k26c43 = k26c34 * exp( -ee* (nu(3,6)-nu32_0200-nu(4,4)) / tt )*1./1.
	1017	! k26d43 = k26d34 * exp( -ee* (nu(3,5)- nu(3,1) -nu(4,4)) / tt )*2./2.
	1018	!
	1019	! k26a24p = k26a24 * exp( -ee* (nu(2,8)-nu22_1000-nu(4,4)) / tt )*1./1.
	1020	! k26c24p = k26c24 * exp( -ee* (nu(2,6)-nu22_0200-nu(4,4)) / tt )*1./1.
	1021	! k26d24p = k26d24 * exp( -ee* (nu(2,5)- nu(2,1) -nu(4,4)) / tt )*2./2.
	1022
	1023
	1024	! this is taken as vv4 in table 4, paper i (in the inverse direction)
	1025	c k41 : co(v) + co2 ---> co(v-1) + co2(001) + de
	1026	! k41_v v=1,2,3,4
	1027	!
	1028	! de = -205.9 cm-1 when v=1
	1029	! de = -232.9 cm-1 when v=2
	1030	! de = -258.6 cm-1 when v=3
	1031	! de = -285.0 cm-1 when v=4
	1032
	1033	k41p_taylor = 1.56d-11 * exp( -30.12/tt**0.333333 ) ! [ s-1 cm+3 ]
	1034	k41p_shved = 7.5d7/sqrt(tt) ! [ s-1 atm-1 ]
	1035	k41p_shved = k41p_shved * 1.38d-16/1013250. * tt ! [ s-1 cm+3 ]
	1036	k41p_starr_hannock = 6.27d3 ! [ s-1 torr-1 ]
	1037
	1038	if (iopt41.eq.1) then
	1039	k41p_1 = k41p_starr_hannock *
	1040	@ 760.1.38d-16/1013250. tt ! [ s-1 cm+3 ]
	1041	elseif (iopt41.eq.2) then
	1042	k41p_1 = 1.6d-12 * exp( -1169/tt + 77601/tt**2.d0 )
	1043	end if
	1044	k41p_1 = k41p_1 * rf41
	1045	k41_1 = k41p_1 * exp(dble( -ee * 205.9/tt ))
	1046	k41_2 = k41_1
	1047	k41p_2 = k41_2 * exp(dble( -ee * (-232.9)/tt ))
	1048	k41_3 = k41_1
	1049	k41p_3 = k41_3 * exp(dble( -ee * (-258.6)/tt ))
	1050	k41_4 = k41_1
	1051	k41p_4 = k41_4 * exp(dble( -ee * (-285.0)/tt ))
	1052
	1053	!k41p_1 = k41p_1 * 1.d-6
	1054	!k41p_2 = k41p_2 * 1.d-6
	1055
	1056	c k41iso : 63(v) + co2 ---> 63(v-1) + co2(001) + de
	1057	! k41iso_v v=1,2,3
	1058	! de = -253 cm-1 when v=1
	1059	! de = -278 cm-1 when v=2
	1060	! de = -303 cm-1 when v=3
	1061
	1062	k41iso_1 = k41_1
	1063	k41iso_1p = k41iso_1 * exp(dble( -ee * (-253.)/tt ))
	1064	k41iso_2 = k41iso_1
	1065	k41iso_2p = k41iso_2 * exp(dble( -ee * (-278.)/tt ))
	1066	k41iso_3 = k41iso_1
	1067	k41iso_3p = k41iso_3 * exp(dble( -ee * (-303.)/tt ))
	1068
	1069
	1070
	1071	c k42 : co(v) + co ---> co(v-1) + co(1) + de=-26.481 si v=2 K42
	1072	! -52.8940 si v=3 k42b
	1073	! -79.2402 si v=4 k42c
	1074	! tomado de stepanova & shved (ellos de powell, 1975), ver pg .. l5
	1075	! solo para v=2 :
	1076
	1077	k42 = 2.89d-10 * (1./sqrt(tt) + 67.4/tt*1.5)
	1078	@ exp(dble(24.78/tt))
	1079	k42 = k42 * rf42
	1080	k42b = k42
	1081	k42c = k42
	1082	k42p = k42 * exp(dble( -ee * (-26.481)/tt ))
	1083	k42bp = k42b * exp(dble( -ee * (-52.894)/tt ))
	1084	k42cp = k42c * exp(dble( -ee * (-79.24)/tt ))
	1085
	1086	c k42iso : 63(v) + 63 ---> 63(v-1) + 63(1) + de=-25.31 si v=2 K42iso
	1087	! -50.57 si v=3 k42isob
	1088	! tomado de stepanova & shved (ellos de powell, 1975), ver pg .. l5
	1089	! solo para v=2 :
	1090
	1091	k42iso = k42
	1092	k42isop = k42iso * exp(dble( -ee * (-25.31)/tt ))
	1093	k42isob = k42
	1094	k42isobp = k42isob * exp(dble( -ee * (-50.57)/tt ))
	1095
	1096
	1097	c k43 : co(v) + o3p ---> co(v-1) + o3p + de=2143
	1098	! tomado de lewittes et. al, 1978 para v=1
	1099
	1100	if (iopt43.eq.1) then
	1101	tt1 = tt - 300.
	1102	k43 = 2.85d-14 * exp( dble( 9.5e-3tt1 + 1.11e-4tt1**2. ) )
	1103	elseif (iopt43.eq.2) then
	1104	k43 = 1.4d-5 * exp( -10957.d0 / tt + 1.486d6 / tt**2.d0 )
	1105	if ( tt.lt.265.0 ) k43 = 2.3d-14
	1106	end if
	1107	k43 = k43 * rf43
	1108	k43p = k43 * exp( -ee * dble(2143.3 / tt) )
	1109
	1110	c k43iso : co63(v) + o3p ---> co63(v-1) + o3p + de=2096
	1111	! Por similitud con el anterior
	1112
	1113	k43iso = k43
	1114	k43isop = k43iso * exp( -ee * dble(2096. / tt) )
	1115
	1116
	1117	c k44 : co62(v) + co63 ---> co62(v-1) + co63(1) + de
	1118	! basado en Lopez-Valverde et al para el caso v=1, solo usamos este
	1119	! k44x x = a --- v=1 de= 147.33
	1120	! b --- v=2 de= 20.7241
	1121	! c --- v=3 de= -5.7
	1122	! d --- v=4 de= -32.0361
	1123
	1124	k44a = 2.0d-12 * rf44 ! Solo vamos a usar este, no los b,c,d
	1125	k44b = k44a
	1126	k44c = k44a
	1127	k44d = k44a
	1128
	1129	de = 147.33
	1130	k44ap = k44a * exp(dble( -ee * de/tt ))
	1131	de = 20.7241
	1132	k44bp = k44b * exp(dble( -ee * de/tt ))
	1133	de = -5.7
	1134	k44cp = k44c * exp(dble( -ee * de/tt ))
	1135	de = -32.0361
	1136	k44dp = k44d * exp(dble( -ee * de/tt ))
	1137
	1138
	1139	co2(hcl) + co2 --> co2 + co2 + de(hcl)
	1140	! este rate tambien lo usamos para los high combination levels (para tra
	1141	! al lte. cualquier valor vale, supongo. es k_vthcl
	1142
	1143	k_vthcl = 3.3d-15 ! similar al valor pequenho del vt2
	1144	k_vthcl = k_vthcl * rf_hcl
	1145
	1146	return
	1147	end

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

Download in other formats:

Original Format