Context Navigation

nlte_calc.F @ 623

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

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

Download in other formats:

Original Format