Context Navigation

calchim.F90 @ 937

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

Mars GCM:

updated high atmosphere photochemistry (jthermcalc.F, param_v4.h, iono.h, paramfoto_compact.F, param_read.F , thermosphere.F).
minor change in calchim.F90 (to not use maxloc(zycol, dim = 2) function which seems to be a problem for g95) .
minor bug fix in perosat.F; set tendency on pdqscloud for h2o2 to zero if none is computed.
in "moldiff.F", changed "tridag" to "tridag_sp", "LUBKSB" to "LUBKSB_SP" and "LUDCMP" to "LUDCMP_SP" to avoid possible conflicts with same routines defined in "moldiff_red.F". Added use of automatic-sized array in "tridag" and "tridag_sp" local array "gam" (to avoid using an a priori oversized local array).

FGG+JYC+EM

File size: 25.6 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
[658]	97	integer :: iloc(1) ! index of major species
[618]	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	! find index of chemical tracers to use
[635]	174	! Listed here are all tracers that can go into photochemistry
[618]	175	nbq = 0 ! to count number of tracers
[635]	176	! Species ALWAYS present if photochem=.T. or thermochem=.T.
[618]	177	i_co2 = igcm_co2
	178	if (i_co2 == 0) then
	179	write(,) "calchim: Error; no CO2 tracer !!!"
	180	stop
[38]	181	else
[618]	182	nbq = nbq + 1
	183	niq(nbq) = i_co2
	184	end if
	185	i_co = igcm_co
	186	if (i_co == 0) then
	187	write(,) "calchim: Error; no CO tracer !!!"
	188	stop
[38]	189	else
[618]	190	nbq = nbq + 1
	191	niq(nbq) = i_co
	192	end if
	193	i_o = igcm_o
	194	if (i_o == 0) then
	195	write(,) "calchim: Error; no O tracer !!!"
	196	stop
[38]	197	else
[618]	198	nbq = nbq + 1
	199	niq(nbq) = i_o
	200	end if
	201	i_o1d = igcm_o1d
	202	if (i_o1d == 0) then
	203	write(,) "calchim: Error; no O1D tracer !!!"
	204	stop
[38]	205	else
[618]	206	nbq = nbq + 1
	207	niq(nbq) = i_o1d
	208	end if
	209	i_o2 = igcm_o2
	210	if (i_o2 == 0) then
	211	write(,) "calchim: Error; no O2 tracer !!!"
	212	stop
[38]	213	else
[618]	214	nbq = nbq + 1
	215	niq(nbq) = i_o2
	216	end if
	217	i_o3 = igcm_o3
	218	if (i_o3 == 0) then
	219	write(,) "calchim: Error; no O3 tracer !!!"
	220	stop
[38]	221	else
[618]	222	nbq = nbq + 1
	223	niq(nbq) = i_o3
	224	end if
	225	i_h = igcm_h
	226	if (i_h == 0) then
	227	write(,) "calchim: Error; no H tracer !!!"
	228	stop
[38]	229	else
[618]	230	nbq = nbq + 1
	231	niq(nbq) = i_h
	232	end if
	233	i_h2 = igcm_h2
	234	if (i_h2 == 0) then
	235	write(,) "calchim: Error; no H2 tracer !!!"
	236	stop
[38]	237	else
[618]	238	nbq = nbq + 1
	239	niq(nbq) = i_h2
	240	end if
	241	i_oh = igcm_oh
	242	if (i_oh == 0) then
	243	write(,) "calchim: Error; no OH tracer !!!"
	244	stop
[38]	245	else
[618]	246	nbq = nbq + 1
	247	niq(nbq) = i_oh
	248	end if
	249	i_ho2 = igcm_ho2
	250	if (i_ho2 == 0) then
	251	write(,) "calchim: Error; no HO2 tracer !!!"
	252	stop
[38]	253	else
[618]	254	nbq = nbq + 1
	255	niq(nbq) = i_ho2
	256	end if
	257	i_h2o2 = igcm_h2o2
	258	if (i_h2o2 == 0) then
	259	write(,) "calchim: Error; no H2O2 tracer !!!"
	260	stop
[38]	261	else
[618]	262	nbq = nbq + 1
	263	niq(nbq) = i_h2o2
	264	end if
	265	i_ch4 = igcm_ch4
	266	if (i_ch4 == 0) then
[635]	267	write(,) "calchim: Error; no CH4 tracer !!!"
[618]	268	write(,) "CH4 will be ignored in the chemistry"
[334]	269	else
[618]	270	nbq = nbq + 1
	271	niq(nbq) = i_ch4
	272	end if
	273	i_n2 = igcm_n2
	274	if (i_n2 == 0) then
	275	write(,) "calchim: Error; no N2 tracer !!!"
	276	stop
[38]	277	else
[618]	278	nbq = nbq + 1
	279	niq(nbq) = i_n2
	280	end if
	281	i_h2o = igcm_h2o_vap
	282	if (i_h2o == 0) then
	283	write(,) "calchim: Error; no water vapor tracer !!!"
	284	stop
[38]	285	else
[618]	286	nbq = nbq + 1
	287	niq(nbq) = i_h2o
	288	end if
[635]	289	!Check tracers needed for thermospheric chemistry
	290	if(thermochem) then
	291	chemthermod=0 !Default: C/O/H chemistry
	292	!Nitrogen chemistry
	293	!NO is used to determine if N chemistry is wanted
	294	!chemthermod=2 -> N chemistry
	295	i_no = igcm_no
	296	if (i_no == 0) then
	297	write(,) "calchim: no NO tracer"
	298	write(,) "C/O/H themosp chemistry only "
	299	else
	300	nbq = nbq + 1
	301	niq(nbq) = i_no
	302	chemthermod=2
	303	write(,) "calchim: NO in traceur.def"
	304	write(,) "Nitrogen chemistry included"
	305	end if
	306	! N
	307	i_n = igcm_n
	308	if(chemthermod == 2) then
	309	if (i_n == 0) then
	310	write(,) "calchim: Error; no N tracer !!!"
	311	write(,) "N is needed if NO is in traceur.def"
	312	stop
	313	else
	314	nbq = nbq + 1
	315	niq(nbq) = i_n
	316	end if
	317	else
	318	if (i_n /= 0) then
	319	write(,) "calchim: Error: N is present, but NO is not!!!"
	320	write(,) "Both must be in traceur.def if N chemistry is wanted"
	321	stop
	322	endif
	323	endif !Of if(chemthermod == 2)
	324	! NO2
	325	i_no2 = igcm_no2
	326	if(chemthermod == 2) then
	327	if (i_no2 == 0) then
	328	write(,) "calchim: Error; no NO2 tracer !!!"
	329	write(,) "NO2 is needed if NO is in traceur.def"
	330	stop
	331	else
	332	nbq = nbq + 1
	333	niq(nbq) = i_no2
	334	end if
	335	else
	336	if (i_no2 /= 0) then
	337	write(,) "calchim: Error: N is present, but NO is not!!!"
	338	write(,) "Both must be in traceur.def if N chemistry is wanted"
	339	stop
	340	endif
	341	endif !Of if(chemthermod == 2)
	342	! N(2D)
	343	if(chemthermod == 2) then
	344	i_n2d = igcm_n2d
	345	if (i_n2d == 0) then
	346	write(,) "calchim: Error; no N2D !!!"
	347	write(,) "N2D is needed if NO is in traceur.def"
	348	stop
	349	else
	350	nbq = nbq + 1
	351	niq(nbq) = i_n2d
	352	end if
	353	else
	354	if (i_n2d /= 0) then
	355	write(,) "calchim: Error: N2D is present, but NO is not!!!"
	356	write(,) "Both must be in traceur.def if N chemistry wanted"
	357	stop
	358	endif
	359	endif !Of if(chemthermod == 2)
	360	! Ions
	361	! O2+ is used to determine if ion chemistry is needed
	362	! chemthermod=3 -> ion chemistry
	363	i_o2plus = igcm_o2plus
	364	if(chemthermod == 2) then
	365	if (i_o2plus == 0) then
	366	write(,) "calchim: no O2+ tracer; no ion chemistry"
	367	else
	368	nbq = nbq + 1
	369	niq(nbq) = i_o2plus
	370	chemthermod = 3
	371	write(,) "calchim: O2+ in traceur.def"
	372	write(,) "Ion chemistry included"
	373	end if
	374	else
	375	if (i_o2plus /= 0) then
	376	write(,) "calchim: O2+ is present, but NO is not!!!"
	377	write(,) "Both must be in traceur.def if ionosphere wanted"
	378	stop
	379	endif
	380	endif !Of if(chemthermod == 2)
	381	! CO2+
	382	i_co2plus = igcm_co2plus
	383	if(chemthermod == 3) then
	384	if (i_co2plus == 0) then
	385	write(,) "calchim: Error; no CO2+ tracer !!!"
	386	write(,) "CO2+ is needed if O2+ is in traceur.def"
	387	stop
	388	else
	389	nbq = nbq + 1
	390	niq(nbq) = i_co2plus
	391	end if
	392	else
	393	if (i_co2plus /= 0) then
	394	write(,) "calchim: Error: CO2+ is present, but O2+ is not!!!"
	395	write(,) "Both must be in traceur.def if ionosphere wanted"
	396	stop
	397	endif
	398	endif !Of if(chemthermod == 3)
	399	! O+
	400	i_oplus = igcm_oplus
	401	if(chemthermod == 3) then
	402	if (i_oplus == 0) then
	403	write(,) "calchim: Error; no O+ tracer !!!"
	404	write(,) "O+ is needed if O2+ is in traceur.def"
	405	stop
	406	else
	407	nbq = nbq + 1
	408	niq(nbq) = i_oplus
	409	end if
	410	else
	411	if (i_oplus /= 0) then
	412	write(,) "calchim: Error: O+ is present, but O2+ is not!!!"
	413	write(,) "Both must be in traceur.def if ionosphere wanted"
	414	stop
	415	endif
	416	endif !Of if (chemthermod == 3)
	417	! CO+
	418	i_coplus = igcm_coplus
	419	if(chemthermod == 3) then
	420	if (i_coplus == 0) then
	421	write(,) "calchim: Error; no CO+ tracer !!!"
	422	write(,) "CO+ is needed if O2+ is in traceur.def"
	423	stop
	424	else
	425	nbq = nbq + 1
	426	niq(nbq) = i_coplus
	427	end if
	428	else
	429	if (i_coplus /= 0) then
	430	write(,) "calchim: Error: CO+ is present, but O2+ is not!!!"
	431	write(,) " Both must be in traceur.def if ionosphere wanted"
	432	stop
	433	endif
	434	endif ! Of if (chemthermod == 3)
	435	! C+
	436	i_cplus = igcm_cplus
	437	if(chemthermod == 3) then
	438	if (i_cplus == 0) then
	439	write(,) "calchim: Error; no C+ tracer !!!"
	440	write(,) "C+ is needed if O2+ is in traceur.def"
	441	stop
	442	else
	443	nbq = nbq + 1
	444	niq(nbq) = i_cplus
	445	end if
	446	else
	447	if (i_cplus /= 0) then
	448	write(,) "calchim: Error; C+ is present, but O2+ is not!!!"
	449	write(,) "Both must be in traceur.def if ionosphere wanted"
	450	stop
	451	endif
	452	endif ! Of if (chemthermod == 3)
	453	! N+
	454	i_nplus = igcm_nplus
	455	if(chemthermod == 3) then
	456	if (i_nplus == 0) then
	457	write(,) "calchim: Error; no N+ tracer !!!"
	458	write(,) "N+ is needed if O2+ is in traceur.def"
	459	stop
	460	else
	461	nbq = nbq + 1
	462	niq(nbq) = i_nplus
	463	end if
	464	else
	465	if (i_nplus /= 0) then
	466	write(,) "calchim: Error: N+ is present, but O2+ is not!!!"
	467	write(,) "Both must be in traceur.def if ionosphere wanted"
	468	stop
	469	endif
	470	endif !Of if (chemthermod == 3)
	471	! NO+
	472	i_noplus = igcm_noplus
	473	if(chemthermod == 3) then
	474	if (i_noplus == 0) then
	475	write(,) "calchim: Error; no NO+ tracer !!!"
	476	write(,) "NO+ is needed if O2+ is in traceur.def"
	477	stop
	478	else
	479	nbq = nbq + 1
	480	niq(nbq) = i_noplus
	481	end if
	482	else
	483	if (i_noplus /= 0) then
	484	write(,) "calchim: Error: NO+ is present, but O2+ is not!!!"
	485	write(,) "Both must be in traceur.def if ionosphere wanted"
	486	stop
	487	endif
	488	endif !Of if (chemthermod == 3)
	489	! N2+
	490	i_n2plus = igcm_n2plus
	491	if (chemthermod == 3) then
	492	if (i_n2plus == 0) then
	493	write(,) "calchim: Error; no N2+ tracer !!!"
	494	write(,) "N2+ is needed if O2+ is in traceur.def"
	495	stop
	496	else
	497	nbq = nbq + 1
	498	niq(nbq) = i_n2plus
	499	end if
	500	else
	501	if (i_n2plus /= 0) then
	502	write(,) "calchim: Error: N2+ is present, but O2+ is not!!!"
	503	write(,) "Both must be in traceur.def if ionosphere wanted"
	504	stop
	505	endif
	506	endif !Of if (chemthermod == 3)
	507	!H+
	508	i_hplus = igcm_hplus
	509	if (chemthermod == 3) then
	510	if (i_hplus == 0) then
	511	write(,) "calchim: Error; no H+ tracer !!!"
	512	write(,) "H+ is needed if O2+ is in traceur.def"
	513	stop
	514	else
	515	nbq = nbq + 1
	516	niq(nbq) = i_hplus
	517	end if
	518	else
	519	if (i_hplus /= 0) then
	520	write(,) "calchim: Error: H+ is present, but O2+ is not!!!"
	521	write(,) "Both must be in traceur.def if ionosphere wanted"
	522	stop
	523	endif
	524	endif !Of if (chemthermod == 3)
	525	! HCO2+
	526	i_hco2plus = igcm_hco2plus
	527	if(chemthermod == 3) then
	528	if (i_hco2plus == 0) then
	529	write(,) "calchim: Error; no HCO2+ tracer !!!"
	530	write(,) "HCO2+ is needed if O2+ is in traceur.def"
	531	stop
	532	else
	533	nbq = nbq + 1
	534	niq(nbq) = i_hco2plus
	535	end if
	536	else
	537	if (i_hco2plus /= 0) then
	538	write(,) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
	539	write(,) "Both must be in traceur.def if ionosphere wanted"
	540	stop
	541	endif
	542	endif !Of if(chemthermod == 3)
	543	!e-
	544	i_elec = igcm_elec
	545	if(chemthermod == 3) then
	546	if (i_elec == 0) then
	547	write(,) "calchim: Error; no e- tracer !!!"
	548	write(,) "e- is needed if O2+ is in traceur.def"
	549	stop
	550	else
	551	nbq = nbq + 1
	552	niq(nbq) = i_elec
	553	end if
	554	else
	555	if(i_elec /= 0) then
	556	write(,) "calchim: Error: e- is present, but O2+ is not!!!"
	557	write(,) "Both must be in traceur.def if ionosphere wanted"
	558	stop
	559	endif
	560	endif !Of if (chemthermod == 3)
	561	endif !Of thermochem
[38]	562
[334]	563	write(,) 'calchim: found nbq = ',nbq,' tracers'
[635]	564
[38]	565	firstcall = .false.
	566	end if ! if (firstcall)
	567
[618]	568	! Initializations
[38]	569
[618]	570	zycol(:,:) = 0.
[635]	571	dqchim(:,:,:) = 0.
	572	dqschim(:,:) = 0.
[618]	573
[334]	574	! latvl1= 22.27
	575	! lonvl1= -47.94
[618]	576	! ig_vl1= 1+ int( (1.5-(latvl1-90.)jjm/180.) -2 )iim + &
	577	! int(1.5+(lonvl1+180)*iim/360.)
	578
	579	!=======================================================================
	580	! loop over grid
	581	!=======================================================================
	582
[38]	583	do ig = 1,ngridmx
[635]	584
[38]	585	foundswitch = 0
	586	do l = 1,nlayermx
[334]	587	do i = 1,nbq
[618]	588	iq = niq(i) ! get tracer index
	589	zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
	590	zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
[334]	591	end do
[618]	592	zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
	593	zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
	594	zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
[38]	595	zpress(l) = pplay(ig,l)/100.
	596	ztemp(l) = zt(ig,l)
	597	zdens(l) = zpress(l)/(kb1.e4ztemp(l))
	598	zlocal(l) = zzlay(ig,l)/1000.
[618]	599
	600	! surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
	601
[459]	602	surfdust1d(l) = surfdust(ig,l)*1.e-2
	603	surfice1d(l) = surfice(ig,l)*1.e-2
[618]	604
	605	! search for switch index between regions
	606
[38]	607	if (photochem .and. thermochem) then
[635]	608	if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
[38]	609	lswitch = l
[618]	610	foundswitch = 1
[38]	611	end if
	612	end if
[618]	613	if (.not. photochem) then
[38]	614	lswitch = 22
	615	end if
[618]	616	if (.not. thermochem) then
	617	lswitch = min(50,nlayermx+1)
[38]	618	end if
[618]	619
[38]	620	end do ! of do l=1,nlayermx
[618]	621
[38]	622	szacol = acos(mu0(ig))*180./pi
[618]	623	taucol = tauref(ig)*(700./610.) ! provisoire en attente de nouveau jmars
	624
	625	!=======================================================================
	626	! call chemical subroutines
	627	!=======================================================================
	628
	629	! chemistry in lower atmosphere
	630
[334]	631	if (photochem) then
[618]	632	call photochemistry(lswitch,zycol,szacol,ptimestep, &
	633	zpress,ztemp,zdens,dist_sol, &
	634	surfdust1d,surfice1d,jo3,taucol)
	635
	636	! ozone photolysis, for output
	637
	638	do l = 1,nlayermx
	639	jo3_3d(ig,l) = jo3(l)
	640	end do
	641
	642	! condensation of h2o2
	643
	644	call perosat(ig,ptimestep,pplev,pplay, &
	645	ztemp,zycol,dqcloud,dqscloud)
[334]	646	end if
[618]	647
	648	! chemistry in upper atmosphere
[635]	649
[334]	650	if (thermochem) then
[635]	651	call chemthermos(ig,lswitch,chemthermod,zycol,ztemp,zdens, &
	652	zpress,zlocal,szacol,ptimestep,zday)
[334]	653	end if
[618]	654
	655	! dry deposition
	656
[334]	657	if (depos) then
[618]	658	call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev,&
	659	zu, zv, zt, zycol, ptimestep, co2ice)
[334]	660	end if
[618]	661	!=======================================================================
	662	! tendencies
	663	!=======================================================================
	664
	665	! index of the most abundant species at each level
	666
[658]	667	! major(:) = maxloc(zycol, dim = 2)
[618]	668
	669	! tendency for the most abundant species = - sum of others
	670	do l = 1,nlayermx
[658]	671	iloc=maxloc(zycol(l,:))
	672	iqmax=iloc(1)
[618]	673	do i = 1,nbq
	674	iq = niq(i) ! get tracer index
	675	if (iq /= iqmax) then
	676	dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l) &
	677	- zq(ig,l,iq))/ptimestep
	678	dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) &
	679	- dqchim(ig,l,iq)
	680	end if
[38]	681	end do
[618]	682	end do ! of do l = 1,nlayermx
	683
	684	!=======================================================================
	685	! end of loop over grid
	686	!=======================================================================
	687
[38]	688	end do ! of do ig=1,ngridmx
[618]	689
	690	!=======================================================================
	691	! write outputs
	692	!=======================================================================
	693
[38]	694	! value of parameter 'output' to trigger writting of outputs
	695	! is set above at the declaration of the variable.
	696
[618]	697	if (photochem .and. output) then
	698	if (ngridmx > 1) then
	699	call writediagfi(ngridmx,'jo3','j o3->o1d', &
	700	's-1',3,jo3_3d(1,1))
[476]	701	if (callstats) then
[618]	702	call wstats(ngridmx,'jo3','j o3->o1d', &
	703	's-1',3,jo3_3d(1,1))
[476]	704	endif
[38]	705	end if ! of if (ngridmx.gt.1)
[618]	706	end if ! of if (output)
	707
	708	return
[635]	709	end
	710

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

Download in other formats:

Original Format