Context Navigation

calchim_asis.F90 @ 1888

Last change on this file since 1888 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.8 KB

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

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

Download in other formats:

Original Format