Context Navigation

nlte_calc.F @ 706

Last change on this file since 706 was 695, checked in by acolaitis, 12 years ago
GCM: modifications to allow for nesting compilation. Transparent to GCM user. MESOSCALE: handling of functions and modules for nesting compilation. Full nested configuration with new physics is now fully compiling.
File size: 172.7 KB

Rev	Line
[498]	1	c***********************************************************************
	2	c mzescape.f
	3	c***********************************************************************
	4	c
	5	c program for calculating atmospheric escape functions, from
	6	c a calculation of transmittances and derivatives of these ones
	7
[690]	8	subroutine mzescape(ig,taustar,tauinf,tauii,
	9	& ib,isot, iirw,iimu)
[498]	10
	11	c jul 2011 malv+fgg adapted to LMD-MGCM
	12	c nov 99 malv adapt mztf to compute taustar (pg.23b-ma
	13	c nov 98 malv allow for overlaping in the lorentz line
	14	c jan 98 malv version for mz1d. based on curtis/mztf.for
	15	c 17-jul-96 mlp&crs change the calculation of mr.
	16	c evitar: divide por cero. anhadiendo: ff
	17	c oct-92 malv correct s(t) dependence for all histogr bands
	18	c june-92 malv proper lower levels for laser bands
	19	c may-92 malv new temperature dependence for laser bands
	20	c @ 991 malv boxing for the averaged absorber amount and t
	21	c ? malv extension up to 200 km altitude in mars
	22	c 13-nov-86 mlp include the temperature weighted to match
	23	c the eqw in the strong doppler limit.
	24	c***********************************************************************
	25
	26	implicit none
	27
	28	include 'nlte_paramdef.h'
	29	include 'nlte_commons.h'
	30
	31
	32	c arguments
	33	integer ig ! ADDED FOR TRACEBACK
	34	real*8 taustar(nl) ! o
	35	real*8 tauinf(nl) ! o
	36	real*8 tauii(nl) ! o
	37	integer ib ! i
	38	integer isot ! i
	39	integer iirw ! i
	40	integer iimu ! i
	41
	42
	43	c local variables and constants
	44	integer i, in, ir, im, k,j
	45	integer nmu
	46	parameter (nmu = 8)
	47	! real*8 tauinf(nl)
	48	real*8 con(nzy), coninf
	49	real*8 c1, c2, ccc
	50	real*8 t1, t2
	51	real*8 p1, p2
	52	real*8 mr1, mr2
	53	real*8 st1, st2
	54	real*8 c1box(70), c2box(70)
	55	real*8 ff ! to avoid too small numbers
	56	real*8 tvtbs(nzy)
	57	real*8 st, beta, ts, eqwmu
	58	real*8 mu(nmu), amu(nmu)
	59	real*8 zld(nl), zyd(nzy)
	60	real*8 correc
	61	real deltanux ! width of vib-rot band (cm-1)
	62	character isotcode*2
	63	real*8 maximum
	64	real*8 csL, psL, Desp, wsL ! for Strong Lorentz limit
	65
	66	c formats
	67	111 format(a1)
	68	112 format(a2)
	69	101 format(i1)
	70	202 format(i2)
	71	180 format(a80)
	72	181 format(a80)
	73	c***********************************************************************
	74
	75	c some needed values
	76	! rl=sqrt(log(2.d0))
	77	! pi2 = 3.14159265358989d0
	78	beta = 1.8d0
	79	! imrco = 0.9865
	80
	81	c esto es para que las subroutines de mztfsub calculen we
	82	c de la forma apropiada para mztf, no para fot
	83	icls=icls_mztf
	84
	85	c codigos para filenames
	86	! if (isot .eq. 1) isotcode = '26'
	87	! if (isot .eq. 2) isotcode = '28'
	88	! if (isot .eq. 3) isotcode = '36'
	89	! if (isot .eq. 4) isotcode = '27'
	90	! if (isot .eq. 5) isotcode = '62'
	91	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	92	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	93	!! encode(2,101,ibcode1)ib
	94	! write ( ibcode1, 101) ib
	95	! else
	96	!! encode(2,202,ibcode2)ib
	97	! write (ibcode2, 202) ib
	98	! endif
	99	! write (*,'( 30h calculating curtis matrix : ,2x,
	100	! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot
	101
	102	c integration in angle !!!!!!!!!!!!!!!!!!!!
	103
	104	c------- diffusivity approx.
	105	if (iimu.eq.1) then
	106	! write (,) ' diffusivity approx. beta = ',beta
	107	mu(1) = 1.0d0
	108	amu(1)= 1.0d0
	109	c-------data for 8 points integration
	110	elseif (iimu.eq.4) then
	111	write (,)' 4 points for the gauss-legendre angle quadrature.'
	112	mu(1)=(1.0d0+0.339981043584856)/2.0d0
	113	mu(2)=(1.0d0-0.339981043584856)/2.0d0
	114	mu(3)=(1.0d0+0.861136311594053)/2.0d0
	115	mu(4)=(1.0d0-0.861136311594053)/2.0d0
	116	amu(1)=0.652145154862546
	117	amu(2)=amu(1)
	118	amu(3)=0.347854845137454
	119	amu(4)=amu(3)
	120	beta=1.0d0
	121	c-------data for 8 points integration
	122	elseif(iimu.eq.8) then
	123	write (,)' 8 points for the gauss-legendre angle quadrature.'
	124	mu(1)=(1.0d0+0.183434642495650)/2.0d0
	125	mu(2)=(1.0d0-0.183434642495650)/2.0d0
	126	mu(3)=(1.0d0+0.525532409916329)/2.0d0
	127	mu(4)=(1.0d0-0.525532409916329)/2.0d0
	128	mu(5)=(1.0d0+0.796666477413627)/2.0d0
	129	mu(6)=(1.0d0-0.796666477413627)/2.0d0
	130	mu(7)=(1.0d0+0.960289856497536)/2.0d0
	131	mu(8)=(1.0d0-0.960289856497536)/2.0d0
	132	amu(1)=0.362683783378362
	133	amu(2)=amu(1)
	134	amu(3)=0.313706645877887
	135	amu(4)=amu(3)
	136	amu(5)=0.222381034453374
	137	amu(6)=amu(5)
	138	amu(7)=0.101228536290376
	139	amu(8)=amu(7)
	140	beta=1.0d0
	141	end if
	142	c!!!!!!!!!!!!!!!!!!!!!!!
	143
	144	ccc
	145	ccc determine abundances included in the absorber amount
	146	ccc
	147
	148	c first, set up the grid ready for interpolation.
	149	do i=1,nzy
	150	zyd(i) = dble(zy(i))
	151	enddo
	152	do i=1,nl
	153	zld(i) = dble(zl(i))
	154	enddo
	155
	156	c vibr. temp of the bending mode :
[695]	157	if (isot.eq.1) call interdp (tvtbs,zyd,nzy,v626t1,zld,nl,1)
	158	if (isot.eq.2) call interdp (tvtbs,zyd,nzy,v628t1,zld,nl,1)
	159	if (isot.eq.3) call interdp (tvtbs,zyd,nzy,v636t1,zld,nl,1)
	160	if (isot.eq.4) call interdp (tvtbs,zyd,nzy,v627t1,zld,nl,1)
[498]	161
	162	c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
	163	c por similitud a la que se hace en cza.for
	164
	165	do i=1,nzy
	166	if (isot.eq.5) then
	167	con(i) = dble( coy(i) * imrco )
	168	else
	169	con(i) = dble( co2y(i) * imr(isot) )
	170	correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )
	171	con(i) = con(i) * ( 1.d0 - correc )
	172	endif
	173	c-----------------------------------------------------------------------
	174	c mlp & cristina. 17 july 1996
	175	c change the calculation of mr. it is used for calculating partial press
	176	c alpha = alpha(self,co2)pp +alpha(n2)(pt-pp)
	177	c for an isotope, if mr is obtained by co2*imr(iso)/nt we are considerin
	178	c collisions with other co2 isotopes (including the major one, 626)
	179	c as if they were with n2. assuming mr as co2/nt, we consider collisions
	180	c of type 628-626 as of 626-626 instead of as 626-n2.
	181	c mrx(i)=con(i)/ntx(i) ! old malv
	182
	183	! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs
	184
	185	c jan 98:
	186	c esta modif de mlp implica anular el correc (deberia revisar esto)
	187	mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98
	188
	189	c-----------------------------------------------------------------------
	190
	191	end do
	192
	193	! como beta y 1.d5 son comunes a todas las weighted absorber amounts,
	194	! los simplificamos:
	195	! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )
	196	coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
	197
	198	! write (,) ' coninf =', coninf
	199
	200	ccc
	201	ccc temp dependence of the band strength and
	202	ccc nlte correction factor for the absorber amount
	203	ccc
	204	call mztf_correccion ( coninf, con, ib, isot, 0 )
	205
	206	ccc
	207	ccc reads histogrammed spectral data (strength for lte and vmr=1)
	208	ccc
	209	!hfile1 = dirspec//'hi'//dn ! Ya no distinguimos entre d/n
	210	!! hfile1 = dirspec//'hid' ! (see why in his.for)
	211	! hfile1='hid'
	212	!! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
	213	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
	214
	215	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	216	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	217	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat'
	218	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat'
	219	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat'
	220	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat'
	221	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
	222	! else
	223	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat'
	224	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat'
	225	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat'
	226	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat'
	227	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat'
	228	! endif
	229	!write (,) ' /MZESCAPE/ hisfile: ', hisfile
	230
	231	! the argument to rhist is to make this compatible with mztf_comp.f,
	232	! which is a useful modification of mztf.f (to change strengths of bands
	233	! call rhist (1.0)
	234	if(ib.eq.1) then
	235	if(isot.eq.1) then !Case 1
	236	mm=mm_c1
	237	nbox=nbox_c1
	238	tmin=tmin_c1
	239	tmax=tmax_c1
	240	do i=1,nbox_max
	241	no(i)=no_c1(i)
	242	dist(i)=dist_c1(i)
	243	do j=1,nhist
	244	sk1(j,i)=sk1_c1(j,i)
	245	xls1(j,i)=xls1_c1(j,i)
	246	xln1(j,i)=xln1_c1(j,i)
	247	xld1(j,i)=xld1_c1(j,i)
	248	enddo
	249	enddo
	250	do j=1,nhist
	251	thist(j)=thist_c1(j)
	252	enddo
	253	else if(isot.eq.2) then !Case 2
	254	mm=mm_c2
	255	nbox=nbox_c2
	256	tmin=tmin_c2
	257	tmax=tmax_c2
	258	do i=1,nbox_max
	259	no(i)=no_c2(i)
	260	dist(i)=dist_c2(i)
	261	do j=1,nhist
	262	sk1(j,i)=sk1_c2(j,i)
	263	xls1(j,i)=xls1_c2(j,i)
	264	xln1(j,i)=xln1_c2(j,i)
	265	xld1(j,i)=xld1_c2(j,i)
	266	enddo
	267	enddo
	268	do j=1,nhist
	269	thist(j)=thist_c2(j)
	270	enddo
	271	else if(isot.eq.3) then !Case 3
	272	mm=mm_c3
	273	nbox=nbox_c3
	274	tmin=tmin_c3
	275	tmax=tmax_c3
	276	do i=1,nbox_max
	277	no(i)=no_c3(i)
	278	dist(i)=dist_c3(i)
	279	do j=1,nhist
	280	sk1(j,i)=sk1_c3(j,i)
	281	xls1(j,i)=xls1_c3(j,i)
	282	xln1(j,i)=xln1_c3(j,i)
	283	xld1(j,i)=xld1_c3(j,i)
	284	enddo
	285	enddo
	286	do j=1,nhist
	287	thist(j)=thist_c3(j)
	288	enddo
	289	else if(isot.eq.4) then !Case 4
	290	mm=mm_c4
	291	nbox=nbox_c4
	292	tmin=tmin_c4
	293	tmax=tmax_c4
	294	do i=1,nbox_max
	295	no(i)=no_c4(i)
	296	dist(i)=dist_c4(i)
	297	do j=1,nhist
	298	sk1(j,i)=sk1_c4(j,i)
	299	xls1(j,i)=xls1_c4(j,i)
	300	xln1(j,i)=xln1_c4(j,i)
	301	xld1(j,i)=xld1_c4(j,i)
	302	enddo
	303	enddo
	304	do j=1,nhist
	305	thist(j)=thist_c4(j)
	306	enddo
	307	else
	308	write(,)'isot must be 2,3 or 4 for ib=1!!'
	309	write(,)'stop at mzescape/312'
	310	stop
	311	endif
	312	else if (ib.eq.2) then
	313	if(isot.eq.1) then !Case 5
	314	mm=mm_c5
	315	nbox=nbox_c5
	316	tmin=tmin_c5
	317	tmax=tmax_c5
	318	do i=1,nbox_max
	319	no(i)=no_c5(i)
	320	dist(i)=dist_c5(i)
	321	do j=1,nhist
	322	sk1(j,i)=sk1_c5(j,i)
	323	xls1(j,i)=xls1_c5(j,i)
	324	xln1(j,i)=xln1_c5(j,i)
	325	xld1(j,i)=xld1_c5(j,i)
	326	enddo
	327	enddo
	328	do j=1,nhist
	329	thist(j)=thist_c5(j)
	330	enddo
	331	else
	332	write(,)'isot must be 1 for ib=2!!'
	333	write(,)'stop at mzescape/336'
	334	stop
	335	endif
	336	else if (ib.eq.3) then
	337	if(isot.eq.1) then !Case 6
	338	mm=mm_c6
	339	nbox=nbox_c6
	340	tmin=tmin_c6
	341	tmax=tmax_c6
	342	do i=1,nbox_max
	343	no(i)=no_c6(i)
	344	dist(i)=dist_c6(i)
	345	do j=1,nhist
	346	sk1(j,i)=sk1_c6(j,i)
	347	xls1(j,i)=xls1_c6(j,i)
	348	xln1(j,i)=xln1_c6(j,i)
	349	xld1(j,i)=xld1_c6(j,i)
	350	enddo
	351	enddo
	352	do j=1,nhist
	353	thist(j)=thist_c6(j)
	354	enddo
	355	else
	356	write(,)'isot must be 1 for ib=3!!'
	357	write(,)'stop at mzescape/360'
	358	stop
	359	endif
	360	else if (ib.eq.4) then
	361	if(isot.eq.1) then !Case 7
	362	mm=mm_c7
	363	nbox=nbox_c7
	364	tmin=tmin_c7
	365	tmax=tmax_c7
	366	do i=1,nbox_max
	367	no(i)=no_c7(i)
	368	dist(i)=dist_c7(i)
	369	do j=1,nhist
	370	sk1(j,i)=sk1_c7(j,i)
	371	xls1(j,i)=xls1_c7(j,i)
	372	xln1(j,i)=xln1_c7(j,i)
	373	xld1(j,i)=xld1_c7(j,i)
	374	enddo
	375	enddo
	376	do j=1,nhist
	377	thist(j)=thist_c7(j)
	378	enddo
	379	else
	380	write(,)'isot must be 1 for ib=4!!'
	381	write(,)'stop at mzescape/384'
	382	stop
	383	endif
	384	else
	385	write(,)'ib must be 1,2,3 or 4!!'
	386	write(,)'stop at mzescape/389'
	387	endif
	388	if (isot.ne.5) deltanux = deltanu(isot,ib)
	389	if (isot.eq.5) deltanux = deltanuco
	390
	391	c******
	392	c****** calculation of tauinf(nl)
	393	c******
	394	call initial
	395
	396	ff=1.0e10
	397
	398	do i=nl,1,-1
	399
	400	if(i.eq.nl)then
	401
	402	call intz (zl(i),c2,p2,mr2,t2, con)
	403	do kr=1,nbox
	404	ta(kr)=t2
	405	end do
	406	! write (,) ' i, t2 =', i, t2
	407	call interstrength (st2,t2,ka,ta)
	408	aa = p2 * coninf * mr2 * (st2 * ff)
	409	bb = p2 * coninf * st2
	410	cc = coninf * st2
	411	dd = t2 * coninf * st2
	412	do kr=1,nbox
	413	ccbox(kr) = coninf * ka(kr)
	414	ddbox(kr) = t2 * ccbox(kr)
	415	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
	416	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	417	end do
	418	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
	419	c2 = c2 * st2 * dble(deltaz)
	420
	421	else
	422	call intz (zl(i),c1,p1,mr1,t1, con)
	423	do kr=1,nbox
	424	ta(kr)=t1
	425	end do
	426	! write (,) ' i, t1 =', i, t1
	427	call interstrength (st1,t1,ka,ta)
	428	do kr=1,nbox
	429	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
	430	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	431	end do
	432	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
	433	c1 = c1 * st1 * dble(deltaz)
	434	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	435	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	436	cc = cc + ( c1 + c2 ) / 2.d0
	437	ccc = ( c1 + c2 ) / 2.d0
	438	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	439	do kr=1,nbox
	440	ccbox(kr) = ccbox(kr) +
	441	@ ( c1box(kr) + c2box(kr) )/2.d0
	442	ddbox(kr) = ddbox(kr) +
	443	@ ( t1c1box(kr)+t2c2box(kr) )/2.d0
	444	end do
	445
	446	mr2 = mr1
	447	c2=c1
	448	do kr=1,nbox
	449	c2box(kr) = c1box(kr)
	450	end do
	451	t2=t1
	452	p2=p1
	453	end if
	454
	455	pt = bb / cc
	456	pp = aa / (cc*ff)
	457
	458	! ta=dd/cc
	459	! tdop = ta
	460	ts = dd/cc
	461	do kr=1,nbox
	462	ta(kr) = ddbox(kr) / ccbox(kr)
	463	end do
	464	! write (,) ' i, ts =', i, ts
	465	call interstrength(st,ts,ka,ta)
	466	! call intershape(alsa,alna,alda,tdop)
	467	call intershape(alsa,alna,alda,ta)
	468
	469	* ua = cc/st
	470
	471	c next loop calculates the eqw for an especified path ua,pp,pt,ta
	472
	473	eqwmu = 0.0d0
	474	do im = 1,iimu
	475	eqw=0.0d0
	476	do kr=1,nbox
	477	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	478	if(ua(kr).lt.0.)write(,)'mzescape/480',ua(kr),
	479	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
	480
	481	call findw (ig,iirw, 0, csL,psL, Desp, wsL)
	482	if ( i_supersat .eq. 0 ) then
	483	eqw=eqw+no(kr)*w
	484	else
	485	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	486	endif
	487	end do
	488	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	489	end do
	490
	491	! tauinf(i) = exp( - eqwmu / dble(deltanux) )
	492	tauinf(i) = 1.d0 - eqwmu / dble(deltanux)
	493	if (tauinf(i).lt.0.d0) tauinf(i) = 0.0d0
	494
	495	if (i.eq.nl) then
	496	taustar(i) = 0.0d0
	497	else
	498	taustar(i) = dble(deltanux) * (tauinf(i+1)-tauinf(i))
	499	! ~ / ( beta * cc * 1.d5 )
	500	~ / ( beta * ccc * 1.d5 )
	501	endif
	502
	503	end do ! i continue
	504
	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
	514	call intz (zl(in), c1,p1,mr1,t1, con)
	515	do kr=1,nbox
	516	ta(kr) = t1
	517	end do
	518	call interstrength (st1,t1,ka,ta)
	519	do kr=1,nbox
	520	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	521	end do
	522	c1 = c1 * st1 * dble(deltaz)
	523
	524	call intz (zl(in+1), c2,p2,mr2,t2, con)
	525	do kr=1,nbox
	526	ta(kr) = t2
	527	end do
	528	call interstrength (st2,t2,ka,ta)
	529	do kr=1,nbox
	530	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	531	end do
	532	c2 = c2 * st2 * dble(deltaz)
	533
	534	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	535	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	536	cc = cc + ( c1 + c2 ) / 2.d0
	537	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	538	do kr=1,nbox
	539	ccbox(kr) = ccbox(kr) + (c1box(kr)+c2box(kr))/2.d0
	540	ddbox(kr) = ddbox(kr) + (t1c1box(kr)+t2c2box(kr))/2.d0
	541	end do
	542
	543	mr1=mr2
	544	t1=t2
	545	c1=c2
	546	p1=p2
	547	do kr=1,nbox
	548	c1box(kr) = c2box(kr)
	549	end do
	550	pt = bb / cc
	551	pp = aa / (cc * ff)
	552	ts = dd/cc
	553	do kr=1,nbox
	554	ta(kr) = ddbox(kr) / ccbox(kr)
	555	end do
	556	call interstrength(st,ts,ka,ta)
	557	call intershape(alsa,alna,alda,ta)
	558
	559	eqwmu = 0.0d0
	560	do im = 1,iimu
	561	eqw=0.0d0
	562	do kr=1,nbox
	563	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	564	if(ua(kr).lt.0.)write(,)'mzescape/566',ua(kr),
	565	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
	566
	567	call findw (ig,iirw, 0, csL,psL, Desp, wsL)
	568	if ( i_supersat .eq. 0 ) then
	569	eqw=eqw+no(kr)*w
	570	else
	571	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	572	endif
	573	end do
	574	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	575	end do
	576
	577	tauii(in) = exp( - eqwmu / dble(deltanux) )
	578	!write (,) 'i,tauii=',in,tauii(in)
	579
	580	1 continue
	581	tauii(nl) = 1.0d0
	582
	583
	584	c end
	585	return
	586	end
	587
	588
	589
	590	c***********************************************************************
	591	c mzescape_normaliz.f
	592	c***********************************************************************
	593	c
	594	c program for correcting some strange values and for normalizing
	595	c the atmospheric escape functions computed by mzescape_15um.f
	596	c possibilities according to istyle (see mzescape_15um.f).
	597	c
	598
	599	subroutine mzescape_normaliz ( taustar, istyle )
	600
	601
	602	c dic 99 malv first version
	603	c jul 2011 malv+fgg Adapted to LMD-MGCM
	604	c***********************************************************************
	605
	606	implicit none
	607	include 'nlte_paramdef.h'
	608	include 'nlte_commons.h'
	609
	610
	611	c arguments
	612	real*8 taustar(nl) ! o
	613	integer istyle ! i
	614
	615	c local variables and constants
	616	integer i, imaximum
	617	real*8 maximum
	618
	619	c***********************************************************************
	620
	621	!
	622	! correcting strange values at top, eliminating local maxima, etc...
	623	!
	624	taustar(nl) = taustar(nl-1)
	625
	626	if ( istyle .eq. 1 ) then
	627	imaximum = nl
	628	maximum = taustar(nl)
	629	do i=1,nl-1
	630	if (taustar(i).gt.maximum) taustar(i) = taustar(nl)
	631	enddo
	632	elseif ( istyle .eq. 2 ) then
	633	imaximum = nl
	634	maximum = taustar(nl)
	635	do i=nl-1,1,-1
	636	if (taustar(i).gt.maximum) then
	637	maximum = taustar(i)
	638	imaximum = i
	639	endif
	640	enddo
	641	do i=imaximum,nl
	642	if (taustar(i).lt.maximum) taustar(i) = maximum
	643	enddo
	644	endif
	645
	646	!
	647	! normalizing
	648	!
	649	do i=1,nl
	650	taustar(i) = taustar(i) / maximum
	651	enddo
	652
	653
	654	c end
	655	return
	656	end
	657
	658
	659
	660	c***********************************************************************
	661	c mzescape_fb.f
	662	c***********************************************************************
	663	subroutine mzescape_fb(ig)
	664
	665	c computes the escape functions of the most important 15um bands
	666	c this calls mzescape ( taustar,tauinf,tauii, ib,isot, iirw,iimu
	667
	668	c nov 99 malv based on cm15um_fb.f
	669	c jul 2011 malv+fgg adapted to LMD-MGCM
	670	c***********************************************************************
	671
	672	implicit none
	673
	674	include 'nlte_paramdef.h'
	675	include 'nlte_commons.h'
	676
	677	c local variables
	678	integer i, ib, ik, istyle
	679	integer ig !ADDED FOR TRACEBACK
	680	real*8 tau_factor
	681	real*8 aux(nl), aux2(nl), aux3(nl)
	682
	683	c***********************************************************************
	684
[695]	685	call mzescape (ig,taustar21,tauinf210,tauii210,1,2
	686	& ,irw_mztf,imu)
	687	call mzescape (ig,taustar31,tauinf310,tauii310,1,3
	688	& ,irw_mztf,imu)
	689	call mzescape (ig,taustar41,tauinf410,tauii410,1,4
	690	& ,irw_mztf,imu)
[498]	691
	692	istyle = 2
	693	call mzescape_normaliz ( taustar21, istyle )
	694	call mzescape_normaliz ( taustar31, istyle )
	695	call mzescape_normaliz ( taustar41, istyle )
	696
	697
	698	c end
	699	return
	700	end
	701
	702
	703
	704	c***********************************************************************
	705	c mzescape_fh.f
	706	c***********************************************************************
	707	subroutine mzescape_fh(ig)
	708
	709	c jul 2011 malv+fgg
	710	c***********************************************************************
	711
	712	implicit none
	713
	714	include 'nlte_paramdef.h'
	715	include 'nlte_commons.h'
	716
	717	c local variables
	718	integer i, ib, ik, istyle
	719	integer ig ! ADDED FOR TRACEBACK
	720	real*8 tau_factor
	721	real*8 aux(nl), aux2(nl), aux3(nl)
	722
	723	c***********************************************************************
	724
	725	call zero4v( aux, taustar12,tauinf121,tauii121, nl)
	726	do ik=1,3
	727	ib=ik+1
	728	call mzescape ( ig,aux,aux2,aux3, ib, 1,irw_mztf,imu )
	729	tau_factor = 1.d0
	730	if (ik.eq.1) tau_factor = dble(667.75/618.03)
	731	if (ik.eq.3) tau_factor = dble(667.75/720.806)
	732	do i=1,nl
	733	taustar12(i) = taustar12(i) + aux(i) * tau_factor
	734	tauinf121(i) = tauinf121(i) + aux2(i) * tau_factor
	735	tauii121(i) = tauii121(i) + aux3(i) * tau_factor
	736	enddo
	737	enddo
	738
	739	istyle = 2
	740	call mzescape_normaliz ( taustar12, istyle )
	741
	742
	743
	744	c end
	745	return
	746	end
	747
	748
	749
	750
	751
	752	c***********************************************************************
	753	c mztud.f
	754	c***********************************************************************
	755
	756	subroutine mztud ( ig,cf,cfup,cfdw,vc,taugr, ib,isot,
	757	@ iirw,iimu,itauout,icfout,itableout )
	758
	759	c program for calculating atmospheric transmittances
	760	c to be used in the calculation of curtis matrix coefficients
	761	c i*out = 1 output of data
	762	c i*out = 0 no output
	763	c itableout = 30 output de toda la C.M. y el VC y las poblaciones de los
	764	c estados 626(020), esta opcion nueva se añade porque
	765	c itableout=1 saca o bien solamente de 5 en 5 capas
	766	c o bien los elementos de C.M. desde una cierta capa
	767	c (consultese elimin_mz1d.f que es quien lo hace); lo
	768	c de las poblaciones (020) lo hace mztf_correcion.f
	769
	770	c jul 2011 malv+fgg Adapted to LMD-MGCM
	771	c jan 07 malv Add new vertical fine grid zy, similar to zx
	772	c sep-oct 01 malv update for fluxes for hb and fb, adapt to Linux
	773	c nov 98 mavl allow for overlaping in the lorentz line
	774	c jan 98 malv version for mz1d. based on curtis/mztf.for
	775	c 17-jul-96 mlp&crs change the calculation of mr.
	776	c evitar: divide por cero. anhadiendo: ff
	777	c oct-92 malv correct s(t) dependence for all histogr bands
	778	c june-92 malv proper lower levels for laser bands
	779	c may-92 malv new temperature dependence for laser bands
	780	c @ 991 malv boxing for the averaged absorber amount and t
	781	c ? malv extension up to 200 km altitude in mars
	782	c 13-nov-86 mlp include the temperature weighted to match
	783	c the eqw in the strong doppler limit.
	784	c***********************************************************************
	785
	786	implicit none
	787
	788	include 'nlte_paramdef.h'
	789	include 'nlte_commons.h'
	790
	791	c arguments
	792	integer ig !ADDED FOR TRACEBACK
	793	real*8 cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o
	794	real*8 vc(nl), taugr(nl) ! o
	795	integer ib ! i
	796	integer isot ! i
	797	integer iirw ! i
	798	integer iimu ! i
	799	integer itauout ! i
	800	integer icfout ! i
	801	integer itableout ! i
	802
	803	c local variables and constants
	804	integer i, in, ir, im, k,j
	805	integer nmu
	806	parameter (nmu = 8)
	807	real*8 tau(nl,nl)
	808	real*8 tauinf(nl)
	809	real*8 con(nzy), coninf
	810	real*8 c1, c2
	811	real*8 t1, t2
	812	real*8 p1, p2
	813	real*8 mr1, mr2
	814	real*8 st1, st2
	815	real*8 c1box(70), c2box(70)
	816	real*8 ff ! to avoid too small numbers
	817	real*8 tvtbs(nzy)
	818	real*8 st, beta, ts, eqwmu
	819	real*8 mu(nmu), amu(nmu)
	820	real*8 zld(nl), zyd(nzy)
	821	real*8 correc
	822	real deltanux ! width of vib-rot band (cm-1)
	823	character isotcode*2
	824	integer idummy
	825	real*8 Desp,wsL
	826
	827	c formats
	828	111 format(a1)
	829	112 format(a2)
	830	101 format(i1)
	831	202 format(i2)
	832	180 format(a80)
	833	181 format(a80)
	834	c***********************************************************************
	835
	836	c some needed values
	837	! rl=sqrt(log(2.d0))
	838	! pi2 = 3.14159265358989d0
	839	beta = 1.8d0
	840	! beta = 1.0d0
	841	idummy = 0
	842	Desp = 0.0d0
	843	wsL = 0.0d0
	844
	845	! write (,) ' MZTUD/ iirw = ', iirw
	846
	847
	848	c esto es para que las subroutines de mztfsub calculen we
	849	c de la forma apropiada para mztf, no para fot
	850	icls=icls_mztf
	851
	852	c codigos para filenames
	853	! if (isot .eq. 1) isotcode = '26'
	854	! if (isot .eq. 2) isotcode = '28'
	855	! if (isot .eq. 3) isotcode = '36'
	856	! if (isot .eq. 4) isotcode = '27'
	857	! if (isot .eq. 5) isotcode = '62'
	858	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	859	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	860	! write (ibcode1,101) ib
	861	! else
	862	! write (ibcode2,202) ib
	863	! endif
	864	! write (*,'( 30h calculating curtis matrix : ,2x,
	865	! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot
	866
	867	c integration in angle !!!!!!!!!!!!!!!!!!!!
	868	c------- diffusivity approx.
	869	if (iimu.eq.1) then
	870	! write (,) ' diffusivity approx. beta = ',beta
	871	mu(1) = 1.0d0
	872	amu(1)= 1.0d0
	873	c-------data for 8 points integration
	874	elseif (iimu.eq.4) then
	875	write (,)' 4 points for the gauss-legendre angle quadrature.'
	876	mu(1)=(1.0d0+0.339981043584856)/2.0d0
	877	mu(2)=(1.0d0-0.339981043584856)/2.0d0
	878	mu(3)=(1.0d0+0.861136311594053)/2.0d0
	879	mu(4)=(1.0d0-0.861136311594053)/2.0d0
	880	amu(1)=0.652145154862546
	881	amu(2)=amu(1)
	882	amu(3)=0.347854845137454
	883	amu(4)=amu(3)
	884	beta=1.0d0
	885	c-------data for 8 points integration
	886	elseif(iimu.eq.8) then
	887	write (,)' 8 points for the gauss-legendre angle quadrature.'
	888	mu(1)=(1.0d0+0.183434642495650)/2.0d0
	889	mu(2)=(1.0d0-0.183434642495650)/2.0d0
	890	mu(3)=(1.0d0+0.525532409916329)/2.0d0
	891	mu(4)=(1.0d0-0.525532409916329)/2.0d0
	892	mu(5)=(1.0d0+0.796666477413627)/2.0d0
	893	mu(6)=(1.0d0-0.796666477413627)/2.0d0
	894	mu(7)=(1.0d0+0.960289856497536)/2.0d0
	895	mu(8)=(1.0d0-0.960289856497536)/2.0d0
	896	amu(1)=0.362683783378362
	897	amu(2)=amu(1)
	898	amu(3)=0.313706645877887
	899	amu(4)=amu(3)
	900	amu(5)=0.222381034453374
	901	amu(6)=amu(5)
	902	amu(7)=0.101228536290376
	903	amu(8)=amu(7)
	904	beta=1.0d0
	905	end if
	906	c!!!!!!!!!!!!!!!!!!!!!!!
	907
	908	ccc
	909	ccc determine abundances included in the absorber amount
	910	ccc
	911
	912	c first, set up the grid ready for interpolation.
	913	do i=1,nzy
	914	zyd(i) = dble(zy(i))
	915	enddo
	916	do i=1,nl
	917	zld(i) = dble(zl(i))
	918	enddo
	919	c vibr. temp of the bending mode :
[695]	920	if (isot.eq.1) call interdp(tvtbs,zyd,nzy, v626t1,zld,nl,1)
	921	if (isot.eq.2) call interdp(tvtbs,zyd,nzy, v628t1,zld,nl,1)
	922	if (isot.eq.3) call interdp(tvtbs,zyd,nzy, v636t1,zld,nl,1)
	923	if (isot.eq.4) call interdp(tvtbs,zyd,nzy, v627t1,zld,nl,1)
[498]	924	!if (isot.eq.5) call interdp ( tvtbs,zxd,nz, vcot1,zld,nl, 1 )
	925
	926	c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
	927	c por similitud a la que se hace en cza.for ; esto solo se hace para CO2
	928	!write (,) 'imr(isot) = ', isot, imr(isot)
	929	do i=1,nzy
	930	if (isot.eq.5) then
	931	con(i) = dble( coy(i) * imrco )
	932	else
	933	con(i) = dble( co2y(i) * imr(isot) )
	934	correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )
	935	con(i) = con(i) * ( 1.d0 - correc )
	936	! write (,) ' iz, correc, co2y(i), con(i) =',
	937	! @ i,correc,co2y(i),con(i)
	938	endif
	939
	940	!-----------------------------------------------------------------
	941	! mlp & cristina. 17 july 1996 change the calculation of mr.
	942	! it is used for calculating partial press
	943	! alpha = alpha(self,co2)pp +alpha(n2)(pt-pp)
	944	! for an isotope, if mr is obtained by
	945	! co2*imr(iso)/nt
	946	! we are considerin collisions with other co2 isotopes
	947	! (including the major one, 626) as if they were with n2.
	948	! assuming mr as co2/nt, we consider collisions
	949	! of type 628-626 as of 626-626 instead of as 626-n2.
	950	! mrx(i)=con(i)/ntx(i) ! old malv
	951	! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs
	952
	953	! jan 98:
	954	! esta modif de mlp implica anular el correc (deberia revisar esto)
	955
	956	mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98
	957
	958	!-----------------------------------------------------------------
	959
	960	end do
	961
	962	! como beta y 1.d5 son comunes a todas las weighted absorber amounts,
	963	! los simplificamos:
	964	! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )
	965	!write (,) ' con(nz), con(nz-1) =', con(nz), con(nz-1)
	966	coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
	967	!write (,) ' coninf =', coninf
	968
	969	ccc
	970	ccc temp dependence of the band strength and
	971	ccc nlte correction factor for the absorber amount
	972	ccc
	973	call mztf_correccion ( coninf, con, ib, isot, itableout )
	974	ccc
	975	ccc reads histogrammed spectral data (strength for lte and vmr=1)
	976	ccc
	977	!hfile1 = dirspec//'hi'//dn !Ya no hacemos distincion d/n en esto
	978	! hfile1 = dirspec//'hid' !(see why in his.for)
	979	! hfile1='hid'
	980	!! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
	981	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
	982
	983	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	984	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	985	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat'
	986	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat'
	987	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat'
	988	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat'
	989	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
	990	! else
	991	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat'
	992	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat'
	993	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat'
	994	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat'
	995	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat'
	996	! endif
	997	if(ib.eq.1) then
	998	if(isot.eq.1) then !Case 1
	999	mm=mm_c1
	1000	nbox=nbox_c1
	1001	tmin=tmin_c1
	1002	tmax=tmax_c1
	1003	do i=1,nbox_max
	1004	no(i)=no_c1(i)
	1005	dist(i)=dist_c1(i)
	1006	do j=1,nhist
	1007	sk1(j,i)=sk1_c1(j,i)
	1008	xls1(j,i)=xls1_c1(j,i)
	1009	xln1(j,i)=xln1_c1(j,i)
	1010	xld1(j,i)=xld1_c1(j,i)
	1011	enddo
	1012	enddo
	1013	do j=1,nhist
	1014	thist(j)=thist_c1(j)
	1015	enddo
	1016	else if(isot.eq.2) then !Case 2
	1017	mm=mm_c2
	1018	nbox=nbox_c2
	1019	tmin=tmin_c2
	1020	tmax=tmax_c2
	1021	do i=1,nbox_max
	1022	no(i)=no_c2(i)
	1023	dist(i)=dist_c2(i)
	1024	do j=1,nhist
	1025	sk1(j,i)=sk1_c2(j,i)
	1026	xls1(j,i)=xls1_c2(j,i)
	1027	xln1(j,i)=xln1_c2(j,i)
	1028	xld1(j,i)=xld1_c2(j,i)
	1029	enddo
	1030	enddo
	1031	do j=1,nhist
	1032	thist(j)=thist_c2(j)
	1033	enddo
	1034	else if(isot.eq.3) then !Case 3
	1035	mm=mm_c3
	1036	nbox=nbox_c3
	1037	tmin=tmin_c3
	1038	tmax=tmax_c3
	1039	do i=1,nbox_max
	1040	no(i)=no_c3(i)
	1041	dist(i)=dist_c3(i)
	1042	do j=1,nhist
	1043	sk1(j,i)=sk1_c3(j,i)
	1044	xls1(j,i)=xls1_c3(j,i)
	1045	xln1(j,i)=xln1_c3(j,i)
	1046	xld1(j,i)=xld1_c3(j,i)
	1047	enddo
	1048	enddo
	1049	do j=1,nhist
	1050	thist(j)=thist_c3(j)
	1051	enddo
	1052	else if(isot.eq.4) then !Case 4
	1053	mm=mm_c4
	1054	nbox=nbox_c4
	1055	tmin=tmin_c4
	1056	tmax=tmax_c4
	1057	do i=1,nbox_max
	1058	no(i)=no_c4(i)
	1059	dist(i)=dist_c4(i)
	1060	do j=1,nhist
	1061	sk1(j,i)=sk1_c4(j,i)
	1062	xls1(j,i)=xls1_c4(j,i)
	1063	xln1(j,i)=xln1_c4(j,i)
	1064	xld1(j,i)=xld1_c4(j,i)
	1065	enddo
	1066	enddo
	1067	do j=1,nhist
	1068	thist(j)=thist_c4(j)
	1069	enddo
	1070	else
	1071	write(,)'isot must be 2,3 or 4 for ib=1!!'
	1072	write(,)'stop at mztud/324'
	1073	stop
	1074	endif
	1075	else if (ib.eq.2) then
	1076	if(isot.eq.1) then !Case 5
	1077	mm=mm_c5
	1078	nbox=nbox_c5
	1079	tmin=tmin_c5
	1080	tmax=tmax_c5
	1081	do i=1,nbox_max
	1082	no(i)=no_c5(i)
	1083	dist(i)=dist_c5(i)
	1084	do j=1,nhist
	1085	sk1(j,i)=sk1_c5(j,i)
	1086	xls1(j,i)=xls1_c5(j,i)
	1087	xln1(j,i)=xln1_c5(j,i)
	1088	xld1(j,i)=xld1_c5(j,i)
	1089	enddo
	1090	enddo
	1091	do j=1,nhist
	1092	thist(j)=thist_c5(j)
	1093	enddo
	1094	else
	1095	write(,)'isot must be 1 for ib=2!!'
	1096	write(,)'stop at mztud/348'
	1097	stop
	1098	endif
	1099	else if (ib.eq.3) then
	1100	if(isot.eq.1) then !Case 6
	1101	mm=mm_c6
	1102	nbox=nbox_c6
	1103	tmin=tmin_c6
	1104	tmax=tmax_c6
	1105	do i=1,nbox_max
	1106	no(i)=no_c6(i)
	1107	dist(i)=dist_c6(i)
	1108	do j=1,nhist
	1109	sk1(j,i)=sk1_c6(j,i)
	1110	xls1(j,i)=xls1_c6(j,i)
	1111	xln1(j,i)=xln1_c6(j,i)
	1112	xld1(j,i)=xld1_c6(j,i)
	1113	enddo
	1114	enddo
	1115	do j=1,nhist
	1116	thist(j)=thist_c6(j)
	1117	enddo
	1118	else
	1119	write(,)'isot must be 1 for ib=3!!'
	1120	write(,)'stop at mztud/372'
	1121	stop
	1122	endif
	1123	else if (ib.eq.4) then
	1124	if(isot.eq.1) then !Case 7
	1125	mm=mm_c7
	1126	nbox=nbox_c7
	1127	tmin=tmin_c7
	1128	tmax=tmax_c7
	1129	do i=1,nbox_max
	1130	no(i)=no_c7(i)
	1131	dist(i)=dist_c7(i)
	1132	do j=1,nhist
	1133	sk1(j,i)=sk1_c7(j,i)
	1134	xls1(j,i)=xls1_c7(j,i)
	1135	xln1(j,i)=xln1_c7(j,i)
	1136	xld1(j,i)=xld1_c7(j,i)
	1137	enddo
	1138	enddo
	1139	do j=1,nhist
	1140	thist(j)=thist_c7(j)
	1141	enddo
	1142	else
	1143	write(,)'isot must be 1 for ib=4!!'
	1144	write(,)'stop at mztud/396'
	1145	stop
	1146	endif
	1147	else
	1148	write(,)'ib must be 1,2,3 or 4!!'
	1149	write(,)'stop at mztud/401'
	1150	endif
	1151
	1152
	1153
	1154
	1155	! write (,) 'hisfile: ', hisfile
	1156	! the argument to rhist is to make this compatible with mztf_comp.f,
	1157	! which is a useful modification of mztf.f (to change strengths of bands
	1158	! call rhist (1.0)
	1159	if (isot.ne.5) deltanux = deltanu(isot,ib)
	1160	if (isot.eq.5) deltanux = deltanuco
	1161
	1162	c******
	1163	c****** calculation of tauinf(nl)
	1164	c******
	1165	call initial
	1166	ff=1.0e10
	1167
	1168	do i=nl,1,-1
	1169
	1170	if(i.eq.nl)then
	1171
	1172	call intz (zl(i),c2,p2,mr2,t2, con)
	1173	do kr=1,nbox
	1174	ta(kr)=t2
	1175	end do
	1176	! write (,) ' i, t2 =', i, t2
	1177	call interstrength (st2,t2,ka,ta)
	1178	aa = p2 * coninf * mr2 * (st2 * ff)
	1179	bb = p2 * coninf * st2
	1180	cc = coninf * st2
	1181	dd = t2 * coninf * st2
	1182	do kr=1,nbox
	1183	ccbox(kr) = coninf * ka(kr)
	1184	ddbox(kr) = t2 * ccbox(kr)
	1185	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
	1186	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	1187	end do
	1188	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
	1189	c2 = c2 * st2 * dble(deltaz)
	1190
	1191	else
	1192	call intz (zl(i),c1,p1,mr1,t1, con)
	1193	do kr=1,nbox
	1194	ta(kr)=t1
	1195	end do
	1196	! write (,) ' i, t1 =', i, t1
	1197	call interstrength (st1,t1,ka,ta)
	1198	do kr=1,nbox
	1199	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
	1200	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	1201	end do
	1202	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
	1203	c1 = c1 * st1 * dble(deltaz)
	1204	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	1205	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	1206	cc = cc + ( c1 + c2 ) / 2.d0
	1207	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	1208	do kr=1,nbox
	1209	ccbox(kr) = ccbox(kr) +
	1210	@ ( c1box(kr) + c2box(kr) )/2.d0
	1211	ddbox(kr) = ddbox(kr) +
	1212	@ ( t1c1box(kr)+t2c2box(kr) )/2.d0
	1213	end do
	1214
	1215	mr2 = mr1
	1216	c2=c1
	1217	do kr=1,nbox
	1218	c2box(kr) = c1box(kr)
	1219	end do
	1220	t2=t1
	1221	p2=p1
	1222	end if
	1223
	1224	pt = bb / cc
	1225	pp = aa / (cc*ff)
	1226
	1227	! ta=dd/cc
	1228	! tdop = ta
	1229	ts = dd/cc
	1230	do kr=1,nbox
	1231	ta(kr) = ddbox(kr) / ccbox(kr)
	1232	end do
	1233	! write (,) ' i, ts =', i, ts
	1234	call interstrength(st,ts,ka,ta)
	1235	! call intershape(alsa,alna,alda,tdop)
	1236	call intershape(alsa,alna,alda,ta)
	1237	* ua = cc/st
	1238
	1239	c next loop calculates the eqw for an especified path uapp,pt,ta
	1240
	1241	eqwmu = 0.0d0
	1242	do im = 1,iimu
	1243	eqw=0.0d0
	1244	do kr=1,nbox
	1245	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	1246	if(ua(kr).lt.0.)write(,)'mztud/504',ua(kr),ccbox(kr),
	1247	$ ka(kr),beta,mu(im),kr,im,i,nl
	1248	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
	1249	if ( i_supersat .eq. 0 ) then
	1250	eqw=eqw+no(kr)*w
	1251	else
	1252	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	1253	endif
	1254	end do
	1255	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	1256	end do
	1257
	1258	tauinf(i) = exp( - eqwmu / dble(deltanux) )
	1259
	1260	end do
	1261	! if ( isot.eq.1 .and. ib.eq.2 ) then
	1262	! write (,) ' tauinf(nl) = ', tauinf(nl)
	1263	! write (,) ' tauinf(1) = ', tauinf(1)
	1264	! endif
	1265
	1266	c******
	1267	c****** calculation of tau(in,ir) for n<=r
	1268	c******
	1269
	1270	do 1 in=1,nl-1
	1271	call initial
	1272	call intz (zl(in), c1,p1,mr1,t1, con)
	1273	do kr=1,nbox
	1274	ta(kr) = t1
	1275	end do
	1276	call interstrength (st1,t1,ka,ta)
	1277	do kr=1,nbox
	1278	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
	1279	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	1280	end do
	1281	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
	1282	c1 = c1 * st1 * dble(deltaz)
	1283
	1284	do 2 ir=in,nl-1
	1285
	1286	if (ir.eq.in) then
	1287	tau(in,ir) = 1.d0
	1288	goto 2
	1289	end if
	1290
	1291	call intz (zl(ir), c2,p2,mr2,t2, con)
	1292	do kr=1,nbox
	1293	ta(kr) = t2
	1294	end do
	1295	call interstrength (st2,t2,ka,ta)
	1296	do kr=1,nbox
	1297	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
	1298	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	1299	end do
	1300	! c2 = c2 * st2 * beta * dble(deltaz) * 1.e5
	1301	c2 = c2 * st2 * dble(deltaz)
	1302
	1303	c aa = aa + ( p1mr1c1 + p2mr2c2 ) / 2.d0
	1304	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	1305	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	1306	cc = cc + ( c1 + c2 ) / 2.d0
	1307	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	1308	do kr=1,nbox
	1309	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
	1310	ddbox(kr) = ddbox(kr) +
	1311	$ ( t1c1box(kr) + t2c2box(kr) ) /2.d0
	1312	end do
	1313
	1314	mr1=mr2
	1315	t1=t2
	1316	c1=c2
	1317	p1=p2
	1318	do kr=1,nbox
	1319	c1box(kr) = c2box(kr)
	1320	end do
	1321
	1322	pt = bb / cc
	1323	pp = aa / (cc * ff)
	1324
	1325	* ta=dd/cc
	1326	* tdop = ta
	1327	ts = dd/cc
	1328	do kr=1,nbox
	1329	ta(kr) = ddbox(kr) / ccbox(kr)
	1330	end do
	1331	call interstrength(st,ts,ka,ta)
	1332	call intershape(alsa,alna,alda,ta)
	1333	* ua = cc/st
	1334
	1335	c next loop calculates the eqw for an especified path ua,pp,pt,ta
	1336
	1337	eqwmu = 0.0d0
	1338	do im = 1,iimu
	1339	eqw=0.0d0
	1340	do kr=1,nbox
	1341	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	1342
	1343	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
	1344	if ( i_supersat .eq. 0 ) then
	1345	eqw=eqw+no(kr)*w
	1346	else
	1347	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	1348	endif
	1349	end do
	1350	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	1351	end do
	1352
	1353	tau(in,ir) = exp( - eqwmu / dble(deltanux) )
	1354
	1355	2 continue
	1356
	1357	1 continue
	1358	! if ( isot.eq.1 .and. ib.eq.2 ) then
	1359	! write (,) ' tau(1,) , =1,20 '
	1360	! write (,) ( sngl(tau(1,k)), k=1,20 )
	1361	! endif
	1362
	1363
	1364	c**********
	1365	c********** calculation of tau(in,ir) for n>r
	1366	c**********
	1367
	1368	in=nl
	1369
	1370	call initial
	1371	call intz (zl(in), c1,p1,mr1,t1, con)
	1372	do kr=1,nbox
	1373	ta(kr) = t1
	1374	end do
	1375	call interstrength (st1,t1,ka,ta)
	1376	do kr=1,nbox
	1377	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
	1378	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	1379	end do
	1380	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
	1381	c1 = c1 * st1 * dble(deltaz)
	1382
	1383	do 4 ir=in-1,1,-1
	1384
	1385	call intz (zl(ir), c2,p2,mr2,t2, con)
	1386	do kr=1,nbox
	1387	ta(kr) = t2
	1388	end do
	1389	call interstrength (st2,t2,ka,ta)
	1390	do kr=1,nbox
	1391	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
	1392	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	1393	end do
	1394	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
	1395	c2 = c2 * st2 * dble(deltaz)
	1396
	1397	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	1398	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	1399	cc = cc + ( c1 + c2 ) / 2.d0
	1400	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	1401	do kr=1,nbox
	1402	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
	1403	ddbox(kr) = ddbox(kr) +
	1404	$ ( t1c1box(kr) + t2c2box(kr) ) /2.d0
	1405	end do
	1406
	1407	mr1=mr2
	1408	c1=c2
	1409	t1=t2
	1410	p1=p2
	1411	do kr=1,nbox
	1412	c1box(kr) = c2box(kr)
	1413	end do
	1414
	1415	pt = bb / cc
	1416	pp = aa / (cc * ff)
	1417	ts = dd / cc
	1418	do kr=1,nbox
	1419	ta(kr) = ddbox(kr) / ccbox(kr)
	1420	end do
	1421	call interstrength (st,ts,ka,ta)
	1422	call intershape (alsa,alna,alda,ta)
	1423
	1424	* ua = cc/st
	1425
	1426	c next loop calculates the eqw for an especified path ua,pp,pt,ta
	1427
	1428	eqwmu = 0.0d0
	1429	do im = 1,iimu
	1430	eqw=0.0d0
	1431	do kr=1,nbox
	1432	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	1433	if(ua(kr).lt.0.)write(,)'mztud/691',ua(kr),ccbox(kr),
	1434	$ ka(kr),beta,mu(im),kr,im,i,nl
	1435
	1436	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
	1437	if ( i_supersat .eq. 0 ) then
	1438	eqw=eqw+no(kr)*w
	1439	else
	1440	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	1441	endif
	1442	end do
	1443	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	1444	end do
	1445
	1446	tau(in,ir) = exp( - eqwmu / dble(deltanux) )
	1447
	1448	4 continue
	1449
	1450	c
	1451	c due to the simmetry of the transmittances
	1452	c
	1453	do in=nl-1,2,-1
	1454	do ir=in-1,1,-1
	1455	tau(in,ir) = tau(ir,in)
	1456	end do
	1457	end do
	1458
	1459
	1460	ccc
	1461	ccc writing out transmittances
	1462	ccc
	1463	if (itauout.eq.1) then
	1464
	1465	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	1466	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	1467	! open( 1, file=
	1468	! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat',
	1469	! @ access='sequential', form='unformatted' )
	1470	! else
	1471	! open( 1, file=
	1472	! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat',
	1473	! @ access='sequential', form='unformatted' )
	1474	! endif
	1475
	1476	! write(1) dummy
	1477	! write(1)' format: (tauinf(n),(tau(n,r),r=1,nl),n=1,nl)'
	1478	! do in=1,nl
	1479	! write (1) tauinf(in), ( tau(in,ir), ir=1,nl )
	1480	! end do
	1481	! close(unit=1)
	1482
	1483	elseif (itauout.eq.2) then
	1484
	1485	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	1486	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	1487	! open( 1, file=
	1488	! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat')
	1489	! else
	1490	! open( 1, file=
	1491	! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat')
	1492	! endif
	1493
	1494	! !write(1,*) dummy
	1495	! !write(1,*) 'tij for curtis matrix calculations '
	1496	! !write(1,*)' cira mars model atmosphere '
	1497	! !write(1,*)' beta= ',beta,'deltanu= ',deltanux
	1498	! write(1,*) nl
	1499	! write(1,*)
	1500	! @ ' format: (tauinf(in),(tau(in,ir),ir=1,nl),in=1,nl)'
	1501
	1502	! do in=1,nl
	1503	! write (1,*) tauinf(in)
	1504	! write (1,*) (tau(in,ir), ir=1,nl)
	1505	! end do
	1506	! close(unit=1)
	1507
	1508	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	1509	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	1510	! write (*,'(1x, 31htransmitances written out in: ,a22)')
	1511	! @ 'taul'//isotcode//dn//ibcode1
	1512	! else
	1513	! write (*,'(1x, 31htransmitances written out in: ,a22)')
	1514	! @ 'taul'//isotcode//dn//ibcode2
	1515	! endif
	1516
	1517	end if
	1518
	1519	c cleaning of transmittances
	1520	! call elimin_tau(tau,tauinf,nl,nan,itableout,nw,dummy,
	1521	! @ isotcode,dn,ibcode2)
	1522
	1523	c construction of the curtis matrix
	1524
	1525	call mzcud ( tauinf,tau, cf,cfup,cfdw, vc,taugr,
	1526	@ ib,isot,icfout,itableout )
	1527
	1528	c end
	1529	return
	1530	end
	1531
	1532
	1533
	1534
	1535
	1536	c***********************************************************************
	1537	c mzcud.f
	1538	c***********************************************************************
	1539
	1540	subroutine mzcud( tauinf,tau, c,cup,cdw,vc,taugr,
	1541	@ ib,isot,icfout,itableout )
	1542
	1543	c old times mlp first version of mzcf
	1544	c a.k.murphy method to avoid extrapolation in the curtis matrix
	1545	c feb-89 malv AKM method to avoid extrapolation in C.M.
	1546	c 25-sept-96 cristina dejar las matrices en doble precision
	1547	c jan 98 malv version para mz1d
	1548	c oct 01 malv update version for fluxes for hb and fb
	1549	c jul 2011 malv+fgg Adapted to LMD-MGCM
	1550	c***********************************************************************
	1551
	1552	implicit none
	1553
	1554	include 'comcstfi.h'
	1555	include 'nlte_paramdef.h'
	1556	include 'nlte_commons.h'
	1557
	1558	c arguments
	1559	real*8 c(nl,nl), cup(nl,nl), cdw(nl,nl) ! o
	1560	real*8 vc(nl), taugr(nl) ! o
	1561	real*8 tau(nl,nl) ! i
	1562	real*8 tauinf(nl) ! i
	1563	integer ib ! i
	1564	integer isot ! i
	1565	integer icfout, itableout ! i
	1566
	1567	c external
	1568	external bandid
	1569	character*2 bandid
	1570
	1571	c local variables
	1572	integer i, in, ir, iw, itblout
	1573	real*8 cfup(nl,nl), cfdw(nl,nl)
	1574	real*8 a(nl,nl), cf(nl,nl)
	1575	character isotcode2, bcode2
	1576
	1577	c formats
	1578	101 format(i1)
	1579	202 format(i2)
	1580	180 format(a80)
	1581	181 format(a80)
	1582	c***********************************************************************
	1583
	1584	if (isot.eq.1) isotcode = '26'
	1585	if (isot.eq.2) isotcode = '28'
	1586	if (isot.eq.3) isotcode = '36'
	1587	if (isot.eq.4) isotcode = '27'
	1588	if (isot.eq.5) isotcode = 'co'
	1589	bcode = bandid( ib )
	1590
	1591	! write (,) ' '
	1592
	1593	do in=1,nl
	1594
	1595	do ir=1,nl
	1596
	1597	cf(in,ir) = 0.0d0
	1598	cfup(in,ir) = 0.0d0
	1599	cfdw(in,ir) = 0.0d0
	1600	c(in,ir) = 0.0d0
	1601	cup(in,ir) = 0.0d0
	1602	cdw(in,ir) = 0.0d0
	1603	a(in,ir) = 0.0d0
	1604
	1605	end do
	1606
	1607	vc(in) = 0.0d0
	1608	taugr(in) = 0.0d0
	1609
	1610	end do
	1611
	1612
	1613	c the next lines are a reduced and equivalent way of calculating
	1614	c the c(in,ir) elements for n=2,nl1 and r=1,nl
	1615
	1616
	1617	c do in=2,nl1
	1618	c do ir=1,nl
	1619	c if(ir.eq.1)then
	1620	c c(in,ir)=tau(in-1,1)-tau(in+1,1)
	1621	c elseif(ir.eq.nl)then
	1622	c c(in,ir)=tau(in+1,nl1)-tauinf(in+1)-tau(in-1,nl1)+tauinf(in-1)
	1623	c else
	1624	c c(in,ir)=tau(in+1,ir-1)-tau(in+1,ir)-tau(in-1,ir-1)+tau(in-1,ir)
	1625	c end if
	1626	c c(in,ir)=c(in,ir)pideltanu(ib)/(2.deltaz1.0e5)
	1627	c end do
	1628	c end do
	1629	c go to 1000
	1630
	1631	c calculation of the matrix cfup(nl,nl)
	1632
	1633	cfup(1,1) = 1.d0 - tau(1,1)
	1634
	1635	do in=2,nl
	1636	do ir=1,in
	1637
	1638	if (ir.eq.1) then
	1639	cfup(in,ir) = tau(in,ir) - tau(in,1)
	1640	elseif (ir.eq.in) then
	1641	cfup(in,ir) = 1.d0 - tau(in,ir-1)
	1642	else
	1643	cfup(in,ir) = tau(in,ir) - tau(in,ir-1)
	1644	end if
	1645
	1646	end do
	1647	end do
	1648
	1649	! contribution to upwards fluxes from bb at bottom :
	1650	do in=1,nl
	1651	taugr(in) = tau(in,1)
	1652	enddo
	1653
	1654	c calculation of the matrix cfdw(nl,nl)
	1655
	1656	cfdw(nl,nl) = 1.d0 - tauinf(nl)
	1657
	1658	do in=1,nl-1
	1659	do ir=in,nl
	1660
	1661	if (ir.eq.in) then
	1662	cfdw(in,ir) = 1.d0 - tau(in,ir)
	1663	elseif (ir.eq.nl) then
	1664	cfdw(in,ir) = tau(in,ir-1) - tauinf(in)
	1665	else
	1666	cfdw(in,ir) = tau(in,ir-1) - tau(in,ir)
	1667	end if
	1668
	1669	end do
	1670	end do
	1671
	1672
	1673	c calculation of the matrix cf(nl,nl)
	1674
	1675	do in=1,nl
	1676	do ir=1,nl
	1677
	1678	if (ir.eq.1) then
	1679	! version con l_bb(tg) = l_bb(t(1))=j(1) (see also vc below)
	1680	! cf(in,ir) = tau(in,ir)
	1681	! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see also vc below)
	1682	cf(in,ir) = tau(in,ir) - tau(in,1)
	1683	elseif (ir.eq.nl) then
	1684	cf(in,ir) = tauinf(in) - tau(in,ir-1)
	1685	else
	1686	cf(in,ir) = tau(in,ir) - tau(in,ir-1)
	1687	end if
	1688
	1689	end do
	1690	end do
	1691
	1692
	1693	c definition of the a(nl,nl) matrix
	1694
	1695	do in=2,nl-1
	1696	do ir=1,nl
	1697	if (ir.eq.in+1) a(in,ir) = -1.d0
	1698	if (ir.eq.in-1) a(in,ir) = +1.d0
	1699	a(in,ir) = a(in,ir) / ( 2.d0deltaz1.d5 )
	1700	end do
	1701	end do
	1702	! this is not needed anymore in the akm scheme
	1703	! a(1,1) = +3.d0
	1704	! a(1,2) = -4.d0
	1705	! a(1,3) = +1.d0
	1706	! a(nl,nl) = -3.d0
	1707	! a(nl,nl1) = +4.d0
	1708	! a(nl,nl2) = -1.d0
	1709
	1710	c calculation of the final curtis matrix ("reduced" by murphy's method)
	1711
	1712	if (isot.ne.5) then
	1713	do in=1,nl
	1714	do ir=1,nl
	1715	cf(in,ir) = cf(in,ir) * pi*deltanu(isot,ib)
	1716	cfup(in,ir) = cfup(in,ir) * pi*deltanu(isot,ib)
	1717	cfdw(in,ir) = cfdw(in,ir) * pi*deltanu(isot,ib)
	1718	end do
	1719	taugr(in) = taugr(in) * pi*deltanu(isot,ib)
	1720	end do
	1721	else
	1722	do in=1,nl
	1723	do ir=1,nl
	1724	cf(in,ir) = cf(in,ir) * pi*deltanuco
	1725	enddo
	1726	taugr(in) = taugr(in) * pi*deltanuco
	1727	enddo
	1728	endif
	1729
	1730	do in=2,nl-1
	1731
	1732	do ir=1,nl
	1733
	1734	do i=1,nl
	1735	! only c contains the matrix a. matrixes cup,cdw dont because
	1736	! these two will be used for flux calculations, not
	1737	! only for flux divergencies
	1738
	1739	c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir)
	1740	! from this matrix we will extract (see below) the
	1741	! nl2 x nl2 "core" for the "reduced" final curtis matrix.
	1742
	1743	end do
	1744	cup(in,ir) = cfup(in,ir)
	1745	cdw(in,ir) = cfdw(in,ir)
	1746
	1747	end do
	1748	! version con l_bb(tg) = l_bb(t(1))=j(1) (see cf above)
	1749	!vc(in) = c(in,1)
	1750	! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see cf above)
	1751	if (isot.ne.5) then
	1752	vc(in) = pideltanu(isot,ib)/( 2.d0deltaz1.d5 )
	1753	@ ( tau(in-1,1) - tau(in+1,1) )
	1754	else
	1755	vc(in) = pideltanuco/( 2.d0deltaz1.d5 )
	1756	@ ( tau(in-1,1) - tau(in+1,1) )
	1757	endif
	1758
	1759	end do
	1760
	1761	5 continue
	1762
	1763	! write (,) 'mztf/1/ c(2,*) =', (c(2,i), i=1,nl)
	1764
	1765	! call elimin_dibuja(c,nl,itableout)
	1766
	1767	c ventana del smoothing de c es nw=3 y de vc es 5 (puesto en lisa):
	1768	c subroutine elimin_mz4(c,vc,ilayer,nl,nan,iw, nw)
	1769
	1770	iw = nan
	1771	if (isot.eq.4) iw = 5 ! eliminates values < 1.d-19
	1772	if (itableout.eq.30) then
	1773	itblout = 0
	1774	else
	1775	itblout = itableout
	1776	endif
	1777	call elimin_mz1d (c,vc,0,iw,itblout,nw)
	1778
	1779	! upper boundary condition
	1780	! j'(nl) = j'(nl1) ==> j(nl) = 2j(nl1) - j(nl2) ==>
	1781	do in=2,nl-1
	1782	c(in,nl-2) = c(in,nl-2) - c(in,nl)
	1783	c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl)
	1784	cup(in,nl-2) = cup(in,nl-2) - cup(in,nl)
	1785	cup(in,nl-1) = cup(in,nl-1) + 2.d0*cup(in,nl)
	1786	cdw(in,nl-2) = cdw(in,nl-2) - cdw(in,nl)
	1787	cdw(in,nl-1) = cdw(in,nl-1) + 2.d0*cdw(in,nl)
	1788	end do
	1789	! j(nl) = j(nl1) ==>
	1790	! do in=2,nl1
	1791	! c(in,nl1) = c(in,nl1) + c(in,nl)
	1792	! end do
	1793
	1794	! 1000 continue
	1795
	1796
	1797	if (icfout.eq.1) then
	1798
	1799	! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then
	1800	! codmatrx = codmatrx_fb
	1801	! else
	1802	! codmatrx = codmatrx_hot
	1803	! end if
	1804	! if (ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	1805	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	1806	! ibcode2 = '0'//ibcode1
	1807	! else
	1808	! write ( ibcode2, 202) ib
	1809	! endif
	1810
	1811	! open ( 1, access='sequential', form='unformatted', file=
	1812	! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat')
	1813	! open ( 2, access='sequential', form='unformatted', file=
	1814	! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat')
	1815	! open ( 3, access='sequential', form='unformatted', file=
	1816	! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat')
	1817	! open ( 4, access='sequential', form='unformatted', file=
	1818	! @ dircurtis//'cflgr'//isotcode//dn//ibcode2//codmatrx//'.dat')
	1819
	1820	! write(1) dummy
	1821	! write(1) ' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)'
	1822	! do in=2,nl-1
	1823	! write(1) vc(in), (c(in,ir) , ir=2,nl-1 )
	1824	!! write (,) in, vc(in)
	1825	! end do
	1826
	1827	! write(2) dummy
	1828	! write(2)' format: (cfup(n,r),r=1,nl), n=1,nl)'
	1829	! do in=1,nl
	1830	! write(2) ( cup(in,ir) , ir=1,nl )
	1831	! end do
	1832
	1833	! write(3) dummy
	1834	! write(3)' format: (cfdw(n,r),r=1,nl), n=1,nl)'
	1835	! do in=1,nl
	1836	! write(3) (cdw(in,ir) , ir=1,nl )
	1837	! end do
	1838
	1839	! write(4) dummy
	1840	! write(4)' format: (taugr(n), n=1,nl)'
	1841	! do in=1,nl
	1842	! write(4) (taugr(in), ir=1,nl )
	1843	! end do
	1844	! !write (,) ' Last value in file: ', taugr(nl)
	1845
	1846	! write (*,'(1x,30hcurtis matrix written out in: ,a,a,a,a)' )
	1847	! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat',
	1848	! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat',
	1849	! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat',
	1850	! @ dircurtis//'cflgr'//isotcode//dn//ibcode2//codmatrx//'.dat'
	1851
	1852	! close (1)
	1853	! close (2)
	1854	! close (3)
	1855	! close (4)
	1856
	1857	else
	1858
	1859	! write (,) ' no curtis matrix output file ', char(10)
	1860
	1861	end if
	1862
	1863	if (itableout.eq.30) then ! Force output of C.M. in ascii file
	1864
	1865	! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then
	1866	! codmatrx = codmatrx_fb
	1867	! else
	1868	! codmatrx = codmatrx_hot
	1869	! end if
	1870	! if (ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	1871	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	1872	! ibcode2 = '0'//ibcode1
	1873	! else
	1874	! write ( ibcode2, 202) ib
	1875	! endif
	1876
	1877	! open (10, file=
	1878	! & dircurtis//'table'//isotcode//dn//ibcode2//codmatrx//'.dat')
	1879	! write(10,*) nl, ' = number of layers '
	1880	! write(10,*) ' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)'
	1881	! do in=2,nl-1
	1882	! write(10,*) vc(in), (c(in,ir) , ir=2,nl-1 )
	1883	! enddo
	1884	! close (10)
	1885	endif
	1886
	1887	c end
	1888	return
	1889	end
	1890
	1891
	1892
	1893
	1894
	1895	c***********************************************************************
	1896	c mztvc
	1897	c***********************************************************************
	1898
	1899	subroutine mztvc ( ig,vc, ib,isot,
	1900	@ iirw,iimu,itauout,icfout,itableout )
	1901
	1902	c jul 2011 malv+fgg
	1903	c***********************************************************************
	1904
	1905	implicit none
	1906
	1907	include 'comcstfi.h'
	1908	include 'nlte_paramdef.h'
	1909	include 'nlte_commons.h'
	1910
	1911	c arguments
	1912	integer ig ! ADDED FOR TRACEBACK
	1913	real*8 cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o
	1914	real*8 vc(nl), taugr(nl) ! o
	1915	integer ib ! i
	1916	integer isot ! i
	1917	integer iirw ! i
	1918	integer iimu ! i
	1919	integer itauout ! i
	1920	integer icfout ! i
	1921	integer itableout ! i
	1922
	1923	c local variables and constants
	1924	integer i, in, ir, im, k ,j
	1925	integer nmu
	1926	parameter (nmu = 8)
	1927	real*8 tau(nl,nl)
	1928	real*8 tauinf(nl)
	1929	real*8 con(nzy), coninf
	1930	real*8 c1, c2
	1931	real*8 t1, t2
	1932	real*8 p1, p2
	1933	real*8 mr1, mr2
	1934	real*8 st1, st2
	1935	real*8 c1box(70), c2box(70)
	1936	real*8 ff ! to avoid too small numbers
	1937	real*8 tvtbs(nzy)
	1938	real*8 st, beta, ts, eqwmu
	1939	real*8 mu(nmu), amu(nmu)
	1940	real*8 zld(nl), zyd(nzy)
	1941	real*8 correc
	1942	real deltanux ! width of vib-rot band (cm-1)
	1943	character isotcode*2
	1944	integer idummy
	1945	real*8 Desp,wsL
	1946
	1947	c formats
	1948	111 format(a1)
	1949	112 format(a2)
	1950	101 format(i1)
	1951	202 format(i2)
	1952	180 format(a80)
	1953	181 format(a80)
	1954	c***********************************************************************
	1955
	1956	c some needed values
	1957	! rl=sqrt(log(2.d0))
	1958	! pi2 = 3.14159265358989d0
	1959	beta = 1.8d0
	1960	! beta = 1.0d0
	1961	idummy = 0
	1962	Desp = 0.0d0
	1963	wsL = 0.0d0
	1964
	1965	!write (,) ' MZTUD/ iirw = ', iirw
	1966
	1967
	1968	c esto es para que las subroutines de mztfsub calculen we
	1969	c de la forma apropiada para mztf, no para fot
	1970	icls=icls_mztf
	1971
	1972	c codigos para filenames
	1973	! if (isot .eq. 1) isotcode = '26'
	1974	! if (isot .eq. 2) isotcode = '28'
	1975	! if (isot .eq. 3) isotcode = '36'
	1976	! if (isot .eq. 4) isotcode = '27'
	1977	! if (isot .eq. 5) isotcode = '62'
	1978	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	1979	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	1980	! write (ibcode1,101) ib
	1981	! else
	1982	! write (ibcode2,202) ib
	1983	! endif
	1984	! write (*,'( 30h calculating curtis matrix : ,2x,
	1985	! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot
	1986
	1987	c integration in angle !!!!!!!!!!!!!!!!!!!!
	1988
	1989	c------- diffusivity approx.
	1990	if (iimu.eq.1) then
	1991	! write (,) ' diffusivity approx. beta = ',beta
	1992	mu(1) = 1.0d0
	1993	amu(1)= 1.0d0
	1994	c-------data for 8 points integration
	1995	elseif (iimu.eq.4) then
	1996	write (,)' 4 points for the gauss-legendre angle quadrature.'
	1997	mu(1)=(1.0d0+0.339981043584856)/2.0d0
	1998	mu(2)=(1.0d0-0.339981043584856)/2.0d0
	1999	mu(3)=(1.0d0+0.861136311594053)/2.0d0
	2000	mu(4)=(1.0d0-0.861136311594053)/2.0d0
	2001	amu(1)=0.652145154862546
	2002	amu(2)=amu(1)
	2003	amu(3)=0.347854845137454
	2004	amu(4)=amu(3)
	2005	beta=1.0d0
	2006	c-------data for 8 points integration
	2007	elseif(iimu.eq.8) then
	2008	write (,)' 8 points for the gauss-legendre angle quadrature.'
	2009	mu(1)=(1.0d0+0.183434642495650)/2.0d0
	2010	mu(2)=(1.0d0-0.183434642495650)/2.0d0
	2011	mu(3)=(1.0d0+0.525532409916329)/2.0d0
	2012	mu(4)=(1.0d0-0.525532409916329)/2.0d0
	2013	mu(5)=(1.0d0+0.796666477413627)/2.0d0
	2014	mu(6)=(1.0d0-0.796666477413627)/2.0d0
	2015	mu(7)=(1.0d0+0.960289856497536)/2.0d0
	2016	mu(8)=(1.0d0-0.960289856497536)/2.0d0
	2017	amu(1)=0.362683783378362
	2018	amu(2)=amu(1)
	2019	amu(3)=0.313706645877887
	2020	amu(4)=amu(3)
	2021	amu(5)=0.222381034453374
	2022	amu(6)=amu(5)
	2023	amu(7)=0.101228536290376
	2024	amu(8)=amu(7)
	2025	beta=1.0d0
	2026	end if
	2027	c!!!!!!!!!!!!!!!!!!!!!!!
	2028
	2029	ccc
	2030	ccc determine abundances included in the absorber amount
	2031	ccc
	2032
	2033	c first, set up the grid ready for interpolation.
	2034	do i=1,nzy
	2035	zyd(i) = dble(zy(i))
	2036	enddo
	2037	do i=1,nl
	2038	zld(i) = dble(zl(i))
	2039	enddo
	2040
	2041	c vibr. temp of the bending mode :
[695]	2042	if (isot.eq.1) call interdp(tvtbs,zyd,nzy, v626t1,zld,nl,1)
	2043	if (isot.eq.2) call interdp(tvtbs,zyd,nzy, v628t1,zld,nl,1)
	2044	if (isot.eq.3) call interdp(tvtbs,zyd,nzy, v636t1,zld,nl,1)
	2045	if (isot.eq.4) call interdp(tvtbs,zyd,nzy, v627t1,zld,nl,1)
[498]	2046	!if (isot.eq.5) call interdp ( tvtbs,zxd,nz, vcot1,zld,nl, 1 )
	2047
	2048	c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
	2049	c por similitud a la que se hace en cza.for ; esto solo se hace para CO2
	2050
	2051	!write (,) 'imr(isot) = ', isot, imr(isot)
	2052	do i=1,nzy
	2053	if (isot.eq.5) then
	2054	con(i) = dble( coy(i) * imrco )
	2055	else
	2056	con(i) = dble( co2y(i) * imr(isot) )
	2057	correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )
	2058	con(i) = con(i) * ( 1.d0 - correc )
	2059	! write (,) ' iz, correc, co2y(i), con(i) =',
	2060	! @ i,correc,co2y(i),con(i)
	2061	endif
	2062
	2063	!-----------------------------------------------------------------
	2064	! mlp & cristina. 17 july 1996 change the calculation of mr.
	2065	! it is used for calculating partial press
	2066	! alpha = alpha(self,co2)pp +alpha(n2)(pt-pp)
	2067	! for an isotope, if mr is obtained by
	2068	! co2*imr(iso)/nt
	2069	! we are considerin collisions with other co2 isotopes
	2070	! (including the major one, 626) as if they were with n2.
	2071	! assuming mr as co2/nt, we consider collisions
	2072	! of type 628-626 as of 626-626 instead of as 626-n2.
	2073	! mrx(i)=con(i)/ntx(i) ! old malv
	2074	! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs
	2075
	2076	! jan 98:
	2077	! esta modif de mlp implica anular el correc (deberia revisar esto)
	2078
	2079	mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98
	2080
	2081	!-----------------------------------------------------------------
	2082
	2083	end do
	2084
	2085	! como beta y 1.d5 son comunes a todas las weighted absorber amounts,
	2086	! los simplificamos:
	2087	! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )
	2088	!write (,) ' con(nz), con(nz-1) =', con(nz), con(nz-1)
	2089	coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
	2090	!write (,) ' coninf =', coninf
	2091
	2092	ccc
	2093	ccc temp dependence of the band strength and
	2094	ccc nlte correction factor for the absorber amount
	2095	ccc
	2096	call mztf_correccion ( coninf, con, ib, isot, itableout )
	2097
	2098	ccc
	2099	ccc reads histogrammed spectral data (strength for lte and vmr=1)
	2100	ccc
	2101	!hfile1 = dirspec//'hi'//dn !Ya no hacemos distincion d/n en esto
	2102	!! hfile1 = dirspec//'hid' !(see why in his.for)
	2103	! hfile1='hid'
	2104	!! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
	2105	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
	2106
	2107	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	2108	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	2109	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat'
	2110	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat'
	2111	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat'
	2112	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat'
	2113	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
	2114	! else
	2115	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat'
	2116	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat'
	2117	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat'
	2118	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat'
	2119	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat'
	2120	! endif
	2121	! write (,) 'hisfile: ', hisfile
	2122
	2123	! the argument to rhist is to make this compatible with mztf_comp.f,
	2124	! which is a useful modification of mztf.f (to change strengths of bands
	2125	! call rhist (1.0)
	2126	if(ib.eq.1) then
	2127	if(isot.eq.1) then !Case 1
	2128	mm=mm_c1
	2129	nbox=nbox_c1
	2130	tmin=tmin_c1
	2131	tmax=tmax_c1
	2132	do i=1,nbox_max
	2133	no(i)=no_c1(i)
	2134	dist(i)=dist_c1(i)
	2135	do j=1,nhist
	2136	sk1(j,i)=sk1_c1(j,i)
	2137	xls1(j,i)=xls1_c1(j,i)
	2138	xln1(j,i)=xln1_c1(j,i)
	2139	xld1(j,i)=xld1_c1(j,i)
	2140	enddo
	2141	enddo
	2142	do j=1,nhist
	2143	thist(j)=thist_c1(j)
	2144	enddo
	2145	else if(isot.eq.2) then !Case 2
	2146	mm=mm_c2
	2147	nbox=nbox_c2
	2148	tmin=tmin_c2
	2149	tmax=tmax_c2
	2150	do i=1,nbox_max
	2151	no(i)=no_c2(i)
	2152	dist(i)=dist_c2(i)
	2153	do j=1,nhist
	2154	sk1(j,i)=sk1_c2(j,i)
	2155	xls1(j,i)=xls1_c2(j,i)
	2156	xln1(j,i)=xln1_c2(j,i)
	2157	xld1(j,i)=xld1_c2(j,i)
	2158	enddo
	2159	enddo
	2160	do j=1,nhist
	2161	thist(j)=thist_c2(j)
	2162	enddo
	2163	else if(isot.eq.3) then !Case 3
	2164	mm=mm_c3
	2165	nbox=nbox_c3
	2166	tmin=tmin_c3
	2167	tmax=tmax_c3
	2168	do i=1,nbox_max
	2169	no(i)=no_c3(i)
	2170	dist(i)=dist_c3(i)
	2171	do j=1,nhist
	2172	sk1(j,i)=sk1_c3(j,i)
	2173	xls1(j,i)=xls1_c3(j,i)
	2174	xln1(j,i)=xln1_c3(j,i)
	2175	xld1(j,i)=xld1_c3(j,i)
	2176	enddo
	2177	enddo
	2178	do j=1,nhist
	2179	thist(j)=thist_c3(j)
	2180	enddo
	2181	else if(isot.eq.4) then !Case 4
	2182	mm=mm_c4
	2183	nbox=nbox_c4
	2184	tmin=tmin_c4
	2185	tmax=tmax_c4
	2186	do i=1,nbox_max
	2187	no(i)=no_c4(i)
	2188	dist(i)=dist_c4(i)
	2189	do j=1,nhist
	2190	sk1(j,i)=sk1_c4(j,i)
	2191	xls1(j,i)=xls1_c4(j,i)
	2192	xln1(j,i)=xln1_c4(j,i)
	2193	xld1(j,i)=xld1_c4(j,i)
	2194	enddo
	2195	enddo
	2196	do j=1,nhist
	2197	thist(j)=thist_c4(j)
	2198	enddo
	2199	else
	2200	write(,)'isot must be 2,3 or 4 for ib=1!!'
	2201	write(,)'stop at mztvc/310'
	2202	stop
	2203	endif
	2204	else if (ib.eq.2) then
	2205	if(isot.eq.1) then !Case 5
	2206	mm=mm_c5
	2207	nbox=nbox_c5
	2208	tmin=tmin_c5
	2209	tmax=tmax_c5
	2210	do i=1,nbox_max
	2211	no(i)=no_c5(i)
	2212	dist(i)=dist_c5(i)
	2213	do j=1,nhist
	2214	sk1(j,i)=sk1_c5(j,i)
	2215	xls1(j,i)=xls1_c5(j,i)
	2216	xln1(j,i)=xln1_c5(j,i)
	2217	xld1(j,i)=xld1_c5(j,i)
	2218	enddo
	2219	enddo
	2220	do j=1,nhist
	2221	thist(j)=thist_c5(j)
	2222	enddo
	2223	else
	2224	write(,)'isot must be 1 for ib=2!!'
	2225	write(,)'stop at mztvc/334'
	2226	stop
	2227	endif
	2228	else if (ib.eq.3) then
	2229	if(isot.eq.1) then !Case 6
	2230	mm=mm_c6
	2231	nbox=nbox_c6
	2232	tmin=tmin_c6
	2233	tmax=tmax_c6
	2234	do i=1,nbox_max
	2235	no(i)=no_c6(i)
	2236	dist(i)=dist_c6(i)
	2237	do j=1,nhist
	2238	sk1(j,i)=sk1_c6(j,i)
	2239	xls1(j,i)=xls1_c6(j,i)
	2240	xln1(j,i)=xln1_c6(j,i)
	2241	xld1(j,i)=xld1_c6(j,i)
	2242	enddo
	2243	enddo
	2244	do j=1,nhist
	2245	thist(j)=thist_c6(j)
	2246	enddo
	2247	else
	2248	write(,)'isot must be 1 for ib=3!!'
	2249	write(,)'stop at mztvc/358'
	2250	stop
	2251	endif
	2252	else if (ib.eq.4) then
	2253	if(isot.eq.1) then !Case 7
	2254	mm=mm_c7
	2255	nbox=nbox_c7
	2256	tmin=tmin_c7
	2257	tmax=tmax_c7
	2258	do i=1,nbox_max
	2259	no(i)=no_c7(i)
	2260	dist(i)=dist_c7(i)
	2261	do j=1,nhist
	2262	sk1(j,i)=sk1_c7(j,i)
	2263	xls1(j,i)=xls1_c7(j,i)
	2264	xln1(j,i)=xln1_c7(j,i)
	2265	xld1(j,i)=xld1_c7(j,i)
	2266	enddo
	2267	enddo
	2268	do j=1,nhist
	2269	thist(j)=thist_c7(j)
	2270	enddo
	2271	else
	2272	write(,)'isot must be 1 for ib=4!!'
	2273	write(,)'stop at mztvc/382'
	2274	stop
	2275	endif
	2276	else
	2277	write(,)'ib must be 1,2,3 or 4!!'
	2278	write(,)'stop at mztvc/387'
	2279	endif
	2280
	2281
	2282	c******
	2283	c****** calculation of tau(1,ir) for 1<=r
	2284	c******
	2285	call initial
	2286
	2287	ff=1.0e10
	2288
	2289	in=1
	2290
	2291	tau(in,1) = 1.d0
	2292
	2293	call initial
	2294	call intz (zl(in), c1,p1,mr1,t1, con)
	2295	do kr=1,nbox
	2296	ta(kr) = t1
	2297	end do
	2298	call interstrength (st1,t1,ka,ta)
	2299	do kr=1,nbox
	2300	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	2301	end do
	2302	c1 = c1 * st1 * dble(deltaz)
	2303
	2304	do 2 ir=2,nl
	2305
	2306	call intz (zl(ir), c2,p2,mr2,t2, con)
	2307	do kr=1,nbox
	2308	ta(kr) = t2
	2309	end do
	2310	call interstrength (st2,t2,ka,ta)
	2311	do kr=1,nbox
	2312	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	2313	end do
	2314	c2 = c2 * st2 * dble(deltaz)
	2315
	2316	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	2317	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	2318	cc = cc + ( c1 + c2 ) / 2.d0
	2319	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	2320	do kr=1,nbox
	2321	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
	2322	ddbox(kr) = ddbox(kr) +
	2323	$ ( t1c1box(kr) + t2c2box(kr) ) /2.d0
	2324	end do
	2325
	2326	mr1=mr2
	2327	t1=t2
	2328	c1=c2
	2329	p1=p2
	2330	do kr=1,nbox
	2331	c1box(kr) = c2box(kr)
	2332	end do
	2333
	2334	pt = bb / cc
	2335	pp = aa / (cc * ff)
	2336
	2337	ts = dd/cc
	2338	do kr=1,nbox
	2339	ta(kr) = ddbox(kr) / ccbox(kr)
	2340	end do
	2341	call interstrength(st,ts,ka,ta)
	2342	call intershape(alsa,alna,alda,ta)
	2343
	2344
	2345	eqwmu = 0.0d0
	2346	do im = 1,iimu
	2347	eqw=0.0d0
	2348	do kr=1,nbox
	2349	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	2350	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
	2351	if ( i_supersat .eq. 0 ) then
	2352	eqw=eqw+no(kr)*w
	2353	else
	2354	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	2355	endif
	2356	end do
	2357	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	2358	end do
	2359
	2360	tau(in,ir) = exp( - eqwmu / dble(deltanu(isot,ib)) )
	2361
	2362	2 continue
	2363
	2364
	2365
	2366	c
	2367	c due to the simmetry of the transmittances
	2368	c
	2369	do in=nl,2,-1
	2370	tau(in,1) = tau(1,in)
	2371	end do
	2372
	2373	vc(1) = 0.0d0
	2374	vc(nl) = 0.0d0
	2375	do in=2,nl-1 ! poner aqui nl-1 luego
	2376	vc(in) = pideltanu(isot,ib)/( 2.d0deltaz1.d5 )
	2377	@ ( tau(in-1,1) - tau(in+1,1) )
	2378	end do
	2379
	2380
	2381	c end
	2382	return
	2383	end
	2384
	2385
	2386
	2387
	2388
	2389	c***********************************************************************
	2390	c mztvc_626fh.F
	2391	c***********************************************************************
	2392
	2393	subroutine mztvc_626fh(ig)
	2394
	2395	c jul 2011 malv+fgg
	2396	c***********************************************************************
	2397
	2398	implicit none
	2399
	2400	!!!!!!!!!!!!!!!!!!!!!!!
	2401	! common variables & constants
	2402
	2403	include 'nlte_paramdef.h'
	2404	include 'nlte_commons.h'
	2405
	2406	!!!!!!!!!!!!!!!!!!!!!!!
	2407	! arguments
	2408
	2409	integer ig ! ADDED FOR TRACEBACK
	2410
	2411	!!!!!!!!!!!!!!!!!!!!!!!
	2412	! local variables
	2413
	2414	real*4 cdummy(nl,nl), csngl(nl,nl)
	2415
	2416	real*8 cax1(nl,nl), cax2(nl,nl), cax3(nl,nl)
	2417	real*8 v1(nl), v2(nl), v3(nl), cm_factor, vc_factor
	2418
	2419	integer itauout,icfout,itableout, interpol,ismooth, isngldble
	2420	integer i,j,ik,ist,isot,ib,itt
	2421
	2422	!character bandcode*2
	2423	character isotcode*2
	2424	!character codmatrx_hot*5
	2425
	2426	!!!!!!!!!!!!!!!!!!!!!!!
	2427	! external functions
	2428
	2429	external bandid
	2430	character*2 bandid
	2431
	2432	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2433	! subroutines called:
	2434	! mz4sub, dmzout, readc_mz4, readcupdw, mztf
	2435
	2436	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2437	! formatos
	2438	132 format(i2)
	2439
	2440	************************************************************************
	2441	************************************************************************
	2442
	2443	isngldble = 1 ! =1 --> dble precission
	2444
	2445	fileroot = 'cfl'
	2446
	2447	ist = 1
	2448	isot = 26
	2449	write (isotcode,132) isot
	2450
	2451	call zerov( vc121, nl )
	2452
	2453	do 11, ik=1,3
	2454
	2455	ib=ik+1
	2456
	2457	call mztvc (ig,v1, ib, 1, irw_mztf, imu, 0,0,0 )
	2458
	2459	do i=1,nl
	2460
	2461	if(ik.eq.1)then
	2462	vc_factor = dble(667.75/618.03)
	2463	elseif(ik.eq.2)then
	2464	vc_factor = 1.d0
	2465	elseif(ik.eq.3)then
	2466	vc_factor = dble(667.75/720.806)
	2467	end if
	2468
	2469	vc121(i) = vc121(i) + v1(i) * vc_factor
	2470
	2471	end do
	2472
	2473	11 continue
	2474
	2475
	2476	return
	2477	end
	2478
	2479
	2480
	2481
	2482
	2483	c***********************************************************************
	2484	c mztf_correccion
	2485	c***********************************************************************
	2486
	2487	subroutine mztf_correccion (coninf, con, ib, isot, icurt_pop)
	2488
	2489	c including the dependence of the absort. coeff. on temp., vibr. temp.,
	2490	c function, etc.., when neccessary. imr is already corrected in his.for
	2491	c we follow pg.39b-43a (l5):
	2492	c tvt1 is the vibr temp of the upper level
	2493	c tvt is the vibr temp of the transition itself
	2494	c tvtbs is the vibr temp of the bending mode (used in qv)
	2495	c for fundamental bands, they are not used at the moment.
	2496	c for the 15 fh and sh bands, only tvt0 is used at the moment.
	2497	c for the laser band, all of them are used following pg. 41a -l5- :
	2498	c we need s(z) and we can read s(tk) from the histogram (also called
	2499	c what we have to calculate now is the factor s(z)/s(tk) or following
	2500	c l5 notebook notation, s_nlte/s_lte.
	2501	c s_nlte/s_lte = xfactor = xlower * xqv * xes
	2502
	2503	c icurt_pop = 30 -> Output of populations of the 0200,0220,1000 states
	2504	c = otro -> no output of these populations
	2505
	2506	c oct 92 malv
	2507	c jan 98 malv version for mz1d
	2508	c jul 2011 malv+fgg adapted to LMD-MGCM
	2509	c***********************************************************************
	2510
	2511	implicit none
	2512
	2513	include 'nlte_paramdef.h'
	2514	include 'nlte_commons.h'
	2515
	2516	c arguments
	2517	integer ib, isot
	2518	integer icurt_pop ! output of Fermi states population
	2519	real*8 con(nzy), coninf
	2520
	2521	! local variables
	2522	integer i
	2523	real*8 tvt0(nzy),tvt1(nzy),tvtbs(nzy), zld(nl),zyd(nzy)
	2524	real xalfa, xbeta, xtv1000, xtv0200, xtv0220, xfactor
	2525	real xqv, xnu_trans, xtv_trans, xes, xlower
	2526	c***********************************************************************
	2527
	2528	xfactor = 1.0
	2529
	2530	do i=1,nzy
	2531	zyd(i) = dble(zy(i))
	2532	enddo
	2533	do i=1,nl
	2534	zld(i) = dble( zl(i) )
	2535	end do
	2536
	2537	! tvtbs is the bending mode of the molecule. used in xqv.
[695]	2538	if (isot.eq.1) call interdp (tvtbs,zyd,nzy,v626t1,zld,nl,1)
	2539	if (isot.eq.2) call interdp (tvtbs,zyd,nzy,v628t1,zld,nl,1)
	2540	if (isot.eq.3) call interdp (tvtbs,zyd,nzy,v636t1,zld,nl,1)
	2541	if (isot.eq.4) call interdp (tvtbs,zyd,nzy,v627t1,zld,nl,1)
	2542	if (isot.eq.5) call interdp (tvtbs,zyd,nzy,vcot1,zld,nl,1)
[498]	2543
	2544	! tvt0 is the lower level of the transition. used in xlower.
	2545	if (ib.eq.2 .or. ib.eq.3 .or. ib.eq.4 .or. ib.eq.15) then
[695]	2546	if (isot.eq.1) call interdp(tvt0,zyd,nzy,v626t1,zld,nl,1)
	2547	if (isot.eq.2) call interdp(tvt0,zyd,nzy,v628t1,zld,nl,1)
	2548	if (isot.eq.3) call interdp(tvt0,zyd,nzy,v636t1,zld,nl,1)
	2549	if (isot.eq.4) call interdp(tvt0,zyd,nzy,v627t1,zld,nl,1)
[498]	2550	elseif (ib.eq.6 .or. ib.eq.8 .or. ib.eq.10
	2551	@ .or. ib.eq.13 .or. ib.eq.14
	2552	@ .or. ib.eq.17 .or. ib.eq.19 .or. ib.eq.20) then
[695]	2553	if (isot.eq.1) call interdp(tvt0,zyd,nzy,v626t2,zld,nl,1)
	2554	if (isot.eq.2) call interdp(tvt0,zyd,nzy,v628t2,zld,nl,1)
	2555	if (isot.eq.3) call interdp(tvt0,zyd,nzy,v636t2,zld,nl,1)
[498]	2556	if (isot.eq.4) then
	2557	call interdp ( tvt0,zyd,nzy, v627t2,zld,nl, 1 )
	2558	endif
	2559	else
	2560	do i=1,nzy
	2561	tvt0(i) = dble( ty(i) )
	2562	end do
	2563	end if
	2564
	2565	c tvt is the vt of the transition. used in xes.
	2566	c since xes=1.0 except for the laser bands, tvt is only needed for them
	2567	c but it is actually calculated from the tv of the upper and lower level
	2568	c of the transition. hence, only tvt1 remains to be read for the laser b
	2569	c tvt1 is the upper level of the transition.
	2570	if (ib.eq.13 .or. ib.eq.14) then
[695]	2571	if (isot.eq.1) call interdp(tvt1,zyd,nzy,v626t4,zld,nl,1)
	2572	if (isot.eq.2) call interdp(tvt1,zyd,nzy,v628t4,zld,nl,1)
	2573	if (isot.eq.3) call interdp(tvt1,zyd,nzy,v636t4,zld,nl,1)
	2574	if (isot.eq.4) call interdp(tvt1,zyd,nzy,v627t4,zld,nl,1)
[498]	2575	end if
	2576
	2577	c here we weight the absorber amount by a factor which compensate the l
	2578	c value of the strength read from hitran. we use that factor in order t
	2579	c correct the product s*m when we later multiply those two variables.
	2580
	2581	! if ( isot.eq.1 .and. icurt_pop.eq.30 ) then
	2582	! open (30, file='020populations.dat')
	2583	! write (30,*) ' z tv(020) tv0200 tv0220 tv1000 '
	2584	! endif
	2585
	2586	do i=1,nzy
	2587
	2588	if (isot.eq.1) then
	2589
	2590	!!! vt of the 3 levels in (020) (see pag. 36a-sn1 for this)
[695]	2591	xalfa = 1.d0/2.d0exp(dble(-ee(nu12_1000-nu(1,2))/ty(i)))
	2592	xbeta = 1.d0/2.d0exp(dble(-ee(nu12_0200-nu(1,2))/ty(i)))
[498]	2593	xtv0200 = dble( - ee * nu12_0200 ) /
	2594	@ ( log( xbeta/(1.d0+xalfa+xbeta) ) -
	2595	@ dble(ee*nu(1,2))/tvt0(i) )
	2596	xtv0220 = dble( - ee * nu(1,2) ) /
	2597	@ ( log( 1.d0/(1.d0+xalfa+xbeta) ) -
	2598	@ dble(ee*nu(1,2))/tvt0(i) )
	2599	xtv1000 = dble( - ee * nu12_1000 ) /
	2600	@ ( log( xalfa/(1.d0+xalfa+xbeta) ) -
	2601	@ dble(ee*nu(1,2))/tvt0(i) )
	2602	!!! correccion 8-Nov-04 (see pag.9b-Marte4-)
	2603	xtv0200 = dble( - ee * nu12_0200 /
	2604	@ (log(4.xbeta/(1.d0+xalfa+xbeta))-eenu(1,2)/tvt0(i)))
	2605	xtv0220 = dble( - ee * nu(1,2) /
	2606	@ ( log(2./(1.d0+xalfa+xbeta)) - ee*nu(1,2)/tvt0(i) ) )
	2607	xtv1000 = dble( - ee * nu12_1000 /
[695]	2608	@ (log(4.xalfa/(1.d0+xalfa+xbeta))-eenu(1,2)/tvt0(i)))
[498]	2609
	2610	! if ( icurt_pop.eq.30 ) then
	2611	! write (30,'( 1x,f7.2, 3x,f8.3, 2x,3(1x,f8.3) )')
	2612	! @ zx(i), tvt0(i), xtv0200, xtv0220, xtv1000
	2613	! endif
	2614
	2615	!!! xlower and xes for the band
	2616	if (ib.eq.19) then
	2617	xlower = exp( dble(eeelow(isot,ib))
	2618	@ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )
	2619	xes = 1.0d0
	2620	elseif (ib.eq.17) then
	2621	xlower = exp( dble(eeelow(isot,ib))
	2622	@ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )
	2623	xes = 1.0d0
	2624	elseif (ib.eq.20) then
	2625	xlower = exp( dble(eeelow(isot,ib))
	2626	@ ( 1.d0/dble(ty(i))-1.d0/xtv0220 ) )
	2627	xes = 1.0d0
	2628	elseif (ib.eq.14) then
	2629	xlower = exp( dble(eenu12_1000)
	2630	@ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )
	2631	xnu_trans = dble( nu(1,4)-nu12_1000 )
	2632	xtv_trans = xnu_trans / dble(nu(1,4)/tvt1(i)-
	2633	@ nu12_1000/xtv1000)
	2634	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
	2635	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
	2636	elseif (ib.eq.13) then
	2637	xlower = exp( dble(eenu12_0200)
	2638	@ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )
	2639	xnu_trans = dble(nu(1,4)-nu12_0200)
	2640	xtv_trans = xnu_trans / dble(nu(1,4)/tvt1(i)-
	2641	@ nu12_0200/xtv0200)
	2642	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
	2643	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
	2644	else
	2645	xlower = exp( dble(eeelow(isot,ib))
	2646	@ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )
	2647	xes = 1.0d0
	2648	end if
	2649	xqv = (1.d0-exp( dble(-ee*667.3801/tvtbs(i)) )) /
	2650	@ (1.d0-exp( dble(-ee*667.3801/ty(i)) ))
	2651	xfactor = xlower * xqv*2.d0 xes
	2652
	2653	elseif (isot.eq.2) then
	2654
	2655	xalfa = 1.d0/2.d0* exp( dble(-ee*(nu22_1000-nu(2,2))/
	2656	@ ty(i)) )
	2657	xbeta = 1.d0/2.d0* exp( dble(-ee*(nu22_0200-nu(2,2))/
	2658	@ ty(i)) )
	2659	xtv0200 = dble( - ee * nu22_0200 ) /
	2660	@ ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(2,2))/
	2661	@ tvt0(i) )
	2662	xtv1000 = dble( - ee * nu22_1000 ) /
	2663	@ ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(2,2))/
	2664	@ tvt0(i) )
	2665
	2666	if (ib.eq.14) then
	2667	xlower = exp( dble(eenu22_1000)
	2668	@ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )
	2669	xnu_trans = dble(nu(2,4)-nu22_1000)
	2670	xtv_trans = xnu_trans / dble(nu(2,4)/tvt1(i)-nu22_1000/
	2671	@ xtv1000)
	2672	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
	2673	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
	2674	elseif (ib.eq.13) then
	2675	xlower = exp( dble(eenu22_0200)
	2676	@ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )
	2677	xnu_trans = dble( nu(2,4)-nu22_0200 )
	2678	xtv_trans = xnu_trans / dble(nu(2,4)/tvt1(i)-nu22_0200/
	2679	@ xtv0200)
	2680	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
	2681	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
	2682	else
	2683	xlower = exp( dble(eeelow(isot,ib))
	2684	@ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )
	2685	xes = 1.0d0
	2686	end if
	2687	xqv = (1.d0-exp( dble(-ee*662.3734/tvtbs(i)) )) /
	2688	@ (1.d0-exp( dble(-ee*662.3734/ty(i)) ))
	2689	xfactor = xlower * xqv*2.d0 xes
	2690
	2691	elseif (isot.eq.3) then
	2692
	2693	xalfa = 1.d0/2.d0* exp( dble(-ee*(nu32_1000-nu(3,2))/
	2694	@ ty(i)) )
	2695	xbeta = 1.d0/2.d0* exp( dble(-ee*(nu32_0200-nu(3,2))/
	2696	@ ty(i)) )
	2697	xtv0200 = dble( - ee * nu32_0200 ) /
	2698	@ ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(3,2))/
	2699	@ tvt0(i) )
	2700	xtv1000 = dble( - ee * nu32_1000 ) /
	2701	@ ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(3,2))/
	2702	@ tvt0(i) )
	2703
	2704	if (ib.eq.14) then
	2705	xlower = exp( dble(eenu32_1000)
	2706	@ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )
	2707	xnu_trans = dble( nu(3,4)-nu32_1000 )
	2708	xtv_trans = xnu_trans / dble(nu(3,4)/tvt1(i)-nu32_1000/
	2709	@ xtv1000)
	2710	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
	2711	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
	2712	elseif (ib.eq.13) then
	2713	xlower = exp( dble(eenu32_0200)
	2714	@ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )
	2715	xnu_trans = dble( nu(3,4)-nu32_0200 )
	2716	xtv_trans = xnu_trans / dble(nu(3,4)/tvt1(i)-nu32_0200/
	2717	@ xtv0200)
	2718	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
	2719	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
	2720	else
	2721	xlower = exp( dble(eeelow(isot,ib))
	2722	@ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )
	2723	xes = 1.0d0
	2724	end if
	2725	xqv = (1.d0-exp( dble(-ee*648.4784/tvtbs(i)) )) /
	2726	@ (1.d0-exp( dble(-ee*648.4784/ty(i)) ))
	2727	xfactor = xlower * xqv*2.d0 xes
	2728
	2729	elseif (isot.eq.4) then
	2730
	2731	xalfa = 1.d0/2.d0* exp( dble(-ee*(nu42_1000-nu(4,2))/
	2732	@ ty(i)) )
	2733	xbeta = 1.d0/2.d0* exp( dble(-ee*(nu42_0200-nu(4,2))/
	2734	@ ty(i)) )
	2735	xtv0200 = dble( - ee * nu42_0200 ) /
	2736	@ ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(4,2))/
	2737	@ tvt0(i) )
	2738	xtv1000 = dble( - ee * nu42_1000 ) /
	2739	@ ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(4,2))/
	2740	@ tvt0(i) )
	2741
	2742	if (ib.eq.14) then
	2743	xlower = exp( dble(eenu42_1000)
	2744	@ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )
	2745	xnu_trans = dble( nu(4,4)-nu42_1000 )
	2746	xtv_trans = xnu_trans / dble(nu(4,4)/tvt1(i)-nu42_1000/
	2747	@ xtv1000)
	2748	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
	2749	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
	2750	elseif (ib.eq.13) then
	2751	xlower = exp( dble(eenu42_0200)
	2752	$ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )
	2753	xnu_trans = dble( nu(4,4)-nu42_0200 )
	2754	xtv_trans = xnu_trans / dble(nu(4,4)/tvt1(i)-nu42_0200/
	2755	@ xtv0200)
	2756	xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /
	2757	@ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))
	2758	else
	2759	xlower = exp( dble(eeelow(isot,ib))
	2760	@ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )
	2761	xes = 1.0d0
	2762	end if
	2763	xqv = (1.d0-exp( dble(-ee*664.7289/tvtbs(i)) )) /
	2764	@ (1.d0-exp( dble(-ee*664.7289/ty(i)) ))
	2765	xfactor = xlower * xqv*2.d0 xes
	2766
	2767	elseif (isot.eq.5 .and. ib.eq.1) then
	2768
	2769	xlower = 1.d0
	2770	xes = 1.0d0
	2771	xqv = (1.d0-exp( dble(-ee*nuco_10/tvtbs(i)) )) /
	2772	@ (1.d0-exp( dble(-ee*nuco_10/ty(i)) ))
	2773	xfactor = xlower * xqv * xes
	2774
	2775	end if
	2776
	2777	con(i) = con(i) * xfactor
	2778	if (i.eq.nzy) coninf = coninf * xfactor
	2779
	2780	end do
	2781
	2782	! if ( isot.eq.1 .and. icurt_pop.eq.30 ) then
	2783	! close (30)
	2784	! endif
	2785
	2786	return
	2787	end
	2788
	2789
	2790
	2791
	2792
	2793	c***********************************************************************
	2794	c mztf.f
	2795	c***********************************************************************
	2796	c
	2797	c program for calculating atmospheric transmittances
	2798	c to be used in the calculation of curtis matrix coefficients
	2799
	2800	subroutine mztf ( ig,cf,cfup,cfdw,vc,taugr, ib,isot,
	2801	@ iirw,iimu,itauout,icfout,itableout )
	2802
	2803	c i*out = 1 output of data
	2804	c i*out = 0 no output
	2805	c
	2806	c jul 2011 malv+fgg adapted to LMD-MGCM
	2807	c nov 98 mavl allow for overlaping in the lorentz line
	2808	c jan 98 malv version for mz1d. based on curtis/mztf.for
	2809	c 17-jul-96 mlp&crs change the calculation of mr.
	2810	c evitar: divide por cero. anhadiendo: ff
	2811	c oct-92 malv correct s(t) dependence for all histogr bands
	2812	c june-92 malv proper lower levels for laser bands
	2813	c may-92 malv new temperature dependence for laser bands
	2814	c @ 991 malv boxing for the averaged absorber amount and t
	2815	c ? malv extension up to 200 km altitude in mars
	2816	c 13-nov-86 mlp include the temperature weighted to match
	2817	c the eqw in the strong doppler limit.
	2818	c***********************************************************************
	2819
	2820	implicit none
	2821
	2822	include 'nlte_paramdef.h'
	2823	include 'nlte_commons.h'
	2824
	2825
	2826	c arguments
	2827	integer ig !ADDED FOR TRACEBACK
	2828	real*8 cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o.
	2829	real*8 vc(nl), taugr(nl) ! o
	2830	integer ib ! i
	2831	integer isot ! i
	2832	integer iirw ! i
	2833	integer iimu ! i
	2834	integer itauout ! i
	2835	integer icfout ! i
	2836	integer itableout ! i
	2837
	2838	c local variables and constants
	2839	integer i, in, ir, im, k ,j
	2840	integer nmu
	2841	parameter (nmu = 8)
	2842	real*8 tau(nl,nl)
	2843	real*8 tauinf(nl)
	2844	real*8 con(nzy), coninf
	2845	real*8 c1, c2
	2846	real*8 t1, t2
	2847	real*8 p1, p2
	2848	real*8 mr1, mr2
	2849	real*8 st1, st2
	2850	real*8 c1box(70), c2box(70)
	2851	real*8 ff ! to avoid too small numbers
	2852	real*8 tvtbs(nzy)
	2853	real*8 st, beta, ts, eqwmu
	2854	real*8 mu(nmu), amu(nmu)
	2855	real*8 zld(nl), zyd(nzy)
	2856	real*8 correc
	2857	real deltanux ! width of vib-rot band (cm-1)
	2858	! character isotcode*2
	2859	integer idummy
	2860	real*8 Desp,wsL
	2861
	2862	c formats
	2863	! 111 format(a1)
	2864	! 112 format(a2)
	2865	101 format(i1)
	2866	202 format(i2)
	2867	! 180 format(a80)
	2868	! 181 format(a80)
	2869	c***********************************************************************
	2870
	2871	c some needed values
	2872	! rl=sqrt(log(2.d0))
	2873	! pi2 = 3.14159265358989d0
	2874	beta = 1.8d0
	2875	idummy = 0
	2876	Desp = 0.d0
	2877	wsL = 0.d0
	2878
	2879	c esto es para que las subroutines de mztfsub calculen we
	2880	c de la forma apropiada para mztf, no para fot
	2881	icls=icls_mztf
	2882
	2883	c codigos para filenames
	2884	! if (isot .eq. 1) isotcode = '26'
	2885	! if (isot .eq. 2) isotcode = '28'
	2886	! if (isot .eq. 3) isotcode = '36'
	2887	! if (isot .eq. 4) isotcode = '27'
	2888	! if (isot .eq. 5) isotcode = '62'
	2889	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	2890	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	2891	! write (ibcode1,101) ib
	2892	! else
	2893	! write (ibcode2,202) ib
	2894	! endif
	2895	! write (*,'( 30h calculating curtis matrix : ,2x,
	2896	! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot
	2897
	2898	c integration in angle !!!!!!!!!!!!!!!!!!!!
	2899
	2900	c------- diffusivity approx.
	2901	if (iimu.eq.1) then
	2902	! write (,) ' diffusivity approx. beta = ',beta
	2903	mu(1) = 1.0d0
	2904	amu(1)= 1.0d0
	2905	c-------data for 8 points integration
	2906	elseif (iimu.eq.4) then
	2907	write (,)' 4 points for the gauss-legendre angle quadrature.'
	2908	mu(1)=(1.0d0+0.339981043584856)/2.0d0
	2909	mu(2)=(1.0d0-0.339981043584856)/2.0d0
	2910	mu(3)=(1.0d0+0.861136311594053)/2.0d0
	2911	mu(4)=(1.0d0-0.861136311594053)/2.0d0
	2912	amu(1)=0.652145154862546
	2913	amu(2)=amu(1)
	2914	amu(3)=0.347854845137454
	2915	amu(4)=amu(3)
	2916	beta=1.0d0
	2917	c-------data for 8 points integration
	2918	elseif(iimu.eq.8) then
	2919	write (,)' 8 points for the gauss-legendre angle quadrature.'
	2920	mu(1)=(1.0d0+0.183434642495650)/2.0d0
	2921	mu(2)=(1.0d0-0.183434642495650)/2.0d0
	2922	mu(3)=(1.0d0+0.525532409916329)/2.0d0
	2923	mu(4)=(1.0d0-0.525532409916329)/2.0d0
	2924	mu(5)=(1.0d0+0.796666477413627)/2.0d0
	2925	mu(6)=(1.0d0-0.796666477413627)/2.0d0
	2926	mu(7)=(1.0d0+0.960289856497536)/2.0d0
	2927	mu(8)=(1.0d0-0.960289856497536)/2.0d0
	2928	amu(1)=0.362683783378362
	2929	amu(2)=amu(1)
	2930	amu(3)=0.313706645877887
	2931	amu(4)=amu(3)
	2932	amu(5)=0.222381034453374
	2933	amu(6)=amu(5)
	2934	amu(7)=0.101228536290376
	2935	amu(8)=amu(7)
	2936	beta=1.0d0
	2937	end if
	2938	c!!!!!!!!!!!!!!!!!!!!!!!
	2939
	2940	ccc
	2941	ccc determine abundances included in the absorber amount
	2942	ccc
	2943
	2944	c first, set up the grid ready for interpolation.
	2945	do i=1,nzy
	2946	zyd(i) = dble(zy(i))
	2947	enddo
	2948	do i=1,nl
	2949	zld(i) = dble(zl(i))
	2950	enddo
	2951
	2952
	2953	c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
	2954	c por similitud a la que se hace en cza.for
	2955
	2956	do i=1,nzy
	2957	if (isot.eq.5) then
	2958	con(i) = dble( coy(i) * imrco )
	2959	else
	2960	con(i) = dble( co2y(i) * imr(isot) )
	2961	c vibr. temp of the bending mode :
	2962	if(isot.eq.1) call interdp(tvtbs,zyd,nzy,v626t1,zld,nl,1)
	2963	if(isot.eq.2) call interdp(tvtbs,zyd,nzy,v628t1,zld,nl,1)
	2964	if(isot.eq.3) call interdp(tvtbs,zyd,nzy,v636t1,zld,nl,1)
	2965	if(isot.eq.4) call interdp(tvtbs,zyd,nzy,v627t1,zld,nl,1)
	2966	correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )
	2967	con(i) = con(i) * ( 1.d0 - correc )
	2968	endif
	2969	c-----------------------------------------------------------------------
	2970	c mlp & cristina. 17 july 1996
	2971	c change the calculation of mr. it is used for calculating partial press
	2972	c alpha = alpha(self,co2)pp +alpha(n2)(pt-pp)
	2973	c for an isotope, if mr is obtained by co2*imr(iso)/nt we are considerin
	2974	c collisions with other co2 isotopes (including the major one, 626)
	2975	c as if they were with n2. assuming mr as co2/nt, we consider collisions
	2976	c of type 628-626 as of 626-626 instead of as 626-n2.
	2977	c mrx(i)=con(i)/ntx(i) ! old malv
	2978
	2979	! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs
	2980
	2981	c jan 98:
	2982	c esta modif de mlp implica anular el correc (deberia revisar esto)
	2983	mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98
	2984
	2985	c-----------------------------------------------------------------------
	2986
	2987	end do
	2988
	2989	! como beta y 1.d5 son comunes a todas las weighted absorber amounts,
	2990	! los simplificamos:
	2991	! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )
	2992	coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
	2993
	2994	! write (,) ' coninf =', coninf
	2995
	2996	ccc
	2997	ccc temp dependence of the band strength and
	2998	ccc nlte correction factor for the absorber amount
	2999	ccc
	3000	call mztf_correccion ( coninf, con, ib, isot, itableout )
	3001
	3002	ccc
	3003	ccc reads histogrammed spectral data (strength for lte and vmr=1)
	3004	ccc
	3005	!hfile1 = dirspec//'hi'//dn ! ya no distinguimos entre d/n
	3006	!! hfile1 = dirspec//'hid' ! (see why in his.for)
	3007	! hfile='hid'
	3008	!! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
	3009	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
	3010	!
	3011	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	3012	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	3013	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat'
	3014	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat'
	3015	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat'
	3016	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat'
	3017	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
	3018	! else
	3019	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat'
	3020	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat'
	3021	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat'
	3022	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat'
	3023	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat'
	3024	! endif
	3025	! write (,) 'hisfile: ', hisfile
	3026
	3027	! the argument to rhist is to make this compatible with mztf_comp.f,
	3028	! which is a useful modification of mztf.f (to change strengths of bands
	3029	! call rhist (1.0)
	3030	if(ib.eq.1) then
	3031	if(isot.eq.1) then !Case 1
	3032	mm=mm_c1
	3033	nbox=nbox_c1
	3034	tmin=tmin_c1
	3035	tmax=tmax_c1
	3036	do i=1,nbox_max
	3037	no(i)=no_c1(i)
	3038	dist(i)=dist_c1(i)
	3039	do j=1,nhist
	3040	sk1(j,i)=sk1_c1(j,i)
	3041	xls1(j,i)=xls1_c1(j,i)
	3042	xln1(j,i)=xln1_c1(j,i)
	3043	xld1(j,i)=xld1_c1(j,i)
	3044	enddo
	3045	enddo
	3046	do j=1,nhist
	3047	thist(j)=thist_c1(j)
	3048	enddo
	3049	else if(isot.eq.2) then !Case 2
	3050	mm=mm_c2
	3051	nbox=nbox_c2
	3052	tmin=tmin_c2
	3053	tmax=tmax_c2
	3054	do i=1,nbox_max
	3055	no(i)=no_c2(i)
	3056	dist(i)=dist_c2(i)
	3057	do j=1,nhist
	3058	sk1(j,i)=sk1_c2(j,i)
	3059	xls1(j,i)=xls1_c2(j,i)
	3060	xln1(j,i)=xln1_c2(j,i)
	3061	xld1(j,i)=xld1_c2(j,i)
	3062	enddo
	3063	enddo
	3064	do j=1,nhist
	3065	thist(j)=thist_c2(j)
	3066	enddo
	3067	else if(isot.eq.3) then !Case 3
	3068	mm=mm_c3
	3069	nbox=nbox_c3
	3070	tmin=tmin_c3
	3071	tmax=tmax_c3
	3072	do i=1,nbox_max
	3073	no(i)=no_c3(i)
	3074	dist(i)=dist_c3(i)
	3075	do j=1,nhist
	3076	sk1(j,i)=sk1_c3(j,i)
	3077	xls1(j,i)=xls1_c3(j,i)
	3078	xln1(j,i)=xln1_c3(j,i)
	3079	xld1(j,i)=xld1_c3(j,i)
	3080	enddo
	3081	enddo
	3082	do j=1,nhist
	3083	thist(j)=thist_c3(j)
	3084	enddo
	3085	else if(isot.eq.4) then !Case 4
	3086	mm=mm_c4
	3087	nbox=nbox_c4
	3088	tmin=tmin_c4
	3089	tmax=tmax_c4
	3090	do i=1,nbox_max
	3091	no(i)=no_c4(i)
	3092	dist(i)=dist_c4(i)
	3093	do j=1,nhist
	3094	sk1(j,i)=sk1_c4(j,i)
	3095	xls1(j,i)=xls1_c4(j,i)
	3096	xln1(j,i)=xln1_c4(j,i)
	3097	xld1(j,i)=xld1_c4(j,i)
	3098	enddo
	3099	enddo
	3100	do j=1,nhist
	3101	thist(j)=thist_c4(j)
	3102	enddo
	3103	else
	3104	write(,)'isot must be 2,3 or 4 for ib=1!!'
	3105	write(,)'stop at mztf_overlap/317'
	3106	stop
	3107	endif
	3108	else if (ib.eq.2) then
	3109	if(isot.eq.1) then !Case 5
	3110	mm=mm_c5
	3111	nbox=nbox_c5
	3112	tmin=tmin_c5
	3113	tmax=tmax_c5
	3114	do i=1,nbox_max
	3115	no(i)=no_c5(i)
	3116	dist(i)=dist_c5(i)
	3117	do j=1,nhist
	3118	sk1(j,i)=sk1_c5(j,i)
	3119	xls1(j,i)=xls1_c5(j,i)
	3120	xln1(j,i)=xln1_c5(j,i)
	3121	xld1(j,i)=xld1_c5(j,i)
	3122	enddo
	3123	enddo
	3124	do j=1,nhist
	3125	thist(j)=thist_c5(j)
	3126	enddo
	3127	else
	3128	write(,)'isot must be 1 for ib=2!!'
	3129	write(,)'stop at mztf_overlap/341'
	3130	stop
	3131	endif
	3132	else if (ib.eq.3) then
	3133	if(isot.eq.1) then !Case 6
	3134	mm=mm_c6
	3135	nbox=nbox_c6
	3136	tmin=tmin_c6
	3137	tmax=tmax_c6
	3138	do i=1,nbox_max
	3139	no(i)=no_c6(i)
	3140	dist(i)=dist_c6(i)
	3141	do j=1,nhist
	3142	sk1(j,i)=sk1_c6(j,i)
	3143	xls1(j,i)=xls1_c6(j,i)
	3144	xln1(j,i)=xln1_c6(j,i)
	3145	xld1(j,i)=xld1_c6(j,i)
	3146	enddo
	3147	enddo
	3148	do j=1,nhist
	3149	thist(j)=thist_c6(j)
	3150	enddo
	3151	else
	3152	write(,)'isot must be 1 for ib=3!!'
	3153	write(,)'stop at mztf_overlap/365'
	3154	stop
	3155	endif
	3156	else if (ib.eq.4) then
	3157	if(isot.eq.1) then !Case 7
	3158	mm=mm_c7
	3159	nbox=nbox_c7
	3160	tmin=tmin_c7
	3161	tmax=tmax_c7
	3162	do i=1,nbox_max
	3163	no(i)=no_c7(i)
	3164	dist(i)=dist_c7(i)
	3165	do j=1,nhist
	3166	sk1(j,i)=sk1_c7(j,i)
	3167	xls1(j,i)=xls1_c7(j,i)
	3168	xln1(j,i)=xln1_c7(j,i)
	3169	xld1(j,i)=xld1_c7(j,i)
	3170	enddo
	3171	enddo
	3172	do j=1,nhist
	3173	thist(j)=thist_c7(j)
	3174	enddo
	3175	else
	3176	write(,)'isot must be 1 for ib=4!!'
	3177	write(,)'stop at mztf_overlap/389'
	3178	stop
	3179	endif
	3180	else
	3181	write(,)'ib must be 1,2,3 or 4!!'
	3182	write(,)'stop at mztf_overlap/394'
	3183	endif
	3184
	3185	if (isot.ne.5) deltanux = deltanu(isot,ib)
	3186	if (isot.eq.5) deltanux = deltanuco
	3187
	3188	c******
	3189	c****** calculation of tauinf(nl)
	3190	c******
	3191	call initial
	3192
	3193	ff=1.0e10
	3194
	3195	do i=nl,1,-1
	3196
	3197	if(i.eq.nl)then
	3198
	3199	call intz (zl(i),c2,p2,mr2,t2, con)
	3200	do kr=1,nbox
	3201	ta(kr)=t2
	3202	end do
	3203	! write (,) ' i, t2 =', i, t2
	3204	call interstrength (st2,t2,ka,ta)
	3205	aa = p2 * coninf * mr2 * (st2 * ff)
	3206	bb = p2 * coninf * st2
	3207	cc = coninf * st2
	3208	dd = t2 * coninf * st2
	3209	do kr=1,nbox
	3210	ccbox(kr) = coninf * ka(kr)
	3211	ddbox(kr) = t2 * ccbox(kr)
	3212	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
	3213	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	3214	end do
	3215	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
	3216	c2 = c2 * st2 * dble(deltaz)
	3217
	3218	else
	3219	call intz (zl(i),c1,p1,mr1,t1, con)
	3220	do kr=1,nbox
	3221	ta(kr)=t1
	3222	end do
	3223	! write (,) ' i, t1 =', i, t1
	3224	call interstrength (st1,t1,ka,ta)
	3225	do kr=1,nbox
	3226	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
	3227	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	3228	end do
	3229	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
	3230	c1 = c1 * st1 * dble(deltaz)
	3231	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	3232	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	3233	cc = cc + ( c1 + c2 ) / 2.d0
	3234	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	3235	do kr=1,nbox
	3236	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) )/2.d0
	3237	ddbox(kr) = ddbox(kr) + (t1c1box(kr)+t2c2box(kr))/2.d0
	3238	end do
	3239
	3240	mr2 = mr1
	3241	c2=c1
	3242	do kr=1,nbox
	3243	c2box(kr) = c1box(kr)
	3244	end do
	3245	t2=t1
	3246	p2=p1
	3247	end if
	3248
	3249	pt = bb / cc
	3250	pp = aa / (cc*ff)
	3251
	3252	! ta=dd/cc
	3253	! tdop = ta
	3254	ts = dd/cc
	3255	do kr=1,nbox
	3256	ta(kr) = ddbox(kr) / ccbox(kr)
	3257	end do
	3258	! write (,) ' i, ts =', i, ts
	3259	call interstrength(st,ts,ka,ta)
	3260	! call intershape(alsa,alna,alda,tdop)
	3261	call intershape(alsa,alna,alda,ta)
	3262
	3263	* ua = cc/st
	3264
	3265	c next loop calculates the eqw for an especified path ua,pp,pt,ta
	3266
	3267	eqwmu = 0.0d0
	3268	do im = 1,iimu
	3269	eqw=0.0d0
	3270	do kr=1,nbox
	3271	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	3272	if(ua(kr).lt.0.)write(,)'mztf_overlap/483',ua(kr),
	3273	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
	3274
	3275	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
	3276	if ( i_supersat .eq. 0 ) then
	3277	eqw=eqw+no(kr)*w
	3278	else
	3279	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	3280	endif
	3281	end do
	3282	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	3283	end do
	3284
	3285	tauinf(i) = exp( - eqwmu / dble(deltanux) )
	3286
	3287	end do ! i continue
	3288
	3289	! if ( isot.eq.1 .and. ib.eq.2 ) then
	3290	! write (,) ' tauinf(nl) = ', tauinf(nl)
	3291	! write (,) ' tauinf(1) = ', tauinf(1)
	3292	! endif
	3293
	3294	c******
	3295	c****** calculation of tau(in,ir) for n<=r
	3296	c******
	3297
	3298	do 1 in=1,nl-1
	3299
	3300	call initial
	3301	call intz (zl(in), c1,p1,mr1,t1, con)
	3302	do kr=1,nbox
	3303	ta(kr) = t1
	3304	end do
	3305	call interstrength (st1,t1,ka,ta)
	3306	do kr=1,nbox
	3307	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
	3308	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	3309	end do
	3310	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
	3311	c1 = c1 * st1 * dble(deltaz)
	3312
	3313	do 2 ir=in,nl-1
	3314
	3315	if (ir.eq.in) then
	3316	tau(in,ir) = 1.d0
	3317	goto 2
	3318	end if
	3319
	3320	call intz (zl(ir), c2,p2,mr2,t2, con)
	3321	do kr=1,nbox
	3322	ta(kr) = t2
	3323	end do
	3324	call interstrength (st2,t2,ka,ta)
	3325	do kr=1,nbox
	3326	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
	3327	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	3328	end do
	3329	! c2 = c2 * st2 * beta * dble(deltaz) * 1.e5
	3330	c2 = c2 * st2 * dble(deltaz)
	3331
	3332	c aa = aa + ( p1mr1c1 + p2mr2c2 ) / 2.d0
	3333	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	3334	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	3335	cc = cc + ( c1 + c2 ) / 2.d0
	3336	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	3337	do kr=1,nbox
	3338	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
	3339	ddbox(kr) = ddbox(kr) +
	3340	$ ( t1c1box(kr) + t2c2box(kr) ) /2.d0
	3341	end do
	3342
	3343	mr1=mr2
	3344	t1=t2
	3345	c1=c2
	3346	p1=p2
	3347	do kr=1,nbox
	3348	c1box(kr) = c2box(kr)
	3349	end do
	3350
	3351	pt = bb / cc
	3352	pp = aa / (cc * ff)
	3353
	3354	* ta=dd/cc
	3355	* tdop = ta
	3356	ts = dd/cc
	3357	do kr=1,nbox
	3358	ta(kr) = ddbox(kr) / ccbox(kr)
	3359	end do
	3360	call interstrength(st,ts,ka,ta)
	3361	call intershape(alsa,alna,alda,ta)
	3362	* ua = cc/st
	3363
	3364	c next loop calculates the eqw for an especified path ua,pp,pt,ta
	3365
	3366	eqwmu = 0.0d0
	3367	do im = 1,iimu
	3368	eqw=0.0d0
	3369	do kr=1,nbox
	3370	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	3371	if(ua(kr).lt.0.)write(,)'mztf_overlap/581',ua(kr),
	3372	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
	3373
	3374	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
	3375	if ( i_supersat .eq. 0 ) then
	3376	eqw=eqw+no(kr)*w
	3377	else
	3378	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	3379	endif
	3380	end do
	3381	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	3382	end do
	3383
	3384	tau(in,ir) = exp( - eqwmu / dble(deltanux) )
	3385
	3386	2 continue
	3387
	3388	1 continue
	3389
	3390	! if ( isot.eq.1 .and. ib.eq.2 ) then
	3391	! write (,) ' tau(1,) , =1,20 '
	3392	! write (,) ( sngl(tau(1,k)), k=1,20 )
	3393	! endif
	3394
	3395
	3396	c**********
	3397	c********** calculation of tau(in,ir) for n>r
	3398	c**********
	3399
	3400	in=nl
	3401
	3402	call initial
	3403	call intz (zl(in), c1,p1,mr1,t1, con)
	3404	do kr=1,nbox
	3405	ta(kr) = t1
	3406	end do
	3407	call interstrength (st1,t1,ka,ta)
	3408	do kr=1,nbox
	3409	! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5
	3410	c1box(kr) = c1 * ka(kr) * dble(deltaz)
	3411	end do
	3412	! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5
	3413	c1 = c1 * st1 * dble(deltaz)
	3414
	3415	do 4 ir=in-1,1,-1
	3416
	3417	call intz (zl(ir), c2,p2,mr2,t2, con)
	3418	do kr=1,nbox
	3419	ta(kr) = t2
	3420	end do
	3421	call interstrength (st2,t2,ka,ta)
	3422	do kr=1,nbox
	3423	! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5
	3424	c2box(kr) = c2 * ka(kr) * dble(deltaz)
	3425	end do
	3426	! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5
	3427	c2 = c2 * st2 * dble(deltaz)
	3428
	3429	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
	3430	bb = bb + ( p1c1 + p2c2 ) / 2.d0
	3431	cc = cc + ( c1 + c2 ) / 2.d0
	3432	dd = dd + ( t1c1 + t2c2 ) / 2.d0
	3433	do kr=1,nbox
	3434	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
	3435	ddbox(kr) = ddbox(kr) + ( t1c1box(kr) + t2c2box(kr) )/2.d0
	3436	end do
	3437
	3438	mr1=mr2
	3439	c1=c2
	3440	t1=t2
	3441	p1=p2
	3442	do kr=1,nbox
	3443	c1box(kr) = c2box(kr)
	3444	end do
	3445
	3446	pt = bb / cc
	3447	pp = aa / (cc * ff)
	3448	ts = dd / cc
	3449	do kr=1,nbox
	3450	ta(kr) = ddbox(kr) / ccbox(kr)
	3451	end do
	3452	call interstrength (st,ts,ka,ta)
	3453	call intershape (alsa,alna,alda,ta)
	3454
	3455	* ua = cc/st
	3456
	3457	c next loop calculates the eqw for an especified path ua,pp,pt,ta
	3458
	3459	eqwmu = 0.0d0
	3460	do im = 1,iimu
	3461	eqw=0.0d0
	3462	do kr=1,nbox
	3463	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
	3464	if(ua(kr).lt.0.)write(,)'mztf_overlap/674',ua(kr),
	3465	$ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
	3466
	3467	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
	3468	if ( i_supersat .eq. 0 ) then
	3469	eqw=eqw+no(kr)*w
	3470	else
	3471	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
	3472	endif
	3473	end do
	3474	eqwmu = eqwmu + eqw * mu(im)*amu(im)
	3475	end do
	3476
	3477	tau(in,ir) = exp( - eqwmu / dble(deltanux) )
	3478
	3479	4 continue
	3480
	3481	c
	3482	c due to the simmetry of the transmittances
	3483	c
	3484	do in=nl-1,2,-1
	3485	do ir=in-1,1,-1
	3486	tau(in,ir) = tau(ir,in)
	3487	end do
	3488	end do
	3489
	3490
	3491	ccc
	3492	ccc writing out transmittances
	3493	ccc
	3494	if (itauout.eq.1) then
	3495
	3496	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	3497	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	3498	! open( 1, file=
	3499	! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat',
	3500	! @ access='sequential', form='unformatted' )
	3501	! else
	3502	! open( 1, file=
	3503	! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat',
	3504	! @ access='sequential', form='unformatted' )
	3505	! endif
	3506
	3507	! write(1) dummy
	3508	! write(1)' format: (tauinf(n),(tau(n,r),r=1,nl),n=1,nl)'
	3509	! do in=1,nl
	3510	! write (1) tauinf(in), ( tau(in,ir), ir=1,nl )
	3511	! end do
	3512	! close(unit=1)
	3513
	3514	elseif (itauout.eq.2) then
	3515
	3516	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	3517	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	3518	! open( 1, file=
	3519	! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat')
	3520	! else
	3521	! open( 1, file=
	3522	! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat')
	3523	! endif
	3524
	3525	! !write(1,*) dummy
	3526	! !write(1,*) 'tij for curtis matrix calculations '
	3527	! !write(1,*)' cira mars model atmosphere '
	3528	! write(1,*)' beta= ',beta,'deltanu= ',deltanux
	3529	! write(1,*)' number of elements (in,ir)= ',nl,nl
	3530	! write(1,*)' format: (tauinf(in),(tau(in,ir),ir=1,nl),in=1,nl)'
	3531
	3532	! do in=1,nl
	3533	! write (1,*) tauinf(in)
	3534	! do ir=1,nl
	3535	! write(1,*) tau(in,ir)
	3536	! end do
	3537	! end do
	3538	! close(unit=1)
	3539
	3540	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	3541	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	3542	! write (*,'(1x, 31htransmitances written out in: ,a22)')
	3543	! @ 'taul'//isotcode//dn//ibcode1
	3544	! else
	3545	! write (*,'(1x, 31htransmitances written out in: ,a22)')
	3546	! @ 'taul'//isotcode//dn//ibcode2
	3547	! endif
	3548
	3549	end if
	3550
	3551	c cleaning of transmittances
	3552	! call elimin_tau(tau,tauinf,nl,nan,itableout,nw,dummy,
	3553	! @ isotcode,dn,ibcode2)
	3554
	3555	c construction of the curtis matrix
	3556
	3557	call mzcf ( tauinf,tau, cf,cfup,cfdw, vc,taugr,
	3558	@ ib,isot,icfout,itableout )
	3559
	3560
	3561	c end
	3562	return
	3563	end
	3564
	3565
	3566
	3567
	3568	c***********************************************************************
	3569	c mzcf
	3570	c***********************************************************************
	3571
	3572	subroutine mzcf( tauinf,tau, c,cup,cdw,vc,taugr,
	3573	@ ib,isot,icfout,itableout )
	3574
	3575	c a.k.murphy method to avoid extrapolation in the curtis matrix
	3576	c feb-89 m. angel granada
	3577	c 25-sept-96 cristina dejar las matrices en doble precision
	3578	c jan 98 malv version para mz1d
	3579	c jul 2011 malv+fgg adapted to LMD-MGCM
	3580	c***********************************************************************
	3581
	3582	implicit none
	3583
	3584	include 'comcstfi.h'
	3585	include 'nlte_paramdef.h'
	3586	include 'nlte_commons.h'
	3587
	3588	c arguments
	3589	real*8 c(nl,nl), cup(nl,nl), cdw(nl,nl) ! o
	3590	real*8 vc(nl), taugr(nl) ! o
	3591	real*8 tau(nl,nl) ! i
	3592	real*8 tauinf(nl) ! i
	3593	integer ib ! i
	3594	integer isot ! i
	3595	integer icfout, itableout ! i
	3596
	3597	c external
	3598	external bandid
	3599	character*2 bandid
	3600
	3601	c local variables
	3602	integer i, in, ir, iw
	3603	real*8 cfup(nl,nl), cfdw(nl,nl)
	3604	real*8 a(nl,nl), cf(nl,nl)
	3605	character isotcode2, bcode2
	3606
	3607	c formats
	3608	101 format(i1)
	3609	202 format(i2)
	3610	180 format(a80)
	3611	181 format(a80)
	3612	c***********************************************************************
	3613
	3614	if (isot.eq.1) isotcode = '26'
	3615	if (isot.eq.2) isotcode = '28'
	3616	if (isot.eq.3) isotcode = '36'
	3617	if (isot.eq.4) isotcode = '27'
	3618	if (isot.eq.5) isotcode = 'co'
	3619	bcode = bandid( ib )
	3620
	3621	! write (,) ' '
	3622
	3623	do in=1,nl
	3624
	3625	do ir=1,nl
	3626
	3627	cf(in,ir) = 0.0d0
	3628	cfup(in,ir) = 0.0d0
	3629	cfdw(in,ir) = 0.0d0
	3630	c(in,ir) = 0.0d0
	3631	cup(in,ir) = 0.0d0
	3632	cdw(in,ir) = 0.0d0
	3633	a(in,ir) = 0.0d0
	3634
	3635	end do
	3636
	3637	vc(in) = 0.0d0
	3638	taugr(in) = 0.0d0
	3639
	3640	end do
	3641
	3642
	3643	c the next lines are a reduced and equivalent way of calculating
	3644	c the c(in,ir) elements for n=2,nl1 and r=1,nl
	3645
	3646
	3647	c do in=2,nl1
	3648	c do ir=1,nl
	3649	c if(ir.eq.1)then
	3650	c c(in,ir)=tau(in-1,1)-tau(in+1,1)
	3651	c elseif(ir.eq.nl)then
	3652	c c(in,ir)=tau(in+1,nl1)-tauinf(in+1)-tau(in-1,nl1)+tauinf(in-1)
	3653	c else
	3654	c c(in,ir)=tau(in+1,ir-1)-tau(in+1,ir)-tau(in-1,ir-1)+tau(in-1,ir)
	3655	c end if
	3656	c c(in,ir)=c(in,ir)pideltanu(ib)/(2.deltaz1.0e5)
	3657	c end do
	3658	c end do
	3659	c go to 1000
	3660
	3661	c calculation of the matrix cfup(nl,nl)
	3662
	3663	cfup(1,1) = 1.d0 - tau(1,1)
	3664
	3665	do in=2,nl
	3666	do ir=1,in
	3667
	3668	if (ir.eq.1) then
	3669	cfup(in,ir) = tau(in,ir) - tau(in,1)
	3670	elseif (ir.eq.in) then
	3671	cfup(in,ir) = 1.d0 - tau(in,ir-1)
	3672	else
	3673	cfup(in,ir) = tau(in,ir) - tau(in,ir-1)
	3674	end if
	3675
	3676	end do
	3677	end do
	3678
	3679	! contribution to upwards fluxes from bb at bottom :
	3680	do in=1,nl
	3681	taugr(in) = tau(in,1)
	3682	enddo
	3683
	3684	c calculation of the matrix cfdw(nl,nl)
	3685
	3686	cfdw(nl,nl) = 1.d0 - tauinf(nl)
	3687
	3688	do in=1,nl-1
	3689	do ir=in,nl
	3690
	3691	if (ir.eq.in) then
	3692	cfdw(in,ir) = 1.d0 - tau(in,ir)
	3693	elseif (ir.eq.nl) then
	3694	cfdw(in,ir) = tau(in,ir-1) - tauinf(in)
	3695	else
	3696	cfdw(in,ir) = tau(in,ir-1) - tau(in,ir)
	3697	end if
	3698
	3699	end do
	3700	end do
	3701
	3702
	3703	c calculation of the matrix cf(nl,nl)
	3704
	3705	do in=1,nl
	3706	do ir=1,nl
	3707
	3708	if (ir.eq.1) then
	3709	! version con l_bb(tg) = l_bb(t(1))=j(1) (see also vc below)
	3710	! cf(in,ir) = tau(in,ir)
	3711	! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see also vc below)
	3712	cf(in,ir) = tau(in,ir) - tau(in,1)
	3713	elseif (ir.eq.nl) then
	3714	cf(in,ir) = tauinf(in) - tau(in,ir-1)
	3715	else
	3716	cf(in,ir) = tau(in,ir) - tau(in,ir-1)
	3717	end if
	3718
	3719	end do
	3720	end do
	3721
	3722
	3723	c definition of the a(nl,nl) matrix
	3724
	3725	do in=2,nl-1
	3726	do ir=1,nl
	3727	if (ir.eq.in+1) a(in,ir) = -1.d0
	3728	if (ir.eq.in-1) a(in,ir) = +1.d0
	3729	a(in,ir) = a(in,ir) / ( 2.d0deltaz1.d5 )
	3730	end do
	3731	end do
	3732	! this is not needed anymore in the akm scheme
	3733	! a(1,1) = +3.d0
	3734	! a(1,2) = -4.d0
	3735	! a(1,3) = +1.d0
	3736	! a(nl,nl) = -3.d0
	3737	! a(nl,nl1) = +4.d0
	3738	! a(nl,nl2) = -1.d0
	3739
	3740	c calculation of the final curtis matrix ("reduced" by murphy's method)
	3741
	3742	if (isot.ne.5) then
	3743	do in=1,nl
	3744	do ir=1,nl
	3745	cf(in,ir) = cf(in,ir) * pi*deltanu(isot,ib)
	3746	cfup(in,ir) = cfup(in,ir) * pi*deltanu(isot,ib)
	3747	cfdw(in,ir) = cfdw(in,ir) * pi*deltanu(isot,ib)
	3748	end do
	3749	taugr(in) = taugr(in) * pi*deltanu(isot,ib)
	3750	end do
	3751	else
	3752	do in=1,nl
	3753	do ir=1,nl
	3754	cf(in,ir) = cf(in,ir) * pi*deltanuco
	3755	enddo
	3756	taugr(in) = taugr(in) * pi*deltanuco
	3757	enddo
	3758	endif
	3759
	3760	do in=2,nl-1
	3761
	3762	do ir=1,nl
	3763
	3764	do i=1,nl
	3765	! only c contains the matrix a. matrixes cup,cdw dont because
	3766	! these two will be used for flux calculations, not
	3767	! only for flux divergencies
	3768
	3769	c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir)
	3770	! from this matrix we will extract (see below) the
	3771	! nl2 x nl2 "core" for the "reduced" final curtis matrix.
	3772
	3773	end do
	3774	cup(in,ir) = cfup(in,ir)
	3775	cdw(in,ir) = cfdw(in,ir)
	3776
	3777	end do
	3778	! version con l_bb(tg) = l_bb(t(1))=j(1) (see cf above)
	3779	!vc(in) = c(in,1)
	3780	! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see cf above)
	3781	vc(in) = pideltanu(isot,ib)/( 2.d0deltaz1.d5 )
	3782	@ ( tau(in-1,1) - tau(in+1,1) )
	3783
	3784	end do
	3785
	3786	5 continue
	3787
	3788	! write (,) 'mztf/1/ c(2,*) =', (c(2,i), i=1,nl)
	3789
	3790	! call elimin_dibuja(c,nl,itableout)
	3791
	3792	c ventana del smoothing de c es nw=3 y de vc es 5 (puesto en lisa):
	3793	c subroutine elimin_mz4(c,vc,ilayer,nl,nan,iw, nw)
	3794
	3795	iw = nan
	3796	if (isot.eq.4) iw = 5
	3797	call elimin_mz1d (c,vc,0,iw,itableout,nw)
	3798
	3799	! upper boundary condition
	3800	! j'(nl) = j'(nl1) ==> j(nl) = 2j(nl1) - j(nl2) ==>
	3801	do in=2,nl-1
	3802	c(in,nl-2) = c(in,nl-2) - c(in,nl)
	3803	c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl)
	3804	cup(in,nl-2) = cup(in,nl-2) - cup(in,nl)
	3805	cup(in,nl-1) = cup(in,nl-1) + 2.d0*cup(in,nl)
	3806	cdw(in,nl-2) = cdw(in,nl-2) - cdw(in,nl)
	3807	cdw(in,nl-1) = cdw(in,nl-1) + 2.d0*cdw(in,nl)
	3808	end do
	3809	! j(nl) = j(nl1) ==>
	3810	! do in=2,nl1
	3811	! c(in,nl1) = c(in,nl1) + c(in,nl)
	3812	! end do
	3813
	3814	! 1000 continue
	3815
	3816	if (icfout.eq.1) then
	3817
	3818	! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then
	3819	! codmatrx = codmatrx_fb
	3820	! else
	3821	! codmatrx = codmatrx_hot
	3822	! end if
	3823
	3824	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	3825	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	3826
	3827	! open ( 1, access='sequential', form='unformatted', file=
	3828	! @ dircurtis//'cfl'//isotcode//dn//ibcode1//codmatrx//'.dat')
	3829	! open ( 2, access='sequential', form='unformatted', file=
	3830	! @ dircurtis//'cflup'//isotcode//dn//ibcode1//codmatrx//'.dat')
	3831	! open ( 3, access='sequential', form='unformatted', file=
	3832	! @ dircurtis//'cfldw'//isotcode//dn//ibcode1//codmatrx//'.dat')
	3833
	3834	! else
	3835
	3836	! open ( 1, access='sequential', form='unformatted', file=
	3837	! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat')
	3838	! open ( 2, access='sequential', form='unformatted', file=
	3839	! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat')
	3840	! open ( 3, access='sequential', form='unformatted', file=
	3841	! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat')
	3842
	3843	! endif
	3844
	3845	! write(1) dummy
	3846	! write(1)' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)'
	3847	! do in=2,nl-1
	3848	! write(1) vc(in), (c(in,ir) , ir=2,nl-1 )
	3849	! es mas importante la precision que ocupar mucho espacio asi que
	3850	! escribiremos las matrices en doble precision y por tanto en
	3851	! [lib]readc_mz4.for no hay que reconvertirlas a doble precision.
	3852	! ch is stored in single prec. to save storage space.
	3853	! end do
	3854
	3855	! write(2) dummy
	3856	! write(2)' format: (cfup(n,r),r=1,nl), n=1,nl)'
	3857	! do in=1,nl
	3858	! write(2) ( cup(in,ir) , ir=1,nl )
	3859	! end do
	3860
	3861	! write(3) dummy
	3862	! write(3)' format: (cfdw(n,r),r=1,nl), n=1,nl)'
	3863	! do in=1,nl
	3864	! write(3) (cdw(in,ir) , ir=1,nl )
	3865	! end do
	3866
	3867	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
	3868	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
	3869	! write (*,'(1x,30hcurtis matrix written out in: ,a50)' )
	3870	! @ dircurtis//'cfl'//isotcode//dn//ibcode1//codmatrx//'.dat'
	3871	! else
	3872	! write (*,'(1x,30hcurtis matrix written out in: ,a50)' )
	3873	! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat'
	3874	! endif
	3875
	3876	else
	3877
	3878	! write (,) ' no curtis matrix output file ', char(10)
	3879
	3880	end if
	3881
	3882
	3883	c end
	3884	return
	3885	end
	3886
	3887
	3888
	3889
	3890
	3891	c***********************************************************************
	3892	c cm15um_hb_simple
	3893	c***********************************************************************
	3894
	3895	subroutine cm15um_hb_simple (ig,icurt)
	3896
	3897	c computing the curtix matrixes for the 15 um hot bands
	3898	c (las de las bandas fudnamentales las calcula cm15um_fb)
	3899
	3900	c jan 98 malv version de mod3/cm_15um.f para mz1d
	3901	c jul 2011 malv+fgg adapted to LMD-MGCM
	3902	c***********************************************************************
	3903
	3904	implicit none
	3905
	3906	!!!!!!!!!!!!!!!!!!!!!!!
	3907	! common variables & constants
	3908
	3909	include 'nlte_paramdef.h'
	3910	include 'nlte_commons.h'
	3911
	3912	!!!!!!!!!!!!!!!!!!!!!!!
	3913	! arguments
	3914
	3915	integer ig ! ADDED FOR TRACEBACK
	3916	integer icurt ! icurt=0,1,2
	3917	! new calculations? (see caa.f heads)
	3918
	3919	!!!!!!!!!!!!!!!!!!!!!!!
	3920	! local variables
	3921
	3922	real*4 cdummy(nl,nl), csngl(nl,nl)
	3923
	3924	real*8 cax1(nl,nl), cax2(nl,nl), cax3(nl,nl)
	3925	real*8 v1(nl), v2(nl), v3(nl), cm_factor, vc_factor
	3926
	3927	integer itauout,icfout,itableout, interpol,ismooth, isngldble
	3928	integer i,j,ik,ist,isot,ib,itt
	3929
	3930	!character bandcode*2
	3931	character isotcode*2
	3932	character codmatrx_hot*5
	3933
	3934	!!!!!!!!!!!!!!!!!!!!!!!
	3935	! external functions
	3936
	3937	external bandid
	3938	character*2 bandid
	3939
	3940	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3941	! subroutines called:
	3942	! mz4sub, dmzout, readc_mz4, readcupdw, mztf
	3943
	3944	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3945	! formatos
	3946	132 format(i2)
	3947
	3948	************************************************************************
	3949	************************************************************************
	3950
	3951	call zerom (c121,nl)
	3952
	3953	call zerov (vc121,nl)
	3954
	3955	call zerom (cup121,nl)
	3956	call zerom (cdw121,nl)
	3957
	3958	call zerov (taugr121,nl)
	3959
	3960
	3961	itauout = 0 ! =1 --> with output of tau
	3962	icfout = 0 ! =1 --> with output of cf
	3963	itableout = 0 ! =1 --> with output of table of taus
	3964	isngldble = 1 ! =1 --> dble precission
	3965
	3966	codmatrx_hot=' '
	3967	if (icurt.eq.2) then
	3968	icfout=1
	3969	elseif (icurt.eq.0) then
	3970	write (*,'(a,a$)')
	3971	@ ' hot bands. code for old matrixes (5 chars): '
	3972	read (*,'(a)') codmatrx_hot
	3973	endif
	3974
	3975	fileroot = 'cfl'
	3976
	3977	! ====================== curtis matrixes for fh bands ==================
	3978
	3979
	3980	! una piedra en el camino ...
	3981	! write (,) ' cm15um_hb/1 '
	3982
	3983	ccc
	3984	if ( input_cza.ge.1 ) then
	3985	ccc
	3986
	3987	if (icurt.eq.2) then
	3988	write (*,'(a,a$)')
	3989	@ ' new calculation of curt. mat. for fh bands.',
	3990	@ ' code for new matrixes : '
	3991	read (*,'(a)') codmatrx_hot
	3992	elseif (icurt.eq.0) then
	3993	write (*,'(a,a$)')
	3994	@ ' reading in curt. mat. for fh bands.',
	3995	@ ' code for old matrixes : '
	3996	read (*,'(a)') codmatrx_hot
	3997	else
	3998	! write (*,'(a)')
	3999	! @ ' new calculation of curt. mat. for fh bands.'
	4000	end if
	4001
	4002	! fh bands for the 626 isotope ================================-
	4003
	4004	ist = 1
	4005	isot = 26
	4006	! encode (2,132,isotcode) isot
	4007	write (isotcode,132) isot
	4008
	4009	do 11, ik=1,3
	4010
	4011	ib=ik+1
	4012
	4013	if (icurt.gt.0) then
	4014	call zero3m (cax1,cax2,cax3,nl)
	4015	! una piedra en el camino ...
	4016	!write (,) ' cm15um_hb/11 '
	4017	!write (,) ' ib, ist, irw, imu =', ib, ist, irw_mztf, imu
[695]	4018	call mztf(ig,cax1,cax2,cax3,v1,v2,ib,ist,irw_mztf,imu,
[498]	4019	@ itauout,icfout,itableout)
	4020	! else
	4021	! bandcode = bandid(ib)
	4022	! filend=isotcode//dn//bandcode//codmatrx_hot
	4023	!! write (,) char(9), fileroot//filend
	4024	! call zero3m (cax1,cax2,cax3,nl)
	4025	! call readcud_mz1d ( cax1,cax2,cax3, v1, v2,
	4026	! @ fileroot,filend, csngl, nl,nan,0,isngldble)
	4027	end if
	4028
	4029	c calculating the total c121(n,r) matrix for the first hot band
	4030	do i=1,nl
	4031
	4032	if ( ib .eq. 4 ) then
	4033	! write (,) ' '
	4034	! write (,) i, ' ib,ist, altura :', ib,ist, zl(i)
	4035	endif
	4036
	4037	! if ( v1(i) .le. 1.d-99 ) v1(i) = 0.0d0
	4038	! if ( v2(i) .le. 1.d-99 ) v2(i) = 0.0d0
	4039
	4040
	4041	if(ik.eq.1)then
	4042	cm_factor = (dble(618.03/667.75))*2.d0
	4043	@ exp( dble(ee*(667.75-618.03)/t(i)) )
	4044	vc_factor = dble(667.75/618.03)
	4045	elseif(ik.eq.2)then
	4046	cm_factor = 1.d0
	4047	vc_factor = 1.d0
	4048	elseif(ik.eq.3)then
	4049	cm_factor = ( dble(720.806/667.75) )*2.d0
	4050	@ exp( dble(ee*(667.75-720.806)/t(i)) )
	4051	vc_factor = dble(667.75/720.806)
	4052	end if
	4053	do j=1,nl
	4054	! if (cax1(i,j) .le. 1.d-99 ) cax1(i,j) = 0.0d0
	4055	! if (cax2(i,j) .le. 1.d-99 ) cax2(i,j) = 0.0d0
	4056	! if (cax3(i,j) .le. 1.d-99 ) cax3(i,j) = 0.0d0
	4057	c121(i,j) = c121(i,j) + cax1(i,j) * cm_factor
	4058	cup121(i,j) = cup121(i,j) + cax2(i,j) * cm_factor
	4059	cdw121(i,j) = cdw121(i,j) + cax3(i,j) * cm_factor
	4060	end do
	4061
	4062	! write (,) ' i =', i
	4063	! write (,) ' vc_factor =', vc_factor
	4064	! write (,) ' v1 =', v1(i)
	4065	! write (,) ' v2 =', v2(i)
	4066	! write (,) vc121(i), taugr121(i)
	4067	! write (,) v1(i) * vc_factor
	4068	! write (,) vc121(i) + v1(i) * vc_factor
	4069
	4070	vc121(i) = vc121(i) + v1(i) * vc_factor
	4071
	4072
	4073	! write (,) v2(i) * vc_factor
	4074	! write (,) taugr121(i) + v2(i) * vc_factor
	4075
	4076	taugr121(i) = taugr121(i) + v2(i) * vc_factor
	4077
	4078	end do
	4079	11 continue
	4080
	4081	ccc
	4082	end if
	4083	ccc
	4084
	4085
	4086	return
	4087	end
	4088
	4089
	4090
	4091
	4092

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

Download in other formats:

Original Format