Context Navigation

jthermcalc_e107.F @ 3567

Last change on this file since 3567 was 2836, checked in by abierjon, 2 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

Rev	Line
[1310]	1	c**********************************************************************
	2
	3	subroutine jthermcalc_e107
[2464]	4	$ (ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit)
[1310]	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**********************************************************************
[2464]	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
[1310]	21	implicit none
	22
[2464]	23	include "clesphys.h"
[1310]	24
[2464]	25
[1310]	26	c input and output variables
	27
[2464]	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
[1310]	35
[2464]	36	! NB: no output variable! (computed jfotsout() is in module param_v4_h)
[1310]	37
	38	c local parameters and variables
	39
[2464]	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)
[1310]	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)
[2464]	82	real*8 wp(nlayer),wm(nlayer)
	83	real*8 auxcolinp(nlayer)
	84	integer auxind(nlayer)
[1310]	85	integer auxi
	86	integer ind
[2464]	87	real*8 cortemp(nlayer)
[1310]	88
	89	real*8 limdown !limits for interpolation
	90	real*8 limup ! ""
	91
[2836]	92	!!!ATTENTION. Here ix_co2 has to have the same value than in euvheat.F90
[1310]	93	!!!If the value is changed there, if has to be changed also here !!!!
[2836]	94	integer,parameter :: ix_co2=1
[1310]	95
[2464]	96	character*20 modname
	97	character*80 abort_message
[1310]	98
	99	c***********************PROGRAM STARTS*****************************
	100
[2464]	101	modname = 'jthermcalc_e107'
	102
[1310]	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
[2464]	117	do i=1,nlayer
[1310]	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,
[2836]	125	$ co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx,
	126	$ n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
[1310]	127
	128	!Auxiliar column to include the temperature dependence
	129	!of CO2 cross section
[2464]	130	coltemp(nlayer)=co2colx(nlayer)*abs(t2(nlayer)-t0(nlayer))
	131	do i=nlayer-1,1,-1
[1310]	132	coltemp(i)=!coltemp(i+1)+ PQ SE ELIMINA? REVISAR
[2836]	133	$ ( rm(i,ix_co2) + rm(i+1,ix_co2) ) * 0.5
[1310]	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)
[2464]	138	do i=1,nlayer
[1310]	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
[2464]	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)
[1310]	152
	153	!Photoabsorption coefficients at TOA as a function of E10.7
	154	do j=1,nabs
	155	do indexint=1,ninter
[2464]	156	jfotsout(indexint,j,nlayer)=coefit0(indexint,j)+
[1310]	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
[2464]	178	do i=1,nlayer
	179	auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint) +
[1310]	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
[2464]	213	$ (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
	214	do i=1,nlayer
[1310]	215	ind=auxind(i)
[2464]	216	auxi=nlayer-i+1
[1310]	217	!CO2 interpolated coefficient
[2464]	218	jfotsout(indexint,1,auxi) = jfotsout(indexint,1,nlayer) *
[1310]	219	$ (wm(i)auxjco2(ind+1) + wp(i)auxjco2(ind))
	220	!O2 interpolated coefficient
[2464]	221	jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
[1310]	222	$ (wm(i)auxjo2(ind+1) + wp(i)auxjo2(ind))
	223	!O3p interpolated coefficient
[2464]	224	jfotsout(indexint,3,auxi) = jfotsout(indexint,3,nlayer) *
[1310]	225	$ (wm(i)auxjo3p(ind+1) + wp(i)auxjo3p(ind))
	226	!H2 interpolated coefficient
[2809]	227	jfotsout(indexint,5,auxi) = jfotsout(indexint,5,nlayer) *
[1310]	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
[2464]	233	do i=1,nlayer
[1310]	234	ind=auxind(i)
[2464]	235	jfotsout(indexint,9,nlayer-i+1) =
	236	$ jfotsout(indexint,9,nlayer) *
[1310]	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
[2464]	255	do i=1,nlayer
	256	auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint)+
[1310]	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
[2464]	302	call interfast(wm,wp,auxind,auxcolinp,nlayer,
[1310]	303	$ auxcoltab,nz2,limdown,limup)
[2464]	304	do i=1,nlayer
[1310]	305	ind=auxind(i)
[2464]	306	auxi = nlayer-i+1
[1310]	307	!O2 interpolated coefficient
	308	jfotsout(indexint,2,auxi) =
[2464]	309	$ jfotsout(indexint,2,nlayer) *
[1310]	310	$ (wm(i)auxjo2(ind+1) + wp(i)auxjo2(ind))
	311	!O3p interpolated coefficient
	312	jfotsout(indexint,3,auxi) =
[2464]	313	$ jfotsout(indexint,3,nlayer) *
[1310]	314	$ (wm(i)auxjo3p(ind+1) + wp(i)auxjo3p(ind))
	315	!CO2 interpolated coefficient
	316	jfotsout(indexint,1,auxi) =
[2464]	317	$ jfotsout(indexint,1,nlayer) *
[1310]	318	$ (wm(i)auxjco2(ind+1) + wp(i)auxjco2(ind))
	319	!H2 interpolated coefficient
	320	jfotsout(indexint,5,auxi) =
[2464]	321	$ jfotsout(indexint,5,nlayer) *
[1310]	322	$ (wm(i)auxjh2(ind+1) + wp(i)auxjh2(ind))
	323	!N2 interpolated coefficient
	324	jfotsout(indexint,8,auxi) =
[2464]	325	$ jfotsout(indexint,8,nlayer) *
[1310]	326	$ (wm(i)auxjn2(ind+1) + wp(i)auxjn2(ind))
	327	!CO interpolated coefficient
	328	jfotsout(indexint,11,auxi) =
[2464]	329	$ jfotsout(indexint,11,nlayer) *
[1310]	330	$ (wm(i)auxjco(ind+1) + wp(i)auxjco(ind))
	331	!H interpolated coefficient
	332	jfotsout(indexint,12,auxi) =
[2464]	333	$ jfotsout(indexint,12,nlayer) *
[1310]	334	$ (wm(i)auxjh(ind+1) + wp(i)auxjh(ind))
	335	enddo
	336	!Only if chemthermod.ge.2
	337	if(chemthermod.ge.2) then
[2464]	338	do i=1,nlayer
[1310]	339	ind=auxind(i)
[2464]	340	auxi = nlayer-i+1
[1310]	341	!N interpolated coefficient
	342	jfotsout(indexint,9,auxi) =
[2464]	343	$ jfotsout(indexint,9,nlayer) *
[1310]	344	$ (wm(i)auxjn(ind+1) + wp(i)auxjn(ind))
	345	!NO interpolated coefficient
	346	jfotsout(indexint,10,auxi)=
[2464]	347	$ jfotsout(indexint,10,nlayer) *
[1310]	348	$ (wm(i)auxjno(ind+1) + wp(i)auxjno(ind))
	349	!NO2 interpolated coefficient
	350	jfotsout(indexint,13,auxi)=
[2464]	351	$ jfotsout(indexint,13,nlayer) *
[1310]	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
[2464]	370	do i=1,nlayer
	371	auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint)+
[1310]	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
[2464]	417	$ (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
	418	do i=1,nlayer
[1310]	419	ind=auxind(i)
[2464]	420	auxi = nlayer-i+1
[1310]	421	!O2 interpolated coefficient
[2464]	422	jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
[1310]	423	$ (wm(i)auxjo2(ind+1) + wp(i)auxjo2(ind))
	424	!CO2 interpolated coefficient
[2464]	425	jfotsout(indexint,1,auxi) = jfotsout(indexint,1,nlayer) *
[1310]	426	$ (wm(i)auxjco2(ind+1) + wp(i)auxjco2(ind))
	427	!O3p interpolated coefficient
[2464]	428	jfotsout(indexint,3,auxi) = jfotsout(indexint,3,nlayer) *
[1310]	429	$ (wm(i)auxjo3p(ind+1) + wp(i)auxjo3p(ind))
	430	!N2 interpolated coefficient
[2464]	431	jfotsout(indexint,8,auxi) = jfotsout(indexint,8,nlayer) *
[1310]	432	$ (wm(i)auxjn2(ind+1) + wp(i)auxjn2(ind))
	433	!CO interpolated coefficient
	434	jfotsout(indexint,11,auxi) =
[2464]	435	$ jfotsout(indexint,11,nlayer) *
[1310]	436	$ (wm(i)auxjco(ind+1) + wp(i)auxjco(ind))
	437	!H interpolated coefficient
	438	jfotsout(indexint,12,auxi) =
[2464]	439	$ jfotsout(indexint,12,nlayer) *
[1310]	440	$ (wm(i)auxjh(ind+1) + wp(i)auxjh(ind))
	441	enddo
	442	!Only if chemthermod.ge.2
	443	if(chemthermod.ge.2) then
[2464]	444	do i=1,nlayer
[1310]	445	ind=auxind(i)
[2464]	446	auxi = nlayer-i+1
[1310]	447	!N interpolated coefficient
	448	jfotsout(indexint,9,auxi) =
[2464]	449	$ jfotsout(indexint,9,nlayer) *
[1310]	450	$ (wm(i)auxjn(ind+1) + wp(i)auxjn(ind))
	451	!NO interpolated coefficient
	452	jfotsout(indexint,10,auxi) =
[2464]	453	$ jfotsout(indexint,10,nlayer) *
[1310]	454	$ (wm(i)auxjno(ind+1) + wp(i)auxjno(ind))
	455	!NO2 interpolated coefficient
	456	jfotsout(indexint,13,auxi) =
[2464]	457	$ jfotsout(indexint,13,nlayer) *
[1310]	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
[2464]	473	do i=1,nlayer
	474	auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
[1310]	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
[2464]	507	$ (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
[1310]	508	!Correction to include T variation of CO2 cross section
	509	if(indexint.eq.24) then
[2464]	510	do i=1,nlayer
	511	auxi = nlayer-i+1
[1310]	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
[2464]	522	do i=1,nlayer
[1310]	523	cortemp(i)=1.
	524	enddo
	525	end if
[2464]	526	do i=1,nlayer
[1310]	527	ind=auxind(i)
[2464]	528	auxi = nlayer-i+1
[1310]	529	!O2 interpolated coefficient
	530	jfotsout(indexint,2,auxi) =
[2464]	531	$ jfotsout(indexint,2,nlayer) *
[1310]	532	$ (wm(i)auxjo2(ind+1) + wp(i)auxjo2(ind)) *
	533	$ cortemp(i)
	534	!CO2 interpolated coefficient
	535	jfotsout(indexint,1,auxi) =
[2464]	536	$ jfotsout(indexint,1,nlayer) *
[1310]	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) =
[2464]	545	$ jfotsout(indexint,8,nlayer) *
[1310]	546	$ (wm(i)auxjn2(ind+1) + wp(i)auxjn2(ind)) *
	547	$ cortemp(i)
	548	!CO interpolated coefficient
	549	jfotsout(indexint,11,auxi) =
[2464]	550	$ jfotsout(indexint,11,nlayer) *
[1310]	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
[2464]	556	do i=1,nlayer
[1310]	557	ind=auxind(i)
[2464]	558	auxi = nlayer-i+1
[1310]	559	!NO interpolated coefficient
	560	jfotsout(indexint,10,auxi)=
[2464]	561	$ jfotsout(indexint,10,nlayer) *
[1310]	562	$ (wm(i)auxjno(ind+1) + wp(i)auxjno(ind)) *
	563	$ cortemp(i)
	564	!NO2 interpolated coefficient
	565	jfotsout(indexint,13,auxi)=
[2464]	566	$ jfotsout(indexint,13,nlayer) *
[1310]	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
[2464]	585	do i=1,nlayer
	586	auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
[1310]	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
[2464]	620	$ (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
	621	do i=1,nlayer
[1310]	622	ind=auxind(i)
[2464]	623	auxi = nlayer-i+1
[1310]	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) =
[2464]	634	$ jfotsout(indexint,1,nlayer) *
[1310]	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) =
[2464]	641	$ jfotsout(indexint,2,nlayer) *
[1310]	642	$ (wm(i)auxjo2(ind+1) + wp(i)auxjo2(ind)) *
	643	$ cortemp(i)
	644	!H2O interpolated coefficient
	645	jfotsout(indexint,4,auxi) =
[2464]	646	$ jfotsout(indexint,4,nlayer) *
[1310]	647	$ (wm(i)auxjh2o(ind+1) + wp(i)auxjh2o(ind)) *
	648	$ cortemp(i)
	649	!H2O2 interpolated coefficient
	650	jfotsout(indexint,6,auxi) =
[2464]	651	$ jfotsout(indexint,6,nlayer) *
[1310]	652	$ (wm(i)auxjh2o2(ind+1) + wp(i)auxjh2o2(ind)) *
	653	$ cortemp(i)
	654	!CO interpolated coefficient
	655	jfotsout(indexint,11,auxi) =
[2464]	656	$ jfotsout(indexint,11,nlayer) *
[1310]	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
[2464]	662	do i=1,nlayer
[1310]	663	ind=auxind(i)
[2464]	664	auxi = nlayer-i+1
[1310]	665	!NO interpolated coefficient
	666	jfotsout(indexint,10,auxi)=
[2464]	667	$ jfotsout(indexint,10,nlayer) *
[1310]	668	$ (wm(i)auxjno(ind+1) + wp(i)auxjno(ind)) *
	669	$ cortemp(i)
	670	!NO2 interpolated coefficient
	671	jfotsout(indexint,13,auxi)=
[2464]	672	$ jfotsout(indexint,13,nlayer) *
[1310]	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
[2464]	693	do i=1,nlayer
	694	auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
[1310]	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
[2464]	727	$ (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
	728	do i=1,nlayer
[1310]	729	ind=auxind(i)
[2464]	730	auxi = nlayer-i+1
[1310]	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) =
[2464]	741	$ jfotsout(indexint,1,nlayer) *
[1310]	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) =
[2464]	748	$ jfotsout(indexint,2,nlayer) *
[1310]	749	$ (wm(i)auxjo2(ind+1) + wp(i)auxjo2(ind)) *
	750	$ cortemp(i)
	751	!H2O interpolated coefficient
	752	jfotsout(indexint,4,auxi) =
[2464]	753	$ jfotsout(indexint,4,nlayer) *
[1310]	754	$ (wm(i)auxjh2o(ind+1) + wp(i)auxjh2o(ind)) *
	755	$ cortemp(i)
	756	!H2O2 interpolated coefficient
	757	jfotsout(indexint,6,auxi) =
[2464]	758	$ jfotsout(indexint,6,nlayer) *
[1310]	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
[2464]	764	do i=1,nlayer
[1310]	765	ind=auxind(i)
[2464]	766	auxi = nlayer-i+1
[1310]	767	!NO interpolated coefficient
	768	jfotsout(indexint,10,auxi)=
[2464]	769	$ jfotsout(indexint,10,nlayer) *
[1310]	770	$ (wm(i)auxjno(ind+1) +wp(i)auxjno(ind)) *
	771	$ cortemp(i)
	772	!NO2 interpolated coefficient
	773	jfotsout(indexint,13,auxi)=
[2809]	774	$ jfotsout(indexint,13,nlayer) *
[1310]	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
[2464]	795	do i=1,nlayer
	796	auxcolinp(nlayer-i+1) =co2colx(i) + o2colx(i) + h2o2colx(i) +
[1310]	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
[2464]	824	$ (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
	825	do i=1,nlayer
[1310]	826	ind=auxind(i)
[2464]	827	auxi = nlayer-i+1
[1310]	828	!Correction to include T variation of CO2 cross section
[2464]	829	if(sigma(indexint,nlayer-i+1)alfa(indexint,auxi)
[1310]	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) =
[2464]	838	$ jfotsout(indexint,1,nlayer) *
[1310]	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) =
[2464]	845	$ jfotsout(indexint,2,nlayer) *
[1310]	846	$ (wm(i)auxjo2(ind+1) + wp(i)auxjo2(ind)) *
	847	$ cortemp(i)
	848	!H2O2 interpolated coefficient
	849	jfotsout(indexint,6,auxi) =
[2464]	850	$ jfotsout(indexint,6,nlayer) *
[1310]	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
[2464]	856	do i=1,nlayer
	857	auxi = nlayer-i+1
[1310]	858	ind=auxind(i)
	859	!NO interpolated coefficient
	860	jfotsout(indexint,10,auxi) =
[2464]	861	$ jfotsout(indexint,10,nlayer) *
[1310]	862	$ (wm(i)auxjno(ind+1) + wp(i)auxjno(ind)) *
	863	$ cortemp(i)
	864	!NO2 interpolated coefficient
	865	jfotsout(indexint,13,auxi) =
[2464]	866	$ jfotsout(indexint,13,nlayer) *
[1310]	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
[2464]	885	do i=1,nlayer
	886	auxcolinp(nlayer-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
[1310]	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
[2464]	908	$ (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
	909	do i=1,nlayer
[1310]	910	ind=auxind(i)
[2464]	911	auxi = nlayer-i+1
[1310]	912	!O2 interpolated coefficient
[2464]	913	jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
[1310]	914	$ (wm(i)auxjo2(ind+1) + wp(i)auxjo2(ind))
	915	!H2O2 interpolated coefficient
[2464]	916	jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayer) *
[1310]	917	$ (wm(i)auxjh2o2(ind+1) + wp(i)auxjh2o2(ind))
	918	enddo
	919	!Only if chemthermod.ge.2
	920	if(chemthermod.ge.2) then
[2464]	921	do i=1,nlayer
[1310]	922	ind=auxind(i)
	923	!NO2 interpolated coefficient
[2464]	924	jfotsout(indexint,13,nlayer-i+1) =
	925	$ jfotsout(indexint,13,nlayer) *
[1310]	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
[2464]	943	do i=1,nlayer
	944	auxcolinp(nlayer-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
[1310]	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
[2464]	969	$ (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
	970	do i=1,nlayer
[1310]	971	ind=auxind(i)
[2464]	972	auxi = nlayer-i+1
[1310]	973	!O2 interpolated coefficient
[2464]	974	jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
[1310]	975	$ (wm(i)auxjo2(ind+1) + wp(i)auxjo2(ind))
	976	!H2O2 interpolated coefficient
[2464]	977	jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayer) *
[1310]	978	$ (wm(i)auxjh2o2(ind+1) + wp(i)auxjh2o2(ind))
	979	!O3 interpolated coefficient
[2464]	980	jfotsout(indexint,7,auxi) = jfotsout(indexint,7,nlayer) *
[1310]	981	$ (wm(i)auxjo3(ind+1) + wp(i)auxjo3(ind))
	982	enddo
	983	!Only if chemthermod.ge.2
	984	if(chemthermod.ge.2) then
[2464]	985	do i=1,nlayer
[1310]	986	ind=auxind(i)
	987	!NO2 interpolated coefficient
[2464]	988	jfotsout(indexint,13,nlayer-i+1) =
	989	$ jfotsout(indexint,13,nlayer) *
[1310]	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
[2464]	1007	do i=1,nlayer
	1008	auxcolinp(nlayer-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
[1310]	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
[2464]	1030	$ (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
	1031	do i=1,nlayer
[1310]	1032	ind=auxind(i)
[2464]	1033	auxi = nlayer-i+1
[1310]	1034	!H2O2 interpolated coefficient
[2464]	1035	jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayer) *
[1310]	1036	$ (wm(i)auxjh2o2(ind+1) + wp(i)auxjh2o2(ind))
	1037	!O3 interpolated coefficient
[2464]	1038	jfotsout(indexint,7,auxi) = jfotsout(indexint,7,nlayer) *
[1310]	1039	$ (wm(i)auxjo3(ind+1) + wp(i)auxjo3(ind))
	1040	enddo
	1041	if(chemthermod.ge.2) then
[2464]	1042	do i=1,nlayer
[1310]	1043	ind=auxind(i)
	1044	!NO2 interpolated coefficient
[2464]	1045	jfotsout(indexint,13,nlayer-i+1) =
	1046	$ jfotsout(indexint,13,nlayer) *
[1310]	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
[2464]	1063	do i=1,nlayer
	1064	auxcolinp(nlayer-i+1) = o3colx(i) + no2colx(i)
[1310]	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
[2464]	1084	$ (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
	1085	do i=1,nlayer
[1310]	1086	ind=auxind(i)
	1087	!O3 interpolated coefficient
[2464]	1088	jfotsout(indexint,7,nlayer-i+1) =
	1089	$ jfotsout(indexint,7,nlayer) *
[1310]	1090	$ (wm(i)auxjo3(ind+1) + wp(i)auxjo3(ind))
	1091	enddo
	1092	!Only if chemthermod.ge.2
	1093	if(chemthermod.ge.2) then
[2464]	1094	do i=1,nlayer
[1310]	1095	ind=auxind(i)
	1096	!NO2 interpolated coefficient
[2464]	1097	jfotsout(indexint,13,nlayer-i+1) =
	1098	$ jfotsout(indexint,13,nlayer) *
[1310]	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
[2464]	1107	c Coefficients are refered to Sun-Mars distance
	1108	c Correction to the Sun-Venus distance (fixed)
[1310]	1109
[2464]	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,
[2836]	1118	$ co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx,
	1119	$ n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
[2464]	1120
[2836]	1121	c Jun 2022 AM readapted to Venus GCM
[2464]	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
[2836]	1140	integer nesptherm !# of species undergoing chemistry, input
[2464]	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
[2836]	1146	real o2colx(klev) !column density of O2(cm^-2), output
[2464]	1147	real o3pcolx(klev) !column density of O(3P)(cm^-2), output
[2836]	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
[2464]	1152	real n2colx(klev) !N2 column density (cm-2), output
[2836]	1153	real ncolx(klev) !N column density (cm-2), output
	1154	real nocolx(klev) !NO column density (cm-2), output
[2464]	1155	real cocolx(klev) !CO column density (cm-2), output
[2836]	1156	real hcolx(klev) !H column density (cm-2), output
	1157	real no2colx(klev) !NO2 column density (cm-2), output
[2464]	1158
[2836]	1159
[2464]	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
[2836]	1167	! density
[2464]	1168	real co2x(klev)
[2836]	1169	real o2x(klev)
[2464]	1170	real o3px(klev)
	1171	real cox(klev)
[2836]	1172	real hx(klev)
	1173	real h2x(klev)
	1174	real h2ox(klev)
	1175	real h2o2x(klev)
	1176	real o3x(klev)
[2464]	1177	real n2x(klev)
	1178	real nx(klev)
[2836]	1179	real nox(klev)
	1180	real no2x(klev)
[2464]	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
[2836]	1202	! integer,parameter :: ix_co2=1
	1203	! integer,parameter :: ix_n2=13
	1204	! integer,parameter :: ix_o=3
	1205	! integer,parameter :: ix_co=4
[2464]	1206
[2836]	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
[2464]	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
[2836]	1237	Ho2 = xx / mmolo2
	1238	Ho3p = xx / mmolo ! Oxygen 3P
	1239	Hh2 = xx / mmolh2
	1240	Hh2o2 = xx / mmolh2o2
	1241	Hh2o = xx / mmolh2o
[2464]	1242
	1243	!Only if O3 chem. required
[2836]	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
[2464]	1255	! first loop in altitude : initialisation
	1256	do i=klev,1,-1
	1257	!Column initialisation
	1258	co2colx(i) = 0.
[2836]	1259	o2colx(i) = 0.
[2464]	1260	o3pcolx(i) = 0.
[2836]	1261	h2colx(i) = 0.
	1262	h2ocolx(i) = 0.
	1263	h2o2colx(i) = 0.
	1264	o3colx(i) = 0.
[2464]	1265	n2colx(i) = 0.
[2836]	1266	ncolx(i) = 0.
	1267	nocolx(i) = 0.
[2464]	1268	cocolx(i) = 0.
[2836]	1269	hcolx(i) = 0.
	1270	no2colx(i) = 0.
[2464]	1271	!--Densities [cm-3]
	1272	co2x(i) = rm(i,ix_co2)
[2836]	1273	o2x(i) = rm(i,ix_o2)
[2464]	1274	o3px(i) = rm(i,ix_o)
[2836]	1275	h2x(i) = rm(i,ix_h2)
	1276	h2ox(i) = rm(i,ix_h2o)
	1277	h2o2x(i) = rm(i,ix_h2o2)
[2464]	1278	cox(i) = rm(i,ix_co)
[2836]	1279	hx(i) = rm(i,ix_h) ! write(,), '--jthermcalc--', co2x(i)
[2464]	1280
	1281	!Only if O3 chem. required
[2836]	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
[2464]	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
[2836]	1300	o2colx(i)=1.e25
[2464]	1301	o3pcolx(i)=1.e25
[2836]	1302	h2colx(i)=1.e25
	1303	h2ocolx(i)=1.e25
	1304	h2o2colx(i)=1.e25
	1305	o3colx(i)=1.e25
[2464]	1306	n2colx(i)=1.e25
	1307	cocolx(i)=1.e25
[2836]	1308	hcolx(i)=1.e25
	1309	ncolx(i)=1.e25
	1310	nocolx(i)=1.e25
	1311	no2colx(i)=1.e25
[2464]	1312	else
	1313	rcmnz = ( radio + iz(klev) ) * 1.e5 ! km --> cm
	1314	rcmmini = ( radio + zmini ) * 1.e5
[2836]	1315	!Column calculation taking into account the geometrical depth
[2464]	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
[2836]	1326	h2o2colx(i)=h2o2colx(i)+
	1327	$ h2o2x(klev)Hh2o2esp(j)*1.e-5
	1328	o2colx(i)=o2colx(i)+o2x(klev)Ho2esp(j)
[2464]	1329	$ *1.e-5
[2836]	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
[2464]	1334	n2colx(i)=n2colx(i)+n2x(klev)Hn2esp(j)
[2836]	1335	$ *1.e-5
	1336	cocolx(i)=cocolx(i)+cox(klev)Hcoesp(j)
[2464]	1337	$ *1.e-5
[2836]	1338	hcolx(i)=hcolx(i)+hx(klev)Hhesp(j)
	1339	$ *1.e-5
[2464]	1340
	1341	!Only if O3 chemistry required
[2836]	1342	if(chemthermod.ge.1) o3colx(i)=
	1343	$ o3colx(i)+o3x(klev)Ho3esp(j)
	1344	$ *1.e-5
[2464]	1345	!Only if N or ion chemistry requested
[2836]	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
[2464]	1354	else if(zenit.gt.60.) then
	1355	espco2 =sqrt((rcmnz+Hco2)2 -rcmmini2) - esp(j)
[2836]	1356	espo2 = sqrt((rcmnz+Ho2)2 -rcmmini2) - esp(j)
[2464]	1357	espo3p = sqrt((rcmnz+Ho3p)2 -rcmmini2)- esp(j)
[2836]	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)
[2464]	1361	espco = sqrt((rcmnz+Hco)2 -rcmmini2) - esp(j)
	1362	espn2 =sqrt((rcmnz+Hn2)2-rcmmini2)-esp(j)
[2836]	1363	esph = sqrt((rcmnz+Hh)2 -rcmmini2) - esp(j)
[2464]	1364
	1365	!Only if O3 chemistry required
[2836]	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
[2464]	1374
	1375	co2colx(i) = co2colx(i) + espco2*co2x(klev)
[2836]	1376	o2colx(i) = o2colx(i) + espo2*o2x(klev)
[2464]	1377	o3pcolx(i) = o3pcolx(i) + espo3p*o3px(klev)
[2836]	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)
[2464]	1382	cocolx(i) = cocolx(i) + espco*cox(klev)
[2836]	1383	hcolx(i) = hcolx(i) + esph*hx(klev)
[2464]	1384
	1385	!Only if O3 chemistry required
[2836]	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
[2464]	1394	endif !Of if zenit.lt.60
	1395	!Other layers
	1396	else
[2836]	1397	co2colx(i) = co2colx(i) +
[2464]	1398	$ esp(j) * (co2x(jj)+co2x(jj+1)) / 2.
[2836]	1399	o2colx(i) = o2colx(i) +
	1400	$ esp(j) * (o2x(jj)+o2x(jj+1)) / 2.
	1401	o3pcolx(i) = o3pcolx(i) +
[2464]	1402	$ esp(j) * (o3px(jj)+o3px(jj+1)) / 2.
[2836]	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.
[2464]	1409	n2colx(i) = n2colx(i) +
	1410	$ esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
[2836]	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.
[2464]	1415
	1416	!Only if O3 chemistry required
[2836]	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
[2464]	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
[1310]	1433	return
	1434
	1435	end
	1436
	1437
[2464]	1438	c**********************************************************************
	1439	c**********************************************************************
[1310]	1440
[2464]	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
[2836]	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
[2464]	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
[1310]	1473
[2464]	1474	c**********************************************************************
	1475	c**********************************************************************
[1310]	1476
[2464]	1477	subroutine espesor_optico_A (ig,capa, szadeg,z,
	1478	@ nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
[1310]	1479
[2464]	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]
[2836]	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
[2464]	1496	real*8 esp(nz3) ! O. layer widths after geometrically
[2836]	1497	! amplified; in [cm] except at TOA
	1498	! where an auxiliary value is used
[2464]	1499	real*8 ilayesp(nz3) ! O. Indexes of layers along ray path
	1500	real*8 szalayesp(nz3) ! O. SZA [deg] " " "
	1501	integer nlayesp
[2836]	1502	! real*8 nlayesp ! O. # layers along ray path at this z
[2464]	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:
[2836]	1539	! 1: First, all layers between z and ztop ("upper part of ray")
[2464]	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
[2836]	1566	! for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately)
[2464]	1567	! We deal with 2 different regions:
[2836]	1568	! 1: First, all layers between z and ztop ("upper part of ray")
[2464]	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]
[2836]	1588	szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
[2464]	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]
[2836]	1637	szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592! [deg]
[2464]	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]
[2836]	1662	szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592! [deg]
[2464]	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]
[2836]	1676	szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592! [deg]
[2464]	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 )
[2836]	1686	$ - sqrt( (radio+diz(j))2 - rkmmini2 )! [km]
	1687	esp(nlayesp) = esp(nlayesp) * 1.d5! [cm]
[2464]	1688	rkmj = radio+diz(j)
[2836]	1689	szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )! [rad]
	1690	szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592! [deg]
[2464]	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