Context Navigation

jthermcalc.F @ 1198

Last change on this file since 1198 was 1119, checked in by emillour, 11 years ago

Mars GCM:

Bug fix: Sun-Mars distance was not correctly taken into account in the solvarmod==1 (daily varying realistic EUV input) case.

FGG

File size: 59.3 KB

Line
1	c**********************************************************************
2
3	subroutine jthermcalc(ig,chemthermod,rm,nesptherm,tx,iz,zenit)
4
5
6	c feb 2002 fgg first version
7	c nov 2002 fgg second version
8	c
9	c modified from paramhr.F
10	c MAC July 2003
11	c**********************************************************************
12
13	implicit none
14
15	c common variables and constants
16	include "dimensions.h"
17	include "dimphys.h"
18	include 'param.h'
19	include 'param_v4.h'
20
21	c input and output variables
22
23	integer ig
24	integer chemthermod
25	integer nesptherm !Number of species considered
26	real rm(nlayermx,nesptherm) !Densities (cm-3)
27	real tx(nlayermx) !temperature
28	real zenit !SZA
29	real iz(nlayermx) !Local altitude
30
31
32	c local parameters and variables
33
34	real co2colx(nlayermx) !column density of CO2 (cm^-2)
35	real o2colx(nlayermx) !column density of O2(cm^-2)
36	real o3pcolx(nlayermx) !column density of O(3P)(cm^-2)
37	real h2colx(nlayermx) !H2 column density (cm-2)
38	real h2ocolx(nlayermx) !H2O column density (cm-2)
39	real h2o2colx(nlayermx) !column density of H2O2(cm^-2)
40	real o3colx(nlayermx) !O3 column density (cm-2)
41	real n2colx(nlayermx) !N2 column density (cm-2)
42	real ncolx(nlayermx) !N column density (cm-2)
43	real nocolx(nlayermx) !NO column density (cm-2)
44	real cocolx(nlayermx) !CO column density (cm-2)
45	real hcolx(nlayermx) !H column density (cm-2)
46	real no2colx(nlayermx) !NO2 column density (cm-2)
47	real t2(nlayermx)
48	real coltemp(nlayermx)
49	real sigma(ninter,nlayermx)
50	real alfa(ninter,nlayermx)
51
52	integer i,j,k,indexint !indexes
53	character dn
54
55
56
57	c variables used in interpolation
58
59	real*8 auxcoltab(nz2)
60	real*8 auxjco2(nz2)
61	real*8 auxjo2(nz2)
62	real*8 auxjo3p(nz2)
63	real*8 auxjh2o(nz2)
64	real*8 auxjh2(nz2)
65	real*8 auxjh2o2(nz2)
66	real*8 auxjo3(nz2)
67	real*8 auxjn2(nz2)
68	real*8 auxjn(nz2)
69	real*8 auxjno(nz2)
70	real*8 auxjco(nz2)
71	real*8 auxjh(nz2)
72	real*8 auxjno2(nz2)
73	real*8 wp(nlayermx),wm(nlayermx)
74	real*8 auxcolinp(nlayermx)
75	integer auxind(nlayermx)
76	integer auxi
77	integer ind
78	real*8 cortemp(nlayermx)
79
80	real*8 limdown !limits for interpolation
81	real*8 limup ! ""
82
83	!!!ATTENTION. Here i_co2 has to have the same value than in chemthermos.F90
84	!!!If the value is changed there, if has to be changed also here !!!!
85	integer,parameter :: i_co2=1
86
87
88	c***********************PROGRAM STARTS*****************************
89
90	if(zenit.gt.140.) then
91	dn='n'
92	else
93	dn='d'
94	end if
95	if(dn.eq.'n') then
96	return
97	endif
98
99	!Initializing the photoabsorption coefficients
100	jfotsout(:,:,:)=0.
101
102	!Auxiliar temperature to take into account the temperature dependence
103	!of CO2 cross section
104	do i=1,nlayermx
105	t2(i)=tx(i)
106	if(t2(i).lt.195.0) t2(i)=195.0
107	if(t2(i).gt.295.0) t2(i)=295.0
108	end do
109
110	!Calculation of column amounts
111	call column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
112	$ co2colx,o2colx,o3pcolx,h2colx,h2ocolx,
113	$ h2o2colx,o3colx,n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
114
115	!Auxiliar column to include the temperature dependence
116	!of CO2 cross section
117	coltemp(nlayermx)=co2colx(nlayermx)*abs(t2(nlayermx)-t0(nlayermx))
118	do i=nlayermx-1,1,-1
119	coltemp(i)=!coltemp(i+1)+ PQ SE ELIMINA? REVISAR
120	$ ( rm(i,i_co2) + rm(i+1,i_co2) ) * 0.5
121	$ * 1e5 * (iz(i+1)-iz(i)) * abs(t2(i)-t0(i))
122	end do
123
124	!Calculation of CO2 cross section at temperature t0(i)
125	do i=1,nlayermx
126	do indexint=24,32
127	sigma(indexint,i)=co2crsc195(indexint-23)
128	alfa(indexint,i)=((co2crsc295(indexint-23)
129	$ /sigma(indexint,i))-1.)/(295.-t0(i))
130	end do
131	end do
132
133	! Interpolation to the tabulated photoabsorption coefficients for each species
134	! in each spectral interval
135
136
137	c auxcolinp-> Actual atmospheric column
138	c auxj*-> Tabulated photoabsorption coefficients
139	c auxcoltab-> Tabulated atmospheric columns
140
141	ccccccccccccccccccccccccccccccc
142	c 0.1,5.0 (int 1)
143	c
144	c Absorption by:
145	c CO2, O2, O, H2, N
146	ccccccccccccccccccccccccccccccc
147
148	c Input atmospheric column
149	indexint=1
150	do i=1,nlayermx
151	auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint) +
152	$ o2colx(i)*crscabsi2(2,indexint) +
153	$ o3pcolx(i)*crscabsi2(3,indexint) +
154	$ h2colx(i)*crscabsi2(5,indexint) +
155	$ ncolx(i)*crscabsi2(9,indexint)
156	end do
157	limdown=1.e-20
158	limup=1.e26
159
160
161	c Interpolations
162
163	do i=1,nz2
164	auxi = nz2-i+1
165	!CO2 tabulated coefficient
166	auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
167	!O2 tabulated coefficient
168	auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
169	!O3p tabulated coefficient
170	auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
171	!H2 tabulated coefficient
172	auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
173	!Tabulated column
174	auxcoltab(i) = c1_16(auxi,indexint)
175	enddo
176	!Only if chemthermod.ge.2
177	!N tabulated coefficient
178	if(chemthermod.ge.2) then
179	do i=1,nz2
180	auxjn(i) = jabsifotsintpar(nz2-i+1,9,indexint)
181	enddo
182	endif
183
184	call interfast
185	$ (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
186	do i=1,nlayermx
187	ind=auxind(i)
188	auxi=nlayermx-i+1
189	!CO2 interpolated coefficient
190	jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
191	$ wp(i)*auxjco2(ind)
192	!O2 interpolated coefficient
193	jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
194	$ wp(i)*auxjo2(ind)
195	!O3p interpolated coefficient
196	jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
197	$ wp(i)*auxjo3p(ind)
198	!H2 interpolated coefficient
199	jfotsout(indexint,5,auxi) = wm(i)*auxjh2(ind+1) +
200	$ wp(i)*auxjh2(ind)
201	enddo
202	!Only if chemthermod.ge.2
203	!N interpolated coefficient
204	if(chemthermod.ge.2) then
205	do i=1,nlayermx
206	ind=auxind(i)
207	jfotsout(indexint,9,nlayermx-i+1) = wm(i)*auxjn(ind+1) +
208	$ wp(i)*auxjn(ind)
209	enddo
210	endif
211
212
213	c End interval 1
214
215
216	ccccccccccccccccccccccccccccccc
217	c 5-80.5nm (int 2-15)
218	c
219	c Absorption by:
220	c CO2, O2, O, H2, N2, N,
221	c NO, CO, H, NO2
222	ccccccccccccccccccccccccccccccc
223
224	c Input atmospheric column
225	do indexint=2,15
226	do i=1,nlayermx
227	auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
228	$ o2colx(i)*crscabsi2(2,indexint)+
229	$ o3pcolx(i)*crscabsi2(3,indexint)+
230	$ h2colx(i)*crscabsi2(5,indexint)+
231	$ n2colx(i)*crscabsi2(8,indexint)+
232	$ ncolx(i)*crscabsi2(9,indexint)+
233	$ nocolx(i)*crscabsi2(10,indexint)+
234	$ cocolx(i)*crscabsi2(11,indexint)+
235	$ hcolx(i)*crscabsi2(12,indexint)+
236	$ no2colx(i)*crscabsi2(13,indexint)
237	end do
238
239	c Interpolations
240
241	do i=1,nz2
242	auxi = nz2-i+1
243	!O2 tabulated coefficient
244	auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
245	!O3p tabulated coefficient
246	auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
247	!CO2 tabulated coefficient
248	auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
249	!H2 tabulated coefficient
250	auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
251	!N2 tabulated coefficient
252	auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
253	!CO tabulated coefficient
254	auxjco(i) = jabsifotsintpar(auxi,11,indexint)
255	!H tabulated coefficient
256	auxjh(i) = jabsifotsintpar(auxi,12,indexint)
257	!tabulated column
258	auxcoltab(i) = c1_16(auxi,indexint)
259	enddo
260	!Only if chemthermod.ge.2
261	if(chemthermod.ge.2) then
262	do i=1,nz2
263	auxi = nz2-i+1
264	!N tabulated coefficient
265	auxjn(i) = jabsifotsintpar(auxi,9,indexint)
266	!NO tabulated coefficient
267	auxjno(i) = jabsifotsintpar(auxi,10,indexint)
268	!NO2 tabulated coefficient
269	auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
270	enddo
271	endif
272
273	call interfast(wm,wp,auxind,auxcolinp,nlayermx,
274	$ auxcoltab,nz2,limdown,limup)
275	do i=1,nlayermx
276	ind=auxind(i)
277	auxi = nlayermx-i+1
278	!O2 interpolated coefficient
279	jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
280	$ wp(i)*auxjo2(ind)
281	!O3p interpolated coefficient
282	jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
283	$ wp(i)*auxjo3p(ind)
284	!CO2 interpolated coefficient
285	jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
286	$ wp(i)*auxjco2(ind)
287	!H2 interpolated coefficient
288	jfotsout(indexint,5,auxi) = wm(i)*auxjh2(ind+1) +
289	$ wp(i)*auxjh2(ind)
290	!N2 interpolated coefficient
291	jfotsout(indexint,8,auxi) = wm(i)*auxjn2(ind+1) +
292	$ wp(i)*auxjn2(ind)
293	!CO interpolated coefficient
294	jfotsout(indexint,11,auxi) = wm(i)*auxjco(ind+1) +
295	$ wp(i)*auxjco(ind)
296	!H interpolated coefficient
297	jfotsout(indexint,12,auxi) = wm(i)*auxjh(ind+1) +
298	$ wp(i)*auxjh(ind)
299	enddo
300	!Only if chemthermod.ge.2
301	if(chemthermod.ge.2) then
302	do i=1,nlayermx
303	ind=auxind(i)
304	auxi = nlayermx-i+1
305	!N interpolated coefficient
306	jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
307	$ wp(i)*auxjn(ind)
308	!NO interpolated coefficient
309	jfotsout(indexint,10,auxi)=wm(i)*auxjno(ind+1) +
310	$ wp(i)*auxjno(ind)
311	!NO2 interpolated coefficient
312	jfotsout(indexint,13,auxi)=wm(i)*auxjno2(ind+1)+
313	$ wp(i)*auxjno2(ind)
314	enddo
315	endif
316	end do
317
318	c End intervals 2-15
319
320
321	ccccccccccccccccccccccccccccccc
322	c 80.6-90.8nm (int16)
323	c
324	c Absorption by:
325	c CO2, O2, O, N2, N, NO,
326	c CO, H, NO2
327	ccccccccccccccccccccccccccccccc
328
329	c Input atmospheric column
330	indexint=16
331	do i=1,nlayermx
332	auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
333	$ o2colx(i)*crscabsi2(2,indexint)+
334	$ o3pcolx(i)*crscabsi2(3,indexint)+
335	$ n2colx(i)*crscabsi2(8,indexint)+
336	$ ncolx(i)*crscabsi2(9,indexint)+
337	$ nocolx(i)*crscabsi2(10,indexint)+
338	$ cocolx(i)*crscabsi2(11,indexint)+
339	$ hcolx(i)*crscabsi2(12,indexint)+
340	$ no2colx(i)*crscabsi2(13,indexint)
341	end do
342
343	c Interpolations
344
345	do i=1,nz2
346	auxi = nz2-i+1
347	!O2 tabulated coefficient
348	auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
349	!CO2 tabulated coefficient
350	auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
351	!O3p tabulated coefficient
352	auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
353	!N2 tabulated coefficient
354	auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
355	!CO tabulated coefficient
356	auxjco(i) = jabsifotsintpar(auxi,11,indexint)
357	!H tabulated coefficient
358	auxjh(i) = jabsifotsintpar(auxi,12,indexint)
359	!NO2 tabulated coefficient
360	auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
361	!Tabulated column
362	auxcoltab(i) = c1_16(auxi,indexint)
363	enddo
364	!Only if chemthermod.ge.2
365	if(chemthermod.ge.2) then
366	do i=1,nz2
367	auxi = nz2-i+1
368	!N tabulated coefficient
369	auxjn(i) = jabsifotsintpar(auxi,9,indexint)
370	!NO tabulated coefficient
371	auxjno(i) = jabsifotsintpar(auxi,10,indexint)
372	!NO2 tabulated coefficient
373	auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
374	enddo
375	endif
376
377	call interfast
378	$ (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
379	do i=1,nlayermx
380	ind=auxind(i)
381	auxi = nlayermx-i+1
382	!O2 interpolated coefficient
383	jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
384	$ wp(i)*auxjo2(ind)
385	!CO2 interpolated coefficient
386	jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
387	$ wp(i)*auxjco2(ind)
388	!O3p interpolated coefficient
389	jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
390	$ wp(i)*auxjo3p(ind)
391	!N2 interpolated coefficient
392	jfotsout(indexint,8,auxi) = wm(i)*auxjn2(ind+1) +
393	$ wp(i)*auxjn2(ind)
394	!CO interpolated coefficient
395	jfotsout(indexint,11,auxi) = wm(i)*auxjco(ind+1) +
396	$ wp(i)*auxjco(ind)
397	!H interpolated coefficient
398	jfotsout(indexint,12,auxi) = wm(i)*auxjh(ind+1) +
399	$ wp(i)*auxjh(ind)
400	enddo
401	!Only if chemthermod.ge.2
402	if(chemthermod.ge.2) then
403	do i=1,nlayermx
404	ind=auxind(i)
405	auxi = nlayermx-i+1
406	!N interpolated coefficient
407	jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
408	$ wp(i)*auxjn(ind)
409	!NO interpolated coefficient
410	jfotsout(indexint,10,auxi) = wm(i)*auxjno(ind+1) +
411	$ wp(i)*auxjno(ind)
412	!NO2 interpolated coefficient
413	jfotsout(indexint,13,auxi) = wm(i)*auxjno2(ind+1) +
414	$ wp(i)*auxjno2(ind)
415	enddo
416	endif
417	c End interval 16
418
419
420	ccccccccccccccccccccccccccccccc
421	c 90.9-119.5nm (int 17-24)
422	c
423	c Absorption by:
424	c CO2, O2, N2, NO, CO, NO2
425	ccccccccccccccccccccccccccccccc
426
427	c Input column
428
429	do i=1,nlayermx
430	auxcolinp(nlayermx-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
431	$ nocolx(i) + cocolx(i) + no2colx(i)
432	end do
433
434	do indexint=17,24
435
436	c Interpolations
437
438	do i=1,nz2
439	auxi = nz2-i+1
440	!CO2 tabulated coefficient
441	auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
442	!O2 tabulated coefficient
443	auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
444	!N2 tabulated coefficient
445	auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
446	!CO tabulated coefficient
447	auxjco(i) = jabsifotsintpar(auxi,11,indexint)
448	!Tabulated column
449	auxcoltab(i) = c17_24(auxi)
450	enddo
451	!Only if chemthermod.ge.2
452	if(chemthermod.ge.2) then
453	do i=1,nz2
454	auxi = nz2-i+1
455	!NO tabulated coefficient
456	auxjno(i) = jabsifotsintpar(auxi,10,indexint)
457	!NO2 tabulated coefficient
458	auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
459	enddo
460	endif
461
462	call interfast
463	$ (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
464	!Correction to include T variation of CO2 cross section
465	if(indexint.eq.24) then
466	do i=1,nlayermx
467	auxi = nlayermx-i+1
468	if(sigma(indexint,auxi)*
469	$ alfa(indexint,auxi)*coltemp(auxi)
470	$ .lt.60.) then
471	cortemp(i)=exp(-sigma(indexint,auxi)*
472	$ alfa(indexint,auxi)*coltemp(auxi))
473	else
474	cortemp(i)=0.
475	end if
476	enddo
477	else
478	do i=1,nlayermx
479	cortemp(i)=1.
480	enddo
481	end if
482	do i=1,nlayermx
483	ind=auxind(i)
484	auxi = nlayermx-i+1
485	!O2 interpolated coefficient
486	jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
487	$ wp(i)auxjo2(ind)) cortemp(i)
488	!CO2 interpolated coefficient
489	jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
490	$ wp(i)auxjco2(ind)) cortemp(i)
491	if(indexint.eq.24) jfotsout(indexint,1,auxi)=
492	$ jfotsout(indexint,1,auxi)*
493	$ (1+alfa(indexint,auxi)*
494	$ (t2(auxi)-t0(auxi)))
495	!N2 interpolated coefficient
496	jfotsout(indexint,8,auxi) = (wm(i)*auxjn2(ind+1) +
497	$ wp(i)auxjn2(ind)) cortemp(i)
498	!CO interpolated coefficient
499	jfotsout(indexint,11,auxi) = (wm(i)*auxjco(ind+1) +
500	$ wp(i)auxjco(ind)) cortemp(i)
501	enddo
502	!Only if chemthermod.ge.2
503	if(chemthermod.ge.2) then
504	do i=1,nlayermx
505	ind=auxind(i)
506	auxi = nlayermx-i+1
507	!NO interpolated coefficient
508	jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
509	$ wp(i)auxjno(ind)) cortemp(i)
510	!NO2 interpolated coefficient
511	jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
512	$ wp(i)auxjno2(ind)) cortemp(i)
513	enddo
514	endif
515	end do
516	c End intervals 17-24
517
518
519	ccccccccccccccccccccccccccccccc
520	c 119.6-167.0nm (int 25-29)
521	c
522	c Absorption by:
523	c CO2, O2, H2O, H2O2, NO,
524	c CO, NO2
525	ccccccccccccccccccccccccccccccc
526
527	c Input atmospheric column
528
529	do i=1,nlayermx
530	auxcolinp(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
531	$ h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
532	end do
533
534	do indexint=25,29
535
536	c Interpolations
537
538	do i=1,nz2
539	auxi = nz2-i+1
540	!CO2 tabulated coefficient
541	auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
542	!O2 tabulated coefficient
543	auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
544	!H2O tabulated coefficient
545	auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
546	!H2O2 tabulated coefficient
547	auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
548	!CO tabulated coefficient
549	auxjco(i) = jabsifotsintpar(auxi,11,indexint)
550	!Tabulated column
551	auxcoltab(i) = c25_29(auxi)
552	enddo
553	!Only if chemthermod.ge.2
554	if(chemthermod.ge.2) then
555	do i=1,nz2
556	auxi = nz2-i+1
557	!NO tabulated coefficient
558	auxjno(i) = jabsifotsintpar(auxi,10,indexint)
559	!NO2 tabulated coefficient
560	auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
561	enddo
562	endif
563	call interfast
564	$ (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
565	do i=1,nlayermx
566	ind=auxind(i)
567	auxi = nlayermx-i+1
568	!Correction to include T variation of CO2 cross section
569	if(sigma(indexint,auxi)alfa(indexint,auxi)
570	$ coltemp(auxi).lt.60.) then
571	cortemp(i)=exp(-sigma(indexint,auxi)*
572	$ alfa(indexint,auxi)*coltemp(auxi))
573	else
574	cortemp(i)=0.
575	end if
576	!CO2 interpolated coefficient
577	jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
578	$ wp(i)auxjco2(ind)) cortemp(i) *
579	$ (1+alfa(indexint,auxi)*
580	$ (t2(auxi)-t0(auxi)))
581	!O2 interpolated coefficient
582	jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
583	$ wp(i)auxjo2(ind)) cortemp(i)
584	!H2O interpolated coefficient
585	jfotsout(indexint,4,auxi) = (wm(i)*auxjh2o(ind+1) +
586	$ wp(i)auxjh2o(ind)) cortemp(i)
587	!H2O2 interpolated coefficient
588	jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
589	$ wp(i)auxjh2o2(ind)) cortemp(i)
590	!CO interpolated coefficient
591	jfotsout(indexint,11,auxi) = (wm(i)*auxjco(ind+1) +
592	$ wp(i)auxjco(ind)) cortemp(i)
593	enddo
594	!Only if chemthermod.ge.2
595	if(chemthermod.ge.2) then
596	do i=1,nlayermx
597	ind=auxind(i)
598	auxi = nlayermx-i+1
599	!NO interpolated coefficient
600	jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
601	$ wp(i)auxjno(ind)) cortemp(i)
602	!NO2 interpolated coefficient
603	jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
604	$ wp(i)auxjno2(ind)) cortemp(i)
605	enddo
606	endif
607
608	end do
609
610	c End intervals 25-29
611
612
613	cccccccccccccccccccccccccccccccc
614	c 167.1-202.5nm (int 30-31)
615	c
616	c Absorption by:
617	c CO2, O2, H2O, H2O2, NO,
618	c NO2
619	cccccccccccccccccccccccccccccccc
620
621	c Input atmospheric column
622
623	do i=1,nlayermx
624	auxcolinp(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
625	$ h2o2colx(i) + nocolx(i) + no2colx(i)
626	end do
627
628	c Interpolation
629
630	do indexint=30,31
631
632	do i=1,nz2
633	auxi = nz2-i+1
634	!CO2 tabulated coefficient
635	auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
636	!O2 tabulated coefficient
637	auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
638	!H2O tabulated coefficient
639	auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
640	!H2O2 tabulated coefficient
641	auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
642	!Tabulated column
643	auxcoltab(i) = c30_31(auxi)
644	enddo
645	!Only if chemthermod.ge.2
646	if(chemthermod.ge.2) then
647	do i=1,nz2
648	auxi = nz2-i+1
649	!NO tabulated coefficient
650	auxjno(i) = jabsifotsintpar(auxi,10,indexint)
651	!NO2 tabulated coefficient
652	auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
653	enddo
654	endif
655
656	call interfast
657	$ (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
658	do i=1,nlayermx
659	ind=auxind(i)
660	auxi = nlayermx-i+1
661	!Correction to include T variation of CO2 cross section
662	if(sigma(indexint,auxi)alfa(indexint,auxi)
663	$ coltemp(auxi).lt.60.) then
664	cortemp(i)=exp(-sigma(indexint,auxi)*
665	$ alfa(indexint,auxi)*coltemp(auxi))
666	else
667	cortemp(i)=0.
668	end if
669	!CO2 interpolated coefficient
670	jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
671	$ wp(i)auxjco2(ind)) cortemp(i) *
672	$ (1+alfa(indexint,auxi)*
673	$ (t2(auxi)-t0(auxi)))
674	!O2 interpolated coefficient
675	jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
676	$ wp(i)auxjo2(ind)) cortemp(i)
677	!H2O interpolated coefficient
678	jfotsout(indexint,4,auxi) = (wm(i)*auxjh2o(ind+1) +
679	$ wp(i)auxjh2o(ind)) cortemp(i)
680	!H2O2 interpolated coefficient
681	jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
682	$ wp(i)auxjh2o2(ind)) cortemp(i)
683	enddo
684	!Only if chemthermod.ge.2
685	if(chemthermod.ge.2) then
686	do i=1,nlayermx
687	ind=auxind(i)
688	auxi = nlayermx-i+1
689	!NO interpolated coefficient
690	jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
691	$ wp(i)auxjno(ind)) cortemp(i)
692	!NO2 interpolated coefficient
693	jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
694	$ wp(i)auxjno2(ind)) cortemp(i)
695	enddo
696	endif
697
698	end do
699
700	c End intervals 30-31
701
702
703	ccccccccccccccccccccccccccccccc
704	c 202.6-210.0nm (int 32)
705	c
706	c Absorption by:
707	c CO2, O2, H2O2, NO, NO2
708	ccccccccccccccccccccccccccccccc
709
710	c Input atmospheric column
711
712	indexint=32
713	do i=1,nlayermx
714	auxcolinp(nlayermx-i+1) =co2colx(i) + o2colx(i) + h2o2colx(i) +
715	$ nocolx(i) + no2colx(i)
716	end do
717
718	c Interpolation
719
720	do i=1,nz2
721	auxi = nz2-i+1
722	!CO2 tabulated coefficient
723	auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
724	!O2 tabulated coefficient
725	auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
726	!H2O2 tabulated coefficient
727	auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
728	!Tabulated column
729	auxcoltab(i) = c32(auxi)
730	enddo
731	!Only if chemthermod.ge.2
732	if(chemthermod.ge.2) then
733	do i=1,nz2
734	auxi = nz2-i+1
735	!NO tabulated coefficient
736	auxjno(i) = jabsifotsintpar(auxi,10,indexint)
737	!NO2 tabulated coefficient
738	auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
739	enddo
740	endif
741	call interfast
742	$ (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
743	do i=1,nlayermx
744	ind=auxind(i)
745	auxi = nlayermx-i+1
746	!Correction to include T variation of CO2 cross section
747	if(sigma(indexint,nlayermx-i+1)alfa(indexint,auxi)
748	$ coltemp(auxi).lt.60.) then
749	cortemp(i)=exp(-sigma(indexint,auxi)*
750	$ alfa(indexint,auxi)*coltemp(auxi))
751	else
752	cortemp(i)=0.
753	end if
754	!CO2 interpolated coefficient
755	jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
756	$ wp(i)auxjco2(ind)) cortemp(i) *
757	$ (1+alfa(indexint,auxi)*
758	$ (t2(auxi)-t0(auxi)))
759	!O2 interpolated coefficient
760	jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
761	$ wp(i)auxjo2(ind)) cortemp(i)
762	!H2O2 interpolated coefficient
763	jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
764	$ wp(i)auxjh2o2(ind)) cortemp(i)
765	enddo
766	!Only if chemthermod.ge.2
767	if(chemthermod.ge.2) then
768	do i=1,nlayermx
769	auxi = nlayermx-i+1
770	ind=auxind(i)
771	!NO interpolated coefficient
772	jfotsout(indexint,10,auxi) = (wm(i)*auxjno(ind+1) +
773	$ wp(i)auxjno(ind)) cortemp(i)
774	!NO2 interpolated coefficient
775	jfotsout(indexint,13,auxi) = (wm(i)*auxjno2(ind+1) +
776	$ wp(i)auxjno2(ind)) cortemp(i)
777	enddo
778	endif
779
780	c End of interval 32
781
782
783	ccccccccccccccccccccccccccccccc
784	c 210.1-231.0nm (int 33)
785	c
786	c Absorption by:
787	c O2, H2O2, NO2
788	ccccccccccccccccccccccccccccccc
789
790	c Input atmospheric column
791
792	indexint=33
793	do i=1,nlayermx
794	auxcolinp(nlayermx-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
795	end do
796
797	c Interpolation
798
799	do i=1,nz2
800	auxi = nz2-i+1
801	!O2 tabulated coefficient
802	auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
803	!H2O2 tabulated coefficient
804	auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
805	!Tabulated column
806	auxcoltab(i) = c33(auxi)
807	enddo
808	!Only if chemthermod.ge.2
809	if(chemthermod.ge.2) then
810	do i=1,nz2
811	!NO2 tabulated coefficient
812	auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
813	enddo
814	endif
815	call interfast
816	$ (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
817	do i=1,nlayermx
818	ind=auxind(i)
819	auxi = nlayermx-i+1
820	!O2 interpolated coefficient
821	jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
822	$ wp(i)*auxjo2(ind)
823	!H2O2 interpolated coefficient
824	jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
825	$ wp(i)*auxjh2o2(ind)
826	enddo
827	!Only if chemthermod.ge.2
828	if(chemthermod.ge.2) then
829	do i=1,nlayermx
830	ind=auxind(i)
831	!NO2 interpolated coefficient
832	jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
833	$ wp(i)*auxjno2(ind)
834	enddo
835	endif
836
837	c End of interval 33
838
839
840	ccccccccccccccccccccccccccccccc
841	c 231.1-240.0nm (int 34)
842	c
843	c Absorption by:
844	c O2, H2O2, O3, NO2
845	ccccccccccccccccccccccccccccccc
846
847	c Input atmospheric column
848
849	indexint=34
850	do i=1,nlayermx
851	auxcolinp(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
852	$ no2colx(i)
853	end do
854
855	c Interpolation
856
857	do i=1,nz2
858	auxi = nz2-i+1
859	!O2 tabulated coefficient
860	auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
861	!H2O2 tabulated coefficient
862	auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
863	!O3 tabulated coefficient
864	auxjo3(i) = jabsifotsintpar(auxi,7,indexint)
865	!Tabulated column
866	auxcoltab(i) = c34(nz2-i+1)
867	enddo
868	!Only if chemthermod.ge.2
869	if(chemthermod.ge.2) then
870	do i=1,nz2
871	!NO2 tabulated coefficient
872	auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
873	enddo
874	endif
875	call interfast
876	$ (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
877	do i=1,nlayermx
878	ind=auxind(i)
879	auxi = nlayermx-i+1
880	!O2 interpolated coefficient
881	jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
882	$ wp(i)*auxjo2(ind)
883	!H2O2 interpolated coefficient
884	jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
885	$ wp(i)*auxjh2o2(ind)
886	!O3 interpolated coefficient
887	jfotsout(indexint,7,auxi) = wm(i)*auxjo3(ind+1) +
888	$ wp(i)*auxjo3(ind)
889	enddo
890	!Only if chemthermod.ge.2
891	if(chemthermod.ge.2) then
892	do i=1,nlayermx
893	ind=auxind(i)
894	!NO2 interpolated coefficient
895	jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
896	$ wp(i)*auxjno2(ind)
897	enddo
898	endif
899
900	c End of interval 34
901
902
903	ccccccccccccccccccccccccccccccc
904	c 240.1-337.7nm (int 35)
905	c
906	c Absorption by:
907	c H2O2, O3, NO2
908	ccccccccccccccccccccccccccccccc
909
910	c Input atmospheric column
911
912	indexint=35
913	do i=1,nlayermx
914	auxcolinp(nlayermx-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
915	end do
916
917	c Interpolation
918
919	do i=1,nz2
920	auxi = nz2-i+1
921	!H2O2 tabulated coefficient
922	auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
923	!O3 tabulated coefficient
924	auxjo3(i) = jabsifotsintpar(auxi,7,indexint)
925	!Tabulated column
926	auxcoltab(i) = c35(auxi)
927	enddo
928	!Only if chemthermod.ge.2
929	if(chemthermod.ge.2) then
930	do i=1,nz2
931	!NO2 tabulated coefficient
932	auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
933	enddo
934	endif
935	call interfast
936	$ (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
937	do i=1,nlayermx
938	ind=auxind(i)
939	auxi = nlayermx-i+1
940	!H2O2 interpolated coefficient
941	jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
942	$ wp(i)*auxjh2o2(ind)
943	!O3 interpolated coefficient
944	jfotsout(indexint,7,auxi) = wm(i)*auxjo3(ind+1) +
945	$ wp(i)*auxjo3(ind)
946	enddo
947	if(chemthermod.ge.2) then
948	do i=1,nlayermx
949	ind=auxind(i)
950	!NO2 interpolated coefficient
951	jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
952	$ wp(i)*auxjno2(ind)
953	enddo
954	endif
955
956	c End of interval 35
957
958	ccccccccccccccccccccccccccccccc
959	c 337.8-800.0 nm (int 36)
960	c
961	c Absorption by:
962	c O3, NO2
963	ccccccccccccccccccccccccccccccc
964
965	c Input atmospheric column
966
967	indexint=36
968	do i=1,nlayermx
969	auxcolinp(nlayermx-i+1) = o3colx(i) + no2colx(i)
970	end do
971
972	c Interpolation
973
974	do i=1,nz2
975	auxi = nz2-i+1
976	!O3 tabulated coefficient
977	auxjo3(i) = jabsifotsintpar(auxi,7,indexint)
978	!Tabulated column
979	auxcoltab(i) = c36(auxi)
980	enddo
981	!Only if chemthermod.ge.2
982	if(chemthermod.ge.2) then
983	do i=1,nz2
984	!NO2 tabulated coefficient
985	auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
986	enddo
987	endif
988	call interfast
989	$ (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
990	do i=1,nlayermx
991	ind=auxind(i)
992	!O3 interpolated coefficient
993	jfotsout(indexint,7,nlayermx-i+1) = wm(i)*auxjo3(ind+1) +
994	$ wp(i)*auxjo3(ind)
995	enddo
996	!Only if chemthermod.ge.2
997	if(chemthermod.ge.2) then
998	do i=1,nlayermx
999	ind=auxind(i)
1000	!NO2 interpolated coefficient
1001	jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
1002	$ wp(i)*auxjno2(ind)
1003	enddo
1004	endif
1005
1006	c End of interval 36
1007
1008	c End of interpolation to obtain photoabsorption rates
1009
1010
1011	return
1012
1013	end
1014
1015
1016
1017	c**********************************************************************
1018	c**********************************************************************
1019
1020	subroutine column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
1021	$ co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx,
1022	$ n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
1023
1024	c nov 2002 fgg first version
1025
1026	c**********************************************************************
1027
1028	use tracer_mod, only: igcm_o, igcm_co2, igcm_o2, igcm_h2,
1029	& igcm_h2o_vap, igcm_h2o2, igcm_co, igcm_h,
1030	& igcm_o3, igcm_n2, igcm_n, igcm_no, igcm_no2,
1031	& mmol
1032	implicit none
1033
1034
1035	c common variables and constants
1036	include "dimensions.h"
1037	include "dimphys.h"
1038	! include "tracer.h"
1039	include 'param.h'
1040	include 'param_v4.h'
1041	include 'callkeys.h'
1042
1043
1044
1045	c local parameters and variables
1046
1047
1048
1049	c input and output variables
1050
1051	integer ig
1052	integer chemthermod
1053	integer nesptherm !# of species undergoing chemistry, input
1054	real rm(nlayermx,nesptherm) !densities (cm-3), input
1055	real tx(nlayermx) !temperature profile, input
1056	real iz(nlayermx+1) !height profile, input
1057	real zenit !SZA, input
1058	real co2colx(nlayermx) !column density of CO2 (cm^-2), output
1059	real o2colx(nlayermx) !column density of O2(cm^-2), output
1060	real o3pcolx(nlayermx) !column density of O(3P)(cm^-2), output
1061	real h2colx(nlayermx) !H2 column density (cm-2), output
1062	real h2ocolx(nlayermx) !H2O column density (cm-2), output
1063	real h2o2colx(nlayermx) !column density of H2O2(cm^-2), output
1064	real o3colx(nlayermx) !O3 column density (cm-2), output
1065	real n2colx(nlayermx) !N2 column density (cm-2), output
1066	real ncolx(nlayermx) !N column density (cm-2), output
1067	real nocolx(nlayermx) !NO column density (cm-2), output
1068	real cocolx(nlayermx) !CO column density (cm-2), output
1069	real hcolx(nlayermx) !H column density (cm-2), output
1070	real no2colx(nlayermx) !NO2 column density (cm-2), output
1071
1072
1073	c local variables
1074
1075	real xx
1076	real grav(nlayermx)
1077	real Hco2,Ho3p,Ho2,Hh2,Hh2o,Hh2o2
1078	real Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2
1079
1080	real co2x(nlayermx)
1081	real o2x(nlayermx)
1082	real o3px(nlayermx)
1083	real cox(nlayermx)
1084	real hx(nlayermx)
1085	real h2x(nlayermx)
1086	real h2ox(nlayermx)
1087	real h2o2x(nlayermx)
1088	real o3x(nlayermx)
1089	real n2x(nlayermx)
1090	real nx(nlayermx)
1091	real nox(nlayermx)
1092	real no2x(nlayermx)
1093
1094	integer i,j,k,icol,indexint !indexes
1095
1096	c variables for optical path calculation
1097
1098	integer nz3
1099	! parameter (nz3=nz*2)
1100
1101	integer jj
1102	real8 esp(nlayermx2)
1103	real8 ilayesp(nlayermx2)
1104	real8 szalayesp(nlayermx2)
1105	integer nlayesp
1106	real*8 zmini
1107	real*8 depth
1108	real*8 espco2, espo2, espo3p, esph2, esph2o, esph2o2,espo3
1109	real*8 espn2,espn,espno,espco,esph,espno2
1110	real*8 rcmnz, rcmmini
1111	real*8 szadeg
1112
1113
1114	! Tracer indexes in the thermospheric chemistry:
1115	!!! ATTENTION. These values have to be identical to those in chemthermos.F90
1116	!!! If the values are changed there, the same has to be done here !!!
1117	integer,parameter :: i_co2=1
1118	integer,parameter :: i_o2=2
1119	integer,parameter :: i_o=3
1120	integer,parameter :: i_co=4
1121	integer,parameter :: i_h=5
1122	integer,parameter :: i_h2=8
1123	integer,parameter :: i_h2o=9
1124	integer,parameter :: i_h2o2=10
1125	integer,parameter :: i_o3=12
1126	integer,parameter :: i_n2=13
1127	integer,parameter :: i_n=14
1128	integer,parameter :: i_no=15
1129	integer,parameter :: i_no2=17
1130
1131
1132	c***********************PROGRAM STARTS*****************************
1133
1134	nz3 = nlayermx*2
1135	do i=1,nlayermx
1136	xx = ( radio + iz(i) ) * 1.e5
1137	grav(i) = gg * masa /(xx**2)
1138	end do
1139
1140	!Scale heights
1141	xx = kboltzman * tx(nlayermx) * n_avog / grav(nlayermx)
1142	Ho3p = xx / mmol(igcm_o)
1143	Hco2 = xx / mmol(igcm_co2)
1144	Ho2 = xx / mmol(igcm_o2)
1145	Hh2 = xx / mmol(igcm_h2)
1146	Hh2o = xx / mmol(igcm_h2o_vap)
1147	Hh2o2 = xx / mmol(igcm_h2o2)
1148	Hco = xx / mmol(igcm_co)
1149	Hh = xx / mmol(igcm_h)
1150	!Only if O3 chem. required
1151	if(chemthermod.ge.1)
1152	$ Ho3 = xx / mmol(igcm_o3)
1153	!Only if N or ion chem.
1154	if(chemthermod.ge.2) then
1155	Hn2 = xx / mmol(igcm_n2)
1156	Hn = xx / mmol(igcm_n)
1157	Hno = xx / mmol(igcm_no)
1158	Hno2 = xx / mmol(igcm_no2)
1159	endif
1160	! first loop in altitude : initialisation
1161	do i=nlayermx,1,-1
1162	!Column initialisation
1163	co2colx(i) = 0.
1164	o2colx(i) = 0.
1165	o3pcolx(i) = 0.
1166	h2colx(i) = 0.
1167	h2ocolx(i) = 0.
1168	h2o2colx(i) = 0.
1169	o3colx(i) = 0.
1170	n2colx(i) = 0.
1171	ncolx(i) = 0.
1172	nocolx(i) = 0.
1173	cocolx(i) = 0.
1174	hcolx(i) = 0.
1175	no2colx(i) = 0.
1176	!Densities
1177	co2x(i) = rm(i,i_co2)
1178	o2x(i) = rm(i,i_o2)
1179	o3px(i) = rm(i,i_o)
1180	h2x(i) = rm(i,i_h2)
1181	h2ox(i) = rm(i,i_h2o)
1182	h2o2x(i) = rm(i,i_h2o2)
1183	cox(i) = rm(i,i_co)
1184	hx(i) = rm(i,i_h)
1185	!Only if O3 chem. required
1186	if(chemthermod.ge.1)
1187	$ o3x(i) = rm(i,i_o3)
1188	!Only if Nitrogen of ion chemistry requested
1189	if(chemthermod.ge.2) then
1190	n2x(i) = rm(i,i_n2)
1191	nx(i) = rm(i,i_n)
1192	nox(i) = rm(i,i_no)
1193	no2x(i) = rm(i,i_no2)
1194	endif
1195	enddo
1196	! second loop in altitude : column calculations
1197	do i=nlayermx,1,-1
1198	!Routine to calculate the geometrical length of each layer
1199	call espesor_optico_A(ig,i,zenit,iz(i),nz3,iz,esp,ilayesp,
1200	$ szalayesp,nlayesp, zmini)
1201	if(ilayesp(nlayesp).eq.-1) then
1202	co2colx(i)=1.e25
1203	o2colx(i)=1.e25
1204	o3pcolx(i)=1.e25
1205	h2colx(i)=1.e25
1206	h2ocolx(i)=1.e25
1207	h2o2colx(i)=1.e25
1208	o3colx(i)=1.e25
1209	n2colx(i)=1.e25
1210	ncolx(i)=1.e25
1211	nocolx(i)=1.e25
1212	cocolx(i)=1.e25
1213	hcolx(i)=1.e25
1214	no2colx(i)=1.e25
1215	else
1216	rcmnz = ( radio + iz(nlayermx) ) * 1.e5
1217	rcmmini = ( radio + zmini ) * 1.e5
1218	!Column calculation taking into account the geometrical depth
1219	!calculated before
1220	do j=1,nlayesp
1221	jj=ilayesp(j)
1222	!Top layer
1223	if(jj.eq.nlayermx) then
1224	if(zenit.le.60.) then
1225	o3pcolx(i)=o3pcolx(i)+o3px(nlayermx)Ho3pesp(j)
1226	$ *1.e-5
1227	co2colx(i)=co2colx(i)+co2x(nlayermx)Hco2esp(j)
1228	$ *1.e-5
1229	h2o2colx(i)=h2o2colx(i)+
1230	$ h2o2x(nlayermx)Hh2o2esp(j)*1.e-5
1231	o2colx(i)=o2colx(i)+o2x(nlayermx)Ho2esp(j)
1232	$ *1.e-5
1233	h2colx(i)=h2colx(i)+h2x(nlayermx)Hh2esp(j)
1234	$ *1.e-5
1235	h2ocolx(i)=h2ocolx(i)+h2ox(nlayermx)Hh2oesp(j)
1236	$ *1.e-5
1237	cocolx(i)=cocolx(i)+cox(nlayermx)Hcoesp(j)
1238	$ *1.e-5
1239	hcolx(i)=hcolx(i)+hx(nlayermx)Hhesp(j)
1240	$ *1.e-5
1241	!Only if O3 chemistry required
1242	if(chemthermod.ge.1) o3colx(i)=
1243	$ o3colx(i)+o3x(nlayermx)Ho3esp(j)
1244	$ *1.e-5
1245	!Only if N or ion chemistry requested
1246	if(chemthermod.ge.2) then
1247	n2colx(i)=n2colx(i)+n2x(nlayermx)Hn2esp(j)
1248	$ *1.e-5
1249	ncolx(i)=ncolx(i)+nx(nlayermx)Hnesp(j)
1250	$ *1.e-5
1251	nocolx(i)=nocolx(i)+nox(nlayermx)Hnoesp(j)
1252	$ *1.e-5
1253	no2colx(i)=no2colx(i)+no2x(nlayermx)Hno2esp(j)
1254	$ *1.e-5
1255	endif
1256	else if(zenit.gt.60.) then
1257	espco2 =sqrt((rcmnz+Hco2)2 -rcmmini2) - esp(j)
1258	espo2 = sqrt((rcmnz+Ho2)2 -rcmmini2) - esp(j)
1259	espo3p = sqrt((rcmnz+Ho3p)2 -rcmmini2)- esp(j)
1260	esph2 = sqrt((rcmnz+Hh2)2 -rcmmini2) - esp(j)
1261	esph2o = sqrt((rcmnz+Hh2o)2 -rcmmini2)- esp(j)
1262	esph2o2= sqrt((rcmnz+Hh2o2)2-rcmmini2)- esp(j)
1263	espco = sqrt((rcmnz+Hco)2 -rcmmini2) - esp(j)
1264	esph = sqrt((rcmnz+Hh)2 -rcmmini2) - esp(j)
1265	!Only if O3 chemistry required
1266	if(chemthermod.ge.1)
1267	$ espo3=sqrt((rcmnz+Ho3)2-rcmmini2)-esp(j)
1268	!Only if N or ion chemistry requested
1269	if(chemthermod.ge.2) then
1270	espn2 =sqrt((rcmnz+Hn2)2-rcmmini2)-esp(j)
1271	espn =sqrt((rcmnz+Hn)2-rcmmini2) - esp(j)
1272	espno =sqrt((rcmnz+Hno)2-rcmmini2) - esp(j)
1273	espno2=sqrt((rcmnz+Hno2)2-rcmmini2)- esp(j)
1274	endif
1275	co2colx(i) = co2colx(i) + espco2*co2x(nlayermx)
1276	o2colx(i) = o2colx(i) + espo2*o2x(nlayermx)
1277	o3pcolx(i) = o3pcolx(i) + espo3p*o3px(nlayermx)
1278	h2colx(i) = h2colx(i) + esph2*h2x(nlayermx)
1279	h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(nlayermx)
1280	h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(nlayermx)
1281	cocolx(i) = cocolx(i) + espco*cox(nlayermx)
1282	hcolx(i) = hcolx(i) + esph*hx(nlayermx)
1283	!Only if O3 chemistry required
1284	if(chemthermod.ge.1)
1285	$ o3colx(i) = o3colx(i) + espo3*o3x(nlayermx)
1286	!Only if N or ion chemistry requested
1287	if(chemthermod.ge.2) then
1288	n2colx(i) = n2colx(i) + espn2*n2x(nlayermx)
1289	ncolx(i) = ncolx(i) + espn*nx(nlayermx)
1290	nocolx(i) = nocolx(i) + espno*nox(nlayermx)
1291	no2colx(i) = no2colx(i) + espno2*no2x(nlayermx)
1292	endif
1293	endif !Of if zenit.lt.60
1294	!Other layers
1295	else
1296	co2colx(i) = co2colx(i) +
1297	$ esp(j) * (co2x(jj)+co2x(jj+1)) / 2.
1298	o2colx(i) = o2colx(i) +
1299	$ esp(j) * (o2x(jj)+o2x(jj+1)) / 2.
1300	o3pcolx(i) = o3pcolx(i) +
1301	$ esp(j) * (o3px(jj)+o3px(jj+1)) / 2.
1302	h2colx(i) = h2colx(i) +
1303	$ esp(j) * (h2x(jj)+h2x(jj+1)) / 2.
1304	h2ocolx(i) = h2ocolx(i) +
1305	$ esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2.
1306	h2o2colx(i) = h2o2colx(i) +
1307	$ esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2.
1308	cocolx(i) = cocolx(i) +
1309	$ esp(j) * (cox(jj)+cox(jj+1)) / 2.
1310	hcolx(i) = hcolx(i) +
1311	$ esp(j) * (hx(jj)+hx(jj+1)) / 2.
1312	!Only if O3 chemistry required
1313	if(chemthermod.ge.1)
1314	$ o3colx(i) = o3colx(i) +
1315	$ esp(j) * (o3x(jj)+o3x(jj+1)) / 2.
1316	!Only if N or ion chemistry requested
1317	if(chemthermod.ge.2) then
1318	n2colx(i) = n2colx(i) +
1319	$ esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
1320	ncolx(i) = ncolx(i) +
1321	$ esp(j) * (nx(jj)+nx(jj+1)) / 2.
1322	nocolx(i) = nocolx(i) +
1323	$ esp(j) * (nox(jj)+nox(jj+1)) / 2.
1324	no2colx(i) = no2colx(i) +
1325	$ esp(j) * (no2x(jj)+no2x(jj+1)) / 2.
1326	endif
1327	endif !Of if jj.eq.nlayermx
1328	end do !Of do j=1,nlayesp
1329	endif !Of ilayesp(nlayesp).eq.-1
1330	enddo !Of do i=nlayermx,1,-1
1331	return
1332
1333
1334	end
1335
1336
1337	c**********************************************************************
1338	c**********************************************************************
1339
1340	subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup)
1341	C
1342	C subroutine to perform linear interpolation in pressure from 1D profile
1343	C escin(nl) sampled on pressure grid pin(nl) to profile
1344	C escout(nlayer) on pressure grid p(nlayer).
1345	C
1346	real*8 wm(nlayer),wp(nlayer),p(nlayer)
1347	integer nm(nlayer)
1348	real*8 pin(nl)
1349	real*8 limup,limdown
1350	integer nl,nlayer,n1,n,np,nini
1351	nini=1
1352	do n1=1,nlayer
1353	if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then
1354	wm(n1) = 0.d0
1355	wp(n1) = 0.d0
1356	else
1357	do n = nini,nl-1
1358	if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
1359	nm(n1)=n
1360	np=n+1
1361	wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n))
1362	wp(n1)=1.d0 - wm(n1)
1363	nini = n
1364	exit
1365	endif
1366	enddo
1367	endif
1368	enddo
1369	return
1370	end
1371
1372
1373	c**********************************************************************
1374	c**********************************************************************
1375
1376	subroutine espesor_optico_A (ig,capa, szadeg,z,
1377	@ nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
1378
1379	c fgg nov 03 Adaptation to Martian model
1380	c malv jul 03 Corrected z grid. Split in alt & frec codes
1381	c fgg feb 03 first version
1382	*************************************************************************
1383
1384	implicit none
1385
1386
1387	c common variables and constants
1388
1389	include "dimensions.h"
1390	include "dimphys.h"
1391	include 'param.h'
1392	include 'param_v4.h'
1393
1394	c arguments
1395
1396	real szadeg ! I. SZA [rad]
1397	real z ! I. altitude of interest [km]
1398	integer nz3,ig ! I. dimension of esp, ylayesp, etc...
1399	! (=2*nlayermx= max# of layers in ray path)
1400	real iz(nlayermx+1) ! I. Altitude of each layer
1401	real*8 esp(nz3) ! O. layer widths after geometrically
1402	! amplified; in [cm] except at TOA
1403	! where an auxiliary value is used
1404	real*8 ilayesp(nz3) ! O. Indexes of layers along ray path
1405	real*8 szalayesp(nz3) ! O. SZA [deg] " " "
1406	integer nlayesp
1407	! real*8 nlayesp ! O. # layers along ray path at this z
1408	real*8 zmini ! O. Minimum altitud of ray path [km]
1409
1410
1411	c local variables and constants
1412
1413	integer j,i,capa
1414	integer jmin ! index of min.altitude along ray path
1415	real*8 szarad ! SZA [deg]
1416	real*8 zz
1417	real*8 diz(nlayermx+1)
1418	real*8 rkmnz ! distance TOA to center of Planet [km]
1419	real*8 rkmmini ! distance zmini to center of P [km]
1420	real*8 rkmj ! intermediate distance to C of P [km]
1421
1422	c external function
1423	external grid_R8 ! Returns index of layer containing the altitude
1424	! of interest, z; for example, if
1425	! zkm(i)=z or zkm(i)<z<zkm(i+1) => grid(z)=i
1426	integer grid_R8
1427
1428	*************************************************************************
1429	szarad = dble(szadeg)*3.141592d0/180.d0
1430	zz=dble(z)
1431	do i=1,nlayermx
1432	diz(i)=dble(iz(i))
1433	enddo
1434	do j=1,nz3
1435	esp(j) = 0.d0
1436	szalayesp(j) = 777.d0
1437	ilayesp(j) = 0
1438	enddo
1439	nlayesp = 0
1440
1441	! First case: szadeg<60
1442	! The optical thickness will be given by 1/cos(sza)
1443	! We deal with 2 different regions:
1444	! 1: First, all layers between z and ztop ("upper part of ray")
1445	! 2: Second, the layer at ztop
1446	if(szadeg.lt.60.d0) then
1447
1448	zmini = zz
1449	if(abs(zz-diz(nlayermx)).lt.1.d-3) goto 1357
1450	! 1st Zone: Upper part of ray
1451	!
1452	do j=grid_R8(zz,diz,nlayermx),nlayermx-1
1453	nlayesp = nlayesp + 1
1454	ilayesp(nlayesp) = j
1455	esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad) ! [km]
1456	esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm]
1457	szalayesp(nlayesp) = szadeg
1458	end do
1459
1460	!
1461	! 2nd Zone: Top layer
1462	1357 continue
1463	nlayesp = nlayesp + 1
1464	ilayesp(nlayesp) = nlayermx
1465	esp(nlayesp) = 1.d0 / cos(szarad) ! aux. non-dimens. factor
1466	szalayesp(nlayesp) = szadeg
1467
1468
1469	! Second case: 60 < szadeg < 90
1470	! The optical thickness is evaluated.
1471	! (the magnitude of the effect of not using cos(sza) is 3.e-5
1472	! for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately)
1473	! We deal with 2 different regions:
1474	! 1: First, all layers between z and ztop ("upper part of ray")
1475	! 2: Second, the layer at ztop ("uppermost layer")
1476	else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then
1477
1478	zmini=(radio+zz)*sin(szarad)-radio
1479	rkmmini = radio + zmini
1480
1481	if(abs(zz-diz(nlayermx)).lt.1.d-4) goto 1470
1482
1483	! 1st Zone: Upper part of ray
1484	!
1485	do j=grid_R8(zz,diz,nlayermx),nlayermx-1
1486	nlayesp = nlayesp + 1
1487	ilayesp(nlayesp) = j
1488	esp(nlayesp) =
1489	# sqrt( (radio+diz(j+1))2 - rkmmini2 ) -
1490	# sqrt( (radio+diz(j))2 - rkmmini2 ) ! [km]
1491	esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm]
1492	rkmj = radio+diz(j)
1493	szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
1494	szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 ! [deg]
1495	end do
1496	1470 continue
1497	! 2nd Zone: Uppermost layer of ray.
1498	!
1499	nlayesp = nlayesp + 1
1500	ilayesp(nlayesp) = nlayermx
1501	rkmnz = radio+diz(nlayermx)
1502	esp(nlayesp) = sqrt( rkmnz2 - rkmmini2 ) ! aux.factor[km]
1503	esp(nlayesp) = esp(nlayesp) * 1.d5 ! aux.f. [cm]
1504	szalayesp(nlayesp) = asin( rkmmini/rkmnz ) ! [rad]
1505	szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
1506
1507
1508	! Third case: szadeg > 90
1509	! The optical thickness is evaluated.
1510	! We deal with 5 different regions:
1511	! 1: all layers between z and ztop ("upper part of ray")
1512	! 2: the layer at ztop ("uppermost layer")
1513	! 3: the lowest layer, at zmini
1514	! 4: the layers increasing from zmini to z (here SZA<90)
1515	! 5: the layers decreasing from z to zmini (here SZA>90)
1516	else if(szadeg.gt.90.d0) then
1517
1518	zmini=(radio+zz)*sin(szarad)-radio
1519	rkmmini = radio + zmini
1520
1521	if(zmini.lt.diz(1)) then ! Can see the sun? No => esp(j)=inft
1522	nlayesp = nlayesp + 1
1523	ilayesp(nlayesp) = - 1 ! Value to mark "no sun on view"
1524	! esp(nlayesp) = 1.e30
1525
1526	else
1527	jmin=grid_R8(zmini,diz,nlayermx)+1
1528
1529
1530	if(abs(zz-diz(nlayermx)).lt.1.d-4) goto 9876
1531
1532	! 1st Zone: Upper part of ray
1533	!
1534	do j=grid_R8(zz,diz,nlayermx),nlayermx-1
1535	nlayesp = nlayesp + 1
1536	ilayesp(nlayesp) = j
1537	esp(nlayesp) =
1538	$ sqrt( (radio+diz(j+1))2 - rkmmini2 ) -
1539	$ sqrt( (radio+diz(j))2 - rkmmini2 ) ! [km]
1540	esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm]
1541	rkmj = radio+diz(j)
1542	szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
1543	szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1544	end do
1545
1546	9876 continue
1547	! 2nd Zone: Uppermost layer of ray.
1548	!
1549	nlayesp = nlayesp + 1
1550	ilayesp(nlayesp) = nlayermx
1551	rkmnz = radio+diz(nlayermx)
1552	esp(nlayesp) = sqrt( rkmnz2 - rkmmini2 ) !aux.factor[km]
1553	esp(nlayesp) = esp(nlayesp) * 1.d5 !aux.f.[cm]
1554	szalayesp(nlayesp) = asin( rkmmini/rkmnz ) ! [rad]
1555	szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1556
1557	! 3er Zone: Lowestmost layer of ray
1558	!
1559	if ( jmin .ge. 2 ) then ! If above the planet's surface
1560	j=jmin-1
1561	nlayesp = nlayesp + 1
1562	ilayesp(nlayesp) = j
1563	esp(nlayesp) = 2. *
1564	$ sqrt( (radio+diz(j+1))2 -rkmmini2 ) ! [km]
1565	esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm]
1566	rkmj = radio+diz(j+1)
1567	szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
1568	szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1569	endif
1570
1571	! 4th zone: Lower part of ray, increasing from zmin to z
1572	! ( layers with SZA < 90 deg )
1573	do j=jmin,grid_R8(zz,diz,nlayermx)-1
1574	nlayesp = nlayesp + 1
1575	ilayesp(nlayesp) = j
1576	esp(nlayesp) =
1577	$ sqrt( (radio+diz(j+1))2 - rkmmini2 )
1578	$ - sqrt( (radio+diz(j))2 - rkmmini2 ) ! [km]
1579	esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm]
1580	rkmj = radio+diz(j)
1581	szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
1582	szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1583	end do
1584
1585	! 5th zone: Lower part of ray, decreasing from z to zmin
1586	! ( layers with SZA > 90 deg )
1587	do j=grid_R8(zz,diz,nlayermx)-1, jmin, -1
1588	nlayesp = nlayesp + 1
1589	ilayesp(nlayesp) = j
1590	esp(nlayesp) =
1591	$ sqrt( (radio+diz(j+1))2 - rkmmini2 )
1592	$ - sqrt( (radio+diz(j))2 - rkmmini2 ) ! [km]
1593	esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm]
1594	rkmj = radio+diz(j)
1595	szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj ) ! [rad]
1596	szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 ! [deg]
1597	end do
1598
1599	end if
1600
1601	end if
1602
1603	return
1604
1605	end
1606
1607
1608
1609	c**********************************************************************
1610	c***********************************************************************
1611
1612	function grid_R8 (z, zgrid, nz)
1613
1614	c Returns the index where z is located within vector zgrid
1615	c The vector zgrid must be monotonously increasing, otherwise program stops.
1616	c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops.
1617	c
1618	c FGG Aug-2004 Correct z.lt.zgrid(i) to .le.
1619	c MALV Jul-2003
1620	c***********************************************************************
1621
1622	implicit none
1623
1624	c Arguments
1625	integer nz
1626	real*8 z
1627	real*8 zgrid(nz)
1628	integer grid_R8
1629
1630	c Local
1631	integer i, nz1, nznew
1632
1633	c*** CODE START
1634
1635	if ( z .lt. zgrid(1) .or. z.gt.zgrid(nz) ) then
1636	write (,) ' GRID/ z outside bounds of zgrid '
1637	write (,) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
1638	stop ' Serious error in GRID.F '
1639	endif
1640	if ( nz .lt. 2 ) then
1641	write (,) ' GRID/ zgrid needs 2 points at least ! '
1642	stop ' Serious error in GRID.F '
1643	endif
1644	if ( zgrid(1) .ge. zgrid(nz) ) then
1645	write (,) ' GRID/ zgrid must increase with index'
1646	stop ' Serious error in GRID.F '
1647	endif
1648
1649	nz1 = 1
1650	nznew = nz/2
1651	if ( z .gt. zgrid(nznew) ) then
1652	nz1 = nznew
1653	nznew = nz
1654	endif
1655	do i=nz1+1,nznew
1656	if ( z. eq. zgrid(i) ) then
1657	grid_R8=i
1658	return
1659	elseif ( z .le. zgrid(i) ) then
1660	grid_R8 = i-1
1661	return
1662	endif
1663	enddo
1664	grid_R8 = nz
1665	return
1666
1667
1668
1669	end
1670
1671
1672
1673	!c***************************************************
1674	!c***************************************************
1675
1676	subroutine flujo(date)
1677
1678
1679	!c fgg nov 2002 first version
1680	!c***************************************************
1681
1682	use comsaison_h, only: dist_sol
1683	implicit none
1684
1685
1686	! common variables and constants
1687	include "dimensions.h"
1688	include "dimphys.h"
1689	! include "comsaison.h"
1690	include 'param.h'
1691	include 'param_v4.h'
1692	include "callkeys.h"
1693
1694
1695	! Arguments
1696
1697	real date
1698
1699
1700	! Local variable and constants
1701
1702	integer i
1703	integer inter
1704	real nada
1705
1706	!c*************************************************
1707
1708	if(date.lt.1985.) date=1985.
1709	if(date.gt.2001.) date=2001.
1710
1711	do i=1,ninter
1712	fluxtop(i)=1.
1713	!Variation of solar flux with 11 years solar cycle
1714	!For more details, see Gonzalez-Galindo et al. 2005
1715	!To be improved in next versions
1716	if(i.le.24.and.solvarmod.eq.0) then
1717	fluxtop(i)=(((ct1(i)+p1(i)*date)/2.)
1718	$ sin(2.3.1416/11.*(date-1985.-3.1416))
1719	$ +(ct2(i)+p2(i)date)+1.)fluxtop(i)
1720	end if
1721	fluxtop(i)=fluxtop(i)(1.52/dist_sol)*2
1722	end do
1723
1724	return
1725	end

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

Download in other formats:

Original Format