Context Navigation

euvheat.F90 @ 2780

Last change on this file since 2780 was 2464, checked in by slebonnois, 5 years ago
SL: major update related to extension to 250 km. New EUV heating as used in Martian GCM, debug of a few routines, adding He, clean up, updating mmean regularly, modification of the order of processes, add the possibility to use Hedin composition in rad transfer
File size: 11.3 KB

Line
1	SUBROUTINE euvheat(nlon, nlev,nqmx, pt,pplev,pplay,zzlay, &
2	mu0,ptimestep,ptime,pq, pdq, pdteuv)
3
4	use chemparam_mod
5	use dimphy
6	use conc, only: rnew, cpnew
7
8	IMPLICIT NONE
9	!=======================================================================
10	! subject:
11	! --------
12	! Computing heating rate due to EUV absorption
13	!
14	! author: MAC 2002
15	! ------
16	!
17	! input:
18	! -----
19	! mu0(klon)
20	! pplay(ngrid,nlayer) pressure at middle of layers (Pa)
21	! zzlay ! altitude at the middle of the layers (m)
22	!
23	! output:
24	! -------
25	!
26	! pdteuv(ngrid,nlayer) Heating rate (K/s)
27	!
28	!=======================================================================
29	!
30	! 0. Declarations :
31	! ------------------
32	!
33	#include "YOMCST.h"
34	#include "clesphys.h"
35	#include "mmol.h"
36	!-----------------------------------------------------------------------
37	! Input/Output
38	! ------------
39
40
41	integer :: nlon
42	integer :: nlev
43	integer :: nqmx
44
45	real :: pt(nlon,nlev)
46	real :: pplev(nlon,nlev+1)
47	real :: pplay(nlon,nlev)
48	real :: zzlay(nlon,nlev)
49	real :: mu0(nlon)
50
51	real :: ptimestep,ptime
52	real :: pq(nlon,nlev,nqmx)
53	real :: pdq(nlon,nlev,nqmx)
54	real :: pdteuv(nlon,nlev)
55	!
56	! Local variables :
57	! -----------------
58
59	integer,save :: nespeuv=17 ! Number of species considered (11, 12 or 17 (with nitrogen))
60	integer,save :: nspeuv_vgcm ! Number of species considered currently considered into VGCM
61
62
63	INTEGER :: l,ig,n
64	integer,save :: euvmod = 0 !0: 4 (main) species 1: O3 chemistry 2: N chemistry, 3: C/O/H
65	real, allocatable, save :: rm(:,:) ! number density (cm-3)
66	real :: zq(nlon,nlev,nqmx) ! local updated tracer quantity
67	real :: zt(nlon,nlev) ! local updated atmospheric temperature
68	real :: zlocal(nlev)
69	real :: zenit
70	real :: jtot(nlev)
71	real :: dens ! Total number density (cm-3)
72	real :: tx(nlev)
73
74	! tracer indexes for the EUV heating:
75	!!! ATTENTION. These values have to be identical to those in hrterm.F90
76	!!! If the values are changed there, the same has to be done here !!!
77
78	integer,parameter :: ix_co2=1
79	integer,parameter :: ix_o=3
80	integer,parameter :: ix_co=4
81	integer,parameter :: ix_n2=13
82
83
84	! Tracer indexes in the GCM:
85	integer,save :: g_co2=0
86	integer,save :: g_o=0
87	integer,save :: g_co=0
88	integer,save :: g_n2=0
89
90	logical,save :: firstcall=.true.
91
92	! Initializations and sanity checks:
93
94
95	if (firstcall) then
96	nspeuv_vgcm=0
97	! ! identify the indexes of the tracers we'll need
98	g_co2=i_co2
99	if (g_co2.eq.0) then
100	write(,) "euvheat: Error; no CO2 tracer !!!"
101	write(,) "CO2 is always needed if calleuv=.true."
102	stop
103	else
104	nspeuv_vgcm=nspeuv_vgcm+1
105	endif
106	g_o=i_o
107	if (g_o.eq.0) then
108	write(,) "euvheat: Error; no O tracer !!!"
109	! write(,) "O is always needed if calleuv=.true."
110	stop
111	else
112	nspeuv_vgcm=nspeuv_vgcm+1
113	endif
114	g_co=i_co
115	if (g_co.eq.0) then
116	write(,) "euvheat: Error; no CO tracer !!!"
117	! write(,) "CO is always needed if calleuv=.true."
118	stop
119	else
120	nspeuv_vgcm=nspeuv_vgcm+1
121	endif
122	! n2
123	g_n2=i_n2
124	if (g_n2.eq.0) then
125	write(,) "euvheat: Error; no N2 tracer !!!"
126	! write(,) "N2 needed if NO is in traceur.def"
127	stop
128	else
129	nspeuv_vgcm=nspeuv_vgcm+1
130	endif
131
132	! g_o2=igcm_o2
133	! if (g_o2.eq.0) then
134	! write(,) "euvheat: Error; no O2 tracer !!!"
135	! write(,) "O2 is always needed if calleuv=.true."
136	! stop
137	! else
138	! nespeuv=nespeuv+1
139	! endif
140	! g_h2=igcm_h2
141	! if (g_h2.eq.0) then
142	! write(,) "euvheat: Error; no H2 tracer !!!"
143	! write(,) "H2 is always needed if calleuv=.true."
144	! stop
145	! else
146	! nespeuv=nespeuv+1
147	! endif
148	! g_oh=igcm_oh
149	! if (g_oh.eq.0) then
150	! write(,) "euvheat: Error; no OH tracer !!!"
151	! write(,) "OH must always be present if thermochem=T"
152	! stop
153	! else
154	! nespeuv=nespeuv+1
155	! endif
156	! g_ho2=igcm_ho2
157	! if (g_ho2.eq.0) then
158	! write(,) "euvheat: Error; no HO2 tracer !!!"
159	! write(,) "HO2 must always be present if thermochem=T"
160	! stop
161	! else
162	! nespeuv=nespeuv+1
163	! endif
164	! g_h2o2=igcm_h2o2
165	! if (g_h2o2.eq.0) then
166	! write(,) "euvheat: Error; no H2O2 tracer !!!"
167	! write(,) "H2O2 is always needed if calleuv=.true."
168	! stop
169	! else
170	! nespeuv=nespeuv+1
171	! endif
172	! g_h2o=igcm_h2o_vap
173	! if (g_h2o.eq.0) then
174	! write(,) "euvheat: Error; no water vapor tracer !!!"
175	! write(,) "H2O is always needed if calleuv=.true."
176	! stop
177	! else
178	! nespeuv=nespeuv+1
179	! endif
180	! g_o1d=igcm_o1d
181	! if (g_o1d.eq.0) then
182	! write(,) "euvheat: Error; no O1D tracer !!!"
183	! write(,) "O1D must always be present if thermochem=T"
184	! stop
185	! else
186	! nespeuv=nespeuv+1
187	! endif
188	! g_h=igcm_h
189	! if (g_h.eq.0) then
190	! write(,) "euvheat: Error; no H tracer !!!"
191	! write(,) "H is always needed if calleuv=.true."
192	! stop
193	! else
194	! nespeuv=nespeuv+1
195	! endif
196
197	! euvmod = 3 !Default: C/O/H chemistry
198	! !Check if O3 is present
199	! g_o3=igcm_o3
200	! if (g_o3.eq.0) then
201	! write(,) "euvheat: Error; no O3 tracer !!!"
202	! write(,) "O3 must be present if calleuv=.true."
203	! stop
204	! else
205	! nespeuv=nespeuv+1
206	! euvmod=1
207	! endif
208
209	!Nitrogen species
210	!NO is used to determine if N chemistry is wanted
211	!euvmod=2 -> N chemistry
212	! g_no=igcm_no
213	! if (g_no.eq.0) then
214	! write(,) "euvheat: no NO tracer"
215	! write(,) "No N species in UV heating"
216	! else if(g_no.ne.0) then
217	! nespeuv=nespeuv+1
218	! euvmod=2
219	! endif
220	! N
221	! g_n=igcm_n
222	! if(euvmod == 2) then
223	! if (g_n.eq.0) then
224	! write(,) "euvheat: Error; no N tracer !!!"
225	! write(,) "N needed if NO is in traceur.def"
226	! stop
227	! else if(g_n.ne.0) then
228	! nespeuv=nespeuv+1
229	! endif
230	! else
231	! if(g_n /= 0) then
232	! write(,) "euvheat: Error: N present, but NO not!!!"
233	! write(,) "Both must be in traceur.def"
234	! stop
235	! endif
236	! endif !Of if(euvmod==2)
237	! !NO2
238	! g_no2=igcm_no2
239	! if(euvmod == 2) then
240	! if (g_no2.eq.0) then
241	! write(,) "euvheat: Error; no NO2 tracer !!!"
242	! write(,) "NO2 needed if NO is in traceur.def"
243	! stop
244	! else if(g_no2.ne.0) then
245	! nespeuv=nespeuv+1
246	! endif
247	! else
248	! if(g_no2 /= 0) then
249	! write(,) "euvheat: Error: NO2 present, but NO not!!!"
250	! write(,) "Both must be in traceur.def"
251	! stop
252	! endif
253	! endif !Of if(euvmod==2)
254	!N2D
255	! g_n2d=igcm_n2d
256	! if(euvmod == 2) then
257	! if (g_n2d.eq.0) then
258	! write(,) "euvheat: Error; no N2D tracer !!!"
259	! write(,) "N2D needed if NO is in traceur.def"
260	! stop
261	! else
262	! nespeuv=nespeuv+1
263	! endif
264	! else
265	! if(g_n2d /= 0) then
266	! write(,) "euvheat: Error: N2D present, but NO not!!!"
267	! write(,) "Both must be in traceur.def"
268	! stop
269	! endif
270	! endif !Of if(euvmod==2)
271
272	!Check if nespeuv is appropriate for the value of euvmod
273	! select case(euvmod)
274	! case(0)
275	! if(nespeuv.ne.11) then
276	! write(,)'euvheat: Wrong number of tracers!'
277	! stop
278	! else
279	! write(,)'euvheat: Computing absorption by',nespeuv, &
280	! ' species'
281	! endif
282	! case(1)
283	! if(nespeuv.ne.12) then
284	! write(,)'euvheat: Wrong number of tracers!',nespeuv
285	! stop
286	! else
287	! write(,)'euvheat: Computing absorption by',nespeuv, &
288	! ' species'
289	! endif
290	! case(2)
291	! if(nespeuv.ne.17) then
292	! write(,)'euvheat: Wrong number of tracers!'
293	! stop
294	! else
295	! write(,)'euvheat: Computing absorption by',nespeuv, &
296	! ' species'
297	! endif
298	! end select
299
300
301	!Allocate density vector
302	allocate(rm(nlev,nespeuv))
303
304	firstcall= .false.
305	endif ! of if (firstcall)
306
307	! write(,), "CHECK n species currently used into VGCM", nspeuv_vgcm
308
309
310	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
311
312	! build local updated values of tracers (if any) and temperature
313
314	do l=1,nlev
315	do ig=1,nlon
316
317	! chemical species
318	zq(ig,l,g_co2)=pq(ig,l,g_co2)
319	zq(ig,l,g_co)=pq(ig,l,g_co)
320	zq(ig,l,g_o)=pq(ig,l,g_o)
321	zq(ig,l,g_n2)=pq(ig,l,g_n2)
322
323	! atmospheric temperature
324	zt(ig,l)=pt(ig,l)
325
326	! write(,), "CHECK update densities L332 euv", zq(ig,l,g_co2)
327
328
329	enddo
330	enddo
331
332	!Solar flux calculation
333
334	do ig=1,nlon
335	zenit=acos(mu0(ig))*180./acos(-1.) !convers from rad to deg
336
337	do l=1,nlev
338
339	! Conversion to number density
340
341	!!! use R specific = R/MolarMass
342
343	dens=pplay(ig,l)/(rnew(ig,l)*zt(ig,l)) / 1.66e-21 ! [g mol-1] [cm-3]
344
345	rm(l,ix_co2) = zq(ig,l,g_co2) * dens / M_tr(g_co2) ! [cm-3]
346	rm(l,ix_o) = zq(ig,l,g_o) * dens / M_tr(g_o)
347	rm(l,ix_co) = zq(ig,l,g_co) * dens / M_tr(g_co)
348	rm(l,ix_n2) = zq(ig,l,g_n2) * dens / M_tr(g_n2)
349
350	! write(,), "CHECK n density", l, rm(l,ix_co2)
351
352
353	enddo
354
355	! zlocal(1)=-log(pplay(ig,1)/pplev(ig,1))
356	! & Rnew(ig,1)zt(ig,1)/g
357	zlocal(1)=zzlay(ig,1)
358	zlocal(1)=zlocal(1)/1000. ! conversion m ---> km
359	tx(1)=zt(ig,1)
360
361	do l=2,nlev
362	tx(l)=zt(ig,l)
363	zlocal(l)=zzlay(ig,l)/1000.
364	enddo
365
366	!Routine to calculate the UV heating
367	call hrtherm (ig,euvmod,rm,nespeuv,tx,zlocal,zenit,jtot)
368
369
370	!Calculates the UV heating from the total photoabsorption coefficient
371	do l=1,nlev
372	! jtot conversion from erg/(scm3) ---> J/(sm3)
373	pdteuv(ig,l)=euveff*jtot(l)/10. &
374	/(cpnew(ig,l)pplay(ig,l)/(rnew(ig,l)zt(ig,l)))
375
376	enddo
377	enddo ! of do ig=1,nlon
378
379	!Deallocations
380	!deallocate(rm)
381
382	return
383	end

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

Download in other formats:

Original Format