Context Navigation

jthermcalc.F @ 1119

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

Mars GCM:

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

FGG

File size: 59.3 KB

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

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

Download in other formats:

Original Format