Context Navigation

euvheat.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: 14.4 KB

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

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

Download in other formats:

Original Format