Context Navigation

jthermcalc.F @ 1974

Last change on this file since 1974 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

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

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

Download in other formats:

Original Format