Context Navigation

calchim_mod.F90 @ 2530

Last change on this file since 2530 was 2461, checked in by emillour, 4 years ago

Mars GCM:

Adding the deuterium chemistry now that the HDO cycle is included.
Chemistry still works as before if deuterium tracers are not present.
Added handling of hdo in molecular diffusion (moldiff_red).

FGG+JYC+EM

File size: 37.9 KB

Rev	Line
[2162]	1	MODULE calchim_mod
	2
	3	IMPLICIT NONE
	4
[2164]	5	INTEGER, SAVE :: ichemistry ! compute chemistry every ichemistry physics step
[2162]	6	REAL,SAVE,ALLOCATABLE :: zdqchim(:,:,:) ! Tendancy on pq due to photochemistry
	7	REAL,SAVE,ALLOCATABLE :: zdqschim(:,:) ! Tendancy on qsurf due to photochemistry
	8
	9	CONTAINS
	10
[2158]	11	subroutine calchim(ngrid,nlayer,nq, &
[1495]	12	ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, &
	13	zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud, &
[2031]	14	dqscloud,tau,co2ice, &
[1495]	15	pu,pdu,pv,pdv,surfdust,surfice)
	16
	17	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, igcm_o2, &
	18	igcm_o3, igcm_h, igcm_h2, igcm_oh, igcm_ho2, &
	19	igcm_h2o2, igcm_ch4, igcm_n2, igcm_h2o_vap, &
	20	igcm_no, igcm_n, igcm_no2, igcm_n2d, &
	21	igcm_o2plus, igcm_co2plus, igcm_oplus, &
	22	igcm_coplus, igcm_cplus, igcm_nplus, &
	23	igcm_noplus, igcm_n2plus, igcm_hplus, &
[2302]	24	igcm_hco2plus, igcm_hcoplus, igcm_h2oplus, &
[2461]	25	igcm_h3oplus, igcm_ohplus, igcm_elec, &
	26	igcm_hdo_vap, igcm_od, igcm_d, igcm_hd, &
	27	igcm_do2, igcm_hdo2, mmol
[1495]	28
	29	use conc_mod, only: mmean ! mean molecular mass of the atmosphere
[2162]	30	use comcstfi_h, only: pi
[2170]	31	use photolysis_mod, only: init_photolysis, nphot
[2158]	32	use iono_h, only: temp_elect
[1495]	33
	34	implicit none
	35
	36	!=======================================================================
	37	!
	38	! subject:
	39	! --------
	40	!
	41	! Prepare the call for the photochemical module, and send back the
	42	! tendencies from photochemistry in the chemical species mass mixing ratios
	43	!
	44	! Arguments:
	45	! ----------
	46	!
	47	! Input:
	48	!
[2031]	49	! ptimestep timestep (s)
[1495]	50	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
	51	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa)
	52	! pt(ngrid,nlayer) Temperature (K)
	53	! pdt(ngrid,nlayer) Temperature tendency (K)
	54	! pu(ngrid,nlayer) u component of the wind (ms-1)
	55	! pdu(ngrid,nlayer) u component tendency (K)
	56	! pv(ngrid,nlayer) v component of the wind (ms-1)
	57	! pdv(ngrid,nlayer) v component tendency (K)
[2031]	58	! dist_sol distance of the sun (AU)
	59	! mu0(ngrid) cos of solar zenith angle (=1 when sun at zenith)
	60	! pq(ngrid,nlayer,nq) advected fields, ie chemical species here
	61	! pdq(ngrid,nlayer,nq) previous tendencies on pq
	62	! tau(ngrid) dust optical depth
	63	! co2ice(ngrid) co2 ice surface layer (kg.m-2)
[1495]	64	! surfdust(ngrid,nlayer) dust surface area (m2/m3)
	65	! surfice(ngrid,nlayer) ice surface area (m2/m3)
	66	!
	67	! Output:
	68	!
[2031]	69	! dqchim(ngrid,nlayer,nq) tendencies on pq due to chemistry
	70	! dqschim(ngrid,nq) tendencies on qsurf
[1495]	71	!
	72	!=======================================================================
	73
[2162]	74	include "callkeys.h"
[1495]	75
	76	! input:
	77
	78	integer,intent(in) :: ngrid ! number of atmospheric columns
	79	integer,intent(in) :: nlayer ! number of atmospheric layers
	80	integer,intent(in) :: nq ! number of tracers
	81	real :: ptimestep
	82	real :: pplay(ngrid,nlayer) ! pressure at the middle of the layers
	83	real :: zzlay(ngrid,nlayer) ! pressure at the middle of the layers
	84	real :: pplev(ngrid,nlayer+1) ! intermediate pressure levels
	85	real :: zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
	86	real :: pt(ngrid,nlayer) ! temperature
	87	real :: pdt(ngrid,nlayer) ! temperature tendency
	88	real :: pu(ngrid,nlayer) ! u component of the wind (m.s-1)
	89	real :: pdu(ngrid,nlayer) ! u component tendency
	90	real :: pv(ngrid,nlayer) ! v component of the wind (m.s-1)
	91	real :: pdv(ngrid,nlayer) ! v component tendency
	92	real :: dist_sol ! distance of the sun (AU)
	93	real :: mu0(ngrid) ! cos of solar zenith angle (=1 when sun at zenith)
	94	real :: pq(ngrid,nlayer,nq) ! tracers mass mixing ratio
	95	real :: pdq(ngrid,nlayer,nq) ! previous tendencies
	96	real :: zday ! date (time since Ls=0, in martian days)
[2031]	97	real :: tau(ngrid) ! dust optical depth
[1495]	98	real :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2)
	99	real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
[2031]	100	real :: surfice(ngrid,nlayer) ! ice surface area (m2/m3)
[1495]	101
	102	! output:
	103
	104	real :: dqchim(ngrid,nlayer,nq) ! tendencies on pq due to chemistry
	105	real :: dqschim(ngrid,nq) ! tendencies on qsurf
	106	real :: dqcloud(ngrid,nlayer,nq) ! tendencies on pq due to condensation
	107	real :: dqscloud(ngrid,nq) ! tendencies on qsurf
	108
	109	! local variables:
	110
	111	integer,save :: nbq ! number of tracers used in the chemistry
	112	integer,allocatable,save :: niq(:) ! array storing the indexes of the tracers
	113	integer :: iloc(1) ! index of major species
	114	integer :: ig,l,i,iq,iqmax
	115	integer :: foundswitch, lswitch
	116	integer,save :: chemthermod
	117
	118	integer,save :: i_co2 = 0
	119	integer,save :: i_co = 0
	120	integer,save :: i_o = 0
	121	integer,save :: i_o1d = 0
	122	integer,save :: i_o2 = 0
	123	integer,save :: i_o3 = 0
	124	integer,save :: i_h = 0
	125	integer,save :: i_h2 = 0
	126	integer,save :: i_oh = 0
	127	integer,save :: i_ho2 = 0
	128	integer,save :: i_h2o2 = 0
	129	integer,save :: i_ch4 = 0
	130	integer,save :: i_n2 = 0
	131	integer,save :: i_h2o = 0
	132	integer,save :: i_n = 0
	133	integer,save :: i_no = 0
	134	integer,save :: i_no2 = 0
	135	integer,save :: i_n2d = 0
	136	integer,save :: i_co2plus=0
	137	integer,save :: i_oplus=0
	138	integer,save :: i_o2plus=0
	139	integer,save :: i_coplus=0
	140	integer,save :: i_cplus=0
	141	integer,save :: i_nplus=0
	142	integer,save :: i_noplus=0
	143	integer,save :: i_n2plus=0
	144	integer,save :: i_hplus=0
	145	integer,save :: i_hco2plus=0
[2284]	146	integer,save :: i_hcoplus=0
[2302]	147	integer,save :: i_h2oplus=0
[2321]	148	integer,save :: i_h3oplus=0
	149	integer,save :: i_ohplus=0
[1495]	150	integer,save :: i_elec=0
[2461]	151	integer,save :: i_hdo=0
	152	integer,save :: i_od=0
	153	integer,save :: i_d=0
	154	integer,save :: i_hd=0
	155	integer,save :: i_do2=0
	156	integer,save :: i_hdo2=0
[1495]	157
	158	integer :: ig_vl1
	159
[2170]	160	integer :: nb_reaction_3_max ! number of quadratic reactions
	161	integer :: nb_reaction_4_max ! number of bimolecular reactions
	162	integer :: nquench ! number of quenching + heterogeneous reactions
	163	integer :: nphotion ! number of photoionizations
[2461]	164	integer :: nb_reaction_4_ion ! quadratic reactions for ionosphere
	165	integer :: nb_reaction_4_deut ! quadratic reactions for deuterium chem
[2170]	166	integer :: nb_phot_max ! total number of photolysis+photoionizations+quenching reactions
	167
	168
[1495]	169	real :: latvl1, lonvl1
	170	real :: zq(ngrid,nlayer,nq) ! pq+pdq*ptimestep before chemistry
	171	! new mole fraction after
	172	real :: zt(ngrid,nlayer) ! temperature
	173	real :: zu(ngrid,nlayer) ! u component of the wind
	174	real :: zv(ngrid,nlayer) ! v component of the wind
[2027]	175	real :: taucol ! dust optical depth at the surface
[2030]	176	real :: kb ! boltzmann constant
[1495]	177
[2170]	178	logical, save :: firstcall = .true.
	179	logical, save :: depos ! switch for dry deposition
	180	logical, save :: ionchem ! switch for ionospheric chemistry
[2461]	181	logical, save :: deutchem ! switch for deuterium chemistry
[2170]	182	logical, save :: jonline ! switch for online photodissociation rates or lookup table
	183	logical, save :: unichim ! only one unified chemical scheme at all
	184	! layers (default), or upper atmospheric
	185	! scheme in the thermosphere
[1495]	186
	187	! for each column of atmosphere:
	188
	189	real :: zpress(nlayer) ! Pressure (mbar)
	190	real :: zdens(nlayer) ! Density (cm-3)
	191	real :: ztemp(nlayer) ! Temperature (K)
[2158]	192	real :: ztelec(nlayer) ! Electronic temperature (K)
[1495]	193	real :: zlocal(nlayer) ! Altitude (km)
	194	real :: zycol(nlayer,nq) ! Composition (mole fractions)
	195	real :: zmmean(nlayer) ! Mean molecular mass (g.mole-1)
	196	real :: szacol ! Solar zenith angle
	197	real :: surfice1d(nlayer) ! Ice surface area (cm2/cm3)
	198	real :: surfdust1d(nlayer) ! Dust surface area (cm2/cm3)
	199	real :: jo3(nlayer) ! Photodissociation rate O3->O1D (s-1)
[2030]	200	real :: jh2o(nlayer) ! Photodissociation rate H2O->H+OH (s-1)
[2170]	201	real :: em_no(nlayer) ! NO nightglow emission rate
	202	real :: em_o2(nlayer) ! O2 nightglow emission rate
[1495]	203
[2170]	204	integer :: iter(nlayer) ! Number of chemical iterations
	205	! within one physical timestep
[1495]	206
	207	! for output:
	208
[2213]	209	logical,save :: output ! to issue calls to writediagfi and stats
[1495]	210	real :: jo3_3d(ngrid,nlayer) ! Photodissociation rate O3->O1D (s-1)
[2030]	211	real :: jh2o_3d(ngrid,nlayer) ! Photodissociation rate H2O->H+OH (s-1)
[1888]	212	real :: emission_no(ngrid,nlayer) !NO emission rate
	213	real :: emission_o2(ngrid,nlayer) !O2 emission rate
[1495]	214	real :: iter_3d(ngrid,nlayer) ! Number of chemical iterations
	215	! within one physical timestep
	216
[2170]	217	!=======================================================================
[2213]	218	! initialization of the chemistry (first call only)
[2170]	219	!=======================================================================
[2158]	220
[2213]	221	if (firstcall) then
[2170]	222
[1495]	223	!=======================================================================
[2213]	224	! main dashboard for the chemistry
[1495]	225	!=======================================================================
	226
[2213]	227	unichim = .true. ! true : unified chemistry ! false : separate models in lower and upper atmosphere
	228	jonline = .true. ! true : on-line calculation of photodissociation rates ! false : lookup table
	229	ionchem = .false. ! switch for ionospheric chemistry
[2461]	230	deutchem= .false. ! switch for deuterium chemistry
[2213]	231	depos = .false. ! switch for dry deposition
[2461]	232	output = .true. ! true : triggers writing of specific outputs (photolysis rates, emission rates, etc)
[1495]	233
[2461]	234	! if (photochem) then
	235	! if (jonline) then
	236	! print*,'calchim: Read UV absorption cross-sections'
	237	! call init_photolysis ! on-line photolysis
	238	! else
	239	! print*,'calchim: Read photolysis lookup table'
	240	! call read_phototable ! off-line photolysis
	241	! end if
	242	! end if
[2158]	243
[2461]	244	! if(.not.unichim) then
[2158]	245	!Read reaction rates from external file if the upper atmospheric
	246	!chemistry is called
[2461]	247	! call chemthermos_readini
	248	! endif
[2158]	249
[1495]	250	! find index of chemical tracers to use
	251	allocate(niq(nq))
	252	! Listed here are all tracers that can go into photochemistry
	253	nbq = 0 ! to count number of tracers
[2158]	254	! Species ALWAYS present if photochem=.T.
[1495]	255	i_co2 = igcm_co2
	256	if (i_co2 == 0) then
	257	write(,) "calchim: Error; no CO2 tracer !!!"
[2302]	258	call abort_physic("calchim","missing co2 tracer",1)
[1495]	259	else
	260	nbq = nbq + 1
	261	niq(nbq) = i_co2
	262	end if
	263	i_co = igcm_co
	264	if (i_co == 0) then
	265	write(,) "calchim: Error; no CO tracer !!!"
[2302]	266	call abort_physic("calchim","missing co tracer",1)
[1495]	267	else
	268	nbq = nbq + 1
	269	niq(nbq) = i_co
	270	end if
	271	i_o = igcm_o
	272	if (i_o == 0) then
	273	write(,) "calchim: Error; no O tracer !!!"
[2302]	274	call abort_physic("calchim","missing o tracer",1)
[1495]	275	else
	276	nbq = nbq + 1
	277	niq(nbq) = i_o
	278	end if
	279	i_o1d = igcm_o1d
	280	if (i_o1d == 0) then
	281	write(,) "calchim: Error; no O1D tracer !!!"
[2302]	282	call abort_physic("calchim","missing o1d tracer",1)
[1495]	283	else
	284	nbq = nbq + 1
	285	niq(nbq) = i_o1d
	286	end if
	287	i_o2 = igcm_o2
	288	if (i_o2 == 0) then
	289	write(,) "calchim: Error; no O2 tracer !!!"
[2302]	290	call abort_physic("calchim","missing o2 tracer",1)
[1495]	291	else
	292	nbq = nbq + 1
	293	niq(nbq) = i_o2
	294	end if
	295	i_o3 = igcm_o3
	296	if (i_o3 == 0) then
	297	write(,) "calchim: Error; no O3 tracer !!!"
[2302]	298	call abort_physic("calchim","missing o3 tracer",1)
[1495]	299	else
	300	nbq = nbq + 1
	301	niq(nbq) = i_o3
	302	end if
	303	i_h = igcm_h
	304	if (i_h == 0) then
	305	write(,) "calchim: Error; no H tracer !!!"
[2302]	306	call abort_physic("calchim","missing h tracer",1)
[1495]	307	else
	308	nbq = nbq + 1
	309	niq(nbq) = i_h
	310	end if
	311	i_h2 = igcm_h2
	312	if (i_h2 == 0) then
	313	write(,) "calchim: Error; no H2 tracer !!!"
[2302]	314	call abort_physic("calchim","missing h2 tracer",1)
[1495]	315	else
	316	nbq = nbq + 1
	317	niq(nbq) = i_h2
	318	end if
	319	i_oh = igcm_oh
	320	if (i_oh == 0) then
	321	write(,) "calchim: Error; no OH tracer !!!"
[2302]	322	call abort_physic("calchim","missing oh tracer",1)
[1495]	323	else
	324	nbq = nbq + 1
	325	niq(nbq) = i_oh
	326	end if
	327	i_ho2 = igcm_ho2
	328	if (i_ho2 == 0) then
	329	write(,) "calchim: Error; no HO2 tracer !!!"
[2302]	330	call abort_physic("calchim","missing ho2 tracer",1)
[1495]	331	else
	332	nbq = nbq + 1
	333	niq(nbq) = i_ho2
	334	end if
	335	i_h2o2 = igcm_h2o2
	336	if (i_h2o2 == 0) then
	337	write(,) "calchim: Error; no H2O2 tracer !!!"
[2302]	338	call abort_physic("calchim","missing h2o2 tracer",1)
[1495]	339	else
	340	nbq = nbq + 1
	341	niq(nbq) = i_h2o2
	342	end if
	343	i_ch4 = igcm_ch4
	344	if (i_ch4 == 0) then
	345	write(,) "calchim: Error; no CH4 tracer !!!"
	346	write(,) "CH4 will be ignored in the chemistry"
	347	else
	348	nbq = nbq + 1
	349	niq(nbq) = i_ch4
	350	end if
	351	i_n2 = igcm_n2
	352	if (i_n2 == 0) then
	353	write(,) "calchim: Error; no N2 tracer !!!"
[2302]	354	call abort_physic("calchim","missing n2 tracer",1)
[1495]	355	else
	356	nbq = nbq + 1
	357	niq(nbq) = i_n2
	358	end if
	359	i_n = igcm_n
	360	if (i_n == 0) then
	361	if (photochem) then
	362	write(,) "calchim: Error; no N tracer !!!"
[2302]	363	call abort_physic("calchim","missing n tracer",1)
[1495]	364	end if
	365	else
	366	nbq = nbq + 1
	367	niq(nbq) = i_n
	368	end if
	369	i_n2d = igcm_n2d
	370	if (i_n2d == 0) then
	371	if (photochem) then
	372	write(,) "calchim: Error; no N2D tracer !!!"
[2302]	373	call abort_physic("calchim","missing n2d tracer",1)
[1495]	374	end if
	375	else
	376	nbq = nbq + 1
	377	niq(nbq) = i_n2d
	378	end if
	379	i_no = igcm_no
	380	if (i_no == 0) then
	381	if (photochem) then
	382	write(,) "calchim: Error; no NO tracer !!!"
[2302]	383	call abort_physic("calchim","missing no tracer",1)
[1495]	384	end if
	385	else
	386	nbq = nbq + 1
	387	niq(nbq) = i_no
	388	end if
	389	i_no2 = igcm_no2
	390	if (i_no2 == 0) then
	391	if (photochem) then
	392	write(,) "calchim: Error; no NO2 tracer !!!"
[2302]	393	call abort_physic("calchim","missing no2 tracer",1)
[1495]	394	end if
	395	else
	396	nbq = nbq + 1
	397	niq(nbq) = i_no2
	398	end if
	399	i_h2o = igcm_h2o_vap
	400	if (i_h2o == 0) then
	401	write(,) "calchim: Error; no water vapor tracer !!!"
[2302]	402	call abort_physic("calchim","missing h2o_vap tracer",1)
[1495]	403	else
	404	nbq = nbq + 1
	405	niq(nbq) = i_h2o
	406	end if
[2461]	407	i_hdo=igcm_hdo_vap
	408	if (i_hdo == 0) then
	409	write(,) "calchim: no HDO tracer !!!"
	410	! call abort_physic("calchim","missing hdo_vap tracer",1)
	411	write(,) "No deuterium chemistry considered"
	412	else
	413	nbq = nbq + 1
	414	niq(nbq) = i_hdo
	415	deutchem = .true.
	416	write(,) "calchim: HDO tracer found in traceur.def"
	417	write(,) "Deuterium chemistry included"
	418	end if
	419	i_od=igcm_od
	420	if(deutchem) then
	421	if (i_od == 0) then
	422	write(,) "calchim: Error, no OD tracer !!!"
	423	write(,) "OD is needed if HDO is in traceur.def"
	424	call abort_physic("calchim","missing od tracer",1)
	425	else
	426	nbq = nbq + 1
	427	niq(nbq) = i_od
	428	end if
	429	else
	430	if (i_oplus /= 0) then
	431	write(,) "calchim: Error: OD is present, but HDO is not!!!"
	432	write(,) "Both must be in traceur.def if deuterium chemistry wanted"
	433	call abort_physic("calchim","missing hdo tracer",1)
	434	endif
	435	endif
	436	i_d=igcm_d
	437	if(deutchem) then
	438	if (i_d == 0) then
	439	write(,) "calchim: Error, no D tracer !!!"
	440	write(,) "D is needed if HDO is in traceur.def"
	441	call abort_physic("calchim","missing d tracer",1)
	442	else
	443	nbq = nbq + 1
	444	niq(nbq) = i_d
	445	end if
	446	else
	447	if (i_d /= 0) then
	448	write(,) "calchim: Error: D is present, but HDO is not!!!"
	449	write(,) "Both must be in traceur.def if deuterium chemistry wanted"
	450	call abort_physic("calchim","missing hdo tracer",1)
	451	endif
	452	endif
	453	i_hd=igcm_hd
	454	if(deutchem) then
	455	if (i_hd == 0) then
	456	write(,) "calchim: Error, no HD tracer !!!"
	457	write(,) "HD is needed if HDO is in traceur.def"
	458	call abort_physic("calchim","missing hd tracer",1)
	459	else
	460	nbq = nbq + 1
	461	niq(nbq) = i_hd
	462	end if
	463	else
	464	if (i_hd /= 0) then
	465	write(,) "calchim: Error: HD is present, but HDO is not!!!"
	466	write(,) "Both must be in traceur.def if deuterium chemistry wanted"
	467	call abort_physic("calchim","missing hd tracer",1)
	468	endif
	469	endif
	470	i_do2=igcm_do2
	471	if(deutchem) then
	472	if (i_do2 == 0) then
	473	write(,) "calchim: Error, no DO2 tracer !!!"
	474	write(,) "DO2 is needed if HDO is in traceur.def"
	475	call abort_physic("calchim","missing do2 tracer",1)
	476	else
	477	nbq = nbq + 1
	478	niq(nbq) = i_do2
	479	end if
	480	else
	481	if (i_do2 /= 0) then
	482	write(,) "calchim: Error: DO2 is present, but HDO is not!!!"
	483	write(,) "Both must be in traceur.def if deuterium chemistry wanted"
	484	call abort_physic("calchim","missing do2 tracer",1)
	485	endif
	486	endif
	487	i_hdo2=igcm_hdo2
	488	if(deutchem) then
	489	if (i_hdo2 == 0) then
	490	write(,) "calchim: Error, no HDO2 tracer !!!"
	491	write(,) "HDO2 is needed if HDO is in traceur.def"
	492	call abort_physic("calchim","missing hdo2 tracer",1)
	493	else
	494	nbq = nbq + 1
	495	niq(nbq) = i_hdo2
	496	end if
	497	else
	498	if (i_hdo2 /= 0) then
	499	write(,) "calchim: Error: HDO2 is present, but HDO is not!!!"
	500	write(,) "Both must be in traceur.def if deuterium chemistry wanted"
	501	call abort_physic("calchim","missing hdo2 tracer",1)
	502	endif
	503	endif
[2158]	504	i_o2plus = igcm_o2plus
	505	if (i_o2plus == 0) then
[2170]	506	write(,) "calchim: no O2+ tracer "
	507	write(,) "Only neutral chemistry"
[2158]	508	else
	509	nbq = nbq + 1
	510	niq(nbq) = i_o2plus
[2170]	511	ionchem = .true.
	512	write(,) "calchim: O2+ tracer found in traceur.def"
	513	write(,) "Ion chemistry included"
[2158]	514	end if
	515	i_co2plus = igcm_co2plus
[2170]	516	if(ionchem) then
	517	if (i_co2plus == 0) then
	518	write(,) "calchim: Error, no CO2+ tracer !!!"
	519	write(,) "CO2+ is needed if O2+ is in traceur.def"
[2302]	520	call abort_physic("calchim","missing co2plus tracer",1)
[2170]	521	else
	522	nbq = nbq + 1
	523	niq(nbq) = i_co2plus
	524	end if
[2158]	525	else
[2170]	526	if (i_co2plus /= 0) then
	527	write(,) "calchim: Error: CO2+ is present, but O2+ is not!!!"
	528	write(,) "Both must be in traceur.def if ion chemistry wanted"
[2302]	529	call abort_physic("calchim","missing o2plus tracer",1)
[2170]	530	endif
	531	endif
[2158]	532	i_oplus=igcm_oplus
[2170]	533	if(ionchem) then
	534	if (i_oplus == 0) then
	535	write(,) "calchim: Error, no O+ tracer !!!"
	536	write(,) "O+ is needed if O2+ is in traceur.def"
[2302]	537	call abort_physic("calchim","missing oplus tracer",1)
[2170]	538	else
	539	nbq = nbq + 1
	540	niq(nbq) = i_oplus
	541	end if
[2158]	542	else
[2170]	543	if (i_oplus /= 0) then
	544	write(,) "calchim: Error: O+ is present, but O2+ is not!!!"
	545	write(,) "Both must be in traceur.def if ion chemistry wanted"
[2302]	546	call abort_physic("calchim","missing o2plus tracer",1)
[2170]	547	endif
	548	endif
[2158]	549	i_noplus=igcm_noplus
[2170]	550	if(ionchem) then
	551	if (i_noplus == 0) then
	552	write(,) "calchim: Error, no NO+ tracer !!!"
	553	write(,) "NO+ is needed if O2+ is in traceur.def"
[2302]	554	call abort_physic("calchim","missing noplus tracer",1)
[2170]	555	else
	556	nbq = nbq + 1
	557	niq(nbq) = i_noplus
	558	end if
[2158]	559	else
[2170]	560	if (i_noplus /= 0) then
	561	write(,) "calchim: Error: NO+ is present, but O2+ is not!!!"
	562	write(,) "Both must be in traceur.def if ion chemistry wanted"
	563	endif
	564	endif
[2158]	565	i_coplus=igcm_coplus
[2170]	566	if(ionchem) then
	567	if (i_coplus == 0) then
	568	write(,) "calchim: Error, no CO+ tracer !!!"
	569	write(,) "CO+ is needed if O2+ is in traceur.def"
[2302]	570	call abort_physic("calchim","missing coplus tracer",1)
[2170]	571	else
	572	nbq = nbq + 1
	573	niq(nbq) = i_coplus
	574	end if
[2158]	575	else
[2170]	576	if (i_coplus /= 0) then
	577	write(,) "calchim: Error: CO+ is present, but O2+ is not!!!"
	578	write(,) "Both must be in traceur.def if ion chemistry wanted"
	579	endif
	580	endif
[2158]	581	i_cplus=igcm_cplus
[2170]	582	if(ionchem) then
	583	if (i_cplus == 0) then
	584	write(,) "calchim: Error, no C+ tracer !!!"
	585	write(,) "C+ is needed if O2+ is in traceur.def"
[2302]	586	call abort_physic("calchim","missing cplus tracer",1)
[2170]	587	else
	588	nbq = nbq + 1
	589	niq(nbq) = i_cplus
	590	end if
[2158]	591	else
[2170]	592	if (i_cplus /= 0) then
	593	write(,) "calchim: Error: C+ is present, but O2+ is not!!!"
	594	write(,) "Both must be in traceur.def if ion chemistry wanted"
	595	endif
	596	endif
[2158]	597	i_n2plus=igcm_n2plus
[2170]	598	if(ionchem) then
	599	if (i_n2plus == 0) then
	600	write(,) "calchim: Error, no N2+ tracer !!!"
	601	write(,) "N2+ is needed if O2+ is in traceur.def"
[2302]	602	call abort_physic("calchim","missing n2plus tracer",1)
[2170]	603	else
	604	nbq = nbq + 1
	605	niq(nbq) = i_n2plus
	606	end if
[2158]	607	else
[2170]	608	if (i_n2plus /= 0) then
	609	write(,) "calchim: Error: N2+ is present, but O2+ is not!!!"
	610	write(,) "Both must be in traceur.def if ion chemistry wanted"
	611	endif
	612	endif
[2158]	613	i_nplus=igcm_nplus
[2170]	614	if(ionchem) then
	615	if (i_nplus == 0) then
	616	write(,) "calchim: Error, no N+ tracer !!!"
	617	write(,) "N+ is needed if O2+ is in traceur.def"
[2302]	618	call abort_physic("calchim","missing nplus tracer",1)
[2170]	619	else
	620	nbq = nbq + 1
	621	niq(nbq) = i_nplus
	622	end if
[2158]	623	else
[2170]	624	if (i_nplus /= 0) then
	625	write(,) "calchim: Error: N+ is present, but O2+ is not!!!"
	626	write(,) "Both must be in traceur.def if ion chemistry wanted"
	627	endif
	628	endif
[2158]	629	i_hplus=igcm_hplus
[2170]	630	if(ionchem) then
	631	if (i_hplus == 0) then
	632	write(,) "calchim: Error, no H+ tracer !!!"
	633	write(,) "H+ is needed if O2+ is in traceur.def"
[2302]	634	call abort_physic("calchim","missing hplus tracer",1)
[2170]	635	else
	636	nbq = nbq + 1
	637	niq(nbq) = i_hplus
	638	end if
[2158]	639	else
[2170]	640	if (i_hplus /= 0) then
	641	write(,) "calchim: Error: H+ is present, but O2+ is not!!!"
	642	write(,) "Both must be in traceur.def if ion chemistry wanted"
	643	endif
	644	endif
[2158]	645	i_hco2plus=igcm_hco2plus
[2170]	646	if(ionchem) then
	647	if (i_hco2plus == 0) then
	648	write(,) "calchim: Error, no HCO2+ tracer !!!"
	649	write(,) "HCO2+ is needed if O2+ is in traceur.def"
[2302]	650	call abort_physic("calchim","missing hco2plus tracer",1)
[2170]	651	else
	652	nbq = nbq + 1
	653	niq(nbq) = i_hco2plus
	654	end if
[2158]	655	else
[2170]	656	if (i_hco2plus /= 0) then
	657	write(,) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
	658	write(,) "Both must be in traceur.def if ion chemistry wanted"
	659	endif
	660	endif
[2284]	661	i_hcoplus=igcm_hcoplus
	662	if(ionchem) then
	663	if (i_hcoplus == 0) then
	664	write(,) "calchim: Error, no HCO+ tracer !!!"
	665	write(,) "HCO+ is needed if O2+ is in traceur.def"
[2302]	666	call abort_physic("calchim","missing hcoplus tracer",1)
[2284]	667	else
	668	nbq = nbq + 1
	669	niq(nbq) = i_hcoplus
	670	end if
	671	else
	672	if (i_hcoplus /= 0) then
	673	write(,) "calchim: Error: HCO+ is present, but O2+ is not!!!"
	674	write(,) "Both must be in traceur.def if ion chemistry wanted"
	675	endif
	676	endif
[2302]	677	i_h2oplus=igcm_h2oplus
	678	if(ionchem) then
	679	if (i_h2oplus == 0) then
	680	write(,) "calchim: Error, no H2O+ tracer !!!"
	681	write(,) "H2O+ is needed if O2+ is in traceur.def"
	682	call abort_physic("calchim","missing h2oplus tracer",1)
	683	else
	684	nbq = nbq + 1
	685	niq(nbq) = i_h2oplus
	686	end if
	687	else
	688	if (i_h2oplus /= 0) then
	689	write(,) "calchim: Error: H2O+ is present, but O2+ is not!!!"
	690	write(,) "Both must be in traceur.def if ion chemistry wanted"
	691	endif
	692	endif
[2321]	693	i_h3oplus=igcm_h3oplus
	694	if(ionchem) then
	695	if (i_h3oplus == 0) then
	696	write(,) "calchim: Error, no H3O+ tracer !!!"
	697	write(,) "H3O+ is needed if O2+ is in traceur.def"
	698	call abort_physic("calchim","missing h3oplus tracer",1)
	699	else
	700	nbq = nbq + 1
	701	niq(nbq) = i_h3oplus
	702	end if
	703	else
	704	if (i_h3oplus /= 0) then
	705	write(,) "calchim: Error: H3O+ is present, but O2+ is not!!!"
	706	write(,) "Both must be in traceur.def if ion chemistry wanted"
	707	endif
	708	endif
	709	i_ohplus=igcm_ohplus
	710	if(ionchem) then
	711	if (i_ohplus == 0) then
	712	write(,) "calchim: Error, no OH+ tracer !!!"
	713	write(,) "OH+ is needed if O2+ is in traceur.def"
	714	call abort_physic("calchim","missing ohplus tracer",1)
	715	else
	716	nbq = nbq + 1
	717	niq(nbq) = i_ohplus
	718	end if
	719	else
	720	if (i_ohplus /= 0) then
	721	write(,) "calchim: Error: OH+ is present, but O2+ is not!!!"
	722	write(,) "Both must be in traceur.def if ion chemistry wanted"
	723	endif
	724	endif
[2158]	725	i_elec = igcm_elec
[2170]	726	if(ionchem) then
	727	if (i_elec == 0) then
	728	write(,) "calchim: Error, no e- tracer !!!"
	729	write(,) "e- is needed if O2+ is in traceur.def"
[2302]	730	call abort_physic("calchim","missing elec tracer",1)
[2170]	731	else
	732	nbq = nbq + 1
	733	niq(nbq) = i_elec
	734	end if
[2158]	735	else
[2170]	736	if (i_elec /= 0) then
	737	write(,) "calchim: Error: e- is present, but O2+ is not!!!"
	738	write(,) "Both must be in traceur.def if ion chemistry wanted"
	739	endif
	740	endif
[1495]	741	write(,) 'calchim: found nbq = ',nbq,' tracers'
[2461]	742
	743	write(,) 'calchim: tracer indices=',niq(:)
	744
	745
	746	if (photochem) then
	747	if (jonline) then
	748	print*,'calchim: Read UV absorption cross-sections'
	749	!Add two photodissociations in deuterium chemistry included
	750	if(deutchem) nphot = nphot + 2
	751	call init_photolysis ! on-line photolysis
	752	else
	753	print*,'calchim: Read photolysis lookup table'
	754	call read_phototable ! off-line photolysis
	755	end if
	756	end if
	757
	758	if(.not.unichim) then
	759	!Read reaction rates from external file if the upper atmospheric
	760	!chemistry is called
	761	call chemthermos_readini
	762	endif
	763
	764	! stop
[1495]	765	firstcall = .false.
	766	end if ! if (firstcall)
	767
	768	! Initializations
	769
	770	zycol(:,:) = 0.
	771	dqchim(:,:,:) = 0.
	772	dqschim(:,:) = 0.
	773
[2030]	774	kb = 1.3806e-23
	775
[1495]	776	! latvl1= 22.27
	777	! lonvl1= -47.94
[2029]	778	! ig_vl1= 1+ int( (1.5-(latvl1-90.)48./180.) -2 )64. + &
	779	! int(1.5+(lonvl1+180)*64./360.)
[1495]	780
	781	!=======================================================================
	782	! loop over grid
	783	!=======================================================================
	784
	785	do ig = 1,ngrid
	786
	787	foundswitch = 0
	788	do l = 1,nlayer
	789	do i = 1,nbq
	790	iq = niq(i) ! get tracer index
	791	zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
	792	zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
	793	end do
	794	zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
	795	zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
	796	zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
	797	zpress(l) = pplay(ig,l)/100.
	798	ztemp(l) = zt(ig,l)
	799	zdens(l) = zpress(l)/(kb1.e4ztemp(l))
	800	zlocal(l) = zzlay(ig,l)/1000.
	801	zmmean(l) = mmean(ig,l)
[2158]	802	ztelec(l) = temp_elect(zlocal(l),ztemp(l),1)
	803	!Electronic temperature. Index 1 -> Viking; Index 2-> MAVEN
[1495]	804
	805	! surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
	806
	807	surfdust1d(l) = surfdust(ig,l)*1.e-2
	808	surfice1d(l) = surfice(ig,l)*1.e-2
	809
[2158]	810	! search for switch index between regions
[2170]	811
	812	if (unichim) then
	813	lswitch = nlayer + 1
[2158]	814	else
[2461]	815	if (foundswitch == 0 .and. pplay(ig,l) < 10.) then
[1495]	816	lswitch = l
	817	foundswitch = 1
	818	end if
[2158]	819	endif
[1495]	820	end do ! of do l=1,nlayer
	821
	822	szacol = acos(mu0(ig))*180./pi
[2031]	823	taucol = tau(ig)
[1495]	824
	825	!=======================================================================
	826	! call chemical subroutines
	827	!=======================================================================
	828
	829	if (photochem) then
[2170]	830	! set number of reactions, depending on ion chemistry or not
[2461]	831	nb_reaction_4_ion = 64
	832	nb_reaction_4_deut = 35
	833
	834	!Default numbers if no ion and no deuterium chemistry included
	835
	836	nb_reaction_4_max = 31 ! set number of bimolecular reactions
	837	nb_reaction_3_max = 6 ! set number of quadratic reactions
	838	nquench = 9 ! set number of quenching + heterogeneous
	839	nphotion = 0 ! set number of photoionizations
	840
[2170]	841	if (ionchem) then
[2461]	842	nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_ion
[2170]	843	nphotion = 18 ! set number of photoionizations
[2461]	844	endif
	845	if(deutchem) then
	846	nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_deut
[2170]	847	end if
	848
	849	! nb_phot_max is the total number of processes that are treated
	850	! numerically as a photolysis:
	851
	852	nb_phot_max = nphot + nphotion + nquench
	853
	854	! call main photochemistry routine
	855
[2461]	856	call photochemistry(nlayer,nq,nbq,ionchem,deutchem, &
	857	nb_reaction_3_max,nb_reaction_4_max,nphot, &
	858	nb_phot_max,nphotion, &
[2433]	859	jonline,ig,lswitch,zycol,szacol,ptimestep, &
[2170]	860	zpress,zlocal,ztemp,ztelec,zdens,zmmean, &
	861	dist_sol,zday,surfdust1d,surfice1d, &
[2321]	862	jo3,jh2o,em_no,em_o2,taucol,iter)
[1495]	863
	864	! ozone photolysis, for output
	865
	866	do l = 1,nlayer
	867	jo3_3d(ig,l) = jo3(l)
[2030]	868	jh2o_3d(ig,l) = jh2o(l)
[1495]	869	iter_3d(ig,l) = iter(l)
	870	end do
[2170]	871
[1495]	872	! condensation of h2o2
	873
	874	call perosat(ngrid, nlayer, nq, &
	875	ig,ptimestep,pplev,pplay, &
	876	ztemp,zycol,dqcloud,dqscloud)
	877
[2170]	878	! case of separate photochemical model in the thermosphere
	879
	880	if (.not.unichim) then
	881	chemthermod = 3 !C/O/H/N/ions
[2158]	882	call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
	883	zdens,zpress,zlocal,szacol,ptimestep,zday,&
	884	em_no,em_o2)
	885	end if
[1495]	886
[2321]	887	do l = 1,nlayer
	888	emission_no(ig,l) = em_no(l)
	889	emission_o2(ig,l) = em_o2(l)
	890	end do
[2170]	891	end if ! photochem
	892
[1495]	893	! dry deposition
	894
	895	if (depos) then
	896	call deposition(ngrid, nlayer, nq, &
	897	ig, ig_vl1, pplay, pplev, zzlay, zzlev, &
	898	zu, zv, zt, zycol, ptimestep, co2ice)
	899	end if
	900
	901	!=======================================================================
	902	! tendencies
	903	!=======================================================================
	904
	905	! index of the most abundant species at each level
	906
	907	! major(:) = maxloc(zycol, dim = 2)
	908
	909	! tendency for the most abundant species = - sum of others
	910
	911	do l = 1,nlayer
	912	iloc=maxloc(zycol(l,:))
	913	iqmax=iloc(1)
	914	do i = 1,nbq
	915	iq = niq(i) ! get tracer index
	916	if (iq /= iqmax) then
	917	dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l) &
	918	- zq(ig,l,iq))/ptimestep
	919	dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) &
	920	- dqchim(ig,l,iq)
	921	end if
	922	end do
	923	end do ! of do l = 1,nlayer
	924
	925	!=======================================================================
	926	! end of loop over grid
	927	!=======================================================================
	928
	929	end do ! of do ig=1,ngrid
	930
	931	!=======================================================================
	932	! write outputs
	933	!=======================================================================
	934
	935	! value of parameter 'output' to trigger writting of outputs
	936	! is set above at the declaration of the variable.
	937
	938	if (photochem .and. output) then
	939	if (ngrid > 1) then
	940	call writediagfi(ngrid,'jo3','j o3->o1d', &
	941	's-1',3,jo3_3d(1,1))
[2030]	942	call writediagfi(ngrid,'jh2o','jh2o', &
	943	's-1',3,jh2o_3d(1,1))
[1495]	944	call writediagfi(ngrid,'iter','iterations', &
	945	' ',3,iter_3d(1,1))
[2170]	946
[2321]	947	! if (.not. unichim) then
[2158]	948	call writediagfi(ngrid,'emission_no', &
	949	'NO nightglow emission rate','cm-3 s-1',3,emission_no)
	950	call writediagfi(ngrid,'emission_o2', &
	951	'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
[2321]	952	! endif
[2042]	953
[1495]	954	if (callstats) then
	955	call wstats(ngrid,'jo3','j o3->o1d', &
	956	's-1',3,jo3_3d(1,1))
[1888]	957	call wstats(ngrid,'emission_no', &
	958	'NO nightglow emission rate','cm-3 s-1',3,emission_no)
	959	call wstats(ngrid,'emission_o2', &
	960	'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
[1495]	961	call wstats(ngrid,'mmean','mean molecular mass', &
	962	'g.mole-1',3,mmean(1,1))
	963	endif
	964	end if ! of if (ngrid.gt.1)
	965	end if ! of if (output)
	966
[2162]	967	end subroutine calchim
	968
	969
	970	subroutine ini_calchim_mod(ngrid,nlayer,nq)
	971
	972	implicit none
	973
	974	integer,intent(in) :: ngrid ! number of atmospheric columns
	975	integer,intent(in) :: nlayer ! number of atmospheric layers
	976	integer,intent(in) :: nq ! number of tracers
	977
	978	allocate(zdqchim(ngrid,nlayer,nq))
	979	zdqchim(:,:,:)=0
	980	allocate(zdqschim(ngrid,nq))
	981	zdqschim(:,:)=0
	982
	983	end subroutine ini_calchim_mod
	984
	985
	986	subroutine end_calchim_mod
	987
	988	implicit none
	989
	990	if (allocated(zdqchim)) deallocate(zdqchim)
	991	if (allocated(zdqschim)) deallocate(zdqschim)
	992
	993	end subroutine end_calchim_mod
	994
	995	END MODULE calchim_mod
	996

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

Download in other formats:

Original Format