Context Navigation

calchim.F90 @ 1036

Last change on this file since 1036 was 1036, checked in by emillour, 11 years ago

Mars GCM: (a first step towards using parallel dynamics)

IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
- advtrac.h (in dyn3d) removed and replaced by module infotrac.F
- tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

File size: 26.3 KB

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

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

Download in other formats:

Original Format