Context Navigation

calchim_mod.F90 @ 3726

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

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

Download in other formats:

Original Format