Context Navigation

euvheat.F90 @ 1521

Last change on this file since 1521 was 1442, checked in by slebonnois, 10 years ago
SL: update of the Venus GCM, + corrections on routines used for newstart/start2archive for Titan and Venus, + some modifications on tools
File size: 11.8 KB

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

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

Download in other formats:

Original Format