Context Navigation

calchim.F90 @ 1448

Last change on this file since 1448 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

File size: 26.5 KB

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

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

Download in other formats:

Original Format