Context Navigation

calchim.F90 @ 629

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

Mars GCM:

Update of the chemistry package:

calchim.F90 :
- upgraded from calchim.F
- can run with or without CH4
- fixed the mass conservation scheme
photochemistry.F : -chemistry timestep is now independant from physics timestep -can run with or without a CH4 tracer
- removed initial tests on species, since these are already done in calchim
removed inichim_readcallphys.F (not used any more)
concentrations.F: adaptated to better handle indexes and molar masses of

tracers. "ncomp" (in chimiedata.h) no longer needed.

inichim_newstart.F90 : rewritten and cleaned (won't be compatible with

old start files where tracer names are q01,...).
Now handles hybrid levels and automaticaly adapts
depending on which tracers are available.

newstart.F : adapted to followup changes in inchim_newstart.F90 and some

cleanup around the initialization of tracer names and surface
values.

FL+EM

File size: 14.5 KB

Line
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
6	implicit none
7
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	!=======================================================================
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
63	! input:
64
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)
85
86	! output:
87
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
92
93	! local variables:
94
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
100
101	integer,save :: i_co2 = 0
102	integer,save :: i_co = 0
103	integer,save :: i_o = 0
104	integer,save :: i_o1d = 0
105	integer,save :: i_o2 = 0
106	integer,save :: i_o3 = 0
107	integer,save :: i_h = 0
108	integer,save :: i_h2 = 0
109	integer,save :: i_oh = 0
110	integer,save :: i_ho2 = 0
111	integer,save :: i_h2o2 = 0
112	integer,save :: i_ch4 = 0
113	integer,save :: i_n2 = 0
114	integer,save :: i_h2o = 0
115
116	integer :: ig_vl1
117
118	real :: latvl1, lonvl1
119	real :: zq(ngridmx,nlayermx,nqmx) ! pq+pdq*ptimestep before chemistry
120	! new mole fraction after
121	real :: zt(ngridmx,nlayermx) ! temperature
122	real :: zu(ngridmx,nlayermx) ! u component of the wind
123	real :: zv(ngridmx,nlayermx) ! v component of the wind
124	real :: taucol ! optical depth at 7 hPa
125
126	logical,save :: firstcall = .true.
127	logical,save :: depos = .false. ! switch for dry deposition
128
129	! for each column of atmosphere:
130
131	real :: zpress(nlayermx) ! Pressure (mbar)
132	real :: zdens(nlayermx) ! Density (cm-3)
133	real :: ztemp(nlayermx) ! Temperature (K)
134	real :: zlocal(nlayermx) ! Altitude (km)
135	real :: zycol(nlayermx,nqmx) ! Composition (mole fractions)
136	real :: szacol ! Solar zenith angle
137	real :: surfice1d(nlayermx) ! Ice surface area (cm2/cm3)
138	real :: surfdust1d(nlayermx) ! Dust surface area (cm2/cm3)
139	real :: jo3(nlayermx) ! Photodissociation rate O3->O1D (s-1)
140
141	! for output:
142
143	logical :: output ! to issue calls to writediagfi and stats
144	parameter (output = .true.)
145	real :: jo3_3d(ngridmx,nlayermx) ! Photodissociation rate O3->O1D (s-1)
146
147	!=======================================================================
148	! initialization of the chemistry (first call only)
149	!=======================================================================
150
151	if (firstcall) then
152
153	if (photochem) then
154	print*,'calchim: Read photolysis lookup table'
155	call read_phototable
156	end if
157
158	! find index of chemical tracers to use
159	nbq = 0 ! to count number of tracers
160	i_co2 = igcm_co2
161	if (i_co2 == 0) then
162	write(,) "calchim: Error; no CO2 tracer !!!"
163	stop
164	else
165	nbq = nbq + 1
166	niq(nbq) = i_co2
167	end if
168	i_co = igcm_co
169	if (i_co == 0) then
170	write(,) "calchim: Error; no CO tracer !!!"
171	stop
172	else
173	nbq = nbq + 1
174	niq(nbq) = i_co
175	end if
176	i_o = igcm_o
177	if (i_o == 0) then
178	write(,) "calchim: Error; no O tracer !!!"
179	stop
180	else
181	nbq = nbq + 1
182	niq(nbq) = i_o
183	end if
184	i_o1d = igcm_o1d
185	if (i_o1d == 0) then
186	write(,) "calchim: Error; no O1D tracer !!!"
187	stop
188	else
189	nbq = nbq + 1
190	niq(nbq) = i_o1d
191	end if
192	i_o2 = igcm_o2
193	if (i_o2 == 0) then
194	write(,) "calchim: Error; no O2 tracer !!!"
195	stop
196	else
197	nbq = nbq + 1
198	niq(nbq) = i_o2
199	end if
200	i_o3 = igcm_o3
201	if (i_o3 == 0) then
202	write(,) "calchim: Error; no O3 tracer !!!"
203	stop
204	else
205	nbq = nbq + 1
206	niq(nbq) = i_o3
207	end if
208	i_h = igcm_h
209	if (i_h == 0) then
210	write(,) "calchim: Error; no H tracer !!!"
211	stop
212	else
213	nbq = nbq + 1
214	niq(nbq) = i_h
215	end if
216	i_h2 = igcm_h2
217	if (i_h2 == 0) then
218	write(,) "calchim: Error; no H2 tracer !!!"
219	stop
220	else
221	nbq = nbq + 1
222	niq(nbq) = i_h2
223	end if
224	i_oh = igcm_oh
225	if (i_oh == 0) then
226	write(,) "calchim: Error; no OH tracer !!!"
227	stop
228	else
229	nbq = nbq + 1
230	niq(nbq) = i_oh
231	end if
232	i_ho2 = igcm_ho2
233	if (i_ho2 == 0) then
234	write(,) "calchim: Error; no HO2 tracer !!!"
235	stop
236	else
237	nbq = nbq + 1
238	niq(nbq) = i_ho2
239	end if
240	i_h2o2 = igcm_h2o2
241	if (i_h2o2 == 0) then
242	write(,) "calchim: Error; no H2O2 tracer !!!"
243	stop
244	else
245	nbq = nbq + 1
246	niq(nbq) = i_h2o2
247	end if
248	i_ch4 = igcm_ch4
249	if (i_ch4 == 0) then
250	write(,) "calchim: no CH4 tracer !!!"
251	write(,) "CH4 will be ignored in the chemistry"
252	else
253	nbq = nbq + 1
254	niq(nbq) = i_ch4
255	end if
256	i_n2 = igcm_n2
257	if (i_n2 == 0) then
258	write(,) "calchim: Error; no N2 tracer !!!"
259	stop
260	else
261	nbq = nbq + 1
262	niq(nbq) = i_n2
263	end if
264	i_h2o = igcm_h2o_vap
265	if (i_h2o == 0) then
266	write(,) "calchim: Error; no water vapor tracer !!!"
267	stop
268	else
269	nbq = nbq + 1
270	niq(nbq) = i_h2o
271	end if
272
273	write(,) 'calchim: found nbq = ',nbq,' tracers'
274
275	firstcall = .false.
276	end if ! if (firstcall)
277
278	! Initializations
279
280	zycol(:,:) = 0.
281	dqchim(:,:,:) = 0
282	dqschim(:,:) = 0
283
284	! latvl1= 22.27
285	! lonvl1= -47.94
286	! ig_vl1= 1+ int( (1.5-(latvl1-90.)jjm/180.) -2 )iim + &
287	! int(1.5+(lonvl1+180)*iim/360.)
288
289	!=======================================================================
290	! loop over grid
291	!=======================================================================
292
293	do ig = 1,ngridmx
294
295	foundswitch = 0
296	do l = 1,nlayermx
297	do i = 1,nbq
298	iq = niq(i) ! get tracer index
299	zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
300	zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
301	end do
302	zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
303	zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
304	zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
305	zpress(l) = pplay(ig,l)/100.
306	ztemp(l) = zt(ig,l)
307	zdens(l) = zpress(l)/(kb1.e4ztemp(l))
308	zlocal(l) = zzlay(ig,l)/1000.
309
310	! surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
311
312	surfdust1d(l) = surfdust(ig,l)*1.e-2
313	surfice1d(l) = surfice(ig,l)*1.e-2
314
315	! search for switch index between regions
316
317	if (photochem .and. thermochem) then
318	if (foundswitch == 0 .and. pplay(ig,l) < 1.e-3) then
319	lswitch = l
320	foundswitch = 1
321	end if
322	end if
323	if (.not. photochem) then
324	lswitch = 22
325	end if
326	if (.not. thermochem) then
327	lswitch = min(50,nlayermx+1)
328	end if
329
330	end do ! of do l=1,nlayermx
331
332	szacol = acos(mu0(ig))*180./pi
333	taucol = tauref(ig)*(700./610.) ! provisoire en attente de nouveau jmars
334
335	!=======================================================================
336	! call chemical subroutines
337	!=======================================================================
338
339	! chemistry in lower atmosphere
340
341	if (photochem) then
342	call photochemistry(lswitch,zycol,szacol,ptimestep, &
343	zpress,ztemp,zdens,dist_sol, &
344	surfdust1d,surfice1d,jo3,taucol)
345
346	! ozone photolysis, for output
347
348	do l = 1,nlayermx
349	jo3_3d(ig,l) = jo3(l)
350	end do
351
352	! condensation of h2o2
353
354	call perosat(ig,ptimestep,pplev,pplay, &
355	ztemp,zycol,dqcloud,dqscloud)
356	end if
357
358	! chemistry in upper atmosphere
359
360	if (thermochem) then
361	call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress, &
362	zlocal,szacol,ptimestep,zday)
363	end if
364
365	! dry deposition
366
367	if (depos) then
368	call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev,&
369	zu, zv, zt, zycol, ptimestep, co2ice)
370	end if
371
372	!=======================================================================
373	! tendencies
374	!=======================================================================
375
376	! index of the most abundant species at each level
377
378	major(:) = maxloc(zycol, dim = 2)
379
380	! tendency for the most abundant species = - sum of others
381
382	do l = 1,nlayermx
383	iqmax = major(l)
384	do i = 1,nbq
385	iq = niq(i) ! get tracer index
386	if (iq /= iqmax) then
387	dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l) &
388	- zq(ig,l,iq))/ptimestep
389	dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) &
390	- dqchim(ig,l,iq)
391	end if
392	end do
393	end do ! of do l = 1,nlayermx
394
395	!=======================================================================
396	! end of loop over grid
397	!=======================================================================
398
399	end do ! of do ig=1,ngridmx
400
401	!=======================================================================
402	! write outputs
403	!=======================================================================
404
405	! value of parameter 'output' to trigger writting of outputs
406	! is set above at the declaration of the variable.
407
408	if (photochem .and. output) then
409	if (ngridmx > 1) then
410	call writediagfi(ngridmx,'jo3','j o3->o1d', &
411	's-1',3,jo3_3d(1,1))
412	if (callstats) then
413	call wstats(ngridmx,'jo3','j o3->o1d', &
414	's-1',3,jo3_3d(1,1))
415	endif
416	end if ! of if (ngridmx.gt.1)
417	end if ! of if (output)
418
419	return
420	end

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

Download in other formats:

Original Format