Context Navigation

calchim_mod.F90 @ 3595

Last change on this file since 3595 was 3466, checked in by emillour, 4 months ago

Mars PCM:
More tidying in aeronomars:

remove unused "inv.F" and remove "dtridgl.F" which is not used here and is a duplicate of the "dtridgl" routine in phymars/swr_toon.F
turn aeronomars routines to modules, for those which aren't in modules yet.

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

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

Download in other formats:

Original Format