Context Navigation

jthermcalc.F @ 2213

Last change on this file since 2213 was 2213, checked in by emillour, 5 years ago
Mars GCM: Bug fix in calchim (only had an impact when computations include the ionosphere), parameters should not be reset at every call. EM
File size: 61.4 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
[2213]	1379	real*8,intent(out) :: wm(nlayer),wp(nlayer) ! interpolation weights
	1380	integer,intent(out) :: nm(nlayer) ! index of nearest point
	1381	real*8,intent(in) :: pin(nl),p(nlayer)
	1382	real*8,intent(in) :: limup,limdown
	1383	integer,intent(in) :: nl,nlayer
	1384	integer :: n1,n,np,nini
[658]	1385	nini=1
	1386	do n1=1,nlayer
[635]	1387	if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then
[658]	1388	wm(n1) = 0.d0
	1389	wp(n1) = 0.d0
[635]	1390	else
[658]	1391	do n = nini,nl-1
[635]	1392	if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
[658]	1393	nm(n1)=n
[635]	1394	np=n+1
[658]	1395	wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n))
	1396	wp(n1)=1.d0 - wm(n1)
	1397	nini = n
	1398	exit
[635]	1399	endif
	1400	enddo
	1401	endif
[658]	1402	enddo
[2213]	1403
[635]	1404	end
	1405
	1406
	1407	c**********************************************************************
	1408	c**********************************************************************
	1409
[1266]	1410	subroutine espesor_optico_A (ig,capa,nlayer, szadeg,z,
[635]	1411	@ nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
	1412
	1413	c fgg nov 03 Adaptation to Martian model
	1414	c malv jul 03 Corrected z grid. Split in alt & frec codes
	1415	c fgg feb 03 first version
	1416	*************************************************************************
	1417
[1266]	1418	use param_v4_h, only: radio
[635]	1419	implicit none
	1420
	1421	c arguments
	1422
	1423	real szadeg ! I. SZA [rad]
	1424	real z ! I. altitude of interest [km]
[1266]	1425	integer nz3,ig,nlayer ! I. dimension of esp, ylayesp, etc...
	1426	! (=2*nlayer= max# of layers in ray path)
	1427	real iz(nlayer+1) ! I. Altitude of each layer
[635]	1428	real*8 esp(nz3) ! O. layer widths after geometrically
	1429	! amplified; in [cm] except at TOA
	1430	! where an auxiliary value is used
	1431	real*8 ilayesp(nz3) ! O. Indexes of layers along ray path
	1432	real*8 szalayesp(nz3) ! O. SZA [deg] " " "
	1433	integer nlayesp
	1434	! real*8 nlayesp ! O. # layers along ray path at this z
	1435	real*8 zmini ! O. Minimum altitud of ray path [km]
	1436
	1437
	1438	c local variables and constants
	1439
	1440	integer j,i,capa
	1441	integer jmin ! index of min.altitude along ray path
	1442	real*8 szarad ! SZA [deg]
	1443	real*8 zz
[1266]	1444	real*8 diz(nlayer+1)
[635]	1445	real*8 rkmnz ! distance TOA to center of Planet [km]
	1446	real*8 rkmmini ! distance zmini to center of P [km]
	1447	real*8 rkmj ! intermediate distance to C of P [km]
	1448	c external function
	1449	external grid_R8 ! Returns index of layer containing the altitude
	1450	! of interest, z; for example, if
	1451	! zkm(i)=z or zkm(i)<z<zkm(i+1) => grid(z)=i
	1452	integer grid_R8
	1453
	1454	*************************************************************************
	1455	szarad = dble(szadeg)*3.141592d0/180.d0
	1456	zz=dble(z)
[1266]	1457	do i=1,nlayer
[635]	1458	diz(i)=dble(iz(i))
	1459	enddo
	1460	do j=1,nz3
	1461	esp(j) = 0.d0
	1462	szalayesp(j) = 777.d0
	1463	ilayesp(j) = 0
	1464	enddo
	1465	nlayesp = 0
	1466
	1467	! First case: szadeg<60
	1468	! The optical thickness will be given by 1/cos(sza)
	1469	! We deal with 2 different regions:
	1470	! 1: First, all layers between z and ztop ("upper part of ray")
	1471	! 2: Second, the layer at ztop
	1472	if(szadeg.lt.60.d0) then
	1473
	1474	zmini = zz
[1266]	1475	if(abs(zz-diz(nlayer)).lt.1.d-3) goto 1357
[635]	1476	! 1st Zone: Upper part of ray
	1477	!
[1266]	1478	do j=grid_R8(zz,diz,nlayer),nlayer-1
[635]	1479	nlayesp = nlayesp + 1
	1480	ilayesp(nlayesp) = j
	1481	esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad) ! [km]
	1482	esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm]
	1483	szalayesp(nlayesp) = szadeg
	1484	end do
	1485
	1486	!
	1487	! 2nd Zone: Top layer
	1488	1357 continue
	1489	nlayesp = nlayesp + 1
[1266]	1490	ilayesp(nlayesp) = nlayer
[635]	1491	esp(nlayesp) = 1.d0 / cos(szarad) ! aux. non-dimens. factor
	1492	szalayesp(nlayesp) = szadeg
	1493
	1494
	1495	! Second case: 60 < szadeg < 90
	1496	! The optical thickness is evaluated.
	1497	! (the magnitude of the effect of not using cos(sza) is 3.e-5
	1498	! for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately)
	1499	! We deal with 2 different regions:
	1500	! 1: First, all layers between z and ztop ("upper part of ray")
	1501	! 2: Second, the layer at ztop ("uppermost layer")
	1502	else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then
	1503
	1504	zmini=(radio+zz)*sin(szarad)-radio
	1505	rkmmini = radio + zmini
	1506
[1266]	1507	if(abs(zz-diz(nlayer)).lt.1.d-4) goto 1470
[635]	1508
	1509	! 1st Zone: Upper part of ray
	1510	!
[1266]	1511	do j=grid_R8(zz,diz,nlayer),nlayer-1
[635]	1512	nlayesp = nlayesp + 1
	1513	ilayesp(nlayesp) = j
	1514	esp(nlayesp) =
	1515	# sqrt( (radio+diz(j+1))2 - rkmmini2 ) -
	1516	# sqrt( (radio+diz(j))2 - rkmmini2 ) ! [km]
	1517	esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm]
	1518	rkmj = radio+diz(j)
	1519	szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
	1520	szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 ! [deg]
	1521	end do
	1522	1470 continue
	1523	! 2nd Zone: Uppermost layer of ray.
	1524	!
	1525	nlayesp = nlayesp + 1
[1266]	1526	ilayesp(nlayesp) = nlayer
	1527	rkmnz = radio+diz(nlayer)
[635]	1528	esp(nlayesp) = sqrt( rkmnz2 - rkmmini2 ) ! aux.factor[km]
	1529	esp(nlayesp) = esp(nlayesp) * 1.d5 ! aux.f. [cm]
	1530	szalayesp(nlayesp) = asin( rkmmini/rkmnz ) ! [rad]
	1531	szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
	1532
	1533
	1534	! Third case: szadeg > 90
	1535	! The optical thickness is evaluated.
	1536	! We deal with 5 different regions:
	1537	! 1: all layers between z and ztop ("upper part of ray")
	1538	! 2: the layer at ztop ("uppermost layer")
	1539	! 3: the lowest layer, at zmini
	1540	! 4: the layers increasing from zmini to z (here SZA<90)
	1541	! 5: the layers decreasing from z to zmini (here SZA>90)
	1542	else if(szadeg.gt.90.d0) then
	1543
	1544	zmini=(radio+zz)*sin(szarad)-radio
[1260]	1545	!zmini should be lower than zz, as SZA<90. However, in situations
	1546	!where SZA is very close to 90, rounding errors can make zmini
	1547	!slightly higher than zz, causing problems in the determination
	1548	!of the jmin index. A correction is implemented in the determination
	1549	!of jmin, some lines below
[635]	1550	rkmmini = radio + zmini
	1551
	1552	if(zmini.lt.diz(1)) then ! Can see the sun? No => esp(j)=inft
	1553	nlayesp = nlayesp + 1
	1554	ilayesp(nlayesp) = - 1 ! Value to mark "no sun on view"
	1555	! esp(nlayesp) = 1.e30
	1556
	1557	else
[1266]	1558	jmin=grid_R8(zmini,diz,nlayer)+1
[1260]	1559	!Correction for possible rounding errors when SZA very close
	1560	!to 90 degrees
[1266]	1561	if(jmin.gt.grid_R8(zz,diz,nlayer)) then
[1260]	1562	write(,)'jthermcalc warning: possible rounding error'
	1563	write(,)'point,sza,layer:',ig,szadeg,capa
[1266]	1564	jmin=grid_R8(zz,diz,nlayer)
[1260]	1565	endif
[635]	1566
[1266]	1567	if(abs(zz-diz(nlayer)).lt.1.d-4) goto 9876
[635]	1568
	1569	! 1st Zone: Upper part of ray
	1570	!
[1266]	1571	do j=grid_R8(zz,diz,nlayer),nlayer-1
[635]	1572	nlayesp = nlayesp + 1
	1573	ilayesp(nlayesp) = j
	1574	esp(nlayesp) =
	1575	$ sqrt( (radio+diz(j+1))2 - rkmmini2 ) -
	1576	$ sqrt( (radio+diz(j))2 - rkmmini2 ) ! [km]
	1577	esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm]
	1578	rkmj = radio+diz(j)
	1579	szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
	1580	szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
	1581	end do
	1582
	1583	9876 continue
	1584	! 2nd Zone: Uppermost layer of ray.
	1585	!
	1586	nlayesp = nlayesp + 1
[1266]	1587	ilayesp(nlayesp) = nlayer
	1588	rkmnz = radio+diz(nlayer)
[635]	1589	esp(nlayesp) = sqrt( rkmnz2 - rkmmini2 ) !aux.factor[km]
	1590	esp(nlayesp) = esp(nlayesp) * 1.d5 !aux.f.[cm]
	1591	szalayesp(nlayesp) = asin( rkmmini/rkmnz ) ! [rad]
	1592	szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
	1593
	1594	! 3er Zone: Lowestmost layer of ray
	1595	!
	1596	if ( jmin .ge. 2 ) then ! If above the planet's surface
	1597	j=jmin-1
	1598	nlayesp = nlayesp + 1
	1599	ilayesp(nlayesp) = j
	1600	esp(nlayesp) = 2. *
	1601	$ sqrt( (radio+diz(j+1))2 -rkmmini2 ) ! [km]
	1602	esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm]
	1603	rkmj = radio+diz(j+1)
	1604	szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
	1605	szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
	1606	endif
	1607
	1608	! 4th zone: Lower part of ray, increasing from zmin to z
	1609	! ( layers with SZA < 90 deg )
[1266]	1610	do j=jmin,grid_R8(zz,diz,nlayer)-1
[635]	1611	nlayesp = nlayesp + 1
	1612	ilayesp(nlayesp) = j
	1613	esp(nlayesp) =
	1614	$ sqrt( (radio+diz(j+1))2 - rkmmini2 )
	1615	$ - sqrt( (radio+diz(j))2 - rkmmini2 ) ! [km]
	1616	esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm]
	1617	rkmj = radio+diz(j)
	1618	szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
	1619	szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
	1620	end do
	1621
	1622	! 5th zone: Lower part of ray, decreasing from z to zmin
	1623	! ( layers with SZA > 90 deg )
[1266]	1624	do j=grid_R8(zz,diz,nlayer)-1, jmin, -1
[635]	1625	nlayesp = nlayesp + 1
	1626	ilayesp(nlayesp) = j
	1627	esp(nlayesp) =
	1628	$ sqrt( (radio+diz(j+1))2 - rkmmini2 )
	1629	$ - sqrt( (radio+diz(j))2 - rkmmini2 ) ! [km]
	1630	esp(nlayesp) = esp(nlayesp) * 1.d5 ! [cm]
	1631	rkmj = radio+diz(j)
	1632	szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj ) ! [rad]
	1633	szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 ! [deg]
	1634	end do
	1635
	1636	end if
	1637
	1638	end if
	1639
[1260]	1640
[635]	1641	return
	1642
	1643	end
	1644
	1645
	1646
	1647	c**********************************************************************
	1648	c***********************************************************************
	1649
	1650	function grid_R8 (z, zgrid, nz)
	1651
	1652	c Returns the index where z is located within vector zgrid
	1653	c The vector zgrid must be monotonously increasing, otherwise program stops.
	1654	c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops.
	1655	c
	1656	c FGG Aug-2004 Correct z.lt.zgrid(i) to .le.
	1657	c MALV Jul-2003
	1658	c***********************************************************************
	1659
	1660	implicit none
	1661
	1662	c Arguments
	1663	integer nz
	1664	real*8 z
	1665	real*8 zgrid(nz)
	1666	integer grid_R8
	1667
	1668	c Local
	1669	integer i, nz1, nznew
	1670
	1671	c*** CODE START
	1672
[1260]	1673	if ( z .lt. zgrid(1) ) then
[635]	1674	write (,) ' GRID/ z outside bounds of zgrid '
	1675	write (,) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
[1260]	1676	z = zgrid(1)
	1677	write(,) 'WARNING: error in grid_r8 (jthermcalc.F)'
	1678	write(,) 'Please check values of z and zgrid above'
[635]	1679	endif
[1260]	1680	if (z .gt. zgrid(nz) ) then
	1681	write (,) ' GRID/ z outside bounds of zgrid '
	1682	write (,) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
	1683	z = zgrid(nz)
	1684	write(,) 'WARNING: error in grid_r8 (jthermcalc.F)'
	1685	write(,) 'Please check values of z and zgrid above'
	1686	endif
[635]	1687	if ( nz .lt. 2 ) then
	1688	write (,) ' GRID/ zgrid needs 2 points at least ! '
	1689	stop ' Serious error in GRID.F '
	1690	endif
	1691	if ( zgrid(1) .ge. zgrid(nz) ) then
	1692	write (,) ' GRID/ zgrid must increase with index'
	1693	stop ' Serious error in GRID.F '
	1694	endif
	1695
	1696	nz1 = 1
	1697	nznew = nz/2
	1698	if ( z .gt. zgrid(nznew) ) then
	1699	nz1 = nznew
	1700	nznew = nz
	1701	endif
	1702	do i=nz1+1,nznew
	1703	if ( z. eq. zgrid(i) ) then
	1704	grid_R8=i
	1705	return
	1706	elseif ( z .le. zgrid(i) ) then
	1707	grid_R8 = i-1
	1708	return
	1709	endif
	1710	enddo
	1711	grid_R8 = nz
	1712	return
	1713
	1714
	1715
	1716	end
	1717
	1718
	1719
	1720	!c***************************************************
	1721	!c***************************************************
	1722
	1723	subroutine flujo(date)
	1724
	1725
	1726	!c fgg nov 2002 first version
	1727	!c***************************************************
	1728
[1047]	1729	use comsaison_h, only: dist_sol
[1266]	1730	use param_v4_h, only: ninter,
	1731	. fluxtop, ct1, ct2, p1, p2
[635]	1732	implicit none
	1733
	1734
	1735	! common variables and constants
[1119]	1736	include "callkeys.h"
[635]	1737
	1738
	1739	! Arguments
	1740
	1741	real date
	1742
	1743
	1744	! Local variable and constants
	1745
	1746	integer i
	1747	integer inter
	1748	real nada
	1749
	1750	!c*************************************************
	1751
	1752	if(date.lt.1985.) date=1985.
	1753	if(date.gt.2001.) date=2001.
	1754
	1755	do i=1,ninter
	1756	fluxtop(i)=1.
	1757	!Variation of solar flux with 11 years solar cycle
	1758	!For more details, see Gonzalez-Galindo et al. 2005
	1759	!To be improved in next versions
[1119]	1760	if(i.le.24.and.solvarmod.eq.0) then
	1761	fluxtop(i)=(((ct1(i)+p1(i)*date)/2.)
	1762	$ sin(2.3.1416/11.*(date-1985.-3.1416))
	1763	$ +(ct2(i)+p2(i)date)+1.)fluxtop(i)
	1764	end if
	1765	fluxtop(i)=fluxtop(i)(1.52/dist_sol)*2
[635]	1766	end do
	1767
	1768	return
	1769	end

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

Download in other formats:

Original Format