Context Navigation

jthermcalc_e107.F @ 3656

Last change on this file since 3656 was 2836, checked in by abierjon, 3 years ago

VENUS GCM:

INCLUDING THE IONOSPHERE CODE IN VENUS GCM

ATTENTION: INCREASED MEMORY DEMAND

NEEDS AT LEAST 135 GB of allocated memory

==============================================================
===> LIST OF APPENDED FILES AND INTERNAL ADDITIONS <==========
==============================================================

NEW MODEL SPECIES

deftank/traceur-chemistry-IONOSPHERE.def

Neutrals: 4 Species: N, NO, NO2, N(2D)

Ions: 15 species:

CO2+, CO+, O+, O2+, H+,

N2+, H2O+, OH+, C+, HCO+,

H3O+, HCO2+, N+, NO+, elec

FROM 36 to 55 chemical species

NEW KEYWORD OF PHYSIQ.DEF

deftank/physiq-96x96x90-chemistry-IONOSPHERE.def

ok_ionchem: keyword supposed to activate ion chemistry.

(be careful that n, no, no2, n2d and the ion species are in the deftank/traceur-chemistry-IONOSPHERE.def)

ok_jonline: keyword supposed to activate the online photochemistry

(be careful that n, no, no2 and n2d are in the traceur-chemistry-IONOSPHERE.def)

==============================================================
===> LIST OF MODIFIED PROGRAMS <=========================
==============================================================

nonoro_gwd_ran_mod.F90

Change EPFLUXMAX value from 5.E-3 to 1.E-3

photochemistry_venus.F90 (krates; photolysis_online; indices;)

Import of keywords: ok_ionchem, tuneupperatm & ok_jonline
addition of ion species in order
Forcing electroneutrality
Update of the reactions a001 and a002 with taking into account the other species
Change of formula for a002 with the formulation of Baulch et al., 1976 (confirmed by Smith and Robertson, 2008)

photolysis_mod.F90

Modification of the rdsolarflux subroutine to include interpolation with Atlas1 and Atlas3 in connection with E10.7

concentration2.F90

Added chemical species n2d, no, no2 and n in the conductivity calculation.

The ions have been excluded because their sum is 10^{5 times less dense than the neutrals and

their thermal conductivity is unknown}

iono_h.F90

Addition of the phdisrate routine
replace the electronic temperature of Mars by that of

origin = 1: Theis et al. 1980 (Venus) with bilinear interpolation altitude/cos(SZA)
origin = 2: Theis et al. 1984 (Venus) with the formula of the electronic temperature at high solar activity

addition of an ion temperature model based on VIRA

cleshphys.h

added ok_ionchem & ok_jonline in COMMON/clesphys_l/

conf_phys.f90

add ok_ionchem & ok_jonline as parameters to read from physiq.def file set to .false. by default

chemparam_mod.F90

Add species in M_tr and corresponding i_X. Set all i_X to zero before reading traceur-chemistry-IONOSPHERE.def
Added Type_tr table to differentiate species: 1 == neutral, 2 == ION, 3 == liquid, 10 == others

euvheat.F90; hrtherm.F; jthermcalc_e107.F; param_read_e107.F

Normalization with Mars

A.M

File size: 61.7 KB

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

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

Download in other formats:

Original Format