Context Navigation

jthermcalc.F @ 1704

Last change on this file since 1704 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

File size: 58.9 KB

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

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

Download in other formats:

Original Format