Context Navigation

photochemistry.F @ 1430

Last change on this file since 1430 was 1430, checked in by flefevre, 10 years ago

aeronomars : preparation de l'implementation du solveur chimique ASIS

1) extension de la table de photodissociation "jmars" a la thermosphere

et restriction aux seules especes necessaires.

2) modernisation de la subroutine d'interpolation dans cette

table de photodissociation.

3) chimiedata.h compatible avec solveur ASIS

ATTENTION!! la nouvelle table de photodissociation

jmars.20140930

doit etre presente dans /datadir sinon crash modele...

File size: 50.0 KB

Rev	Line
[334]	1	c*****************************************************************
	2	c
	3	c Photochemical routine
	4	c
	5	c Author: Franck Lefevre
	6	c ------
	7	c
[1125]	8	c Version: 11/12/2013
[334]	9	c
	10	c*****************************************************************
[1125]	11
	12	subroutine photochemistry(nlayer, nq,
	13	$ lswitch, zycol, sza, ptimestep,
[690]	14	$ press,temp, dens, dist_sol, surfdust1d,
	15	$ surfice1d, jo3, tau)
[1125]	16
[334]	17	implicit none
[1125]	18
[334]	19	#include "chimiedata.h"
	20	#include "callkeys.h"
[1125]	21
[334]	22	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	23	c inputs:
[334]	24	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	25
	26	integer, intent(in) :: nlayer ! number of atmospheric layers
	27	integer, intent(in) :: nq ! number of tracers
	28
	29	real :: sza ! solar zenith angle (deg)
	30	real :: ptimestep ! physics timestep (s)
	31	real :: press(nlayer) ! pressure (hpa)
	32	real :: temp(nlayer) ! temperature (k)
	33	real :: dens(nlayer) ! density (cm-3)
	34	real :: dist_sol ! sun distance (au)
	35	real :: surfdust1d(nlayer) ! dust surface area (cm2/cm3)
	36	real :: surfice1d(nlayer) ! ice surface area (cm2/cm3)
	37	real :: tau ! optical depth at 7 hpa
	38
[334]	39	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	40	c input/output:
[334]	41	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	42
	43	real :: zycol(nlayer,nq) ! chemical species volume mixing ratio
	44
[334]	45	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	46	c output:
	47	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	48
	49	real :: jo3(nlayer) ! photodissociation rate o3 -> o1d
	50
[334]	51	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	52	c local:
	53	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	54
	55	integer, parameter :: nesp = 16 ! number of species in the chemistry
	56
	57	integer :: phychemrat ! ratio physics timestep/chemistry timestep
	58	integer :: istep
	59	integer :: i_co2,i_o3,j_o3_o1d,lev
	60	integer :: lswitch
	61
	62	real :: stimestep ! standard timestep for the chemistry (s)
	63	real :: ctimestep ! real timestep for the chemistry (s)
	64	real :: rm(nlayer,nesp) ! species volume mixing ratio
	65	real :: j(nlayer,nd) ! interpolated photolysis rates (s-1)
	66	real :: rmco2(nlayer) ! co2 mixing ratio
	67	real :: rmo3(nlayer) ! ozone mixing ratio
	68
[334]	69	c reaction rates
[1125]	70
	71	real :: a001(nlayer), a002(nlayer), a003(nlayer)
	72	real :: b001(nlayer), b002(nlayer), b003(nlayer),
	73	$ b004(nlayer), b005(nlayer), b006(nlayer),
	74	$ b007(nlayer), b008(nlayer), b009(nlayer)
	75	real :: c001(nlayer), c002(nlayer), c003(nlayer),
	76	$ c004(nlayer), c005(nlayer), c006(nlayer),
	77	$ c007(nlayer), c008(nlayer), c009(nlayer),
	78	$ c010(nlayer), c011(nlayer), c012(nlayer),
	79	$ c013(nlayer), c014(nlayer), c015(nlayer),
	80	$ c016(nlayer), c017(nlayer), c018(nlayer)
	81	real :: d001(nlayer), d002(nlayer), d003(nlayer)
	82	real :: e001(nlayer), e002(nlayer), e003(nlayer)
	83	real :: h001(nlayer), h002(nlayer), h003(nlayer),
	84	$ h004(nlayer), h005(nlayer)
	85	real :: t001(nlayer), t002(nlayer), t003(nlayer)
	86
[334]	87	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[618]	88	c stimestep : standard timestep for the chemistry (s) c
	89	c ctimestep : real timestep for the chemistry (s) c
	90	c phychemrat : ratio physical/chemical timestep c
[334]	91	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	92
[618]	93	stimestep = 600. ! standard value : 10 mn
[1125]	94
[618]	95	phychemrat = nint(ptimestep/stimestep)
[1125]	96
[408]	97	ctimestep = ptimestep/real(phychemrat)
[1125]	98
[334]	99	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	100	c initialisation of chemical species and families c
	101	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	102
	103	call gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm)
	104
[334]	105	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	106	c compute reaction rates c
	107	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	108
	109	call chemrates(nlayer,
	110	$ lswitch, dens, press, temp,
[334]	111	$ surfdust1d, surfice1d,
	112	$ a001, a002, a003,
	113	$ b001, b002, b003, b004, b005, b006,
	114	$ b007, b008, b009,
	115	$ c001, c002, c003, c004, c005, c006,
	116	$ c007, c008, c009, c010, c011, c012,
	117	$ c013, c014, c015, c016, c017, c018,
	118	$ d001, d002, d003,
	119	$ e001, e002, e003,
	120	$ h001, h002, h003, h004, h005,
	121	$ t001, t002, t003, tau)
[1125]	122
[334]	123	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	124	c interpolation of photolysis rates in the lookup table c
	125	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	126
[334]	127	i_co2 = 1
	128	i_o3 = 6
[1125]	129
[334]	130	do lev = 1,lswitch-1
	131	rmco2(lev) = rm(lev,i_co2)
	132	rmo3(lev) = rm(lev, i_o3)
	133	end do
[1125]	134
[1430]	135	call photolysis(nlayer, lswitch, press, temp, sza, tau, dist_sol,
	136	$ rmco2, rmo3, j)
[1125]	137
[334]	138	j_o3_o1d = 5
[1125]	139
[334]	140	do lev = 1,lswitch-1
	141	jo3(lev) = j(lev,j_o3_o1d)
	142	end do
[1125]	143
[334]	144	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	145	c loop over time within the physical timestep c
	146	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	147
[334]	148	do istep = 1,phychemrat
[1125]	149	call chimie(nlayer, lswitch, nesp, rm, j, dens, ctimestep,
[334]	150	$ press, temp, sza, dist_sol,
	151	$ a001, a002, a003,
	152	$ b001, b002, b003, b004, b005, b006,
	153	$ b007, b008, b009,
	154	$ c001, c002, c003, c004, c005, c006,
	155	$ c007, c008, c009, c010, c011, c012,
	156	$ c013, c014, c015, c016, c017, c018,
	157	$ d001, d002, d003,
	158	$ e001, e002, e003,
	159	$ h001, h002, h003, h004, h005,
	160	$ t001, t002, t003)
	161	end do
[1125]	162
[334]	163	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	164	c save chemical species for the gcm c
	165	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	166
	167	call chimtogcm(nlayer, nq, zycol, lswitch, nesp, rm)
	168
[334]	169	return
	170	end
[1125]	171
[334]	172	c*****************************************************************
[1125]	173
	174	subroutine chimie(nlayer,
	175	$ lswitch, nesp, rm, j, dens, dt,
[334]	176	$ press, t, sza, dist_sol,
	177	$ a001, a002, a003,
	178	$ b001, b002, b003, b004, b005, b006,
	179	$ b007, b008, b009,
	180	$ c001, c002, c003, c004, c005, c006,
	181	$ c007, c008, c009, c010, c011, c012,
	182	$ c013, c014, c015, c016, c017, c018,
	183	$ d001, d002, d003,
	184	$ e001, e002, e003,
	185	$ h001, h002, h003, h004, h005,
	186	$ t001, t002, t003)
[1125]	187
[334]	188	c*****************************************************************
[1125]	189
[334]	190	implicit none
[1125]	191
[334]	192	#include "chimiedata.h"
	193	#include "callkeys.h"
[1125]	194
[334]	195	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	196	c inputs:
	197	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	198
	199	integer, intent(in) :: nlayer ! number of atmospheric layers
	200	integer :: lswitch ! interface level between chemistries
	201	integer :: nesp ! number of species
	202
	203	real :: dens(nlayer) ! density (cm-3)
	204	real :: dt ! chemistry timestep (s)
	205	real :: j(nlayer,nd) ! interpolated photolysis rates (s-1)
	206	real :: press(nlayer) ! pressure (hpa)
	207	real :: t(nlayer) ! temperature (k)
	208	real :: sza ! solar zenith angle (deg)
	209	real :: dist_sol ! sun distance (au)
	210
[334]	211	c reaction rates
[1125]	212
	213	real :: a001(nlayer), a002(nlayer), a003(nlayer)
	214	real :: b001(nlayer), b002(nlayer), b003(nlayer),
	215	$ b004(nlayer), b005(nlayer), b006(nlayer),
	216	$ b007(nlayer), b008(nlayer), b009(nlayer)
	217	real :: c001(nlayer), c002(nlayer), c003(nlayer),
	218	$ c004(nlayer), c005(nlayer), c006(nlayer),
	219	$ c007(nlayer), c008(nlayer), c009(nlayer),
	220	$ c010(nlayer), c011(nlayer), c012(nlayer),
	221	$ c013(nlayer), c014(nlayer), c015(nlayer),
	222	$ c016(nlayer), c017(nlayer), c018(nlayer)
	223	real :: d001(nlayer), d002(nlayer), d003(nlayer)
	224	real :: e001(nlayer), e002(nlayer), e003(nlayer)
	225	real :: h001(nlayer), h002(nlayer), h003(nlayer),
	226	$ h004(nlayer), h005(nlayer)
	227	real :: t001(nlayer), t002(nlayer), t003(nlayer)
	228
[334]	229	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	230	c input/output:
	231	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	232
	233	real :: rm(nlayer,nesp) ! volume mixing ratios
	234
	235	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[334]	236	c local:
	237	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	238
	239	real, parameter :: hetero_ice = 1. ! switch for het. chem. on ice clouds
	240	real, parameter :: hetero_dust = 0. ! switch for het. chem. on dust
	241
	242	integer :: iesp, iter, l
	243	integer :: i_co2, i_co, i_o2, i_h2, i_h2o, i_h2o2, i_hox, i_ox,
	244	$ i_o1d, i_o, i_o3, i_h, i_oh, i_ho2, i_ch4, i_n2
	245	integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d,
	246	$ j_o3_o1d, j_o3_o, j_h2o, j_hdo, j_h2o2,
[1430]	247	$ j_ho2, j_no, j_no2, j_ch4_ch3_h, j_ch4_1ch2_h2,
[1125]	248	$ j_ch4_3ch2_h_h, j_ch4_ch_h2_h, j_ch3o2h,
	249	$ j_ch2o_co, j_ch2o_hco, j_ch3oh, j_c2h6, j_hcl,
	250	$ j_hocl, j_clo, j_so2, j_so, j_h2s, j_so3,
	251	$ j_hno3, j_hno4
	252
	253	integer, parameter :: niter = 5 ! iterations in the chemical scheme
	254
	255	real :: cc0(nlayer,nesp) ! initial number density (cm-3)
	256	real :: cc(nlayer,nesp) ! final number density (cm-3)
	257	real :: nox(nlayer) ! nox number density (cm-3)
	258	real :: no(nlayer) ! no number density (cm-3)
	259	real :: no2(nlayer) ! no2 number density (cm-3)
	260	real :: production(nlayer,nesp) ! production rate
	261	real :: loss(nlayer,nesp) ! loss rate
	262
	263	real :: ro_o3, rh_ho2, roh_ho2, rno2_no
	264
[334]	265	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	266	c tracer numbering in the chemistry
	267	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	268
[334]	269	i_co2 = 1
	270	i_co = 2
	271	i_o = 3
	272	i_o1d = 4
	273	i_o2 = 5
	274	i_o3 = 6
	275	i_h = 7
	276	i_h2 = 8
	277	i_oh = 9
	278	i_ho2 = 10
	279	i_h2o2 = 11
	280	i_ch4 = 12
	281	i_h2o = 13
	282	i_n2 = 14
	283	i_hox = 15
	284	i_ox = 16
[1125]	285
[334]	286	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	287	c numbering of photolysis rates
	288	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	289
[334]	290	j_o2_o = 1 ! o2 + hv -> o + o
	291	j_o2_o1d = 2 ! o2 + hv -> o + o(1d)
	292	j_co2_o = 3 ! co2 + hv -> co + o
	293	j_co2_o1d = 4 ! co2 + hv -> co + o(1d)
	294	j_o3_o1d = 5 ! o3 + hv -> o2 + o(1d)
	295	j_o3_o = 6 ! o3 + hv -> o2 + o
	296	j_h2o = 7 ! h2o + hv -> h + oh
[1430]	297	j_h2o2 = 8 ! h2o2 + hv -> oh + oh
	298	j_ho2 = 9 ! ho2 + hv -> oh + o
	299	j_no = 10 ! no + hv -> n + o
[334]	300	j_no2 = 11 ! no2 + hv -> no + o
[1430]	301	j_hno3 = 12 ! hno3 + hv -> oh + no2
	302	j_hno4 = 13 ! hno4 + hv -> ho2 + no2
[1125]	303
[334]	304	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	305	c volume mixing ratio -> density conversion
	306	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	307
[334]	308	do iesp = 1,nesp
	309	do l = 1,lswitch-1
	310	cc0(l,iesp) = rm(l,iesp)*dens(l)
	311	cc(l,iesp) = cc0(l,iesp)
	312	end do
	313	end do
[1125]	314
[334]	315	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	316	c co2 and nox number densities (cm-3)
	317	c
	318	c nox mixing ratio: 6.e-10 (yung and demore, 1999)
	319	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	320
[334]	321	do l = 1,lswitch-1
	322	nox(l) = 6.e-10*dens(l)
	323	end do
[1125]	324
[334]	325	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	326	c loop over iterations in the chemical scheme
	327	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	328
[334]	329	do iter = 1,niter
[1125]	330
[334]	331	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	332	c nox species
	333	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	334	c no2/no
	335	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	336
[334]	337	do l = 1,lswitch-1
[1125]	338
[334]	339	rno2_no = (d002(l)cc(l,i_o3) + d003(l)cc(l,i_ho2))
	340	$ /(j(l,j_no2) +
	341	$ d001(l)max(cc(l,i_o),1.e-30dens(l)))
[1125]	342
[334]	343	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	344	c no
	345	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	346
[334]	347	no(l) = nox(l)/(1. + rno2_no)
[1125]	348
[334]	349	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	350	c no2
	351	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	352
[334]	353	no2(l) = no(l)*rno2_no
[1125]	354
[334]	355	end do
[1125]	356
[334]	357	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	358	c hox species
	359	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	360	c photochemical equilibrium for oh and ho2
	361	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	362	c h/ho2
	363	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	364
[334]	365	do l = 1,lswitch-1
[1125]	366
[334]	367	rh_ho2 = (c001(l)*cc(l,i_o)
	368	$ + c004(l)*cc(l,i_h)
	369	$ + c005(l)*cc(l,i_h)
	370	$ + c006(l)*cc(l,i_h)
	371	$ + c007(l)*cc(l,i_oh)
	372	$ + 2.c008(l)cc(l,i_ho2)
	373	$ + c015(l)*cc(l,i_o3)
	374	$ + 2.c016(l)cc(l,i_ho2)
	375	$ + d003(l)*no(l) ! ajout 20090401
	376	$ + j(l,j_ho2)
	377	$ + h001(l)*hetero_ice
	378	$ + h003(l)*hetero_dust)
	379	$ /( c011(l)*cc(l,i_o2)
	380	$ + t001(l)*cc(l,i_h2o)
	381	$ /max(cc(l,i_h),dens(l)*1.e-30) ! ajout 20090401
	382	$ )
[1125]	383
[334]	384	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	385	c oh/ho2
	386	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	387
[334]	388	roh_ho2 = (c001(l)*cc(l,i_o)
	389	$ + c003(l)cc(l,i_o3)rh_ho2
	390	$ + 2.c004(l)cc(l,i_h)
	391	$ + 2.c008(l)cc(l,i_ho2)
	392	$ + c015(l)*cc(l,i_o3)
	393	$ + d003(l)*no(l)
	394	$ + j(l,j_ho2)
	395	$ + 2.b002(l)cc(l,i_o1d)*cc(l,i_h2o) ! ajout 20101210
	396	$ /max(cc(l,i_ho2),dens(l)*1.e-30) ! ajout 20101210
	397	$ + b003(l)cc(l,i_o1d)cc(l,i_h2) ! ajout 20101210
	398	$ /max(cc(l,i_ho2),dens(l)*1.e-30) ! ajout 20101210
	399	$ + j(l,j_h2o)*cc(l,i_h2o)
	400	$ /max(cc(l,i_ho2),dens(l)*1.e-30)
	401	$ + t001(l)*cc(l,i_h2o) ! suppression 20090401
	402	$ /max(cc(l,i_ho2),dens(l)*1.e-30) ! suppression 20090401
	403	$ )
	404	$ /(c002(l)*cc(l,i_o)
	405	$ + c007(l)*cc(l,i_ho2)
	406	$ + c009(l)*cc(l,i_h2o2) ! ajout 20090401
	407	$ + 2.c013(l)cc(l,i_oh) ! ajout 20090401
	408	$ + 2.c017(l)cc(l,i_oh) ! ajout 20090401
	409	$ + e001(l)*cc(l,i_co)
	410	$ + h002(l)*hetero_ice)
[1125]	411
[334]	412	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	413	c h
	414	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	415
[334]	416	cc(l,i_h) = cc(l,i_hox)
	417	$ /(1. + (1. + roh_ho2)/rh_ho2)
[1125]	418
[334]	419	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	420	c ho2
	421	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	422
[334]	423	cc(l,i_ho2) = cc(l,i_h)/rh_ho2
[1125]	424
[334]	425	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	426	c oh
	427	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	428
[334]	429	cc(l,i_oh) = cc(l,i_ho2)*roh_ho2
[1125]	430
[334]	431	end do
[1125]	432
[334]	433	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	434	c ox species
	435	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	436	c day:
	437	c - o1d at photochemical equilibrium
	438	c - o3 at photochemical equilibrium
	439	c - continuity equation for ox
	440	c night:
	441	c - o1d = 0
	442	c - continuity equation for o3
	443	c - continuity equation for o
	444	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	445
[334]	446	if (sza .le. 95.) then
[1125]	447
[334]	448	do l = 1,lswitch-1
[1125]	449
[334]	450	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	451	c o(1d)
	452	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	453
[334]	454	cc(l,i_o1d) = (j(l,j_co2_o1d)*cc(l,i_co2)
	455	$ + j(l,j_o2_o1d)*cc(l,i_o2)
	456	$ + j(l,j_o3_o1d)*cc(l,i_o3))
	457	$ /(b001(l)*cc(l,i_co2)
	458	$ + b002(l)*cc(l,i_h2o)
	459	$ + b003(l)*cc(l,i_h2)
	460	$ + b004(l)*cc(l,i_o2)
	461	$ + b005(l)*cc(l,i_o3)
	462	$ + b006(l)*cc(l,i_o3))
[1125]	463
[334]	464	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	465	c o/o3
	466	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	467
[334]	468	ro_o3 = (j(l,j_o3_o1d) + j(l,j_o3_o)
	469	$ + a003(l)*cc(l,i_o)
	470	$ + c003(l)*cc(l,i_h)
	471	$ + c014(l)*cc(l,i_oh)
	472	$ + c015(l)*cc(l,i_ho2)
	473	$ )
	474	$ /(a001(l)*cc(l,i_o2))
[1125]	475
[334]	476	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	477	c o3
	478	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	479
[334]	480	cc(l,i_o3) = cc(l,i_ox)/(1. + ro_o3)
[1125]	481
[334]	482	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	483	c o
	484	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	485
[334]	486	cc(l,i_o) = cc(l,i_o3)*ro_o3
[1125]	487
[334]	488	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	489	c ox = o + o3
	490	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	491
[334]	492	production(l,i_ox) =
	493	$ + j(l,j_co2_o)*cc(l,i_co2)
	494	$ + j(l,j_co2_o1d)*cc(l,i_co2)
	495	$ + j(l,j_ho2)*cc(l,i_ho2)
	496	$ + 2.j(l,j_o2_o)cc(l,i_o2)
	497	$ + 2.j(l,j_o2_o1d)cc(l,i_o2)
	498	$ + c006(l)cc(l,i_h)cc(l,i_ho2)
	499	$ + c013(l)cc(l,i_oh)cc(l,i_oh)
	500	$ + d003(l)cc(l,i_ho2)no(l)
[1125]	501
[334]	502	loss(l,i_ox) = 2.a002(l)cc(l,i_o)*cc(l,i_o)
	503	$ + 2.a003(l)cc(l,i_o)*cc(l,i_o3)
	504	$ + c001(l)cc(l,i_ho2)cc(l,i_o)
	505	$ + c002(l)cc(l,i_oh)cc(l,i_o)
	506	$ + c003(l)cc(l,i_h)cc(l,i_o3)
	507	$ + c012(l)cc(l,i_o)cc(l,i_h2o2)
	508	$ + c014(l)cc(l,i_o3)cc(l,i_oh)
	509	$ + c015(l)cc(l,i_o3)cc(l,i_ho2)
	510	$ + d001(l)cc(l,i_o)no2(l)
	511	$ + e002(l)cc(l,i_o)cc(l,i_co)
[1125]	512
[334]	513	loss(l,i_ox) = loss(l,i_ox)/cc(l,i_ox)
[1125]	514
[334]	515	end do
[1125]	516
[334]	517	else
[1125]	518
[334]	519	do l = 1,lswitch-1
[1125]	520
[334]	521	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	522	c o(1d)
	523	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	524
[334]	525	cc(l,i_o1d) = 0.
[1125]	526
[334]	527	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	528	c o3
	529	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	530
[334]	531	production(l,i_o3) = a001(l)cc(l,i_o2)cc(l,i_o)
[1125]	532
[334]	533	loss(l,i_o3) = a003(l)*cc(l,i_o)
	534	$ + c003(l)*cc(l,i_h)
	535	$ + c014(l)*cc(l,i_oh)
	536	$ + c015(l)*cc(l,i_ho2)
	537
	538	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	539	c o
	540	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	541
[334]	542	production(l,i_o) = c006(l)cc(l,i_h)cc(l,i_ho2)
	543	$ + c013(l)cc(l,i_oh)cc(l,i_oh)
[1125]	544
[334]	545	loss(l,i_o) = a001(l)*cc(l,i_o2)
	546	$ + 2.a002(l)cc(l,i_o)
	547	$ + a003(l)*cc(l,i_o3)
	548	$ + c001(l)*cc(l,i_ho2)
	549	$ + c002(l)*cc(l,i_oh)
	550	$ + c012(l)*cc(l,i_h2o2)
	551	$ + e002(l)*cc(l,i_co)
[1125]	552
[334]	553	end do
	554	end if
[1125]	555
[334]	556	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	557	c other species
	558	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	559
[334]	560	do l = 1,lswitch-1
[1125]	561
[334]	562	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	563	c co2
	564	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	565
[334]	566	production(l,i_co2) = e001(l)cc(l,i_oh)cc(l,i_co)
	567	$ + e002(l)cc(l,i_o)cc(l,i_co)
	568	$ + t002(l)cc(l,i_ch4)16./44. ! ajout 20090401
	569	$ + t003(l)cc(l,i_co2)8./44. ! ajout 20090401
[1125]	570
[334]	571	loss(l,i_co2) = j(l,j_co2_o)
	572	$ + j(l,j_co2_o1d)
	573	$ + t003(l)
[1125]	574
[334]	575	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	576	c co
	577	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	578
[334]	579	production(l,i_co) = j(l,j_co2_o)*cc(l,i_co2)
	580	$ + j(l,j_co2_o1d)*cc(l,i_co2)
	581	$ + t003(l)*cc(l,i_co2)
[1125]	582
[334]	583	loss(l,i_co) = e001(l)*cc(l,i_oh)
	584	$ + e002(l)*cc(l,i_o)
[1125]	585
[334]	586	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	587	c o2
	588	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	589
[334]	590	production(l,i_o2) =
	591	$ j(l,j_o3_o)*cc(l,i_o3)
	592	$ + j(l,j_o3_o1d)*cc(l,i_o3)
	593	$ + a002(l)cc(l,i_o)cc(l,i_o)
	594	$ + 2.a003(l)cc(l,i_o)*cc(l,i_o3)
	595	$ + 2.b005(l)cc(l,i_o1d)*cc(l,i_o3)
	596	$ + b006(l)cc(l,i_o1d)cc(l,i_o3)
	597	$ + c001(l)cc(l,i_o)cc(l,i_ho2)
	598	$ + c002(l)cc(l,i_o)cc(l,i_oh)
	599	$ + c003(l)cc(l,i_h)cc(l,i_o3)
	600	$ + c005(l)cc(l,i_h)cc(l,i_ho2)
	601	$ + c007(l)cc(l,i_oh)cc(l,i_ho2)
	602	$ + c008(l)cc(l,i_ho2)cc(l,i_ho2)
	603	$ + c014(l)cc(l,i_o3)cc(l,i_oh)
	604	$ + 2.c015(l)cc(l,i_o3)*cc(l,i_ho2)
	605	$ + c016(l)cc(l,i_ho2)cc(l,i_ho2)
	606	$ + d001(l)cc(l,i_o)no2(l)
[1125]	607
[334]	608	loss(l,i_o2) = j(l,j_o2_o)
	609	$ + j(l,j_o2_o1d)
	610	$ + a001(l)*cc(l,i_o)
	611	$ + c011(l)*cc(l,i_h)
[1125]	612
[334]	613	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	614	c h2
	615	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	616
[334]	617	production(l,i_h2) = c005(l)cc(l,i_h)cc(l,i_ho2)
	618	$ + c018(l)cc(l,i_h)cc(l,i_h)
[1125]	619
[334]	620	loss(l,i_h2) = b003(l)*cc(l,i_o1d)
	621	$ + c010(l)*cc(l,i_oh)
[1125]	622
[334]	623	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	624	c h2o
	625	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	626
[334]	627	production(l,i_h2o) =
	628	$ c006(l)cc(l,i_h)cc(l,i_ho2)
	629	$ + c007(l)cc(l,i_oh)cc(l,i_ho2)
	630	$ + c009(l)cc(l,i_oh)cc(l,i_h2o2)
	631	$ + c010(l)cc(l,i_oh)cc(l,i_h2)
	632	$ + c013(l)cc(l,i_oh)cc(l,i_oh)
	633	$ + h004(l)cc(l,i_h2o2)hetero_ice
[1125]	634
[334]	635	loss(l,i_h2o) = j(l,j_h2o)
	636	$ + b002(l)*cc(l,i_o1d)
	637	$ + t001(l)
[1125]	638
[334]	639	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	640	c h2o2
	641	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	642
[334]	643	production(l,i_h2o2) =
	644	$ c008(l)cc(l,i_ho2)cc(l,i_ho2)
	645	$ + c016(l)cc(l,i_ho2)cc(l,i_ho2)
	646	$ + c017(l)cc(l,i_oh)cc(l,i_oh)
	647	c $ + 0.5h001(l)cc(l,i_ho2)*hetero_ice
	648	c $ + 0.5h002(l)cc(l,i_oh)*hetero_ice
[1125]	649
[334]	650	loss(l,i_h2o2) = j(l,j_h2o2)
	651	$ + c009(l)*cc(l,i_oh)
	652	$ + c012(l)*cc(l,i_o)
	653	$ + h004(l)*hetero_ice
	654	$ + h005(l)*hetero_dust
[1125]	655
[334]	656	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	657	c hox = h + oh + ho2
	658	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	659
[334]	660	production(l,i_hox) =
	661	$ 2.j(l,j_h2o)cc(l,i_h2o)
	662	$ + 2.j(l,j_h2o2)cc(l,i_h2o2)
	663	$ + 2.b002(l)cc(l,i_o1d)*cc(l,i_h2o)
	664	$ + 2.b003(l)cc(l,i_o1d)*cc(l,i_h2)
	665	$ + 2.c012(l)cc(l,i_o)*cc(l,i_h2o2)
	666	$ + 2.t001(l)cc(l,i_h2o)
[1125]	667
[334]	668	loss(l,i_hox) = 2.c005(l)cc(l,i_h)*cc(l,i_ho2)
	669	$ + 2.c006(l)cc(l,i_h)*cc(l,i_ho2)
	670	$ + 2.c007(l)cc(l,i_oh)*cc(l,i_ho2)
	671	$ + 2.c008(l)cc(l,i_ho2)*cc(l,i_ho2)
	672	$ + 2.c013(l)cc(l,i_oh)*cc(l,i_oh)
	673	$ + 2.c016(l)cc(l,i_ho2)*cc(l,i_ho2)
	674	$ + 2.c017(l)cc(l,i_oh)*cc(l,i_oh)
	675	$ + 2.c018(l)cc(l,i_h)*cc(l,i_h)
	676	$ + h001(l)cc(l,i_ho2)hetero_ice
	677	$ + h002(l)cc(l,i_oh)hetero_ice
	678	$ + h003(l)cc(l,i_ho2)hetero_dust
[1125]	679
[334]	680	loss(l,i_hox) = loss(l,i_hox)/cc(l,i_hox)
[1125]	681
[334]	682	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	683	c ch4
	684	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	685
[334]	686	production(l,i_ch4) = 0.
[1125]	687
[334]	688	loss(l,i_ch4) = j(l,j_ch4_ch3_h)
	689	$ + j(l,j_ch4_1ch2_h2)
	690	$ + j(l,j_ch4_3ch2_h_h)
	691	$ + j(l,j_ch4_ch_h2_h)
	692	$ + b007(l)*cc(l,i_o1d)
	693	$ + b008(l)*cc(l,i_o1d)
	694	$ + b009(l)*cc(l,i_o1d)
	695	$ + e003(l)*cc(l,i_oh)
[1125]	696
[334]	697	end do
[1125]	698
[334]	699	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	700	c update number densities
	701	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	702
[334]	703	c long-lived species
[1125]	704
[334]	705	do l = 1,lswitch-1
	706	cc(l,i_co2) = (cc0(l,i_co2) + production(l,i_co2)*dt)
	707	$ /(1. + loss(l,i_co2)*dt)
	708	cc(l,i_co) = (cc0(l,i_co) + production(l,i_co)*dt)
	709	$ /(1. + loss(l,i_co)*dt)
	710	cc(l,i_o2) = (cc0(l,i_o2) + production(l,i_o2)*dt)
	711	$ /(1. + loss(l,i_o2)*dt)
	712	cc(l,i_h2) = (cc0(l,i_h2) + production(l,i_h2)*dt)
	713	$ /(1. + loss(l,i_h2)*dt)
	714	cc(l,i_h2o2)= (cc0(l,i_h2o2) + production(l,i_h2o2)*dt)
	715	$ /(1. + loss(l,i_h2o2)*dt)
	716	cc(l,i_h2o) = (cc0(l,i_h2o) + production(l,i_h2o)*dt)
	717	$ /(1. + loss(l,i_h2o)*dt)
	718	cc(l,i_hox) = (cc0(l,i_hox) + production(l,i_hox)*dt)
	719	$ /(1. + loss(l,i_hox)*dt)
	720	cc(l,i_ch4) = (cc0(l,i_ch4) + production(l,i_ch4)*dt)
	721	$ /(1. + loss(l,i_ch4)*dt)
	722	end do
[1125]	723
[334]	724	c ox species
[1125]	725
[334]	726	if (sza .le. 95.) then
	727	do l = 1,lswitch-1
	728	cc(l,i_ox) = (cc0(l,i_ox) + production(l,i_ox)*dt)
	729	$ /(1. + loss(l,i_ox)*dt)
	730	end do
	731	else
	732	do l = 1,lswitch-1
	733	cc(l,i_o) = (cc0(l,i_o) + production(l,i_o)*dt)
	734	$ /(1. + loss(l,i_o)*dt)
	735	cc(l,i_o3) = (cc0(l,i_o3) + production(l,i_o3)*dt)
	736	$ /(1. + loss(l,i_o3)*dt)
	737	cc(l,i_ox) = cc(l,i_o) + cc(l,i_o3)
	738	end do
	739	end if
[1125]	740
[334]	741	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	742	c end of loop over iterations
	743	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	744
[334]	745	end do
[1125]	746
[334]	747	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	748	c density -> volume mixing ratio conversion
	749	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	750
[334]	751	do iesp = 1,nesp
	752	do l = 1,lswitch-1
	753	rm(l,iesp) = max(cc(l,iesp)/dens(l), 1.e-30)
	754	end do
	755	end do
[1125]	756
[334]	757	return
	758	end
[1125]	759
[334]	760	c*****************************************************************
[1125]	761
	762	subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm)
	763
[334]	764	c*****************************************************************
[1125]	765
	766	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
[1036]	767	& igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
	768	& igcm_ho2, igcm_h2o2, igcm_n2, igcm_h2o_vap,
	769	& igcm_ch4
[1125]	770
[334]	771	implicit none
[1125]	772
[334]	773	#include "callkeys.h"
[1125]	774
[334]	775	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	776	c inputs:
	777	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	778
	779	integer, intent(in) :: nlayer ! number of atmospheric layers
	780	integer, intent(in) :: nq ! number of tracers
	781	real :: zycol(nlayer,nq) ! species volume mixing ratio in the gcm
	782
	783	integer :: nesp ! number of species in the chemistry
	784	integer :: lswitch ! interface level between chemistries
	785
[334]	786	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	787	c outputs:
	788	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	789
	790	real :: rm(nlayer,nesp) ! species volume mixing ratio
	791
[334]	792	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	793	c local:
	794	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	795
	796	integer :: l, iq
[334]	797
	798	c tracer indexes in the chemistry:
	799
	800	integer,parameter :: i_co2 = 1
	801	integer,parameter :: i_co = 2
	802	integer,parameter :: i_o = 3
	803	integer,parameter :: i_o1d = 4
	804	integer,parameter :: i_o2 = 5
	805	integer,parameter :: i_o3 = 6
	806	integer,parameter :: i_h = 7
	807	integer,parameter :: i_h2 = 8
	808	integer,parameter :: i_oh = 9
	809	integer,parameter :: i_ho2 = 10
	810	integer,parameter :: i_h2o2 = 11
	811	integer,parameter :: i_ch4 = 12
	812	integer,parameter :: i_h2o = 13
	813	integer,parameter :: i_n2 = 14
	814	integer,parameter :: i_hox = 15
	815	integer,parameter :: i_ox = 16
[1125]	816
[334]	817	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	818	c initialise chemical species
	819	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	820
[334]	821	do l = 1,lswitch-1
	822	rm(l,i_co2) = max(zycol(l, igcm_co2), 1.e-30)
	823	rm(l,i_co) = max(zycol(l, igcm_co), 1.e-30)
	824	rm(l,i_o) = max(zycol(l, igcm_o), 1.e-30)
	825	rm(l,i_o1d) = max(zycol(l, igcm_o1d), 1.e-30)
	826	rm(l,i_o2) = max(zycol(l, igcm_o2), 1.e-30)
	827	rm(l,i_o3) = max(zycol(l, igcm_o3), 1.e-30)
	828	rm(l,i_h) = max(zycol(l, igcm_h), 1.e-30)
	829	rm(l,i_h2) = max(zycol(l, igcm_h2), 1.e-30)
	830	rm(l,i_oh) = max(zycol(l, igcm_oh), 1.e-30)
	831	rm(l,i_ho2) = max(zycol(l, igcm_ho2), 1.e-30)
	832	rm(l,i_h2o2) = max(zycol(l, igcm_h2o2), 1.e-30)
	833	rm(l,i_n2) = max(zycol(l, igcm_n2), 1.e-30)
	834	rm(l,i_h2o) = max(zycol(l, igcm_h2o_vap), 1.e-30)
	835	end do
[618]	836
	837	if (igcm_ch4 .eq. 0) then
	838	do l = 1,lswitch-1
	839	rm(l,i_ch4) = 0.
	840	end do
	841	else
	842	do l = 1,lswitch-1
	843	rm(l,i_ch4) = max(zycol(l,igcm_ch4), 1.e-30)
	844	end do
	845	end if
[1125]	846
[334]	847	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	848	c initialise chemical families c
	849	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	850
[334]	851	do l = 1,lswitch-1
	852	rm(l,i_hox) = rm(l,i_h)
	853	$ + rm(l,i_oh)
	854	$ + rm(l,i_ho2)
	855	rm(l,i_ox) = rm(l,i_o)
	856	$ + rm(l,i_o3)
	857	end do
[1125]	858
[334]	859	return
	860	end
[1125]	861
[334]	862	c*****************************************************************
	863	c
[1125]	864	subroutine chimtogcm(nlayer, nq, zycol, lswitch, nesp, rm)
[334]	865	c
	866	c*****************************************************************
	867	c
[1125]	868	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
[1036]	869	& igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
	870	& igcm_ho2, igcm_h2o2, igcm_n2, igcm_h2o_vap,
	871	& igcm_ch4
[1125]	872
[334]	873	implicit none
[1125]	874
[334]	875	#include "callkeys.h"
[1125]	876
[334]	877	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	878	c inputs:
	879	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	880
	881	integer, intent(in) :: nlayer ! number of atmospheric layers
	882	integer, intent(in) :: nq ! number of tracers
	883	integer :: nesp ! number of species in the chemistry
	884	integer :: lswitch ! interface level between chemistries
	885
	886	real :: rm(nlayer,nesp) ! species volume mixing ratio
	887
[334]	888	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	889	c output:
	890	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	891
	892	real :: zycol(nlayer,nq) ! species volume mixing ratio in the gcm
	893
[334]	894	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	895	c local:
	896	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	897
	898	integer l, iq
[334]	899
	900	c tracer indexes in the chemistry:
	901
	902	integer,parameter :: i_co2 = 1
	903	integer,parameter :: i_co = 2
	904	integer,parameter :: i_o = 3
	905	integer,parameter :: i_o1d = 4
	906	integer,parameter :: i_o2 = 5
	907	integer,parameter :: i_o3 = 6
	908	integer,parameter :: i_h = 7
	909	integer,parameter :: i_h2 = 8
	910	integer,parameter :: i_oh = 9
	911	integer,parameter :: i_ho2 = 10
	912	integer,parameter :: i_h2o2 = 11
	913	integer,parameter :: i_ch4 = 12
	914	integer,parameter :: i_h2o = 13
	915	integer,parameter :: i_n2 = 14
	916	integer,parameter :: i_hox = 15
	917	integer,parameter :: i_ox = 16
[1125]	918
[334]	919	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	920	c save mixing ratios for the gcm
	921	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	922
[334]	923	do l = 1,lswitch-1
	924	zycol(l, igcm_co2) = rm(l,i_co2)
	925	zycol(l, igcm_co) = rm(l,i_co)
	926	zycol(l, igcm_o) = rm(l,i_o)
	927	zycol(l, igcm_o1d) = rm(l,i_o1d)
	928	zycol(l, igcm_o2) = rm(l,i_o2)
	929	zycol(l, igcm_o3) = rm(l,i_o3)
	930	zycol(l, igcm_h) = rm(l,i_h)
	931	zycol(l, igcm_h2) = rm(l,i_h2)
	932	zycol(l, igcm_oh) = rm(l,i_oh)
	933	zycol(l, igcm_ho2) = rm(l,i_ho2)
	934	zycol(l, igcm_h2o2) = rm(l,i_h2o2)
	935	zycol(l, igcm_n2) = rm(l,i_n2)
	936	zycol(l, igcm_h2o_vap) = rm(l,i_h2o)
	937	end do
[618]	938
	939	if (igcm_ch4 .ne. 0) then
	940	do l = 1,lswitch-1
	941	zycol(l,igcm_ch4) = rm(l,i_ch4)
	942	end do
	943	end if
[1125]	944
[334]	945	return
	946	end
[1125]	947
[334]	948	c*****************************************************************
[1125]	949
	950	subroutine chemrates(nlayer,
	951	$ lswitch, dens, press, t,
[334]	952	$ surfdust1d, surfice1d,
	953	$ a001, a002, a003,
	954	$ b001, b002, b003, b004, b005, b006,
	955	$ b007, b008, b009,
	956	$ c001, c002, c003, c004, c005, c006,
	957	$ c007, c008, c009, c010, c011, c012,
	958	$ c013, c014, c015, c016, c017, c018,
	959	$ d001, d002, d003,
	960	$ e001, e002, e003,
	961	$ h001, h002, h003, h004, h005,
	962	$ t001, t002, t003, tau)
[1125]	963
[334]	964	c*****************************************************************
[1125]	965
[1226]	966	USE comcstfi_h
[334]	967	implicit none
[1125]	968
[334]	969	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	970	c inputs: c
	971	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	972
[1125]	973	integer, intent(in) :: nlayer ! number of atmospheric layers
	974	integer :: lswitch ! interface level between chemistries
	975
	976	real :: dens(nlayer) ! density (cm-3)
	977	real :: press(nlayer) ! pressure (hpa)
	978	real :: t(nlayer) ! temperature (k)
	979	real :: surfdust1d(nlayer) ! dust surface area (cm^2/cm^3)
	980	real :: surfice1d(nlayer) ! ice surface area (cm^2/cm^3)
	981	real :: tau ! dust opacity at 7 hpa
	982
	983	real, parameter :: tribo = 0. ! switch for triboelectricity
	984
[334]	985	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	986	c outputs: c
	987	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	988
	989	real :: a001(nlayer), a002(nlayer), a003(nlayer)
	990	real :: b001(nlayer), b002(nlayer), b003(nlayer),
	991	$ b004(nlayer), b005(nlayer), b006(nlayer),
	992	$ b007(nlayer), b008(nlayer), b009(nlayer)
	993	real :: c001(nlayer), c002(nlayer), c003(nlayer),
	994	$ c004(nlayer), c005(nlayer), c006(nlayer),
	995	$ c007(nlayer), c008(nlayer), c009(nlayer),
	996	$ c010(nlayer), c011(nlayer), c012(nlayer),
	997	$ c013(nlayer), c014(nlayer), c015(nlayer),
	998	$ c016(nlayer), c017(nlayer), c018(nlayer)
	999	real :: d001(nlayer), d002(nlayer), d003(nlayer)
	1000	real :: e001(nlayer), e002(nlayer), e003(nlayer)
	1001	real :: h001(nlayer), h002(nlayer), h003(nlayer),
	1002	$ h004(nlayer), h005(nlayer)
	1003	real :: t001(nlayer), t002(nlayer), t003(nlayer)
	1004
[334]	1005	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1006	c local: c
	1007	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	1008
	1009	real :: ak0, ak1, rate, rate1, rate2, xpo, xpo1, xpo2
	1010	real :: ef, efmax, lossh2o, lossch4, lossco2
	1011
	1012	integer :: l
	1013	real :: k1a, k1b, k1a0, k1b0, k1ainf
	1014	real :: x, y, fc, fx
	1015
[334]	1016	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1017	c compute reaction rates
	1018	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[1125]	1019
[334]	1020	do l = 1,lswitch-1
[1125]	1021
[334]	1022	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1023	c oxygen compounds
	1024	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1025	c
	1026	ccc a001: o + o2 + co2 -> o3 + co2
	1027	c
	1028	c jpl 2003
	1029	c
[576]	1030	c co2 efficiency as a third body (2.075)
	1031	c from sehested et al., j. geophys. res., 100, 1995.
	1032	c
	1033	a001(l) = 2.075
[334]	1034	$ 6.0e-34(t(l)/300.)*(-2.4)dens(l)
	1035	c
	1036	c mulcahy and williams, 1968
	1037	c
	1038	c a001(l) = 2.68e-33(t(l)/298.)(-2.4)dens(l)
	1039	c
	1040	c nair et al., 1994
	1041	c
	1042	c a001(l) = 1.3e-34exp(724./t(l))dens(l)
	1043	c
	1044	ccc a002: o + o + co2 -> o2 + co2
	1045	c
	1046	c Tsang and Hampson, J. Chem. Phys. Ref. Data, 15, 1087, 1986
	1047	c
	1048	c a002(l) = 2.55.2e-35exp(900./t(l))*dens(l)
	1049	c
	1050	c Campbell and Gray, Chem. Phys. Lett., 18, 607, 1973
	1051	c
	1052	c a002(l) = 1.2e-32(300./t(l))(2.0)dens(l) ! yung expression
	1053	c
	1054	a002(l) = 2.59.46e-34exp(485./t(l))*dens(l) ! nist expression
	1055	c
	1056	c baulch et al., 1976 confirmed by smith and robertson, 2008
	1057	c
	1058	c a002(l) = 2.52.76e-34exp(720./t(l))*dens(l)
	1059	c
	1060	ccc a003: o + o3 -> o2 + o2
	1061	c
	1062	c jpl 2003
	1063	c
	1064	a003(l) = 8.0e-12*exp(-2060./t(l))
	1065	c
	1066	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1067	c reactions with o(1d)
	1068	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1069	c
	1070	ccc b001: o(1d) + co2 -> o + co2
	1071	c
	1072	c jpl 2003
	1073	c
	1074	c b001(l) = 7.4e-11*exp(120./t(l))
	1075	c
	1076	c jpl 2006
	1077	c
	1078	b001(l) = 7.5e-11*exp(115./t(l))
	1079	c
	1080	ccc b002: o(1d) + h2o -> oh + oh
	1081	c
	1082	c jpl 2003
	1083	c
	1084	c b002(l) = 2.2e-10
	1085	c
	1086	c jpl 2006
	1087	c
	1088	b002(l) = 1.63e-10*exp(60./t(l))
	1089	c
	1090	ccc b003: o(1d) + h2 -> oh + h
	1091	c
[408]	1092	c jpl 2011
[334]	1093	c
[408]	1094	b003(l) = 1.2e-10
[334]	1095	c
	1096	ccc b004: o(1d) + o2 -> o + o2
	1097	c
	1098	c jpl 2003
	1099	c
	1100	c b004(l) = 3.2e-11*exp(70./t(l))
	1101	c
	1102	c jpl 2006
	1103	c
	1104	b004(l) = 3.3e-11*exp(55./t(l))
	1105	c
	1106	ccc b005: o(1d) + o3 -> o2 + o2
	1107	c
	1108	c jpl 2003
	1109	c
	1110	b005(l) = 1.2e-10
	1111	c
	1112	ccc b006: o(1d) + o3 -> o2 + o + o
	1113	c
	1114	c jpl 2003
	1115	c
	1116	b006(l) = 1.2e-10
	1117	c
	1118	ccc b007: o(1d) + ch4 -> ch3 + oh
	1119	c
	1120	c jpl 2003
	1121	c
	1122	b007(l) = 1.5e-10*0.75
	1123	c
	1124	ccc b008: o(1d) + ch4 -> ch3o + h
	1125	c
	1126	c jpl 2003
	1127	c
	1128	b008(l) = 1.5e-10*0.20
	1129	c
	1130	ccc b009: o(1d) + ch4 -> ch2o + h2
	1131	c
	1132	c jpl 2003
	1133	c
	1134	b009(l) = 1.5e-10*0.05
	1135	c
	1136	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1137	c hydrogen compounds
	1138	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1139	c
	1140	ccc c001: o + ho2 -> oh + o2
	1141	c
	1142	c jpl 2003
	1143	c
	1144	c001(l) = 3.0e-11*exp(200./t(l))
	1145	c
	1146	ccc c002: o + oh -> o2 + h
	1147	c
[408]	1148	c jpl 2011
[334]	1149	c
[408]	1150	c002(l) = 1.8e-11*exp(180./t(l))
[334]	1151	c
	1152	c robertson and smith, j. chem. phys. a 110, 6673, 2006
	1153	c
	1154	c c002(l) = 11.2e-11t(l)(-0.32)exp(177./t(l))
	1155	c
	1156	ccc c003: h + o3 -> oh + o2
	1157	c
	1158	c jpl 2003
	1159	c
	1160	c003(l) = 1.4e-10*exp(-470./t(l))
	1161	c
	1162	ccc c004: h + ho2 -> oh + oh
	1163	c
	1164	c jpl 2003
	1165	c
	1166	c c004(l) = 8.1e-11*0.90
	1167	c
	1168	c jpl 2006
	1169	c
	1170	c004(l) = 7.2e-11
	1171	c
	1172	ccc c005: h + ho2 -> h2 + o2
	1173	c
	1174	c jpl 2003
	1175	c
	1176	c c005(l) = 8.1e-11*0.08
	1177	c
	1178	c jpl 2006
	1179	c
	1180	c005(l) = 6.9e-12
	1181	c
	1182	ccc c006: h + ho2 -> h2o + o
	1183	c
	1184	c jpl 2003
	1185	c
	1186	c c006(l) = 8.1e-11*0.02
	1187	c
	1188	c jpl 2006
	1189	c
	1190	c006(l) = 1.6e-12
	1191	c
	1192	ccc c007: oh + ho2 -> h2o + o2
	1193	c
	1194	c jpl 2003
	1195	c
	1196	c007(l) = 4.8e-11*exp(250./t(l))
	1197	c
	1198	c jpl 2003 +20% d'apres canty et al., grl, 2006
	1199	c
	1200	c c007(l) = 4.8e-11exp(250./t(l))1.2
	1201	c
	1202	ccc c008: ho2 + ho2 -> h2o2 + o2
	1203	c
	1204	c jpl 2003
	1205	c
	1206	c c008(l) = 2.3e-13*exp(600./t(l))
	1207	c
	1208	c christensen et al., grl, 13, 2002
	1209	c
	1210	c008(l) = 1.5e-12*exp(19./t(l))
	1211	c
	1212	ccc c009: oh + h2o2 -> h2o + ho2
	1213	c
	1214	c jpl 2003
	1215	c
	1216	c c009(l) = 2.9e-12*exp(-160./t(l))
	1217	c
	1218	c jpl 2006
	1219	c
	1220	c009(l) = 1.8e-12
	1221	c
	1222	ccc c010: oh + h2 -> h2o + h
	1223	c
	1224	c jpl 2003
	1225	c
	1226	c c010(l) = 5.5e-12*exp(-2000./t(l))
	1227	c
	1228	c jpl 2006
	1229	c
	1230	c010(l) = 2.8e-12*exp(-1800./t(l))
	1231	c
	1232	ccc c011: h + o2 + co2 -> ho2 + co2
	1233	c
[408]	1234	c jpl 2011
[334]	1235	c
	1236	ak0 = 2.54.4e-32(t(l)/300.)**(-1.3)
[408]	1237	ak1 = 7.5e-11(t(l)/300.)*(0.2)
[334]	1238	c
	1239	rate = (ak0dens(l))/(1. + ak0dens(l)/ak1)
	1240	xpo = 1./(1. + alog10((ak0dens(l))/ak1)*2)
	1241	c011(l) = rate0.6*xpo
	1242	c
	1243	ccc c012: o + h2o2 -> oh + ho2
	1244	c
	1245	c jpl 2003
	1246	c
	1247	c012(l) = 1.4e-12*exp(-2000./t(l))
	1248	c
	1249	ccc c013: oh + oh -> h2o + o
	1250	c
	1251	c jpl 2003
	1252	c
	1253	c c013(l) = 4.2e-12*exp(-240./t(l))
	1254	c
	1255	c jpl 2006
	1256	c
	1257	c013(l) = 1.8e-12
	1258	c
	1259	ccc c014: oh + o3 -> ho2 + o2
	1260	c
	1261	c jpl 2003
	1262	c
	1263	c014(l) = 1.7e-12*exp(-940./t(l))
	1264	c
	1265	c jpl 2000
	1266	c
	1267	c c014(l) = 1.5e-12*exp(-880./t(l))
	1268	c
	1269	c nair et al., 1994 (jpl 1997)
	1270	c
	1271	c c014(l) = 1.6e-12*exp(-940./t(l))
	1272	c
	1273	ccc c015: ho2 + o3 -> oh + o2 + o2
	1274	c
	1275	c jpl 2003
	1276	c
	1277	c015(l) = 1.0e-14*exp(-490./t(l))
	1278	c
	1279	c jpl 2000
	1280	c
	1281	c c015(l) = 2.0e-14*exp(-680./t(l))
	1282	c
	1283	c nair et al., 1994 (jpl 1997)
	1284	c
	1285	c c015(l) = 1.1e-14*exp(-500./t(l))
	1286	c
	1287	ccc c016: ho2 + ho2 + co2 -> h2o2 + o2 + co2
	1288	c
[408]	1289	c jpl 2011
[334]	1290	c
[408]	1291	c016(l) = 2.5*2.1e-33
	1292	$ exp(920./t(l))dens(l)
[334]	1293	c
	1294	ccc c017: oh + oh + co2 -> h2o2 + co2
	1295	c
	1296	c jpl 2003
	1297	c
	1298	ak0 = 2.56.9e-31(t(l)/300.)**(-1.0)
	1299	ak1 = 2.6e-11(t(l)/300.)*(0.0)
	1300	c
	1301	c jpl 1997
	1302	c
	1303	c ak0 = 2.56.2e-31(t(l)/300.)**(-1.0)
	1304	c ak1 = 2.6e-11(t(l)/300.)*(0.0)
	1305	c
	1306	c nair et al., 1994
	1307	c
	1308	c ak0 = 2.57.1e-31(t(l)/300.)**(-0.8)
	1309	c ak1 = 1.5e-11(t(l)/300.)*(0.0)
	1310	c
	1311	rate = (ak0dens(l))/(1. + ak0dens(l)/ak1)
	1312	xpo = 1./(1. + alog10((ak0dens(l))/ak1)*2)
	1313	c017(l) = rate0.6*xpo
	1314	c
	1315	ccc c018: h + h + co2 -> h2 + co2
	1316	c
	1317	c baulch et al., 1992
	1318	c
	1319	c c018(l) = 2.58.85e-33(t(l)/298.)*(-0.6)dens(l)
	1320	c
	1321	c baulch et al., 2005
	1322	c
	1323	c018(l) = 2.51.8e-30(t(l)*(-1.0))dens(l)
	1324	c
	1325	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1326	c nitrogen compounds
	1327	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1328	c
	1329	ccc d001: no2 + o -> no + o2
	1330	c
	1331	c jpl 2003
	1332	c
	1333	c d001(l) = 5.6e-12*exp(180./t(l))
	1334	c
	1335	ccc jpl 2006
	1336	c
	1337	d001(l) = 5.1e-12*exp(210./t(l))
	1338	c
	1339	ccc d002: no + o3 -> no2 + o2
	1340	c
	1341	c jpl 2003
	1342	c
	1343	d002(l) = 3.0e-12*exp(-1500./t(l))
	1344	c
	1345	ccc d003: no + ho2 -> no2 + oh
	1346	c
[408]	1347	c jpl 2011
[334]	1348	c
[408]	1349	d003(l) = 3.3e-12*exp(270./t(l))
[334]	1350	c
	1351	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1352	c carbon compounds
	1353	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1354	c
	1355	ccc e001: oh + co -> co2 + h
	1356	c
	1357	c jpl 2003
	1358	c
	1359	c e001(l) = 1.5e-13(1 + 0.6press(l)/1013.)
	1360	c
	1361	c mccabe et al., grl, 28, 3135, 2001
	1362	c
	1363	c e001(l) = 1.57e-13 + 3.54e-33*dens(l)
	1364	c
	1365	c atkinson et al. 2006
	1366	c
	1367	c e001(l) = 1.44e-13 + 3.43e-33*dens(l)
	1368	c
	1369	c joshi et al., 2006
	1370	c
	1371	k1a0 = 1.342.5dens(l)
	1372	$ 1/(1/(3.62e-26t(l)*(-2.739)exp(-20./t(l)))
	1373	$ + 1/(6.48e-33t(l)(0.14)exp(-57./t(l)))) ! corrige de l'erreur publi
	1374	k1b0 = 1.17e-19t(l)(2.053)exp(139./t(l))
	1375	$ + 9.56e-12t(l)(-0.664)exp(-167./t(l))
	1376	k1ainf = 1.52e-17t(l)(1.858)exp(28.8/t(l))
	1377	$ + 4.78e-8t(l)(-1.851)exp(-318./t(l))
	1378	x = k1a0/(k1ainf - k1b0)
	1379	y = k1b0/(k1ainf - k1b0)
	1380	fc = 0.628exp(-1223./t(l)) + (1. - 0.628)exp(-39./t(l))
	1381	$ + exp(-t(l)/255.)
	1382	fx = fc(1./(1. + (alog(x))2)) ! corrige de l'erreur publi
	1383	k1a = k1a0((1. + y)/(1. + x))fx
	1384	k1b = k1b0(1./(1.+x))fx
	1385	c
	1386	e001(l) = k1a + k1b
	1387	c
	1388	ccc e002: o + co + m -> co2 + m
	1389	c
	1390	c tsang and hampson, 1986.
	1391	c
	1392	e002(l) = 2.56.5e-33exp(-2184./t(l))*dens(l)
	1393	c
	1394	c baulch et al., butterworths, 1976.
	1395	c
	1396	c e002(l) = 1.6e-32exp(-2184./t(l))dens(l)
	1397	c
	1398	ccc e003: ch4 + oh -> ch3 + h2o
	1399	c
	1400	c jpl 2003
	1401	c
	1402	c e003(l) = 2.45e-12*exp(-1775./t(l))
	1403	c
	1404	c jpl 2003, three-parameter expression
	1405	c
	1406	e003(l) = 2.80e-14(t(l)0.667)exp(-1575./t(l))
	1407	c
	1408	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1409	c heterogenous chemistry
	1410	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1411	c
	1412	c k = (surfacevgamma)/4 (s-1)
	1413	c v = 100sqrt(8rt/(pim)) (cm s-1)
	1414	c
	1415	ccc h001: ho2 + ice -> products
	1416	c
	1417	c cooper and abbatt, 1996: gamma = 0.025
	1418	c
	1419	h001(l) = surfice1d(l)
	1420	$ 100.sqrt(8.8.31t(l)/(33.e-3pi))0.025/4.
	1421	c
	1422	c h002: oh + ice -> products
	1423	c
	1424	c cooper and abbatt, 1996: gamma = 0.03
	1425	c
	1426	h002(l) = surfice1d(l)
	1427	$ 100.sqrt(8.8.31t(l)/(17.e-3pi))0.03/4.
	1428	c
	1429	c h003: ho2 + dust -> products
	1430	c
	1431	c jacob, 2000: gamma = 0.2
	1432	c see dereus et al., atm. chem. phys., 2005
	1433	c
	1434	c h003(l) = surfdust1d(l)
	1435	c $ 100.sqrt(8.8.31t(l)/(33.e-3pi))0.2/4.
	1436	h003(l) = 0. ! advised
	1437	c
	1438	ccc h004: h2o2 + ice -> products
	1439	c
	1440	c gamma = 1.e-3 test value
	1441	c
	1442	c h004(l) = surfice1d(l)
	1443	c $ 100.sqrt(8.8.31t(l)/(34.e-3pi))0.001/4.
	1444	h004(l) = 0. ! advised
	1445	c
	1446	c h005: h2o2 + dust -> products
	1447	c
	1448	c gamma = 5.e-4
	1449	c see dereus et al., atm. chem. phys., 2005
	1450	c
	1451	h005(l) = surfdust1d(l)
	1452	$ 100.sqrt(8.8.31t(l)/(34.e-3pi))5.e-4/4.
	1453	h005(l) = 0. ! advised
	1454	c
	1455	end do
	1456	c
	1457	if (tribo .eq. 1.) then
	1458	c
	1459	c electrochemical reactions
	1460	c
	1461	c efmax: maximum electric field (kv.m-1)
	1462	c
	1463	efmax = 23.3
	1464	c
	1465	c ef: actual electric field, scaled by tau.
	1466	c
	1467	c if (tau .ge. 1.) then
	1468	c ef = efmax
	1469	c else
	1470	c ef = 0.
	1471	c end if
	1472	c ef = min(efmax,efmax*tau/1.0)
	1473	c
	1474	ef = (efmax/0.5)tau - (efmax/0.5)0.5
	1475	c
	1476	ef = max(ef, 0.)
	1477	ef = min(ef, efmax)
	1478	c
	1479	ccc t001: h2o + e -> oh + h-
	1480	c
	1481	c lossh2o: fit of oh/h- production rates
	1482	c given by delory et al., astrobiology, 6, 451, 2006
	1483	c
	1484	if (ef .eq. 0.) then
	1485	lossh2o = 0.
	1486	else if (ef .lt. 10.) then
	1487	lossh2o = 0.054136exp(1.0978ef)
	1488	else if (ef .lt. 16.) then
	1489	lossh2o = 64.85exp(0.38894ef)
	1490	else if (ef .le. 20.) then
	1491	lossh2o = 0.2466exp(0.73719ef)
	1492	else
	1493	lossh2o = 2.3269e-8exp(1.546ef)
	1494	end if
	1495	c
	1496	c production rates are given for h2o = 20 prec. microns.
	1497	c t001 is converted to first-order reaction rate
	1498	c assuming h2o number density at the surface = 5e13 mol cm-3
	1499	c
	1500	do l = 1,21 ! 70 km
	1501	t001(l) = lossh2o/5.e13 ! s-1
	1502	end do
	1503	do l = 22,lswitch-1
	1504	t001(l) = 0.
	1505	end do
	1506	c
	1507	ccc t002: ch4 + e -> products
	1508	c
	1509	c lossch4: fit of ch4 loss rates
	1510	c given by farrell et al., grl, 33, 2006
	1511	c
	1512	if (ef .eq. 0.) then
	1513	lossch4 = 0.
	1514	else if (ef .gt. 20.) then
	1515	lossch4 = 1.113e-21exp(1.6065ef)
	1516	else if (ef .gt. 17.5) then
	1517	lossch4 = 1.e-15exp(0.92103ef)
	1518	else if (ef .gt. 14.) then
	1519	lossch4 = 1.e-13exp(0.65788ef)
	1520	else
	1521	lossch4 = 8.9238e-15exp(0.835ef)
	1522	end if
	1523	c
	1524	do l = 1,21 ! 70 km
	1525	t002(l) = lossch4 ! s-1
	1526	end do
	1527	do l = 22,lswitch-1
	1528	t002(l) = 0.
	1529	end do
	1530	c
	1531	ccc t003: co2 + e -> co + o-
	1532	c
	1533	c lossco2: fit of co/o- production rates
	1534	c given by delory et al., astrobiology, 6, 451, 2006
	1535	c
	1536	if (ef .eq. 0.) then
	1537	lossco2 = 0.
	1538	else if (ef .lt. 10.) then
	1539	lossco2 = 22.437exp(1.045ef)
	1540	else if (ef .lt. 16.) then
	1541	lossco2 = 17518.exp(0.37896ef)
	1542	else if (ef .lt. 20.) then
	1543	lossco2 = 54.765exp(0.73946ef)
	1544	else
	1545	lossco2 = 4.911e-6exp(1.5508ef)
	1546	end if
	1547	c
	1548	c production rates are assumed to be given for p = 6 hPa
	1549	c lossco2 is converted to first-order reaction rate
	1550	c assuming co2 number density at the surface = 2e17 mol cm-3
	1551	c
	1552	do l = 1,21 ! 70 km
	1553	t003(l) = lossco2/2.e17 ! s-1
	1554	end do
	1555	do l = 22,lswitch-1
	1556	t003(l) = 0.
	1557	end do
	1558	else
	1559	do l = 1,lswitch-1
	1560	t001(l) = 0.
	1561	t002(l) = 0.
	1562	t003(l) = 0.
	1563	end do
	1564	end if
	1565	c
	1566	return
	1567	end

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

Download in other formats:

Original Format