Context Navigation

jthermcalc.F @ 784

Last change on this file since 784 was 658, checked in by emillour, 13 years ago

Mars GCM:

updated high atmosphere photochemistry (jthermcalc.F, param_v4.h, iono.h, paramfoto_compact.F, param_read.F , thermosphere.F).
minor change in calchim.F90 (to not use maxloc(zycol, dim = 2) function which seems to be a problem for g95) .
minor bug fix in perosat.F; set tendency on pdqscloud for h2o2 to zero if none is computed.
in "moldiff.F", changed "tridag" to "tridag_sp", "LUBKSB" to "LUBKSB_SP" and "LUDCMP" to "LUDCMP_SP" to avoid possible conflicts with same routines defined in "moldiff_red.F". Added use of automatic-sized array in "tridag" and "tridag_sp" local array "gam" (to avoid using an a priori oversized local array).

FGG+JYC+EM

File size: 59.0 KB

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

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

Download in other formats:

Original Format