Context Navigation

calchim.F90 @ 1918

Last change on this file since 1918 was 1888, checked in by emillour, 7 years ago

Mars GCM:
Updated the calculation of the dissociation and ionization
branching ratios, using the values in the Schunk and Nagy book.
New datafiles *branchingratio_schunkandnagy2000_param.dat must be loaded
for the "EUVDAT" subdirectory of the standard "datadir" directory.
Main changes are:

param_v4_h.F90 -> New definition for the O2 ionization branching ratio
param_read_e107.F -> Read the new files containing the S&N branching ratios
paramfoto_compact.F -> Mainly cleaning of the code and the comments.

Also correction of a bug affecting the calculation of CO losses

chemthermos.F90 -> Small modification to add the possibility of including

NO and O2 nightglow rates to the outputs

calchim.F90 and calchim_asis.F90 -> account for change in arguments in

calls to chemthermos
FGG

File size: 27.4 KB

Rev	Line
[1047]	1	subroutine calchim(ngrid,nlayer,nq, &
[1036]	2	ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, &
[618]	3	zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud, &
	4	dqscloud,tauref,co2ice, &
	5	pu,pdu,pv,pdv,surfdust,surfice)
	6
[1036]	7	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, igcm_o2, &
	8	igcm_o3, igcm_h, igcm_h2, igcm_oh, igcm_ho2, &
	9	igcm_h2o2, igcm_ch4, igcm_n2, igcm_h2o_vap, &
	10	igcm_no, igcm_n, igcm_no2, igcm_n2d, &
	11	igcm_o2plus, igcm_co2plus, igcm_oplus, &
	12	igcm_coplus, igcm_cplus, igcm_nplus, &
	13	igcm_noplus, igcm_n2plus, igcm_hplus, &
	14	igcm_hco2plus, igcm_elec, mmol
[1125]	15
[1047]	16	use conc_mod, only: mmean ! mean molecular mass of the atmosphere
[1226]	17	USE comcstfi_h
[1047]	18
[38]	19	implicit none
	20
[618]	21	!=======================================================================
	22	!
	23	! subject:
	24	! --------
	25	!
	26	! Prepare the call for the photochemical module, and send back the
	27	! tendencies from photochemistry in the chemical species mass mixing ratios
	28	!
	29	! Author: Sebastien Lebonnois (08/11/2002)
	30	! -------
	31	! update 12/06/2003 for water ice clouds and compatibility with dust
	32	! update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
	33	! update 03/05/2005 cosmetic changes (Franck Lefevre)
	34	! update sept. 2008 identify tracers by their names (Ehouarn Millour)
	35	! update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
	36	! update 16/03/2012 optimization (Franck Lefevre)
[1125]	37	! update 11/12/2013 optimization (Franck Lefevre)
[618]	38	!
	39	! Arguments:
	40	! ----------
	41	!
	42	! Input:
	43	!
	44	! ptimestep timestep (s)
[1047]	45	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
	46	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa)
	47	! pt(ngrid,nlayer) Temperature (K)
	48	! pdt(ngrid,nlayer) Temperature tendency (K)
	49	! pu(ngrid,nlayer) u component of the wind (ms-1)
	50	! pdu(ngrid,nlayer) u component tendency (K)
	51	! pv(ngrid,nlayer) v component of the wind (ms-1)
	52	! pdv(ngrid,nlayer) v component tendency (K)
[618]	53	! dist_sol distance of the sun (AU)
[1047]	54	! mu0(ngrid) cos of solar zenith angle (=1 when sun at zenith)
	55	! pq(ngrid,nlayer,nq) Advected fields, ie chemical species here
	56	! pdq(ngrid,nlayer,nq) Previous tendencies on pq
	57	! tauref(ngrid) Optical depth at 7 hPa
	58	! co2ice(ngrid) co2 ice surface layer (kg.m-2)
	59	! surfdust(ngrid,nlayer) dust surface area (m2/m3)
	60	! surfice(ngrid,nlayer) ice surface area (m2/m3)
[618]	61	!
	62	! Output:
	63	!
[1047]	64	! dqchim(ngrid,nlayer,nq) ! tendencies on pq due to chemistry
	65	! dqschim(ngrid,nq) ! tendencies on qsurf
[618]	66	!
	67	!=======================================================================
[38]	68
	69	#include "chimiedata.h"
	70	#include "callkeys.h"
	71
[618]	72	! input:
[38]	73
[1125]	74	integer,intent(in) :: ngrid ! number of atmospheric columns
	75	integer,intent(in) :: nlayer ! number of atmospheric layers
	76	integer,intent(in) :: nq ! number of tracers
[618]	77	real :: ptimestep
[1047]	78	real :: pplay(ngrid,nlayer) ! pressure at the middle of the layers
	79	real :: zzlay(ngrid,nlayer) ! pressure at the middle of the layers
	80	real :: pplev(ngrid,nlayer+1) ! intermediate pressure levels
	81	real :: zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
	82	real :: pt(ngrid,nlayer) ! temperature
	83	real :: pdt(ngrid,nlayer) ! temperature tendency
	84	real :: pu(ngrid,nlayer) ! u component of the wind (m.s-1)
	85	real :: pdu(ngrid,nlayer) ! u component tendency
	86	real :: pv(ngrid,nlayer) ! v component of the wind (m.s-1)
	87	real :: pdv(ngrid,nlayer) ! v component tendency
[1125]	88	real :: dist_sol ! distance of the sun (AU)
	89	real :: mu0(ngrid) ! cos of solar zenith angle (=1 when sun at zenith)
	90	real :: pq(ngrid,nlayer,nq) ! tracers mass mixing ratio
	91	real :: pdq(ngrid,nlayer,nq) ! previous tendencies
	92	real :: zday ! date (time since Ls=0, in martian days)
	93	real :: tauref(ngrid) ! optical depth at 7 hPa
	94	real :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2)
[1047]	95	real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
	96	real :: surfice(ngrid,nlayer) ! ice surface area (m2/m3)
[38]	97
[618]	98	! output:
[38]	99
[1125]	100	real :: dqchim(ngrid,nlayer,nq) ! tendencies on pq due to chemistry
[1047]	101	real :: dqschim(ngrid,nq) ! tendencies on qsurf
[1125]	102	real :: dqcloud(ngrid,nlayer,nq) ! tendencies on pq due to condensation
[1047]	103	real :: dqscloud(ngrid,nq) ! tendencies on qsurf
[38]	104
[618]	105	! local variables:
[38]	106
[618]	107	integer,save :: nbq ! number of tracers used in the chemistry
[1036]	108	integer,allocatable,save :: niq(:) ! array storing the indexes of the tracers
[1125]	109	integer :: iloc(1) ! index of major species
[618]	110	integer :: ig,l,i,iq,iqmax
	111	integer :: foundswitch, lswitch
[635]	112	integer,save :: chemthermod
[38]	113
[618]	114	integer,save :: i_co2 = 0
	115	integer,save :: i_co = 0
	116	integer,save :: i_o = 0
	117	integer,save :: i_o1d = 0
	118	integer,save :: i_o2 = 0
	119	integer,save :: i_o3 = 0
	120	integer,save :: i_h = 0
	121	integer,save :: i_h2 = 0
	122	integer,save :: i_oh = 0
	123	integer,save :: i_ho2 = 0
	124	integer,save :: i_h2o2 = 0
	125	integer,save :: i_ch4 = 0
	126	integer,save :: i_n2 = 0
	127	integer,save :: i_h2o = 0
[635]	128	integer,save :: i_n = 0
	129	integer,save :: i_no = 0
	130	integer,save :: i_no2 = 0
	131	integer,save :: i_n2d = 0
	132	integer,save :: i_co2plus=0
	133	integer,save :: i_oplus=0
	134	integer,save :: i_o2plus=0
	135	integer,save :: i_coplus=0
	136	integer,save :: i_cplus=0
	137	integer,save :: i_nplus=0
	138	integer,save :: i_noplus=0
	139	integer,save :: i_n2plus=0
	140	integer,save :: i_hplus=0
	141	integer,save :: i_hco2plus=0
	142	integer,save :: i_elec=0
[334]	143
[618]	144	integer :: ig_vl1
[38]	145
[618]	146	real :: latvl1, lonvl1
[1125]	147	real :: zq(ngrid,nlayer,nq) ! pq+pdq*ptimestep before chemistry
	148	! new mole fraction after
[1047]	149	real :: zt(ngrid,nlayer) ! temperature
	150	real :: zu(ngrid,nlayer) ! u component of the wind
	151	real :: zv(ngrid,nlayer) ! v component of the wind
[1125]	152	real :: taucol ! optical depth at 7 hPa
[38]	153
[618]	154	logical,save :: firstcall = .true.
[1125]	155	logical,save :: depos = .false. ! switch for dry deposition
[38]	156
[618]	157	! for each column of atmosphere:
[334]	158
[1125]	159	real :: zpress(nlayer) ! Pressure (mbar)
	160	real :: zdens(nlayer) ! Density (cm-3)
	161	real :: ztemp(nlayer) ! Temperature (K)
	162	real :: zlocal(nlayer) ! Altitude (km)
	163	real :: zycol(nlayer,nq) ! Composition (mole fractions)
	164	real :: szacol ! Solar zenith angle
	165	real :: surfice1d(nlayer) ! Ice surface area (cm2/cm3)
	166	real :: surfdust1d(nlayer) ! Dust surface area (cm2/cm3)
	167	real :: jo3(nlayer) ! Photodissociation rate O3->O1D (s-1)
[1888]	168	real :: em_no(nlayer) ! NO nightglow emission rate
	169	real :: em_o2(nlayer) ! O2 nightglow emission rate
[618]	170
	171	! for output:
	172
[1125]	173	logical :: output ! to issue calls to writediagfi and stats
[618]	174	parameter (output = .true.)
[1047]	175	real :: jo3_3d(ngrid,nlayer) ! Photodissociation rate O3->O1D (s-1)
[1888]	176	real :: emission_no(ngrid,nlayer) !NO emission rate
	177	real :: emission_o2(ngrid,nlayer) !O2 emission rate
[618]	178
	179	!=======================================================================
	180	! initialization of the chemistry (first call only)
	181	!=======================================================================
	182
[38]	183	if (firstcall) then
[618]	184
[38]	185	if (photochem) then
[334]	186	print*,'calchim: Read photolysis lookup table'
	187	call read_phototable
[38]	188	end if
	189	! find index of chemical tracers to use
[1036]	190	allocate(niq(nq))
[635]	191	! Listed here are all tracers that can go into photochemistry
[618]	192	nbq = 0 ! to count number of tracers
[635]	193	! Species ALWAYS present if photochem=.T. or thermochem=.T.
[618]	194	i_co2 = igcm_co2
	195	if (i_co2 == 0) then
	196	write(,) "calchim: Error; no CO2 tracer !!!"
	197	stop
[38]	198	else
[618]	199	nbq = nbq + 1
	200	niq(nbq) = i_co2
	201	end if
	202	i_co = igcm_co
	203	if (i_co == 0) then
	204	write(,) "calchim: Error; no CO tracer !!!"
	205	stop
[38]	206	else
[618]	207	nbq = nbq + 1
	208	niq(nbq) = i_co
	209	end if
	210	i_o = igcm_o
	211	if (i_o == 0) then
	212	write(,) "calchim: Error; no O tracer !!!"
	213	stop
[38]	214	else
[618]	215	nbq = nbq + 1
	216	niq(nbq) = i_o
	217	end if
	218	i_o1d = igcm_o1d
	219	if (i_o1d == 0) then
	220	write(,) "calchim: Error; no O1D tracer !!!"
	221	stop
[38]	222	else
[618]	223	nbq = nbq + 1
	224	niq(nbq) = i_o1d
	225	end if
	226	i_o2 = igcm_o2
	227	if (i_o2 == 0) then
	228	write(,) "calchim: Error; no O2 tracer !!!"
	229	stop
[38]	230	else
[618]	231	nbq = nbq + 1
	232	niq(nbq) = i_o2
	233	end if
	234	i_o3 = igcm_o3
	235	if (i_o3 == 0) then
	236	write(,) "calchim: Error; no O3 tracer !!!"
	237	stop
[38]	238	else
[618]	239	nbq = nbq + 1
	240	niq(nbq) = i_o3
	241	end if
	242	i_h = igcm_h
	243	if (i_h == 0) then
	244	write(,) "calchim: Error; no H tracer !!!"
	245	stop
[38]	246	else
[618]	247	nbq = nbq + 1
	248	niq(nbq) = i_h
	249	end if
	250	i_h2 = igcm_h2
	251	if (i_h2 == 0) then
	252	write(,) "calchim: Error; no H2 tracer !!!"
	253	stop
[38]	254	else
[618]	255	nbq = nbq + 1
	256	niq(nbq) = i_h2
	257	end if
	258	i_oh = igcm_oh
	259	if (i_oh == 0) then
	260	write(,) "calchim: Error; no OH tracer !!!"
	261	stop
[38]	262	else
[618]	263	nbq = nbq + 1
	264	niq(nbq) = i_oh
	265	end if
	266	i_ho2 = igcm_ho2
	267	if (i_ho2 == 0) then
	268	write(,) "calchim: Error; no HO2 tracer !!!"
	269	stop
[38]	270	else
[618]	271	nbq = nbq + 1
	272	niq(nbq) = i_ho2
	273	end if
	274	i_h2o2 = igcm_h2o2
	275	if (i_h2o2 == 0) then
	276	write(,) "calchim: Error; no H2O2 tracer !!!"
	277	stop
[38]	278	else
[618]	279	nbq = nbq + 1
	280	niq(nbq) = i_h2o2
	281	end if
	282	i_ch4 = igcm_ch4
	283	if (i_ch4 == 0) then
[635]	284	write(,) "calchim: Error; no CH4 tracer !!!"
[618]	285	write(,) "CH4 will be ignored in the chemistry"
[334]	286	else
[618]	287	nbq = nbq + 1
	288	niq(nbq) = i_ch4
	289	end if
	290	i_n2 = igcm_n2
	291	if (i_n2 == 0) then
	292	write(,) "calchim: Error; no N2 tracer !!!"
	293	stop
[38]	294	else
[618]	295	nbq = nbq + 1
	296	niq(nbq) = i_n2
	297	end if
	298	i_h2o = igcm_h2o_vap
	299	if (i_h2o == 0) then
	300	write(,) "calchim: Error; no water vapor tracer !!!"
	301	stop
[38]	302	else
[618]	303	nbq = nbq + 1
	304	niq(nbq) = i_h2o
	305	end if
[635]	306	!Check tracers needed for thermospheric chemistry
	307	if(thermochem) then
	308	chemthermod=0 !Default: C/O/H chemistry
	309	!Nitrogen chemistry
	310	!NO is used to determine if N chemistry is wanted
	311	!chemthermod=2 -> N chemistry
	312	i_no = igcm_no
	313	if (i_no == 0) then
	314	write(,) "calchim: no NO tracer"
	315	write(,) "C/O/H themosp chemistry only "
	316	else
	317	nbq = nbq + 1
	318	niq(nbq) = i_no
	319	chemthermod=2
	320	write(,) "calchim: NO in traceur.def"
	321	write(,) "Nitrogen chemistry included"
	322	end if
	323	! N
	324	i_n = igcm_n
	325	if(chemthermod == 2) then
	326	if (i_n == 0) then
	327	write(,) "calchim: Error; no N tracer !!!"
	328	write(,) "N is needed if NO is in traceur.def"
	329	stop
	330	else
	331	nbq = nbq + 1
	332	niq(nbq) = i_n
	333	end if
	334	else
	335	if (i_n /= 0) then
	336	write(,) "calchim: Error: N is present, but NO is not!!!"
	337	write(,) "Both must be in traceur.def if N chemistry is wanted"
	338	stop
	339	endif
	340	endif !Of if(chemthermod == 2)
	341	! NO2
	342	i_no2 = igcm_no2
	343	if(chemthermod == 2) then
	344	if (i_no2 == 0) then
	345	write(,) "calchim: Error; no NO2 tracer !!!"
	346	write(,) "NO2 is needed if NO is in traceur.def"
	347	stop
	348	else
	349	nbq = nbq + 1
	350	niq(nbq) = i_no2
	351	end if
	352	else
	353	if (i_no2 /= 0) then
	354	write(,) "calchim: Error: N is present, but NO is not!!!"
	355	write(,) "Both must be in traceur.def if N chemistry is wanted"
	356	stop
	357	endif
	358	endif !Of if(chemthermod == 2)
	359	! N(2D)
	360	if(chemthermod == 2) then
	361	i_n2d = igcm_n2d
	362	if (i_n2d == 0) then
	363	write(,) "calchim: Error; no N2D !!!"
	364	write(,) "N2D is needed if NO is in traceur.def"
	365	stop
	366	else
	367	nbq = nbq + 1
	368	niq(nbq) = i_n2d
	369	end if
	370	else
	371	if (i_n2d /= 0) then
	372	write(,) "calchim: Error: N2D is present, but NO is not!!!"
	373	write(,) "Both must be in traceur.def if N chemistry wanted"
	374	stop
	375	endif
	376	endif !Of if(chemthermod == 2)
	377	! Ions
	378	! O2+ is used to determine if ion chemistry is needed
	379	! chemthermod=3 -> ion chemistry
	380	i_o2plus = igcm_o2plus
	381	if(chemthermod == 2) then
	382	if (i_o2plus == 0) then
	383	write(,) "calchim: no O2+ tracer; no ion chemistry"
	384	else
	385	nbq = nbq + 1
	386	niq(nbq) = i_o2plus
	387	chemthermod = 3
	388	write(,) "calchim: O2+ in traceur.def"
	389	write(,) "Ion chemistry included"
	390	end if
	391	else
	392	if (i_o2plus /= 0) then
	393	write(,) "calchim: O2+ is present, but NO is not!!!"
	394	write(,) "Both must be in traceur.def if ionosphere wanted"
	395	stop
	396	endif
	397	endif !Of if(chemthermod == 2)
	398	! CO2+
	399	i_co2plus = igcm_co2plus
	400	if(chemthermod == 3) then
	401	if (i_co2plus == 0) then
	402	write(,) "calchim: Error; no CO2+ tracer !!!"
	403	write(,) "CO2+ is needed if O2+ is in traceur.def"
	404	stop
	405	else
	406	nbq = nbq + 1
	407	niq(nbq) = i_co2plus
	408	end if
	409	else
	410	if (i_co2plus /= 0) then
	411	write(,) "calchim: Error: CO2+ is present, but O2+ is not!!!"
	412	write(,) "Both must be in traceur.def if ionosphere wanted"
	413	stop
	414	endif
	415	endif !Of if(chemthermod == 3)
	416	! O+
	417	i_oplus = igcm_oplus
	418	if(chemthermod == 3) then
	419	if (i_oplus == 0) then
	420	write(,) "calchim: Error; no O+ tracer !!!"
	421	write(,) "O+ is needed if O2+ is in traceur.def"
	422	stop
	423	else
	424	nbq = nbq + 1
	425	niq(nbq) = i_oplus
	426	end if
	427	else
	428	if (i_oplus /= 0) then
	429	write(,) "calchim: Error: O+ is present, but O2+ is not!!!"
	430	write(,) "Both must be in traceur.def if ionosphere wanted"
	431	stop
	432	endif
	433	endif !Of if (chemthermod == 3)
	434	! CO+
	435	i_coplus = igcm_coplus
	436	if(chemthermod == 3) then
	437	if (i_coplus == 0) then
	438	write(,) "calchim: Error; no CO+ tracer !!!"
	439	write(,) "CO+ is needed if O2+ is in traceur.def"
	440	stop
	441	else
	442	nbq = nbq + 1
	443	niq(nbq) = i_coplus
	444	end if
	445	else
	446	if (i_coplus /= 0) then
	447	write(,) "calchim: Error: CO+ is present, but O2+ is not!!!"
	448	write(,) " Both must be in traceur.def if ionosphere wanted"
	449	stop
	450	endif
	451	endif ! Of if (chemthermod == 3)
	452	! C+
	453	i_cplus = igcm_cplus
	454	if(chemthermod == 3) then
	455	if (i_cplus == 0) then
	456	write(,) "calchim: Error; no C+ tracer !!!"
	457	write(,) "C+ is needed if O2+ is in traceur.def"
	458	stop
	459	else
	460	nbq = nbq + 1
	461	niq(nbq) = i_cplus
	462	end if
	463	else
	464	if (i_cplus /= 0) then
	465	write(,) "calchim: Error; C+ is present, but O2+ is not!!!"
	466	write(,) "Both must be in traceur.def if ionosphere wanted"
	467	stop
	468	endif
	469	endif ! Of if (chemthermod == 3)
	470	! N+
	471	i_nplus = igcm_nplus
	472	if(chemthermod == 3) then
	473	if (i_nplus == 0) then
	474	write(,) "calchim: Error; no N+ tracer !!!"
	475	write(,) "N+ is needed if O2+ is in traceur.def"
	476	stop
	477	else
	478	nbq = nbq + 1
	479	niq(nbq) = i_nplus
	480	end if
	481	else
	482	if (i_nplus /= 0) then
	483	write(,) "calchim: Error: N+ is present, but O2+ is not!!!"
	484	write(,) "Both must be in traceur.def if ionosphere wanted"
	485	stop
	486	endif
	487	endif !Of if (chemthermod == 3)
	488	! NO+
	489	i_noplus = igcm_noplus
	490	if(chemthermod == 3) then
	491	if (i_noplus == 0) then
	492	write(,) "calchim: Error; no NO+ tracer !!!"
	493	write(,) "NO+ is needed if O2+ is in traceur.def"
	494	stop
	495	else
	496	nbq = nbq + 1
	497	niq(nbq) = i_noplus
	498	end if
	499	else
	500	if (i_noplus /= 0) then
	501	write(,) "calchim: Error: NO+ is present, but O2+ is not!!!"
	502	write(,) "Both must be in traceur.def if ionosphere wanted"
	503	stop
	504	endif
	505	endif !Of if (chemthermod == 3)
	506	! N2+
	507	i_n2plus = igcm_n2plus
	508	if (chemthermod == 3) then
	509	if (i_n2plus == 0) then
	510	write(,) "calchim: Error; no N2+ tracer !!!"
	511	write(,) "N2+ is needed if O2+ is in traceur.def"
	512	stop
	513	else
	514	nbq = nbq + 1
	515	niq(nbq) = i_n2plus
	516	end if
	517	else
	518	if (i_n2plus /= 0) then
	519	write(,) "calchim: Error: N2+ is present, but O2+ is not!!!"
	520	write(,) "Both must be in traceur.def if ionosphere wanted"
	521	stop
	522	endif
	523	endif !Of if (chemthermod == 3)
	524	!H+
	525	i_hplus = igcm_hplus
	526	if (chemthermod == 3) then
	527	if (i_hplus == 0) then
	528	write(,) "calchim: Error; no H+ tracer !!!"
	529	write(,) "H+ is needed if O2+ is in traceur.def"
	530	stop
	531	else
	532	nbq = nbq + 1
	533	niq(nbq) = i_hplus
	534	end if
	535	else
	536	if (i_hplus /= 0) then
	537	write(,) "calchim: Error: H+ is present, but O2+ is not!!!"
	538	write(,) "Both must be in traceur.def if ionosphere wanted"
	539	stop
	540	endif
	541	endif !Of if (chemthermod == 3)
	542	! HCO2+
	543	i_hco2plus = igcm_hco2plus
	544	if(chemthermod == 3) then
	545	if (i_hco2plus == 0) then
	546	write(,) "calchim: Error; no HCO2+ tracer !!!"
	547	write(,) "HCO2+ is needed if O2+ is in traceur.def"
	548	stop
	549	else
	550	nbq = nbq + 1
	551	niq(nbq) = i_hco2plus
	552	end if
	553	else
	554	if (i_hco2plus /= 0) then
	555	write(,) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
	556	write(,) "Both must be in traceur.def if ionosphere wanted"
	557	stop
	558	endif
	559	endif !Of if(chemthermod == 3)
	560	!e-
	561	i_elec = igcm_elec
	562	if(chemthermod == 3) then
	563	if (i_elec == 0) then
	564	write(,) "calchim: Error; no e- tracer !!!"
	565	write(,) "e- is needed if O2+ is in traceur.def"
	566	stop
	567	else
	568	nbq = nbq + 1
	569	niq(nbq) = i_elec
	570	end if
	571	else
	572	if(i_elec /= 0) then
	573	write(,) "calchim: Error: e- is present, but O2+ is not!!!"
	574	write(,) "Both must be in traceur.def if ionosphere wanted"
	575	stop
	576	endif
	577	endif !Of if (chemthermod == 3)
	578	endif !Of thermochem
[38]	579
[334]	580	write(,) 'calchim: found nbq = ',nbq,' tracers'
[635]	581
[38]	582	firstcall = .false.
	583	end if ! if (firstcall)
	584
[618]	585	! Initializations
[38]	586
[618]	587	zycol(:,:) = 0.
[635]	588	dqchim(:,:,:) = 0.
	589	dqschim(:,:) = 0.
[618]	590
[334]	591	! latvl1= 22.27
	592	! lonvl1= -47.94
[618]	593	! ig_vl1= 1+ int( (1.5-(latvl1-90.)jjm/180.) -2 )iim + &
	594	! int(1.5+(lonvl1+180)*iim/360.)
	595
	596	!=======================================================================
	597	! loop over grid
	598	!=======================================================================
	599
[1047]	600	do ig = 1,ngrid
[635]	601
[38]	602	foundswitch = 0
[1047]	603	do l = 1,nlayer
[334]	604	do i = 1,nbq
[618]	605	iq = niq(i) ! get tracer index
	606	zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
	607	zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
[334]	608	end do
[618]	609	zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
	610	zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
	611	zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
[38]	612	zpress(l) = pplay(ig,l)/100.
	613	ztemp(l) = zt(ig,l)
	614	zdens(l) = zpress(l)/(kb1.e4ztemp(l))
	615	zlocal(l) = zzlay(ig,l)/1000.
[618]	616
	617	! surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
	618
[459]	619	surfdust1d(l) = surfdust(ig,l)*1.e-2
	620	surfice1d(l) = surfice(ig,l)*1.e-2
[618]	621
	622	! search for switch index between regions
	623
[38]	624	if (photochem .and. thermochem) then
[635]	625	if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
[38]	626	lswitch = l
[618]	627	foundswitch = 1
[38]	628	end if
	629	end if
[618]	630	if (.not. photochem) then
[38]	631	lswitch = 22
	632	end if
[618]	633	if (.not. thermochem) then
[1047]	634	lswitch = min(50,nlayer+1)
[38]	635	end if
[618]	636
[1047]	637	end do ! of do l=1,nlayer
[618]	638
[38]	639	szacol = acos(mu0(ig))*180./pi
[618]	640	taucol = tauref(ig)*(700./610.) ! provisoire en attente de nouveau jmars
	641
	642	!=======================================================================
	643	! call chemical subroutines
	644	!=======================================================================
	645
	646	! chemistry in lower atmosphere
	647
[334]	648	if (photochem) then
[1125]	649	call photochemistry(nlayer,nq, &
	650	lswitch,zycol,szacol,ptimestep, &
[618]	651	zpress,ztemp,zdens,dist_sol, &
	652	surfdust1d,surfice1d,jo3,taucol)
	653
	654	! ozone photolysis, for output
	655
[1047]	656	do l = 1,nlayer
[618]	657	jo3_3d(ig,l) = jo3(l)
	658	end do
	659
	660	! condensation of h2o2
	661
[1047]	662	call perosat(ngrid, nlayer, nq, &
	663	ig,ptimestep,pplev,pplay, &
[618]	664	ztemp,zycol,dqcloud,dqscloud)
[334]	665	end if
[618]	666
	667	! chemistry in upper atmosphere
[635]	668
[334]	669	if (thermochem) then
[1266]	670	call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
[1888]	671	zdens,zpress,zlocal,szacol,ptimestep,zday,&
	672	em_no,em_o2)
	673	do l=1,nlayer
	674	emission_no(ig,l)=em_no(l)
	675	emission_o2(ig,l)=em_o2(l)
	676	enddo
[334]	677	end if
[618]	678
	679	! dry deposition
	680
[334]	681	if (depos) then
[1047]	682	call deposition(ngrid, nlayer, nq, &
	683	ig, ig_vl1, pplay, pplev, zzlay, zzlev, &
[618]	684	zu, zv, zt, zycol, ptimestep, co2ice)
[334]	685	end if
[618]	686	!=======================================================================
	687	! tendencies
	688	!=======================================================================
	689
	690	! index of the most abundant species at each level
	691
[658]	692	! major(:) = maxloc(zycol, dim = 2)
[618]	693
	694	! tendency for the most abundant species = - sum of others
[1047]	695	do l = 1,nlayer
[658]	696	iloc=maxloc(zycol(l,:))
	697	iqmax=iloc(1)
[618]	698	do i = 1,nbq
	699	iq = niq(i) ! get tracer index
	700	if (iq /= iqmax) then
	701	dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l) &
	702	- zq(ig,l,iq))/ptimestep
	703	dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) &
	704	- dqchim(ig,l,iq)
	705	end if
[38]	706	end do
[1047]	707	end do ! of do l = 1,nlayer
[618]	708
	709	!=======================================================================
	710	! end of loop over grid
	711	!=======================================================================
	712
[1047]	713	end do ! of do ig=1,ngrid
[618]	714
	715	!=======================================================================
	716	! write outputs
	717	!=======================================================================
	718
[38]	719	! value of parameter 'output' to trigger writting of outputs
	720	! is set above at the declaration of the variable.
	721
[618]	722	if (photochem .and. output) then
[1047]	723	if (ngrid > 1) then
	724	call writediagfi(ngrid,'jo3','j o3->o1d', &
[618]	725	's-1',3,jo3_3d(1,1))
[1888]	726	call writediagfi(ngrid,'emission_no', &
	727	'NO nightglow emission rate','cm-3 s-1',3,emission_no)
	728	call writediagfi(ngrid,'emission_o2', &
	729	'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
	730
[476]	731	if (callstats) then
[1047]	732	call wstats(ngrid,'jo3','j o3->o1d', &
[618]	733	's-1',3,jo3_3d(1,1))
[1888]	734	call wstats(ngrid,'emission_no', &
	735	'NO nightglow emission rate','cm-3 s-1',3,emission_no)
	736	call wstats(ngrid,'emission_o2', &
	737	'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
[476]	738	endif
[1047]	739	end if ! of if (ngrid.gt.1)
[618]	740	end if ! of if (output)
	741
	742	return
[635]	743	end
	744

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

Download in other formats:

Original Format