Context Navigation

calchim.F90 @ 635

Last change on this file since 635 was 635, checked in by emillour, 13 years ago

Mars GCM: Update of the chemistry package, including:

93 reactions are accounted for (instead of 22); tracking 28 species (instead of 11)
computation of photoabsorption using raytracing
improved time stepping in the photochemistry
updated parameters (cross-sections); with this new version input files

are in 'EUV/param_v5' of "datafile" directory.

transition between lower and upper atmosphere chemistry set to 0.1 Pa (calchim.F90)
Lots of code clean-up: removed obsolete files column.F, param_v3.h, flujo.F, phdisrate.F, ch.F, interpfast.F, paramfoto.F, getch.F Converted chemtermos.F -> chemthermos.F90 and euvheat.F -> euvheat.F90. Added paramfoto_compact.F , param_v4.h and iono.h
Upadted surfacearea.F
Cleaned initracer.F and callkeys.h (removed use of obsolete "nqchem" and "oldnames" case when initializing tracers).
Minor correction in "callsedim": compute "rdust" and/or "rice" only when it makes sense.

FGG+FL+EM

File size: 25.7 KB

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

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

Download in other formats:

Original Format