Context Navigation

euvheat.F90 @ 2187

Last change on this file since 2187 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

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

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

Download in other formats:

Original Format