Context Navigation

mztf_overlap.F @ 426

Last change on this file since 426 was 414, checked in by aslmd, 14 years ago

LMDZ.MARS : NEW NLTE MODEL FROM GRANADA AMIGOS

23/11/11 == FGG + MALV

New parameterization of the NLTE 15 micron cooling. The old parameterization is kept as an option, including or not variable atomic oxygen concentration. A new flag is introduced in callphys.def, nltemodel, to select which parameterization wants to be used (new one, old one with variable [O], or old one with fixed [O], see below). Includes many new subroutines and commons in phymars. Some existing routines are also modified:

-physiq.F. Call to the new subroutine NLTE_leedat in first call. Call to nltecool modified to allow for variable atomic oxygen. Depending on the value of nltemodel, the new subroutine NLTEdlvr09_TCOOL is called instead of nltecool.

-inifis.F. Reading of nltemodel is added.

-callkeys.h Declaration of nltemodel is added.

The following lines should be added to callphys.def (ideally after setting callnlte):

# NLTE 15um scheme to use.
# 0-> Old scheme, static oxygen
# 1-> Old scheme, dynamic oxygen
# 2-> New scheme
nltemodel = 2

A new directory, NLTEDAT, has to be included in datagcm.

Improvements into NLTE NIR heating parameterization to take into account variability with O/CO2 ratio and SZA. A new subroutine, NIR_leedat.F, and a new common, NIRdata.h, are included. A new flag, nircorr, is added in callphys.def, to include or not these corrections. The following files are modified:

-nirco2abs.F: nq and pq are added in the arguments. The corrections factors are interpolated to the GCM grid and included in the clculation. A new subroutine, interpnir, is added at the end of the file.

-physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs

-inifis: Reading new flag nircorr.

-callkeys.h: Declaration of nircorr.

The following lines have to be added to callphys.def (ideally after callnirco2):

# NIR NLTE correction for variable SZA and O/CO2?
# matters only if callnirco2=T
# 0-> no correction
# 1-> include correction
nircorr=1

A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm.

Small changes to the molecular diffusion scheme to fix the number of species considered, to avoid problems when compiling with more than 15 tracers (for example, when CH4 is included). Modified subroutines: aeronomars/moldiff.F and aeronomars/moldiffcoeff.F

File size: 24.1 KB

Rev	Line
[414]	1	c***********************************************************************
	2	c mztf.f
	3	c***********************************************************************
	4	c
	5	c program for calculating atmospheric transmittances
	6	c to be used in the calculation of curtis matrix coefficients
	7
	8	subroutine mztf ( ig,cf,cfup,cfdw,vc,taugr, ib,isot,
	9	@ iirw,iimu,itauout,icfout,itableout )
	10
	11	c i*out = 1 output of data
	12	c i*out = 0 no output
	13	c
	14	c jul 2011 malv+fgg adapted to LMD-MGCM
	15	c nov 98 mavl allow for overlaping in the lorentz line
	16	c jan 98 malv version for mz1d. based on curtis/mztf.for
	17	c 17-jul-96 mlp&crs change the calculation of mr.
	18	c evitar: divide por cero. anhadiendo: ff
	19	c oct-92 malv correct s(t) dependence for all histogr bands
	20	c june-92 malv proper lower levels for laser bands
	21	c may-92 malv new temperature dependence for laser bands
	22	c @ 991 malv boxing for the averaged absorber amount and t
	23	c ? malv extension up to 200 km altitude in mars
	24	c 13-nov-86 mlp include the temperature weighted to match
	25	c the eqw in the strong doppler limit.
	26	c***********************************************************************
	27
	28	implicit none
	29
	30	include 'nltedefs.h'
	31	include 'nlte_atm.h'
	32	include 'nlte_data.h'
	33	include 'nlte_curtis.h'
	34	include 'tcr_15um.h'
	35	include 'nlte_results.h'
	36
	37
	38	c arguments
	39	integer ig !ADDED FOR TRACEBACK
	40	real*8 cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o.
	41	real*8 vc(nl), taugr(nl) ! o
	42	integer ib ! i
	43	integer isot ! i
	44	integer iirw ! i
	45	integer iimu ! i
	46	integer itauout ! i
	47	integer icfout ! i
	48	integer itableout ! i
	49
	50	c local variables and constants
	51	integer i, in, ir, im, k ,j
	52	integer nmu
	53	parameter (nmu = 8)
	54	real*8 tau(nl,nl)
	55	real*8 tauinf(nl)
	56	real*8 con(nzy), coninf
	57	real*8 c1, c2
	58	real*8 t1, t2
	59	real*8 p1, p2
	60	real*8 mr1, mr2
	61	real*8 st1, st2
	62	real*8 c1box(70), c2box(70)
	63	real*8 ff ! to avoid too small numbers
	64	real*8 tvtbs(nzy)
	65	real*8 st, beta, ts, eqwmu
	66	real*8 mu(nmu), amu(nmu)
	67	real*8 zld(nl), zyd(nzy)
	68	real*8 correc
	69	real deltanux ! width of vib-rot band (cm-1)
	70	! character isotcode*2
	71	integer idummy
	72	real*8 Desp,wsL
	73
	74	c formats
	75	! 111 format(a1)
	76	! 112 format(a2)
	77	101 format(i1)
	78	202 format(i2)
	79	! 180 format(a80)
	80	! 181 format(a80)
	81	c***********************************************************************
	82
	83	c some needed values
	84	! rl=sqrt(log(2.d0))
	85	! pi2 = 3.14159265358989d0
	86	beta = 1.8d0
	87	idummy = 0
	88	Desp = 0.d0
	89	wsL = 0.d0
	90
	91	c esto es para que las subroutines de mztfsub calculen we
	92	c de la forma apropiada para mztf, no para fot
	93	icls=icls_mztf
	94
	95	c codigos para filenames
	96	! if (isot .eq. 1) isotcode = '26'
	97	! if (isot .eq. 2) isotcode = '28'
	98	! if (isot .eq. 3) isotcode = '36'
	99	! if (isot .eq. 4) isotcode = '27'
	100	! if (isot .eq. 5) isotcode = '62'
	101	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	102	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	103	! write (ibcode1,101) ib
	104	! else
	105	! write (ibcode2,202) ib
	106	! endif
	107	! write (*,'( 30h calculating curtis matrix : ,2x,
	108	! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot
	109
	110	c integration in angle !!!!!!!!!!!!!!!!!!!!
	111
	112	c------- diffusivity approx.
	113	if (iimu.eq.1) then
	114	! write (,) ' diffusivity approx. beta = ',beta
	115	mu(1) = 1.0d0
	116	amu(1)= 1.0d0
	117	c-------data for 8 points integration
	118	elseif (iimu.eq.4) then
	119	write (,)' 4 points for the gauss-legendre angle quadrature.'
	120	mu(1)=(1.0d0+0.339981043584856)/2.0d0
	121	mu(2)=(1.0d0-0.339981043584856)/2.0d0
	122	mu(3)=(1.0d0+0.861136311594053)/2.0d0
	123	mu(4)=(1.0d0-0.861136311594053)/2.0d0
	124	amu(1)=0.652145154862546
	125	amu(2)=amu(1)
	126	amu(3)=0.347854845137454
	127	amu(4)=amu(3)
	128	beta=1.0d0
	129	c-------data for 8 points integration
	130	elseif(iimu.eq.8) then
	131	write (,)' 8 points for the gauss-legendre angle quadrature.'
	132	mu(1)=(1.0d0+0.183434642495650)/2.0d0
	133	mu(2)=(1.0d0-0.183434642495650)/2.0d0
	134	mu(3)=(1.0d0+0.525532409916329)/2.0d0
	135	mu(4)=(1.0d0-0.525532409916329)/2.0d0
	136	mu(5)=(1.0d0+0.796666477413627)/2.0d0
	137	mu(6)=(1.0d0-0.796666477413627)/2.0d0
	138	mu(7)=(1.0d0+0.960289856497536)/2.0d0
	139	mu(8)=(1.0d0-0.960289856497536)/2.0d0
	140	amu(1)=0.362683783378362
	141	amu(2)=amu(1)
	142	amu(3)=0.313706645877887
	143	amu(4)=amu(3)
	144	amu(5)=0.222381034453374
	145	amu(6)=amu(5)
	146	amu(7)=0.101228536290376
	147	amu(8)=amu(7)
	148	beta=1.0d0
	149	end if
	150	c!!!!!!!!!!!!!!!!!!!!!!!
	151
	152	ccc
	153	ccc determine abundances included in the absorber amount
	154	ccc
	155
	156	c first, set up the grid ready for interpolation.
	157	do i=1,nzy
	158	zyd(i) = dble(zy(i))
	159	enddo
	160	do i=1,nl
	161	zld(i) = dble(zl(i))
	162	enddo
	163
	164
	165	c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
	166	c por similitud a la que se hace en cza.for
	167
	168	do i=1,nzy
	169	if (isot.eq.5) then
	170	con(i) = dble( coy(i) * imrco )
	171	else
	172	con(i) = dble( co2y(i) * imr(isot) )
	173	c vibr. temp of the bending mode :
	174	if (isot.eq.1) call interdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
	175	if (isot.eq.2) call interdp ( tvtbs,zyd,nzy, v628t1,zld,nl, 1 )
	176	if (isot.eq.3) call interdp ( tvtbs,zyd,nzy, v636t1,zld,nl, 1 )
	177	if (isot.eq.4) call interdp ( tvtbs,zyd,nzy, v627t1,zld,nl, 1 )
	178	correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )
	179	con(i) = con(i) * ( 1.d0 - correc )
	180	endif
	181	c-----------------------------------------------------------------------
	182	c mlp & cristina. 17 july 1996
	183	c change the calculation of mr. it is used for calculating partial press
	184	c alpha = alpha(self,co2)pp +alpha(n2)(pt-pp)
	185	c for an isotope, if mr is obtained by co2*imr(iso)/nt we are considerin
	186	c collisions with other co2 isotopes (including the major one, 626)
	187	c as if they were with n2. assuming mr as co2/nt, we consider collisions
	188	c of type 628-626 as of 626-626 instead of as 626-n2.
	189	c mrx(i)=con(i)/ntx(i) ! old malv
	190
	191	! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs
	192
	193	c jan 98:
	194	c esta modif de mlp implica anular el correc (deberia revisar esto)
	195	mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98
	196
	197	c-----------------------------------------------------------------------
	198
	199	end do
	200
	201	! como beta y 1.d5 son comunes a todas las weighted absorber amounts,
	202	! los simplificamos:
	203	! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )
	204	coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
	205
	206	! write (,) ' coninf =', coninf
	207
	208	ccc
	209	ccc temp dependence of the band strength and
	210	ccc nlte correction factor for the absorber amount
	211	ccc
	212	call mztf_correccion ( coninf, con, ib, isot, itableout )
	213
	214	ccc
	215	ccc reads histogrammed spectral data (strength for lte and vmr=1)
	216	ccc
	217	!hfile1 = dirspec//'hi'//dn ! ya no distinguimos entre d/n
	218	!! hfile1 = dirspec//'hid' ! (see why in his.for)
	219	! hfile='hid'
	220	!! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
	221	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
	222	!
	223	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	224	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	225	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat'
	226	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat'
	227	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat'
	228	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat'
	229	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
	230	! else
	231	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat'
	232	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat'
	233	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat'
	234	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat'
	235	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat'
	236	! endif
	237	! write (,) 'hisfile: ', hisfile
	238
	239	! the argument to rhist is to make this compatible with mztf_comp.f,
	240	! which is a useful modification of mztf.f (to change strengths of bands
	241	! call rhist (1.0)
	242	if(ib.eq.1) then
	243	if(isot.eq.1) then !Case 1
	244	mm=mm_c1
	245	nbox=nbox_c1
	246	tmin=tmin_c1
	247	tmax=tmax_c1
	248	do i=1,nbox_max
	249	no(i)=no_c1(i)
	250	dist(i)=dist_c1(i)
	251	do j=1,nhist
	252	sk1(j,i)=sk1_c1(j,i)
	253	xls1(j,i)=xls1_c1(j,i)
	254	xln1(j,i)=xln1_c1(j,i)
	255	xld1(j,i)=xld1_c1(j,i)
	256	enddo
	257	enddo
	258	do j=1,nhist
	259	thist(j)=thist_c1(j)
	260	enddo
	261	else if(isot.eq.2) then !Case 2
	262	mm=mm_c2
	263	nbox=nbox_c2
	264	tmin=tmin_c2
	265	tmax=tmax_c2
	266	do i=1,nbox_max
	267	no(i)=no_c2(i)
	268	dist(i)=dist_c2(i)
	269	do j=1,nhist
	270	sk1(j,i)=sk1_c2(j,i)
	271	xls1(j,i)=xls1_c2(j,i)
	272	xln1(j,i)=xln1_c2(j,i)
	273	xld1(j,i)=xld1_c2(j,i)
	274	enddo
	275	enddo
	276	do j=1,nhist
	277	thist(j)=thist_c2(j)
	278	enddo
	279	else if(isot.eq.3) then !Case 3
	280	mm=mm_c3
	281	nbox=nbox_c3
	282	tmin=tmin_c3
	283	tmax=tmax_c3
	284	do i=1,nbox_max
	285	no(i)=no_c3(i)
	286	dist(i)=dist_c3(i)
	287	do j=1,nhist
	288	sk1(j,i)=sk1_c3(j,i)
	289	xls1(j,i)=xls1_c3(j,i)
	290	xln1(j,i)=xln1_c3(j,i)
	291	xld1(j,i)=xld1_c3(j,i)
	292	enddo
	293	enddo
	294	do j=1,nhist
	295	thist(j)=thist_c3(j)
	296	enddo
	297	else if(isot.eq.4) then !Case 4
	298	mm=mm_c4
	299	nbox=nbox_c4
	300	tmin=tmin_c4
	301	tmax=tmax_c4
	302	do i=1,nbox_max
	303	no(i)=no_c4(i)
	304	dist(i)=dist_c4(i)
	305	do j=1,nhist
	306	sk1(j,i)=sk1_c4(j,i)
	307	xls1(j,i)=xls1_c4(j,i)
	308	xln1(j,i)=xln1_c4(j,i)
	309	xld1(j,i)=xld1_c4(j,i)
	310	enddo
	311	enddo
	312	do j=1,nhist
	313	thist(j)=thist_c4(j)
	314	enddo
	315	else
	316	write(,)'isot must be 2,3 or 4 for ib=1!!'
	317	write(,)'stop at mztf_overlap/317'
	318	stop
	319	endif
	320	else if (ib.eq.2) then
	321	if(isot.eq.1) then !Case 5
	322	mm=mm_c5
	323	nbox=nbox_c5
	324	tmin=tmin_c5
	325	tmax=tmax_c5
	326	do i=1,nbox_max
	327	no(i)=no_c5(i)
	328	dist(i)=dist_c5(i)
	329	do j=1,nhist
	330	sk1(j,i)=sk1_c5(j,i)
	331	xls1(j,i)=xls1_c5(j,i)
	332	xln1(j,i)=xln1_c5(j,i)
	333	xld1(j,i)=xld1_c5(j,i)
	334	enddo
	335	enddo
	336	do j=1,nhist
	337	thist(j)=thist_c5(j)
	338	enddo
	339	else
	340	write(,)'isot must be 1 for ib=2!!'
	341	write(,)'stop at mztf_overlap/341'
	342	stop
	343	endif
	344	else if (ib.eq.3) then
	345	if(isot.eq.1) then !Case 6
	346	mm=mm_c6
	347	nbox=nbox_c6
	348	tmin=tmin_c6
	349	tmax=tmax_c6
	350	do i=1,nbox_max
	351	no(i)=no_c6(i)
	352	dist(i)=dist_c6(i)
	353	do j=1,nhist
	354	sk1(j,i)=sk1_c6(j,i)
	355	xls1(j,i)=xls1_c6(j,i)
	356	xln1(j,i)=xln1_c6(j,i)
	357	xld1(j,i)=xld1_c6(j,i)
	358	enddo
	359	enddo
	360	do j=1,nhist
	361	thist(j)=thist_c6(j)
	362	enddo
	363	else
	364	write(,)'isot must be 1 for ib=3!!'
	365	write(,)'stop at mztf_overlap/365'
	366	stop
	367	endif
	368	else if (ib.eq.4) then
	369	if(isot.eq.1) then !Case 7
	370	mm=mm_c7
	371	nbox=nbox_c7
	372	tmin=tmin_c7
	373	tmax=tmax_c7
	374	do i=1,nbox_max
	375	no(i)=no_c7(i)
	376	dist(i)=dist_c7(i)
	377	do j=1,nhist
	378	sk1(j,i)=sk1_c7(j,i)
	379	xls1(j,i)=xls1_c7(j,i)
	380	xln1(j,i)=xln1_c7(j,i)
	381	xld1(j,i)=xld1_c7(j,i)
	382	enddo
	383	enddo
	384	do j=1,nhist
	385	thist(j)=thist_c7(j)
	386	enddo
	387	else
	388	write(,)'isot must be 1 for ib=4!!'
	389	write(,)'stop at mztf_overlap/389'
	390	stop
	391	endif
	392	else
	393	write(,)'ib must be 1,2,3 or 4!!'
	394	write(,)'stop at mztf_overlap/394'
	395	endif
	396
	397	if (isot.ne.5) deltanux = deltanu(isot,ib)
	398	if (isot.eq.5) deltanux = deltanuco
	399
	400	c******
	401	c****** calculation of tauinf(nl)
	402	c******
	403	call initial
	404
	405	ff=1.0e10
	406
	407	do i=nl,1,-1
	408
	409	if(i.eq.nl)then
	410
	411	call intz (zl(i),c2,p2,mr2,t2, con)
	412	do kr=1,nbox
	413	ta(kr)=t2
	414	end do
	415	! write (,) ' i, t2 =', i, t2
	416	call interstrength (st2,t2,ka,ta)
	417	aa = p2 * coninf * mr2 * (st2 * ff)
	418	bb = p2 * coninf * st2
	419	cc = coninf * st2
	420	dd = t2 * coninf * st2
	421	do kr=1,nbox
	422	ccbox(kr) = coninf * ka(kr)
	423	ddbox(kr) = t2 * ccbox(kr)
	424	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
	425	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	426	end do
	427	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
	428	c2 = c2 * st2 * dble(deltaz)
	429
	430	else
	431	call intz (zl(i),c1,p1,mr1,t1, con)
	432	do kr=1,nbox
	433	ta(kr)=t1
	434	end do
	435	! write (,) ' i, t1 =', i, t1
	436	call interstrength (st1,t1,ka,ta)
	437	do kr=1,nbox
	438	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
	439	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	440	end do
	441	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
	442	c1 = c1 * st1 * dble(deltaz)
	443	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	444	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	445	cc = cc + ( c1 + c2 ) / 2.d0
	446	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	447	do kr=1,nbox
	448	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) )/2.d0
	449	ddbox(kr) = ddbox(kr) + ( t1c1box(kr)+t2c2box(kr) )/2.d0
	450	end do
	451
	452	mr2 = mr1
	453	c2=c1
	454	do kr=1,nbox
	455	c2box(kr) = c1box(kr)
	456	end do
	457	t2=t1
	458	p2=p1
	459	end if
	460
	461	pt = bb / cc
	462	pp = aa / (cc*ff)
	463
	464	! ta=dd/cc
	465	! tdop = ta
	466	ts = dd/cc
	467	do kr=1,nbox
	468	ta(kr) = ddbox(kr) / ccbox(kr)
	469	end do
	470	! write (,) ' i, ts =', i, ts
	471	call interstrength(st,ts,ka,ta)
	472	! call intershape(alsa,alna,alda,tdop)
	473	call intershape(alsa,alna,alda,ta)
	474
	475	* ua = cc/st
	476
	477	c next loop calculates the eqw for an especified path ua,pp,pt,ta
	478
	479	eqwmu = 0.0d0
	480	do im = 1,iimu
	481	eqw=0.0d0
	482	do kr=1,nbox
	483	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	484	if(ua(kr).lt.0.)write(,)'mztf_overlap/483',ua(kr),
	485	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
	486
	487	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
	488	if ( i_supersat .eq. 0 ) then
	489	eqw=eqw+no(kr)*w
	490	else
	491	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	492	endif
	493	end do
	494	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	495	end do
	496
	497	tauinf(i) = exp( - eqwmu / dble(deltanux) )
	498
	499	end do ! i continue
	500
	501	! if ( isot.eq.1 .and. ib.eq.2 ) then
	502	! write (,) ' tauinf(nl) = ', tauinf(nl)
	503	! write (,) ' tauinf(1) = ', tauinf(1)
	504	! endif
	505
	506	c******
	507	c****** calculation of tau(in,ir) for n<=r
	508	c******
	509
	510	do 1 in=1,nl-1
	511
	512	call initial
	513	call intz (zl(in), c1,p1,mr1,t1, con)
	514	do kr=1,nbox
	515	ta(kr) = t1
	516	end do
	517	call interstrength (st1,t1,ka,ta)
	518	do kr=1,nbox
	519	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
	520	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	521	end do
	522	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
	523	c1 = c1 * st1 * dble(deltaz)
	524
	525	do 2 ir=in,nl-1
	526
	527	if (ir.eq.in) then
	528	tau(in,ir) = 1.d0
	529	goto 2
	530	end if
	531
	532	call intz (zl(ir), c2,p2,mr2,t2, con)
	533	do kr=1,nbox
	534	ta(kr) = t2
	535	end do
	536	call interstrength (st2,t2,ka,ta)
	537	do kr=1,nbox
	538	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
	539	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	540	end do
	541	! c2 = c2 * st2 * beta * dble(deltaz) * 1.e5
	542	c2 = c2 * st2 * dble(deltaz)
	543
	544	c aa = aa + ( p1mr1c1 + p2mr2c2 ) / 2.d0
	545	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	546	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	547	cc = cc + ( c1 + c2 ) / 2.d0
	548	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	549	do kr=1,nbox
	550	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
	551	ddbox(kr) = ddbox(kr) + ( t1c1box(kr) + t2c2box(kr) ) /2.d0
	552	end do
	553
	554	mr1=mr2
	555	t1=t2
	556	c1=c2
	557	p1=p2
	558	do kr=1,nbox
	559	c1box(kr) = c2box(kr)
	560	end do
	561
	562	pt = bb / cc
	563	pp = aa / (cc * ff)
	564
	565	* ta=dd/cc
	566	* tdop = ta
	567	ts = dd/cc
	568	do kr=1,nbox
	569	ta(kr) = ddbox(kr) / ccbox(kr)
	570	end do
	571	call interstrength(st,ts,ka,ta)
	572	call intershape(alsa,alna,alda,ta)
	573	* ua = cc/st
	574
	575	c next loop calculates the eqw for an especified path ua,pp,pt,ta
	576
	577	eqwmu = 0.0d0
	578	do im = 1,iimu
	579	eqw=0.0d0
	580	do kr=1,nbox
	581	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	582	if(ua(kr).lt.0.)write(,)'mztf_overlap/581',ua(kr),
	583	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
	584
	585	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
	586	if ( i_supersat .eq. 0 ) then
	587	eqw=eqw+no(kr)*w
	588	else
	589	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	590	endif
	591	end do
	592	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	593	end do
	594
	595	tau(in,ir) = exp( - eqwmu / dble(deltanux) )
	596
	597	2 continue
	598
	599	1 continue
	600
	601	! if ( isot.eq.1 .and. ib.eq.2 ) then
	602	! write (,) ' tau(1,) , =1,20 '
	603	! write (,) ( sngl(tau(1,k)), k=1,20 )
	604	! endif
	605
	606
	607	c**********
	608	c********** calculation of tau(in,ir) for n>r
	609	c**********
	610
	611	in=nl
	612
	613	call initial
	614	call intz (zl(in), c1,p1,mr1,t1, con)
	615	do kr=1,nbox
	616	ta(kr) = t1
	617	end do
	618	call interstrength (st1,t1,ka,ta)
	619	do kr=1,nbox
	620	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
	621	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	622	end do
	623	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
	624	c1 = c1 * st1 * dble(deltaz)
	625
	626	do 4 ir=in-1,1,-1
	627
	628	call intz (zl(ir), c2,p2,mr2,t2, con)
	629	do kr=1,nbox
	630	ta(kr) = t2
	631	end do
	632	call interstrength (st2,t2,ka,ta)
	633	do kr=1,nbox
	634	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
	635	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	636	end do
	637	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
	638	c2 = c2 * st2 * dble(deltaz)
	639
	640	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	641	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	642	cc = cc + ( c1 + c2 ) / 2.d0
	643	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	644	do kr=1,nbox
	645	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
	646	ddbox(kr) = ddbox(kr) + ( t1c1box(kr) + t2c2box(kr) ) /2.d0
	647	end do
	648
	649	mr1=mr2
	650	c1=c2
	651	t1=t2
	652	p1=p2
	653	do kr=1,nbox
	654	c1box(kr) = c2box(kr)
	655	end do
	656
	657	pt = bb / cc
	658	pp = aa / (cc * ff)
	659	ts = dd / cc
	660	do kr=1,nbox
	661	ta(kr) = ddbox(kr) / ccbox(kr)
	662	end do
	663	call interstrength (st,ts,ka,ta)
	664	call intershape (alsa,alna,alda,ta)
	665
	666	* ua = cc/st
	667
	668	c next loop calculates the eqw for an especified path ua,pp,pt,ta
	669
	670	eqwmu = 0.0d0
	671	do im = 1,iimu
	672	eqw=0.0d0
	673	do kr=1,nbox
	674	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	675	if(ua(kr).lt.0.)write(,)'mztf_overlap/674',ua(kr),
	676	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
	677
	678	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
	679	if ( i_supersat .eq. 0 ) then
	680	eqw=eqw+no(kr)*w
	681	else
	682	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	683	endif
	684	end do
	685	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	686	end do
	687
	688	tau(in,ir) = exp( - eqwmu / dble(deltanux) )
	689
	690	4 continue
	691
	692	c
	693	c due to the simmetry of the transmittances
	694	c
	695	do in=nl-1,2,-1
	696	do ir=in-1,1,-1
	697	tau(in,ir) = tau(ir,in)
	698	end do
	699	end do
	700
	701
	702	ccc
	703	ccc writing out transmittances
	704	ccc
	705	if (itauout.eq.1) then
	706
	707	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	708	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	709	! open( 1, file=
	710	! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat',
	711	! @ access='sequential', form='unformatted' )
	712	! else
	713	! open( 1, file=
	714	! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat',
	715	! @ access='sequential', form='unformatted' )
	716	! endif
	717
	718	! write(1) dummy
	719	! write(1)' format: (tauinf(n),(tau(n,r),r=1,nl),n=1,nl)'
	720	! do in=1,nl
	721	! write (1) tauinf(in), ( tau(in,ir), ir=1,nl )
	722	! end do
	723	! close(unit=1)
	724
	725	elseif (itauout.eq.2) then
	726
	727	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	728	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	729	! open( 1, file=
	730	! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat')
	731	! else
	732	! open( 1, file=
	733	! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat')
	734	! endif
	735
	736	! !write(1,*) dummy
	737	! !write(1,*) 'tij for curtis matrix calculations '
	738	! !write(1,*)' cira mars model atmosphere '
	739	! write(1,*)' beta= ',beta,'deltanu= ',deltanux
	740	! write(1,*)' number of elements (in,ir)= ',nl,nl
	741	! write(1,*)' format: (tauinf(in),(tau(in,ir),ir=1,nl),in=1,nl)'
	742
	743	! do in=1,nl
	744	! write (1,*) tauinf(in)
	745	! do ir=1,nl
	746	! write(1,*) tau(in,ir)
	747	! end do
	748	! end do
	749	! close(unit=1)
	750
	751	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	752	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	753	! write (*,'(1x, 31htransmitances written out in: ,a22)')
	754	! @ 'taul'//isotcode//dn//ibcode1
	755	! else
	756	! write (*,'(1x, 31htransmitances written out in: ,a22)')
	757	! @ 'taul'//isotcode//dn//ibcode2
	758	! endif
	759
	760	end if
	761
	762	c cleaning of transmittances
	763	! call elimin_tau(tau,tauinf,nl,nan,itableout,nw,dummy,
	764	! @ isotcode,dn,ibcode2)
	765
	766	c construction of the curtis matrix
	767
	768	call mzcf ( tauinf,tau, cf,cfup,cfdw, vc,taugr,
	769	@ ib,isot,icfout,itableout )
	770
	771
	772	c end
	773	return
	774	end

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

Download in other formats:

Original Format