Context Navigation

euvheat.F90 @ 937

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

Mars GCM:

Added possibility to run with a varying EUV cycle following real one. The flag solvarmod=1 triggers this behaviour, with companion flag solvaryear=## , where ## is the Mars Year (from 23 to 30). Setting solvarmod=0 reverts to 'old' behaviour, where there is a constant EUV forcing throughout the run, set by the "solarcondate" flag.
Needs corresponding input data files ("param_v6" subdirectory of "EUV" subdirectory in "datadir").
Added files jthermcalc_e107.F and param_read_e107.F in "aeronomars", modified files euvheat.F90, hrtherm.F, chemthermos.F90, param_v4.h and param_read.F in "aeronomars" and files inifis.F, physiq.F and callkeys.h in "phymars".

FGG

File size: 14.2 KB

Line
1	SUBROUTINE euvheat(pt,pdt,pplev,pplay,zzlay, &
2	mu0,ptimestep,ptime,zday,pq,pdq,pdteuv)
3
4	IMPLICIT NONE
5	!=======================================================================
6	! subject:
7	! --------
8	! Computing heating rate due to EUV absorption
9	!
10	! author: MAC 2002
11	! ------
12	!
13	! input:
14	! -----
15	! mu0(ngridmx)
16	! pplay(ngrid,nlayer) pressure at middle of layers (Pa)
17	!
18	! output:
19	! -------
20	!
21	! pdteuv(ngrid,nlayer) Heating rate (K/s)
22	!
23	!=======================================================================
24	!
25	! 0. Declarations :
26	! ------------------
27	!
28	#include "dimensions.h"
29	#include "dimphys.h"
30	#include "comcstfi.h"
31	#include "callkeys.h"
32	#include "comdiurn.h"
33	#include "param.h"
34	#include "param_v4.h"
35	#include "chimiedata.h"
36	#include "tracer.h"
37	#include "conc.h"
38	!-----------------------------------------------------------------------
39	! Input/Output
40	! ------------
41
42	real :: pt(ngridmx,nlayermx)
43	real :: pdt(ngridmx,nlayermx)
44	real :: pplev(ngridmx,nlayermx+1)
45	real :: pplay(ngridmx,nlayermx)
46	real :: zzlay(ngridmx,nlayermx)
47	real :: mu0(ngridmx)
48	real :: ptimestep,ptime
49	real :: zday
50	real :: pq(ngridmx,nlayermx,nqmx)
51	real :: pdq(ngridmx,nlayermx,nqmx)
52
53	real :: pdteuv(ngridmx,nlayermx)
54	!
55	! Local variables :
56	! -----------------
57	integer,save :: nespeuv ! Number of species considered
58
59	INTEGER :: l,ig,n
60	integer,save :: euvmod
61	real, allocatable :: rm(:,:) ! number density (cm-3)
62	real :: zq(ngridmx,nlayermx,nqmx) ! local updated tracer quantity
63	real :: zt(ngridmx,nlayermx) ! local updated atmospheric temperature
64	real :: zlocal(nlayermx)
65	real :: zenit
66	real :: jtot(nlayermx)
67	real :: dens ! amu/cm-3
68	real :: tx(nlayermx)
69	! real euveff !UV heating efficiency
70
71	! tracer indexes for the EUV heating:
72	!!! ATTENTION. These values have to be identical to those in chemthermos.F90
73	!!! If the values are changed there, the same has to be done here !!!
74	integer,parameter :: i_co2=1
75	integer,parameter :: i_o2=2
76	integer,parameter :: i_o=3
77	integer,parameter :: i_co=4
78	integer,parameter :: i_h=5
79	integer,parameter :: i_oh=6
80	integer,parameter :: i_ho2=7
81	integer,parameter :: i_h2=8
82	integer,parameter :: i_h2o=9
83	integer,parameter :: i_h2o2=10
84	integer,parameter :: i_o1d=11
85	integer,parameter :: i_o3=12
86	integer,parameter :: i_n2=13
87	integer,parameter :: i_n=14
88	integer,parameter :: i_no=15
89	integer,parameter :: i_n2d=16
90	integer,parameter :: i_no2=17
91
92
93	! Tracer indexes in the GCM:
94	integer,save :: g_co2=0
95	integer,save :: g_o=0
96	integer,save :: g_o2=0
97	integer,save :: g_h2=0
98	integer,save :: g_h2o2=0
99	integer,save :: g_h2o=0
100	integer,save :: g_o3=0
101	integer,save :: g_n2=0
102	integer,save :: g_n=0
103	integer,save :: g_no=0
104	integer,save :: g_co=0
105	integer,save :: g_h=0
106	integer,save :: g_no2=0
107	integer,save :: g_oh=0
108	integer,save :: g_ho2=0
109	integer,save :: g_o1d=0
110	integer,save :: g_n2d=0
111
112
113	logical,save :: firstcall=.true.
114
115	! Initializations and sanity checks:
116
117
118	if (firstcall) then
119	nespeuv=0
120	! identify the indexes of the tracers we'll need
121	g_co2=igcm_co2
122	if (g_co2.eq.0) then
123	write(,) "euvheat: Error; no CO2 tracer !!!"
124	write(,) "CO2 is always needed if calleuv=.true."
125	stop
126	else
127	nespeuv=nespeuv+1
128	endif
129	g_o=igcm_o
130	if (g_o.eq.0) then
131	write(,) "euvheat: Error; no O tracer !!!"
132	write(,) "O is always needed if calleuv=.true."
133	stop
134	else
135	nespeuv=nespeuv+1
136	endif
137	g_o2=igcm_o2
138	if (g_o2.eq.0) then
139	write(,) "euvheat: Error; no O2 tracer !!!"
140	write(,) "O2 is always needed if calleuv=.true."
141	stop
142	else
143	nespeuv=nespeuv+1
144	endif
145	g_h2=igcm_h2
146	if (g_h2.eq.0) then
147	write(,) "euvheat: Error; no H2 tracer !!!"
148	write(,) "H2 is always needed if calleuv=.true."
149	stop
150	else
151	nespeuv=nespeuv+1
152	endif
153	g_oh=igcm_oh
154	if (g_oh.eq.0) then
155	write(,) "euvheat: Error; no OH tracer !!!"
156	write(,) "OH must always be present if thermochem=T"
157	stop
158	else
159	nespeuv=nespeuv+1
160	endif
161	g_ho2=igcm_ho2
162	if (g_ho2.eq.0) then
163	write(,) "euvheat: Error; no HO2 tracer !!!"
164	write(,) "HO2 must always be present if thermochem=T"
165	stop
166	else
167	nespeuv=nespeuv+1
168	endif
169	g_h2o2=igcm_h2o2
170	if (g_h2o2.eq.0) then
171	write(,) "euvheat: Error; no H2O2 tracer !!!"
172	write(,) "H2O2 is always needed if calleuv=.true."
173	stop
174	else
175	nespeuv=nespeuv+1
176	endif
177	g_h2o=igcm_h2o_vap
178	if (g_h2o.eq.0) then
179	write(,) "euvheat: Error; no water vapor tracer !!!"
180	write(,) "H2O is always needed if calleuv=.true."
181	stop
182	else
183	nespeuv=nespeuv+1
184	endif
185	g_o1d=igcm_o1d
186	if (g_o1d.eq.0) then
187	write(,) "euvheat: Error; no O1D tracer !!!"
188	write(,) "O1D must always be present if thermochem=T"
189	stop
190	else
191	nespeuv=nespeuv+1
192	endif
193	g_co=igcm_co
194	if (g_co.eq.0) then
195	write(,) "euvheat: Error; no CO tracer !!!"
196	write(,) "CO is always needed if calleuv=.true."
197	stop
198	else
199	nespeuv=nespeuv+1
200	endif
201	g_h=igcm_h
202	if (g_h.eq.0) then
203	write(,) "euvheat: Error; no H tracer !!!"
204	write(,) "H is always needed if calleuv=.true."
205	stop
206	else
207	nespeuv=nespeuv+1
208	endif
209
210	euvmod = 0 !Default: C/O/H chemistry
211	!Check if O3 is present
212	g_o3=igcm_o3
213	if (g_o3.eq.0) then
214	write(,) "euvheat: Error; no O3 tracer !!!"
215	write(,) "O3 must be present if calleuv=.true."
216	stop
217	else
218	nespeuv=nespeuv+1
219	euvmod=1
220	endif
221
222	!Nitrogen species
223	!NO is used to determine if N chemistry is wanted
224	!euvmod=2 -> N chemistry
225	g_no=igcm_no
226	if (g_no.eq.0) then
227	write(,) "euvheat: no NO tracer"
228	write(,) "No N species in UV heating"
229	else if(g_no.ne.0) then
230	nespeuv=nespeuv+1
231	euvmod=2
232	endif
233	! n2
234	g_n2=igcm_n2
235	if(euvmod.eq.2) then
236	if (g_n2.eq.0) then
237	write(,) "euvheat: Error; no N2 tracer !!!"
238	write(,) "N2 needed if NO is in traceur.def"
239	stop
240	else
241	nespeuv=nespeuv+1
242	endif
243	endif ! Of if(euvmod.eq.2)
244	! N
245	g_n=igcm_n
246	if(euvmod == 2) then
247	if (g_n.eq.0) then
248	write(,) "euvheat: Error; no N tracer !!!"
249	write(,) "N needed if NO is in traceur.def"
250	stop
251	else if(g_n.ne.0) then
252	nespeuv=nespeuv+1
253	endif
254	else
255	if(g_n /= 0) then
256	write(,) "euvheat: Error: N present, but NO not!!!"
257	write(,) "Both must be in traceur.def"
258	stop
259	endif
260	endif !Of if(euvmod==2)
261	!NO2
262	g_no2=igcm_no2
263	if(euvmod == 2) then
264	if (g_no2.eq.0) then
265	write(,) "euvheat: Error; no NO2 tracer !!!"
266	write(,) "NO2 needed if NO is in traceur.def"
267	stop
268	else if(g_no2.ne.0) then
269	nespeuv=nespeuv+1
270	endif
271	else
272	if(g_no2 /= 0) then
273	write(,) "euvheat: Error: NO2 present, but NO not!!!"
274	write(,) "Both must be in traceur.def"
275	stop
276	endif
277	endif !Of if(euvmod==2)
278	!N2D
279	g_n2d=igcm_n2d
280	if(euvmod == 2) then
281	if (g_n2d.eq.0) then
282	write(,) "euvheat: Error; no N2D tracer !!!"
283	write(,) "N2D needed if NO is in traceur.def"
284	stop
285	else
286	nespeuv=nespeuv+1
287	endif
288	else
289	if(g_n2d /= 0) then
290	write(,) "euvheat: Error: N2D present, but NO not!!!"
291	write(,) "Both must be in traceur.def"
292	stop
293	endif
294	endif !Of if(euvmod==2)
295
296	!Check if nespeuv is appropriate for the value of euvmod
297	select case(euvmod)
298	case(0)
299	if(nespeuv.ne.11) then
300	write(,)'euvheat: Wrong number of tracers!'
301	stop
302	else
303	write(,)'euvheat: Computing absorption by',nespeuv, &
304	' species'
305	endif
306	case(1)
307	if(nespeuv.ne.12) then
308	write(,)'euvheat: Wrong number of tracers!',nespeuv
309	stop
310	else
311	write(,)'euvheat: Computing absorption by',nespeuv, &
312	' species'
313	endif
314	case(2)
315	if(nespeuv.ne.17) then
316	write(,)'euvheat: Wrong number of tracers!'
317	stop
318	else
319	write(,)'euvheat: Computing absorption by',nespeuv, &
320	' species'
321	endif
322	end select
323
324	firstcall= .false.
325	endif ! of if (firstcall)
326
327	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
328
329	!Number of species if not firstcall
330
331
332	!Allocate density vector
333	allocate(rm(nlayermx,nespeuv))
334	! build local updated values of tracers and temperature
335	do l=1,nlayermx
336	do ig=1,ngridmx
337	! chemical species
338	zq(ig,l,g_co2)=pq(ig,l,g_co2)+pdq(ig,l,g_co2)*ptimestep
339	zq(ig,l,g_o2)=pq(ig,l,g_o2)+pdq(ig,l,g_o2)*ptimestep
340	zq(ig,l,g_o)=pq(ig,l,g_o)+pdq(ig,l,g_o)*ptimestep
341	zq(ig,l,g_h2)=pq(ig,l,g_h2)+pdq(ig,l,g_h2)*ptimestep
342	zq(ig,l,g_h2o2)=pq(ig,l,g_h2o2)+pdq(ig,l,g_h2o2)*ptimestep
343	zq(ig,l,g_h2o)=pq(ig,l,g_h2o)+pdq(ig,l,g_h2o)*ptimestep
344	zq(ig,l,g_n2)=pq(ig,l,g_n2)+pdq(ig,l,g_n2)*ptimestep
345	zq(ig,l,g_co)=pq(ig,l,g_co)+pdq(ig,l,g_co)*ptimestep
346	zq(ig,l,g_h)=pq(ig,l,g_h)+pdq(ig,l,g_h)*ptimestep
347	if(euvmod.ge.1) &
348	zq(ig,l,g_o3)=pq(ig,l,g_o3)+pdq(ig,l,g_o3)*ptimestep
349	if(euvmod.eq.2) then
350	zq(ig,l,g_n)=pq(ig,l,g_n)+pdq(ig,l,g_n)*ptimestep
351	zq(ig,l,g_no)=pq(ig,l,g_no)+pdq(ig,l,g_no)*ptimestep
352	zq(ig,l,g_no2)=pq(ig,l,g_no2)+pdq(ig,l,g_no2)*ptimestep
353	endif
354	if(euvmod.gt.2.or.euvmod.lt.0) then
355	write(,)'euvheat: bad value for euvmod. Stop'
356	stop
357	endif
358	! atmospheric temperature
359	zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
360	enddo
361	enddo
362
363	!Solar flux calculation
364	if(solvarmod.eq.0) call flujo(solarcondate)
365
366	! Not recommended for long runs
367	! (e.g. to build the MCD), as the
368	! solar conditions at the end will
369	! be different to the ones initially
370	! set
371
372	do ig=1,ngridmx
373	zenit=acos(mu0(ig))*180./acos(-1.)
374
375	do l=1,nlayermx
376	!Conversion to number density
377	dens=pplay(ig,l)/(rnew(ig,l)*zt(ig,l)) / 1.66e-21
378	rm(l,i_co2) = zq(ig,l,g_co2) * dens / mmol(g_co2)
379	rm(l,i_o2) = zq(ig,l,g_o2) * dens / mmol(g_o2)
380	rm(l,i_o) = zq(ig,l,g_o) * dens / mmol(g_o)
381	rm(l,i_h2) = zq(ig,l,g_h2) * dens / mmol(g_h2)
382	rm(l,i_h2o) = zq(ig,l,g_h2o) * dens / mmol(g_h2o)
383	rm(l,i_h2o2) = zq(ig,l,g_h2o2) * dens / mmol(g_h2o2)
384	rm(l,i_co) = zq(ig,l,g_co) * dens / mmol(g_co)
385	rm(l,i_h) = zq(ig,l,g_h) * dens / mmol(g_h)
386	!Only if O3, N or ion chemistry requested
387	if(euvmod.ge.1) &
388	rm(l,i_o3) = zq(ig,l,g_o3) * dens / mmol(g_o3)
389	!Only if N or ion chemistry requested
390	if(euvmod.ge.2) then
391	rm(l,i_n2) = zq(ig,l,g_n2) * dens / mmol(g_n2)
392	rm(l,i_n) = zq(ig,l,g_n) * dens / mmol(g_n)
393	rm(l,i_no) = zq(ig,l,g_no) * dens / mmol(g_no)
394	rm(l,i_no2) = zq(ig,l,g_no2) * dens / mmol(g_no2)
395	endif
396	enddo
397
398	! zlocal(1)=-log(pplay(ig,1)/pplev(ig,1))
399	! & Rnew(ig,1)zt(ig,1)/g
400	zlocal(1)=zzlay(ig,1)
401	zlocal(1)=zlocal(1)/1000.
402	tx(1)=zt(ig,1)
403
404	do l=2,nlayermx
405	tx(l)=zt(ig,l)
406	zlocal(l)=zzlay(ig,l)/1000.
407	enddo
408	!Routine to calculate the UV heating
409	call hrtherm (ig,euvmod,rm,nespeuv,tx,zlocal,zenit,zday,jtot)
410
411	! euveff=0.16 !UV heating efficiency. Following Fox et al. ASR 1996
412	!should vary between 19% and 23%. Lower values
413	!(i.e. 16%) can be used to compensate underestimation
414	!of 15-um cooling (see Forget et al. JGR 2009 and
415	!Gonzalez-Galindo et al. JGR 2009) for details
416	!Calculates the UV heating from the total photoabsorption coefficient
417	do l=1,nlayermx
418	pdteuv(ig,l)=euveff*jtot(l)/10. &
419	/(cpnew(ig,l)pplay(ig,l)/(rnew(ig,l)zt(ig,l)))
420	! & (1.52/dist_sol)*2 !The solar flux calculated in
421	!flujo.F is already corrected for
422	!the actual Mars-Sun distance
423	enddo
424	enddo ! of do ig=1,ngridmx
425	!Deallocations
426	deallocate(rm)
427
428	return
429	end

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

Download in other formats:

Original Format