Context Navigation

jthermcalc.F @ 1242

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

Mars GCM and common dynamics:

Common dynamics:

correction in inidissip (only matters in Martian case)
added correction in addfi on theta to account for surface pressure change.

Mars GCM:
Some fixes and updates:

addfi (dyn3d): Add correction on theta when surface pressure changes
vdif_cd (phymars): Correction for coefficients in stable nighttime case
jthermcalc (aeronomars): Fix for some pathological cases (further investigations on the origin of these is needed)

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

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

Download in other formats:

Original Format