Context Navigation

jthermcalc.F @ 2212

Last change on this file since 2212 was 2212, checked in by aslmd, 5 years ago
MESOSCALE: precompiling flags 1) change in physiq_mod to avoid new GCM subgrid-scale slope wind param 2) addition in jthermcalc because abort_gcm is no longer available in mesoscale setting and this is about some lines related to pure testing.
File size: 61.2 KB

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

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

Download in other formats:

Original Format