Context Navigation

calchim_mod.F90 @ 2883

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

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

Download in other formats:

Original Format