Context Navigation

calchim_mod.F90 @ 2573

Last change on this file since 2573 was 2563, checked in by emillour, 4 years ago

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

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

Download in other formats:

Original Format