Index: trunk/LMDZ.MARS/libf/phymars/nlte_calc.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/nlte_calc.F	(revision 695)
+++ trunk/LMDZ.MARS/libf/phymars/nlte_calc.F	(revision 757)
@@ -1,4092 +1,1678 @@
-c***********************************************************************
-c     mzescape.f                                    
-c***********************************************************************
-c                                               
-c     program  for calculating atmospheric escape functions, from           
-c     a calculation of transmittances and derivatives of these ones   
+c**********************************************************************
+c     
+c     Includes the following 1-d model subroutines:
+c     
+c     -MZESC110_dlvr11_03.f
+c     -MZTUD110_dlvr11_03.f
+c     -MZMC121_dlvr11_03.f
+c     -MZTUD121_dlvr11_03.f
+c     -MZESC121_dlvr11_03.f
+c     -MZESC121sub_dlvr11_03.f
+c     -MZTVC121_dlvr11.f
+c     -MZTVC121sub_dlvr11_03.f
+
+
+
+c     *** Old MZESC110_dlvr11_03.f
+
+c**********************************************************************
+
+c***********************************************************************
+      subroutine MZESC110 (nl_cts_real, nzy_cts_real)
+c***********************************************************************
+
+      implicit none
+
+      include 'datafile.h'
+      include 'nlte_paramdef.h'
+      include 'nlte_commons.h'
+
+c     arguments 
+      integer     nl_cts_real, nzy_cts_real ! i
+
+c     old arguments
+      integer         ierr      ! o
+      real*8          varerr    ! o
+
+c     local variables and constants
+      integer 	i, in, ir, iaquiHIST , iaquiZ
+      integer         ib, isot
+      real*8 		argumento
+      real*8 		tauinf(nl_cts)
+      real*8 		con(nzy_cts), coninf
+      real*8 		c1, c2 , ccc
+      real*8 		t1, t2
+      real*8 		p1, p2
+      real*8		mr1, mr2
+      real*8 		st1, st2
+      real*8 		c1box(nbox_max), c2box(nbox_max)
+      real*8		ff      ! to avoid too small numbers
+      real*8 		st, beta, ts
+      real*8  	tyd(nzy_cts)
+      real*8 		correc
+      real*8		deltanudbl, deltazdbl
+      real*8          yy
+
+c     external function
+      external        we_clean
+      real*8          we_clean
+
+c***********************************************************************
+
+c     
+      ib = 1
+      beta = 1.8d5
+      ibcode1 = '1'
+      isot = 1
+      deltanudbl = dble(deltanu(1,1))
+      deltazdbl = dble(deltaz_cts)
+      ff=1.0d10
+
+ccc   
+      do i=1,nzy_cts_real
+         tyd(i) = dble(ty_cts(i))
+         con(i) =  dble( co2y_cts(i) * imr(isot) )
+         correc = 2.d0 * dexp( -ee*dble(elow(isot,2))/tyd(i) )
+         con(i) = con(i) * ( 1.d0 - correc )
+         mr_cts(i) = dble(co2y_cts(i)/nty_cts(i))
+      end do
+      coninf = dble( con(nzy_cts_real) /
+     @     log( con(nzy_cts_real-1) / con(nzy_cts_real) ) )
+                                ! Correccion pequeña para la FB, nos la saltamos
+                                !call mztf_correccion_cts ( coninf, con, ib )
+
+ccc   
+      call gethist_03 ( 1 )
+
+c     
+c     tauinf
+c     
+      call initial
+
+      iaquiHIST = nhist/2
+      iaquiZ = nzy_cts_real - 2
+
+      do i=nl_cts_real,1,-1
+
+         if(i.eq.nl_cts_real)then
+
+            call intzhunt_cts (iaquiZ, zl_cts(i), nzy_cts_real, 
+     @           c2,p2,mr2,t2, con)
+            do kr=1,nbox
+               ta(kr)=t2
+	    end do
+            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
+                                ! Check interpolation errors :
+            if (c2.le.0.0d0) then
+               ierr=15
+               varerr=c2
+               return
+       	    elseif (p2.le.0.0d0) then
+               ierr=16
+               varerr=p2
+               return
+	    elseif (mr2.le.0.0d0) then
+               ierr=17
+               varerr=mr2
+               return
+	    elseif (t2.le.0.0d0) then
+               ierr=18
+               varerr=t2
+               return
+	    elseif (st2.le.0.0d0) then
+               ierr=19
+               varerr=st2
+               return
+	    endif
+                                !
+            aa = p2 * coninf * mr2 * (st2 * ff)
+            cc = coninf * st2
+            dd = t2 * coninf * st2
+            do kr=1,nbox
+               ccbox(kr) = coninf * ka(kr)
+               ddbox(kr) = t2 * ccbox(kr)
+               c2box(kr) = c2 * ka(kr) * deltazdbl
+            end do
+            c2 = c2 * st2 * deltazdbl
+
+         else
+
+            call intzhunt_cts (iaquiZ, zl_cts(i), nzy_cts_real, 
+     @           c1,p1,mr1,t1, con)
+            do kr=1,nbox
+               ta(kr)=t1
+            end do
+            call interstrhunt (iaquiHIST, st1,t1,ka,ta)
+            do kr=1,nbox
+               c1box(kr) = c1 * ka(kr) * deltazdbl
+            end do
+            c1 = c1 * st1 * deltazdbl
+            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
+            cc = cc + ( c1 + c2 ) / 2.d0
+            ccc = ( c1 + c2 ) / 2.d0
+            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
+            do kr=1,nbox
+               ccbox(kr) = ccbox(kr) +
+     @              ( c1box(kr) + c2box(kr) )/2.d0
+               ddbox(kr) = ddbox(kr) +
+     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0
+            end do
+
+            mr2 = mr1
+            c2=c1
+            do kr=1,nbox
+               c2box(kr) = c1box(kr)
+            end do
+            t2=t1
+            p2=p1
+         end if
+
+         pp = aa / (cc*ff)
+
+         ts = dd/cc
+         do kr=1,nbox
+   	    ta(kr) = ddbox(kr) / ccbox(kr)
+         end do
+         call interstrhunt(iaquiHIST, st,ts,ka,ta)
+         call intershphunt(iaquiHIST, alsa,alda,ta)
+
+c     
+         eqw=0.0d0
+         do  kr=1,nbox
+            yy = ccbox(kr) * beta
+            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
+            eqw = eqw + no(kr)*w
+         end do
+
+         argumento = eqw / deltanudbl
+         tauinf(i) = dexp( - argumento )
+
+         if (i.eq.nl_cts_real) then
+            taustar11_cts(i) = 0.0d0
+         else
+            taustar11_cts(i) = deltanudbl * (tauinf(i+1)-tauinf(i))
+     @           / ( beta * ccc )
+         endif
+
+      end do
+
+
+      call mzescape_normaliz_02 ( taustar11_cts, nl_cts_real, 2 )
+
+c     end
+      return
+      end
+
+
+c     *** Old MZTUD110_dlvr11_03.f
+
+c***********************************************************************
+      subroutine MZTUD110( ierr, varerr )
+c***********************************************************************
+
+      implicit none
+
+      include 'datafile.h'
+      include 'nlte_paramdef.h'
+      include 'nlte_commons.h'
+
+
+c     arguments
+      integer         ierr      ! o
+      real*8          varerr    ! o
+
+c     local variables and constants
+      integer 	i, in, ir, iaquiHIST , iaquiZ
+      integer         ib, isot
+      real*8 		tau(nl,nl), argumento
+      real*8 		tauinf(nl)
+      real*8 		con(nzy), coninf
+      real*8 		c1, c2
+      real*8 		t1, t2
+      real*8 		p1, p2
+      real*8		mr1, mr2
+      real*8 		st1, st2
+      real*8 		c1box(nbox_max), c2box(nbox_max)
+      real*8		ff      ! to avoid too small numbers
+      real*8		tvtbs(nzy)
+      real*8 		st, beta, ts
+      real*8  	zld(nl), zyd(nzy), deltazdbl
+      real*8 		correc
+      real*8		deltanudbl
+      real*8          maxtau, yy
+
+c     external function
+      external        we_clean
+      real*8          we_clean
+
+c***********************************************************************
+
+c     
+      ib = 1
+      beta = 1.8d5
+      ibcode1 = '1'
+      isot = 1
+      deltanudbl = dble(deltanu(1,1))
+      deltazdbl = dble(deltaz)
+      ff=1.0d10
+
+ccc   
+      do i=1,nzy
+         zyd(i) = dble(zy(i))
+      enddo
+      do i=1,nl
+         zld(i) = dble(zl(i))
+      enddo
+      call interhuntdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
+      do i=1,nzy
+         con(i) =  dble( co2y(i) * imr(isot) )
+         correc = 2.d0 * dexp( -ee*dble(elow(isot,2))/tvtbs(i) )
+         con(i) = con(i) * ( 1.d0 - correc )
+         mr(i) = dble(co2y(i)/nty(i))
+      end do
+      coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
+      call mztf_correccion ( coninf, con, ib )
+
+ccc   
+      call gethist_03 ( 1 )
+
+c     
+c     tauinf
+c     
+      call initial
+
+      iaquiHIST = nhist/2
+      iaquiZ = nzy - 2
+
+      do i=nl,1,-1
+
+         if(i.eq.nl)then
+
+            call intzhunt (iaquiZ, zl(i),c2,p2,mr2,t2, con)
+            do kr=1,nbox
+               ta(kr)=t2
+            end do
+            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
+                                ! Check interpolation errors :
+            if (c2.le.0.0d0) then
+               ierr=15
+               varerr=c2
+               return
+            elseif (p2.le.0.0d0) then
+               ierr=16
+               varerr=p2
+               return
+            elseif (mr2.le.0.0d0) then
+               ierr=17
+               varerr=mr2
+               return
+            elseif (t2.le.0.0d0) then
+               ierr=18
+               varerr=t2
+               return
+            elseif (st2.le.0.0d0) then
+               ierr=19
+               varerr=st2
+               return
+            endif
+                                !
+            aa = p2 * coninf * mr2 * (st2 * ff)
+            cc = coninf * st2
+            dd = t2 * coninf * st2
+            do kr=1,nbox
+               ccbox(kr) = coninf * ka(kr)
+               ddbox(kr) = t2 * ccbox(kr)
+               c2box(kr) = c2 * ka(kr) * deltazdbl
+            end do
+            c2 = c2 * st2 * deltazdbl
+
+         else
+
+            call intzhunt (iaquiZ, zl(i),c1,p1,mr1,t1, con)
+            do kr=1,nbox
+               ta(kr)=t1
+            end do
+            call interstrhunt (iaquiHIST, st1,t1,ka,ta)
+            do kr=1,nbox
+               c1box(kr) = c1 * ka(kr) * deltazdbl
+            end do
+            c1 = c1 * st1 * deltazdbl
+            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
+            cc = cc + ( c1 + c2 ) / 2.d0
+            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
+            do kr=1,nbox
+               ccbox(kr) = ccbox(kr) +
+     @              ( c1box(kr) + c2box(kr) )/2.d0
+               ddbox(kr) = ddbox(kr) +
+     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0
+            end do
+
+            mr2 = mr1
+            c2=c1
+            do kr=1,nbox
+               c2box(kr) = c1box(kr)
+            end do
+            t2=t1
+            p2=p1
+         end if
+
+         pp = aa / (cc*ff)
+
+         ts = dd/cc
+         do kr=1,nbox
+   	    ta(kr) = ddbox(kr) / ccbox(kr)
+         end do
+         call interstrhunt(iaquiHIST, st,ts,ka,ta)
+         call intershphunt(iaquiHIST, alsa,alda,ta)
+
+c     
+         eqw=0.0d0
+         do  kr=1,nbox
+            yy = ccbox(kr) * beta
+            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
+            eqw = eqw + no(kr)*w
+         end do
+
+         argumento = eqw / deltanudbl
+         tauinf(i) = dexp( - argumento )
+
+      end do
+
+
+c     
+c     tau
+c     
+
+      iaquiZ = 2
+      do 1 in=1,nl-1
+
+         call initial
+         call intzhunt (iaquiZ, zl(in), c1,p1,mr1,t1, con)
+         do kr=1,nbox
+            ta(kr) = t1
+         end do
+         call interstrhunt (iaquiHIST, st1,t1,ka,ta)
+         do kr=1,nbox
+            c1box(kr) = c1 * ka(kr) * deltazdbl
+         end do
+         c1 = c1 * st1 * deltazdbl
+
+         do 2 ir=in,nl-1
+
+            if (ir.eq.in) then
+               tau(in,ir) = 1.d0
+               goto 2
+            end if
+
+            call intzhunt (iaquiZ, zl(ir), c2,p2,mr2,t2, con)
+            do kr=1,nbox
+               ta(kr) = t2
+            end do
+            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
+            do kr=1,nbox
+               c2box(kr) = c2 * ka(kr) * deltazdbl
+            end do
+            c2 = c2 * st2 * deltazdbl
+
+            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
+            cc = cc + ( c1 + c2 ) / 2.d0
+            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
+            do kr=1,nbox
+               ccbox(kr) = ccbox(kr) + 
+     $              ( c1box(kr) + c2box(kr) ) / 2.d0
+               ddbox(kr) = ddbox(kr) + 
+     $              ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
+            end do
+
+            mr1=mr2
+            t1=t2
+            c1=c2
+            p1=p2
+            do kr=1,nbox
+               c1box(kr) = c2box(kr)
+            end do
+
+            pp = aa / (cc * ff)
+
+            ts = dd/cc
+            do kr=1,nbox
+               ta(kr) = ddbox(kr) / ccbox(kr)
+            end do
+            call interstrhunt(iaquiHIST, st,ts,ka,ta)
+            call intershphunt(iaquiHIST, alsa,alda,ta)
+c     
+            eqw=0.0d0
+            do kr=1,nbox
+               yy = ccbox(kr) * beta
+               w = we_clean ( yy, pp, alsa(kr),alda(kr) )
+               eqw = eqw + no(kr)*w
+            end do
+
+            argumento = eqw / deltanudbl
+            tau(in,ir) = exp( - argumento )
+
+
+ 2       continue
+
+ 1    continue
+
+
+c     
+c     tau(in,ir) for n>r
+c     
+
+      in=nl
+
+      call initial
+
+      iaquiZ = nzy - 2
+      call intzhunt (iaquiZ, zl(in), c1,p1,mr1,t1, con)
+      do kr=1,nbox
+         ta(kr) = t1
+      end do
+      call interstrhunt (iaquiHIST,st1,t1,ka,ta)
+      do kr=1,nbox
+         c1box(kr) = c1 * ka(kr) * deltazdbl
+      end do
+      c1 = c1 * st1 * deltazdbl
+
+      do 4 ir=in-1,1,-1
+
+         call intzhunt (iaquiZ, zl(ir), c2,p2,mr2,t2, con)
+         do kr=1,nbox
+            ta(kr) = t2
+         end do
+         call interstrhunt (iaquiHIST, st2,t2,ka,ta)
+         do kr=1,nbox
+            c2box(kr) = c2 * ka(kr) * deltazdbl
+         end do
+         c2 = c2 * st2 * deltazdbl
+
+         aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
+         cc = cc + ( c1 + c2 ) / 2.d0
+         dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
+         do kr=1,nbox
+            ccbox(kr) = ccbox(kr) + 
+     $           ( c1box(kr) + c2box(kr) ) / 2.d0
+            ddbox(kr) = ddbox(kr) + 
+     $           ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
+         end do
+
+         mr1=mr2
+         c1=c2
+         t1=t2
+         p1=p2
+         do kr=1,nbox
+            c1box(kr) = c2box(kr)
+         end do
+
+         pp = aa / (cc * ff)
+         ts = dd / cc
+         do kr=1,nbox
+            ta(kr) = ddbox(kr) / ccbox(kr)
+         end do
+         call interstrhunt (iaquiHIST, st,ts,ka,ta)
+         call intershphunt (iaquiHIST, alsa,alda,ta)
+
+c     
+
+         eqw=0.0d0
+         do kr=1,nbox
+            yy = ccbox(kr) * beta
+            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
+            eqw = eqw + no(kr)*w
+         end do
+
+         argumento = eqw / deltanudbl
+         tau(in,ir) = exp( - argumento )
+
+
+ 4    continue
+
+c     
+c     
+c     
+      do in=nl-1,2,-1
+         do ir=in-1,1,-1
+            tau(in,ir) = tau(ir,in)
+         end do
+      end do
+
+c     
+c     Tracking potential numerical errors
+c     
+      maxtau = 0.0d0
+      do in=nl-1,2,-1
+         do ir=in-1,1,-1
+            maxtau = max( maxtau, tau(in,ir) )
+         end do
+      end do
+      if (maxtau .gt. 1.0d0) then
+         ierr = 13
+         varerr = maxtau
+         return
+      endif
+
+
+c     
+      call MZCUD110 ( tauinf,tau )
+
+c     end
+      return
+      end
+
+
+c     *** Old file MZCUD_dlvr11.f ***
+
+c***********************************************************************
+
+      subroutine MZCUD110 ( tauinf,tau )
+
+c***********************************************************************
+
+      implicit none
+
+      include 'nlte_paramdef.h'
+      include 'nlte_commons.h'
+
+c     arguments
+      real*8 		tau(nl,nl) ! i
+      real*8		tauinf(nl) ! i
+
+
+c     local variables
+      integer 	i, in, ir
+      real*8		a(nl,nl), cf(nl,nl), pideltanu, deltazdp, pi
+
+c***********************************************************************
+
+      pi = 3.141592
+      pideltanu = pi * dble(deltanu(1,1))
+      deltazdp = 2.0d5 * dble(deltaz)
+
+      do in=1,nl
+         do ir=1,nl
+            cf(in,ir) = 0.0d0
+            c110(in,ir) = 0.0d0
+            a(in,ir) = 0.0d0
+         end do
+         vc110(in) = 0.0d0
+      end do
+
+c     
+      do in=1,nl
+         do ir=1,nl
+
+            if (ir.eq.1) then
+               cf(in,ir) = tau(in,ir) - tau(in,1)
+            elseif (ir.eq.nl) then
+               cf(in,ir) = tauinf(in) - tau(in,ir-1)
+            else
+               cf(in,ir) = tau(in,ir) - tau(in,ir-1)
+            end if
+
+         end do
+      end do
+
+c     
+      do in=2,nl-1
+         do ir=1,nl
+            if (ir.eq.in+1) a(in,ir) = -1.d0
+            if (ir.eq.in-1) a(in,ir) = +1.d0
+            a(in,ir) = a(in,ir) / deltazdp
+         end do
+      end do
+
+c     
+      do in=1,nl
+         do ir=1,nl
+	    cf(in,ir) = cf(in,ir) * pideltanu
+         end do
+      end do
+
+      do in=2,nl-1
+         do ir=1,nl
+	    do i=1,nl
+               c110(in,ir) = c110(in,ir) + a(in,i) * cf(i,ir)
+	    end do
+         end do
+      end do
+
+      do in=2,nl-1
+         vc110(in) =  pideltanu/deltazdp *
+     @        ( tau(in-1,1) - tau(in+1,1) )
+      end do
+
+
+c     
+      do in=2,nl-1
+         c110(in,nl-2) = c110(in,nl-2) - c110(in,nl)
+         c110(in,nl-1) = c110(in,nl-1) + 2.d0*c110(in,nl)
+      end do
+
+c     end
+      return
+      end
+
+
+c     *** Old MZMC121_dlvr11_03.f ***
+
+c***********************************************************************
+
+      subroutine MZMC121
+
+c***********************************************************************
+
+      implicit none
+
+                                ! common variables & constants
       
-      subroutine mzescape(ig,taustar,tauinf,tauii,
-     &  ib,isot, iirw,iimu) 
-
-c     jul 2011        malv+fgg   adapted to LMD-MGCM                        
-c     nov 99          malv    adapt mztf to compute taustar (pg.23b-ma
-c     nov 98          malv    allow for overlaping in the lorentz line
-c     jan 98		malv	version for mz1d. based on curtis/mztf.for   
-c     17-jul-96	mlp&crs	change the calculation of mr.     
-c     evitar: divide por cero. anhadiendo: ff    
-c     oct-92		malv 	correct s(t) dependence for all histogr bands
-c     june-92		malv	proper lower levels for laser bands         
-c     may-92		malv 	new temperature dependence for laser bands  
-c     @    991 		malv 	boxing for the averaged absorber amount and t
-c     ?		malv	extension up to 200 km altitude in mars 
-c     13-nov-86	mlp	include the temperature weighted to match 
-c				the eqw in the strong doppler limit.       
-c***********************************************************************
-                                                            
-      implicit none                                  
-                                                            
+      include 'nlte_paramdef.h'
+      include 'nlte_commons.h'
+
+                                ! local variables
+
+      real*8 cax1(nl,nl)
+      real*8 v1(nl), cm_factor, vc_factor
+      real nuaux1, nuaux2, nuaux3
+      real*8 faux2,faux3, daux2,daux3
+
+      integer i,j,ik,ib
+
+************************************************************************
+
+      c121(1:nl,1:nl)=0.d0
+!      call zerom (c121,nl)
+      vc121(1:nl)=0.d0
+!      call zerov (vc121,nl)
+
+      nuaux1 = nu(1,2) - nu(1,1) ! 667.75
+      nuaux2 = nu12_0200-nu(1,1) ! 618.03
+      nuaux3 = nu12_1000-nu(1,1) ! 720.81
+      faux2 = dble(nuaux2/nuaux1)
+      faux3 = dble(nuaux3/nuaux1)
+      daux2 = dble(nuaux1-nuaux2)
+      daux3 = dble(nuaux1-nuaux3)
+
+      do 11, ik=1,3
+
+         ib=ik+1
+         cax1(1:nl,1:nl)=0.d0
+!         call zerom (cax1,nl)
+         call MZTUD121 ( cax1,v1, ib )
+
+         do i=1,nl
+
+	    if(ik.eq.1)then
+               cm_factor = faux2**2.d0 * exp( daux2*ee/dble(t(i)) )
+               vc_factor = 1.d0/faux2
+	    elseif(ik.eq.2)then
+               cm_factor = 1.d0
+               vc_factor = 1.d0
+	    elseif(ik.eq.3)then
+               cm_factor = faux3**2.d0 * exp( daux3*ee/dble(t(i)) )
+               vc_factor = 1.d0 / faux3
+            else
+               write (*,*) ' Error in 626 hot band index  ik =', ik
+               stop ' ik can only be = 2,3,4.   Check needed.'
+	    end if
+	    do j=1,nl
+               c121(i,j) = c121(i,j) + cax1(i,j) * cm_factor
+	    end do
+
+	    vc121(i) = vc121(i) + v1(i) * vc_factor
+
+         end do
+
+ 11   continue
+
+      return
+      end
+
+
+c     *** Old MZTUD121_dlvr11_03.f ***
+
+c***********************************************************************
+      subroutine MZTUD121 ( cf,vc, ib )
+c***********************************************************************
+
+      implicit none
+
+      include 'datafile.h'
       include 'nlte_paramdef.h'
       include 'nlte_commons.h'
       
-                                                            
-c arguments          
-      integer         ig        ! ADDED FOR TRACEBACK
-      real*8 		taustar(nl) ! o                   
-      real*8 		tauinf(nl) ! o                   
-      real*8 		tauii(nl) ! o                   
-      integer		ib      ! i 
-      integer		isot    ! i
-      integer		iirw    ! i
-      integer		iimu    ! i
-                                                            
-                                                            
-c local variables and constants                 
-      integer 	i, in, ir, im, k,j                      
-      integer 	nmu                                   
-      parameter 	(nmu = 8)                          
-!     real*8 		tauinf(nl)                           
-      real*8 		con(nzy), coninf                       
-      real*8 		c1, c2, ccc                           
-      real*8 		t1, t2                                
-      real*8 		p1, p2                                
-      real*8		mr1, mr2                               
-      real*8 		st1, st2                              
-      real*8 		c1box(70), c2box(70)                  
+
+c     arguments
+      real*8  	cf(nl,nl)	! o.
+      real*8		vc(nl)  ! o
+      integer		ib      ! i
+
+
+c     local variables and constants
+      integer 	i, in, ir, iaquiHIST, iaquiZ
+      integer 	nmu
+      parameter 	(nmu = 8)
+      real*8 		tau(nl,nl), argumento, deltazdbl
+      real*8 		tauinf(nl)
+      real*8 		con(nzy), coninf
+      real*8 		c1, c2
+      real*8 		t1, t2
+      real*8 		p1, p2
+      real*8		mr1, mr2
+      real*8 		st1, st2
+      real*8 		c1box(70), c2box(70)
       real*8		ff      ! to avoid too small numbers
-      real*8		tvtbs(nzy)                              
-      real*8 		st, beta, ts, eqwmu                   
-      real*8 		mu(nmu), amu(nmu)                     
-      real*8  	zld(nl), zyd(nzy)                               
-      real*8 		correc                                
-      real 		deltanux ! width of vib-rot band (cm-1)
-      character	isotcode*2	                          
-      real*8          maximum                        
-      real*8          csL, psL, Desp, wsL ! for Strong Lorentz limit
-                                                            
-c formats                                       
- 111  format(a1)                                 
- 112  format(a2)                                 
- 101  format(i1)                                 
- 202  format(i2)                                 
- 180  format(a80)                                
- 181  format(a80)                                
-c***********************************************************************
-                                                            
-c some needed values                            
-!     rl=sqrt(log(2.d0))                             
-!     pi2 = 3.14159265358989d0                       
-      beta = 1.8d0                                   
-!     imrco = 0.9865                                 
-      
-c  esto es para que las subroutines de mztfsub calculen we  
-c  de la forma apropiada para mztf, no para fot 
-      icls=icls_mztf                                 
-                                                            
-c codigos para filenames                        
-!     if (isot .eq. 1)  isotcode = '26'              
-!     if (isot .eq. 2)  isotcode = '28'              
-!     if (isot .eq. 3)  isotcode = '36'              
-!     if (isot .eq. 4)  isotcode = '27'              
-!     if (isot .eq. 5)  isotcode = '62'              
-!     if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5     
-!     @    	.or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then      
-!!		encode(2,101,ibcode1)ib                       
-!     write ( ibcode1, 101) ib                       
-!     else                                           
-!!		encode(2,202,ibcode2)ib
-!     write (ibcode2, 202) ib
-!     endif                                          
-!     write (*,'( 30h calculating curtis matrix :  ,2x,         
-!     @    	8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot         
-                                                            
-c integration in angle !!!!!!!!!!!!!!!!!!!!     
-                                                            
-c------- diffusivity approx.                    
-      if (iimu.eq.1) then                            
-!	  write (*,*)  ' diffusivity approx. beta = ',beta             
-         mu(1) = 1.0d0                                
-         amu(1)= 1.0d0                                
-c-------data for 8 points integration           
-      elseif (iimu.eq.4) then                        
-         write (*,*)' 4 points for the gauss-legendre angle quadrature.'
-         mu(1)=(1.0d0+0.339981043584856)/2.0d0 	      
-         mu(2)=(1.0d0-0.339981043584856)/2.0d0 	      
-         mu(3)=(1.0d0+0.861136311594053)/2.0d0        
-         mu(4)=(1.0d0-0.861136311594053)/2.0d0 	      
-         amu(1)=0.652145154862546 	                   
-         amu(2)=amu(1) 	                              
-         amu(3)=0.347854845137454 	                   
-         amu(4)=amu(3)                                
-         beta=1.0d0                                   
-c-------data for 8 points integration           
-      elseif(iimu.eq.8) then                         
-         write (*,*)' 8 points for the gauss-legendre angle quadrature.'
-         mu(1)=(1.0d0+0.183434642495650)/2.0d0        
-         mu(2)=(1.0d0-0.183434642495650)/2.0d0        
-         mu(3)=(1.0d0+0.525532409916329)/2.0d0        
-         mu(4)=(1.0d0-0.525532409916329)/2.0d0        
-         mu(5)=(1.0d0+0.796666477413627)/2.0d0        
-         mu(6)=(1.0d0-0.796666477413627)/2.0d0        
-         mu(7)=(1.0d0+0.960289856497536)/2.0d0        
-         mu(8)=(1.0d0-0.960289856497536)/2.0d0        
-         amu(1)=0.362683783378362                     
-         amu(2)=amu(1)                                
-         amu(3)=0.313706645877887                     
-         amu(4)=amu(3)                                
-         amu(5)=0.222381034453374                     
-         amu(6)=amu(5)                                
-         amu(7)=0.101228536290376                     
-         amu(8)=amu(7)                                
-         beta=1.0d0                                   
-      end if                                         
-c!!!!!!!!!!!!!!!!!!!!!!!                        
-                                                            
-ccc                                             
-ccc  determine abundances included in the absorber amount   
-ccc                                             
-                                                            
-c first, set up the grid ready for interpolation.           
-      do i=1,nzy                                      
-         zyd(i) = dble(zy(i))                         
-      enddo                                          
-      do i=1,nl                                      
-         zld(i) = dble(zl(i))                         
-      enddo                                          
-                                                            
-c vibr. temp of the bending mode :              
-      if (isot.eq.1) call interdp (tvtbs,zyd,nzy,v626t1,zld,nl,1)  
-      if (isot.eq.2) call interdp (tvtbs,zyd,nzy,v628t1,zld,nl,1)  
-      if (isot.eq.3) call interdp (tvtbs,zyd,nzy,v636t1,zld,nl,1)  
-      if (isot.eq.4) call interdp (tvtbs,zyd,nzy,v627t1,zld,nl,1)  
-      
-c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
-c por similitud a la que se hace en cza.for     
-                                                            
-      do i=1,nzy                                      
-         if (isot.eq.5) then                          
-	    con(i) = dble( coy(i) * imrco )            
-         else                                         
-	    con(i) =  dble( co2y(i) * imr(isot) )      
-	    correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )            
-	    con(i) = con(i) * ( 1.d0 - correc )        
-         endif                                        
-c-----------------------------------------------------------------------
-c mlp & cristina. 17 july 1996                  
-c change the calculation of mr. it is used for calculating partial press
-c alpha = alpha(self,co2)*pp +alpha(n2)*(pt-pp) 
-c for an isotope, if mr is obtained by co2*imr(iso)/nt we are considerin
-c collisions with other co2 isotopes (including the major one, 626)     
-c as if they were with n2. assuming mr as co2/nt, we consider collisions
-c of type 628-626 as of 626-626 instead of as 626-n2.       
-c	  mrx(i)=con(i)/ntx(i) ! old malv             
-                                                            
-!	  mrx(i)= dble(co2x(i)/ntx(i))	! mlp & crs    
-                                                            
-c jan 98:                                       
-c esta modif de mlp implica anular el correc (deberia revisar esto)     
-         mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98  
-                                                            
-c-----------------------------------------------------------------------
-                                                            
-      end do                                         
-                                                            
-! como  beta y 1.d5 son comunes a todas las weighted absorber amounts,  
-! los simplificamos:                            
-!	coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )      
-      coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )      
-                                                            
-!	write (*,*)  ' coninf =', coninf                               
-                                                            
-ccc                                             
-ccc  temp dependence of the band strength and   
-ccc  nlte correction factor for the absorber amount         
-ccc                                             
-      call mztf_correccion ( coninf, con, ib, isot, 0 ) 
-                                                            
-ccc                                             
-ccc reads histogrammed spectral data (strength for lte and vmr=1)       
-ccc                                             
-	!hfile1 = dirspec//'hi'//dn        ! Ya no distinguimos entre d/n
-!!	hfile1 = dirspec//'hid'            ! (see why in his.for)
-!        hfile1='hid'
-!!	if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'       
-!        if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'  
-                                                            
-!	if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5     
-!     @    	.or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then      
-!	   if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat' 
-!	   if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat' 
-!	   if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat' 
-!	   if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat' 
-!	   if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat' 
-!	else                                           
-!	   if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat' 
-!	   if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat' 
-!	   if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat' 
-!	   if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat' 
-!	   if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat' 
-!	endif                                          
-	!write (*,*) ' /MZESCAPE/ hisfile: ', hisfile                         
-                                                            
-! the argument to rhist is to make this compatible with mztf_comp.f,    
-! which is a useful modification of mztf.f (to change strengths of bands
-!	call rhist (1.0)                               
-      if(ib.eq.1) then
-         if(isot.eq.1) then     !Case 1
-            mm=mm_c1
-            nbox=nbox_c1
-            tmin=tmin_c1
-            tmax=tmax_c1
-            do i=1,nbox_max
-               no(i)=no_c1(i)
-               dist(i)=dist_c1(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c1(j,i)
-                  xls1(j,i)=xls1_c1(j,i)
-                  xln1(j,i)=xln1_c1(j,i)
-                  xld1(j,i)=xld1_c1(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c1(j)
-            enddo
-         else if(isot.eq.2) then !Case 2
-            mm=mm_c2
-            nbox=nbox_c2
-            tmin=tmin_c2
-            tmax=tmax_c2
-            do i=1,nbox_max
-               no(i)=no_c2(i)
-               dist(i)=dist_c2(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c2(j,i)
-                  xls1(j,i)=xls1_c2(j,i)
-                  xln1(j,i)=xln1_c2(j,i)
-                  xld1(j,i)=xld1_c2(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c2(j)
-            enddo
-         else if(isot.eq.3) then !Case 3
-            mm=mm_c3
-            nbox=nbox_c3
-            tmin=tmin_c3
-            tmax=tmax_c3
-            do i=1,nbox_max
-               no(i)=no_c3(i)
-               dist(i)=dist_c3(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c3(j,i)
-                  xls1(j,i)=xls1_c3(j,i)
-                  xln1(j,i)=xln1_c3(j,i)
-                  xld1(j,i)=xld1_c3(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c3(j)
-            enddo
-         else if(isot.eq.4) then !Case 4
-            mm=mm_c4
-            nbox=nbox_c4
-            tmin=tmin_c4
-            tmax=tmax_c4
-            do i=1,nbox_max
-               no(i)=no_c4(i)
-               dist(i)=dist_c4(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c4(j,i)
-                  xls1(j,i)=xls1_c4(j,i)
-                  xln1(j,i)=xln1_c4(j,i)
-                  xld1(j,i)=xld1_c4(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c4(j)
-            enddo
-         else
-            write(*,*)'isot must be 2,3 or 4 for ib=1!!'
-            write(*,*)'stop at mzescape/312'
-            stop
-         endif
-      else if (ib.eq.2) then
-         if(isot.eq.1) then	!Case 5
-            mm=mm_c5
-            nbox=nbox_c5
-            tmin=tmin_c5
-            tmax=tmax_c5
-            do i=1,nbox_max
-               no(i)=no_c5(i)
-               dist(i)=dist_c5(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c5(j,i)
-                  xls1(j,i)=xls1_c5(j,i)
-                  xln1(j,i)=xln1_c5(j,i)
-                  xld1(j,i)=xld1_c5(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c5(j)
-            enddo
-         else
-            write(*,*)'isot must be 1 for ib=2!!'
-            write(*,*)'stop at mzescape/336'
-            stop
-         endif
-      else if (ib.eq.3) then
-         if(isot.eq.1) then	!Case 6
-            mm=mm_c6
-            nbox=nbox_c6
-            tmin=tmin_c6
-            tmax=tmax_c6
-            do i=1,nbox_max
-               no(i)=no_c6(i)
-               dist(i)=dist_c6(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c6(j,i)
-                  xls1(j,i)=xls1_c6(j,i)
-                  xln1(j,i)=xln1_c6(j,i)
-                  xld1(j,i)=xld1_c6(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c6(j)
-            enddo
-         else
-            write(*,*)'isot must be 1 for ib=3!!'
-            write(*,*)'stop at mzescape/360'
-            stop
-         endif
-      else if (ib.eq.4) then
-         if(isot.eq.1) then	!Case 7
-            mm=mm_c7
-            nbox=nbox_c7
-            tmin=tmin_c7
-            tmax=tmax_c7
-            do i=1,nbox_max
-               no(i)=no_c7(i)
-               dist(i)=dist_c7(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c7(j,i)
-                  xls1(j,i)=xls1_c7(j,i)
-                  xln1(j,i)=xln1_c7(j,i)
-                  xld1(j,i)=xld1_c7(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c7(j)
-            enddo
-         else
-            write(*,*)'isot must be 1 for ib=4!!'
-            write(*,*)'stop at mzescape/384'
-            stop
-         endif
-      else 
-         write(*,*)'ib must be 1,2,3 or 4!!'
-         write(*,*)'stop at mzescape/389'
-      endif                                                    
-      if (isot.ne.5) deltanux = deltanu(isot,ib)     
-      if (isot.eq.5) deltanux = deltanuco            
-      
-c******                                         
-c****** calculation of tauinf(nl)               
-c******                                         
-      call initial                                   
-                                                            
-      ff=1.0e10                                      
-                                                            
-      do i=nl,1,-1                                   
-                                                            
-         if(i.eq.nl)then                              
-                                                            
-            call intz (zl(i),c2,p2,mr2,t2, con)           
-            do kr=1,nbox                                  
-               ta(kr)=t2                                    
-            end do                                      
-!	write (*,*)  ' i, t2 =', i, t2                                 
-            call interstrength (st2,t2,ka,ta)             
-            aa = p2 * coninf * mr2 * (st2 * ff)           
-            bb = p2 * coninf * st2                        
-            cc = coninf * st2                             
-            dd = t2 * coninf * st2                        
-            do kr=1,nbox                                  
-               ccbox(kr) = coninf * ka(kr)          
-               ddbox(kr) = t2 * ccbox(kr)                  
-!		  c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5   
-               c2box(kr) = c2 * ka(kr) * dble(deltaz)      
-            end do                                        
-!		c2 = c2 * st2 * beta * dble(deltaz) * 1.d5   
-            c2 = c2 * st2 * dble(deltaz)                  
-                                                            
-         else                                         
-            call intz (zl(i),c1,p1,mr1,t1, con)           
-            do kr=1,nbox                                  
-               ta(kr)=t1                                    
-            end do                                      
-!	write (*,*)  ' i, t1 =', i, t1                                 
-            call interstrength (st1,t1,ka,ta)             
-            do kr=1,nbox                                  
-!		  c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5   
-               c1box(kr) = c1 * ka(kr) * dble(deltaz)      
-            end do                                        
-!		c1 = c1 * st1 * beta * dble(deltaz) * 1.d5   
-            c1 = c1 * st1 * dble(deltaz)                  
-            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0       
-            bb = bb + ( p1*c1 + p2*c2 ) / 2.d0            
-            cc = cc + ( c1 + c2 ) / 2.d0                  
-            ccc = ( c1 + c2 ) / 2.d0                      
-            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0            
-            do kr=1,nbox                                  
-               ccbox(kr) = ccbox(kr) + 
-     @              ( c1box(kr) + c2box(kr) )/2.d0       
-               ddbox(kr) = ddbox(kr) + 
-     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0
-            end do                                        
-                                                            
-            mr2 = mr1                                     
-            c2=c1                                         
-            do kr=1,nbox	                                 
-               c2box(kr) = c1box(kr)                       
-            end do                                        
-            t2=t1                                         
-            p2=p1                                         
-         end if                                       
-                                                            
-         pt = bb / cc                                 
-         pp = aa / (cc*ff)                            
-                                                            
-!	  ta=dd/cc                                    
-!	  tdop = ta                                   
-         ts = dd/cc                                   
-         do kr=1,nbox                          
-   	    ta(kr) = ddbox(kr) / ccbox(kr)          
-         end do                                       
-!	write (*,*)  ' i, ts =', i, ts                                 
-         call interstrength(st,ts,ka,ta)              
-!	  call intershape(alsa,alna,alda,tdop)        
-         call intershape(alsa,alna,alda,ta)           
-                                                            
-*	  ua = cc/st                                  
-                                                            
-c  	next loop calculates the eqw for an especified path ua,pp,pt,ta     
-                                                            
-         eqwmu = 0.0d0                                
-         do im = 1,iimu                               
-	    eqw=0.0d0                                  
-            do  kr=1,nbox                       
-               ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)       
-               if(ua(kr).lt.0.)write(*,*)'mzescape/480',ua(kr),
-     $              ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
-
-               call findw (ig,iirw, 0, csL,psL, Desp, wsL)                       
-               if ( i_supersat .eq. 0 ) then                 
-	          eqw=eqw+no(kr)*w                     
-               else                                          
-                  eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )           
-               endif                                         
-	    end do                                     
-	    eqwmu = eqwmu + eqw * mu(im)*amu(im)       
-         end do                                       
-                                                            
-! 	  tauinf(i) = exp( - eqwmu / dble(deltanux) )            
-         tauinf(i) = 1.d0 - eqwmu / dble(deltanux)    
-         if (tauinf(i).lt.0.d0) tauinf(i) = 0.0d0     
-                                                            
-         if (i.eq.nl) then                            
-            taustar(i) = 0.0d0                        
-         else                                         
-            taustar(i) = dble(deltanux) * (tauinf(i+1)-tauinf(i)) 
-!     ~            / ( beta * cc * 1.d5 )       
-     ~           / ( beta * ccc * 1.d5 )       
-         endif                                        
-                                                            
-      end do                    ! i continue                           
-                                                            
-                                                            
-c******                                         
-c****** calculation of tau(in,ir) for n<=r      
-c******                                         
-                                                            
-      do 1 in=1,nl-1                          
-                                                            
-         call initial                          
-                                                            
-         call intz (zl(in), c1,p1,mr1,t1, con) 
-         do kr=1,nbox                          
-            ta(kr) = t1                         
-         end do                                
-         call interstrength (st1,t1,ka,ta)     
-         do kr=1,nbox                          
-            c1box(kr) = c1 * ka(kr) * dble(deltaz)          
-         end do                                
-         c1 = c1 * st1 * dble(deltaz)          
-                                                            
-         call intz (zl(in+1), c2,p2,mr2,t2, con)           
-         do kr=1,nbox                          
-            ta(kr) = t2                         
-         end do                                
-         call interstrength (st2,t2,ka,ta)     
-         do kr=1,nbox                          
-            c2box(kr) = c2 * ka(kr) * dble(deltaz)          
-         end do                                
-         c2 = c2 * st2 * dble(deltaz)          
-                                                            
-         aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0           
-         bb = bb + ( p1*c1 + p2*c2 ) / 2.d0    
-         cc = cc + ( c1 + c2 ) / 2.d0          
-         dd = dd + ( t1*c1 + t2*c2 ) / 2.d0    
-         do kr=1,nbox                          
-            ccbox(kr) = ccbox(kr) + (c1box(kr)+c2box(kr))/2.d0          
-            ddbox(kr) = ddbox(kr) + (t1*c1box(kr)+t2*c2box(kr))/2.d0    
-         end do                                
-                                                            
-         mr1=mr2                               
-         t1=t2                                 
-         c1=c2                                 
-         p1=p2                                 
-         do kr=1,nbox                          
-            c1box(kr) = c2box(kr)               
-         end do                                
-         pt = bb / cc                          
-         pp = aa / (cc * ff)                   
-         ts = dd/cc                            
-         do kr=1,nbox                          
-            ta(kr) = ddbox(kr) / ccbox(kr)      
-         end do                                
-         call interstrength(st,ts,ka,ta)       
-         call intershape(alsa,alna,alda,ta)    
-                                                            
-         eqwmu = 0.0d0                         
-         do im = 1,iimu                        
-            eqw=0.0d0                           
-            do kr=1,nbox      
-               ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
-               if(ua(kr).lt.0.)write(*,*)'mzescape/566',ua(kr),
-     $              ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
-               
-               call findw (ig,iirw, 0, csL,psL, Desp, wsL)      
-               if ( i_supersat .eq. 0 ) then               
-                  eqw=eqw+no(kr)*w                          
-               else                                        
-                  eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )         
-               endif                                       
-            end do                              
-            eqwmu = eqwmu + eqw * mu(im)*amu(im)            
-         end do                                
-                                                            
-          tauii(in) = exp( - eqwmu / dble(deltanux) )                         
-          !write (*,*) 'i,tauii=',in,tauii(in)  
-
- 1     continue                                
-       tauii(nl) = 1.0d0
-                           
-                                                            
-c end                                           
-       return                                         
-       end   
-
-
-
-c***********************************************************************
-c     mzescape_normaliz.f                           
-c***********************************************************************
-c                                               
-c     program  for correcting some strange values and for normalizing       
-c     the atmospheric escape functions computed by mzescape_15um.f    
-c     possibilities according to istyle (see mzescape_15um.f).        
-c                                               
-                                                            
-      subroutine mzescape_normaliz ( taustar, istyle )           
-                                                            
-                                                            
-c     dic 99          malv    first version   
-c     jul 2011 malv+fgg       Adapted to LMD-MGCM
-c***********************************************************************
-                                                            
-      implicit none                                  
-      include 'nlte_paramdef.h'
-      include 'nlte_commons.h'
-                                                            
-                                                            
-c arguments                                     
-      real*8 		taustar(nl) ! o                   
-      integer         istyle    ! i            
-                                                            
-c local variables and constants                 
-      integer 	i, imaximum                           
-      real*8          maximum                        
-                                                            
-c***********************************************************************
-                                                            
-!                                               
-! correcting strange values at top, eliminating local maxima, etc...    
-!                                               
-      taustar(nl) = taustar(nl-1)                    
-                                                            
-      if ( istyle .eq. 1 ) then                      
-         imaximum = nl                                
-         maximum = taustar(nl)                        
-         do i=1,nl-1                                  
-	    if (taustar(i).gt.maximum) taustar(i) = taustar(nl)    
-         enddo                                        
-      elseif ( istyle .eq. 2 ) then                  
-         imaximum = nl                                
-         maximum = taustar(nl)                        
-         do i=nl-1,1,-1                               
-	    if (taustar(i).gt.maximum) then            
-	       maximum = taustar(i)                    
-	       imaximum = i                            
-	    endif                                      
-         enddo                                        
-         do i=imaximum,nl                             
-	    if (taustar(i).lt.maximum) taustar(i) = maximum        
-         enddo                                        
-      endif                                          
-                                                            
-!                                               
-! normalizing                                   
-!                                               
-      do i=1,nl                                      
-         taustar(i) = taustar(i) / maximum            
-      enddo                                          
-                                                            
-                                                            
-c end                                           
-      return                                         
-      end 
-
-
-
-c***********************************************************************
-c     mzescape_fb.f                           
-c***********************************************************************
-      subroutine mzescape_fb(ig)                       
-                                                            
-c     computes the escape functions of the most important 15um bands       
-c     this calls mzescape ( taustar,tauinf,tauii,  ib,isot, iirw,iimu 
-                                                            
-c     nov 99 	malv		based on cm15um_fb.f           
-c     jul 2011 malv+fgg       adapted to LMD-MGCM
-c***********************************************************************
-                                                            
-      implicit none                                  
-                                                            
-      include 'nlte_paramdef.h'
-      include 'nlte_commons.h'
-                                                            
-c local variables                               
-      integer 	i, ib, ik, istyle                     
-      integer         ig        !ADDED FOR TRACEBACK
-      real*8          tau_factor                     
-      real*8          aux(nl), aux2(nl), aux3(nl)    
-                                                            
-c***********************************************************************
-                                                            
-      call mzescape (ig,taustar21,tauinf210,tauii210,1,2
-     & ,irw_mztf,imu) 
-      call mzescape (ig,taustar31,tauinf310,tauii310,1,3
-     & ,irw_mztf,imu)
-      call mzescape (ig,taustar41,tauinf410,tauii410,1,4
-     & ,irw_mztf,imu)
-                                                            
-      istyle = 2                                     
-      call mzescape_normaliz ( taustar21, istyle )   
-      call mzescape_normaliz ( taustar31, istyle )   
-      call mzescape_normaliz ( taustar41, istyle )   
-                                                            
-                                                            
-c end                                           
-      return                                         
-      end 
-
-
-
-c***********************************************************************
-c     mzescape_fh.f  
-c***********************************************************************
-      subroutine mzescape_fh(ig)                     
-              
-c     jul 2011 malv+fgg                                
-c***********************************************************************
-                                                            
-      implicit none                                  
-                                                            
-      include 'nlte_paramdef.h'
-      include 'nlte_commons.h'
-                                                            
-c local variables                               
-      integer 	i, ib, ik, istyle
-      integer         ig        ! ADDED FOR TRACEBACK
-      real*8          tau_factor                     
-      real*8          aux(nl), aux2(nl), aux3(nl)    
-                                                            
-c***********************************************************************
-                                                            
-      call zero4v( aux, taustar12,tauinf121,tauii121, nl)
-      do ik=1,3                               
-         ib=ik+1                               
-         call mzescape ( ig,aux,aux2,aux3, ib, 1,irw_mztf,imu )         
-         tau_factor = 1.d0                     
-         if (ik.eq.1) tau_factor = dble(667.75/618.03)            
-         if (ik.eq.3) tau_factor = dble(667.75/720.806)    
-         do i=1,nl                                    
-            taustar12(i) = taustar12(i) + aux(i) * tau_factor           
-            tauinf121(i) = tauinf121(i) + aux2(i) * tau_factor          
-            tauii121(i) = tauii121(i) + aux3(i) * tau_factor            
-         enddo                                        
-      enddo                                          
-                                                            
-      istyle = 2                                     
-      call mzescape_normaliz ( taustar12, istyle )   
-                                                            
-                                                           
-                                                            
-c end                                           
-      return                                         
-      end 
-
-
-
-
-
-c***********************************************************************
-c     mztud.f                                        
-c***********************************************************************
-
-      subroutine mztud ( ig,cf,cfup,cfdw,vc,taugr, ib,isot,          
-     @     iirw,iimu,itauout,icfout,itableout )   
-            
-c     program  for calculating atmospheric transmittances       
-c     to be used in the calculation of curtis matrix coefficients           
-c     i*out = 1	output of data 
-c     i*out = 0	no output   
-c     itableout = 30  output de toda la C.M. y el VC y las poblaciones de los
-c                         estados 626(020), esta opcion nueva se añade porque 
-c                         itableout=1 saca o bien solamente de 5 en 5 capas 
-c                         o bien los elementos de C.M. desde una cierta capa 
-c                         (consultese elimin_mz1d.f que es quien lo hace); lo
-c                         de las poblaciones (020) lo hace mztf_correcion.f
-
-c     jul 2011        malv+fgg Adapted to LMD-MGCM  
-c     jan 07          malv    Add new vertical fine grid zy, similar to zx
-c     sep-oct 01      malv    update for fluxes for hb and fb, adapt to Linux
-c     nov 98          mavl    allow for overlaping in the lorentz line
-c     jan 98		malv	version for mz1d. based on curtis/mztf.for   
-c     17-jul-96	mlp&crs	change the calculation of mr.     
-c				evitar: divide por cero. anhadiendo: ff    
-c     oct-92		malv 	correct s(t) dependence for all histogr bands
-c     june-92		malv	proper lower levels for laser bands         
-c     may-92		malv 	new temperature dependence for laser bands  
-c     @    991 		malv 	boxing for the averaged absorber amount and t
-c     ?		malv	extension up to 200 km altitude in mars
-c     13-nov-86	mlp	include the temperature weighted to match
-c				the eqw in the strong doppler limit.       
-c***********************************************************************
-            
-      implicit none      
-            
-      include 'nlte_paramdef.h'
-      include 'nlte_commons.h'
-                         
-c arguments             
-      integer         ig        !ADDED FOR TRACEBACK
-      real*8  	cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o    
-      real*8		vc(nl),  taugr(nl) ! o        
-      integer		ib      ! i    
-      integer		isot    ! i  
-      integer		iirw    ! i  
-      integer		iimu    ! i  
-      integer		itauout ! i           
-      integer		icfout  ! i            
-      integer		itableout ! i          
-            
-c local variables and constants     
-      integer 	i, in, ir, im, k,j
-      integer 	nmu           
-      parameter 	(nmu = 8)  
-      real*8 		tau(nl,nl)    
-      real*8 		tauinf(nl)    
-      real*8 		con(nzy), coninf           
-      real*8 		c1, c2        
-      real*8 		t1, t2        
-      real*8 		p1, p2        
-      real*8		mr1, mr2       
-      real*8 		st1, st2      
-      real*8 		c1box(70), c2box(70)      
-      real*8		ff      ! to avoid too small numbers      
-      real*8		tvtbs(nzy)      
-      real*8 		st, beta, ts, eqwmu       
-      real*8 		mu(nmu), amu(nmu)         
+      real*8		tvtbs(nzy)
+      real*8 		st, beta, ts
+
       real*8  	zld(nl), zyd(nzy)
-      real*8 		correc        
-      real 		deltanux ! width of vib-rot band (cm-1)    
-      character	isotcode*2 
-      integer         idummy
-      real*8          Desp,wsL
-       
-c formats   
- 111  format(a1)         
- 112  format(a2)         
- 101  format(i1)         
- 202  format(i2)         
- 180  format(a80)        
- 181  format(a80)        
-c***********************************************************************
-            
-c some needed values    
-!     rl=sqrt(log(2.d0))     
-!     pi2 = 3.14159265358989d0           
-      beta = 1.8d0           
-!     beta = 1.0d0           
-      idummy = 0
-      Desp = 0.0d0
-      wsL = 0.0d0
-      
-!	write (*,*) ' MZTUD/ iirw = ', iirw
-
-
-c  esto es para que las subroutines de mztfsub calculen we  
-c  de la forma apropiada para mztf, no para fot 
-      icls=icls_mztf
-            
-c codigos para filenames            
-!     if (isot .eq. 1)  isotcode = '26'  
-!     if (isot .eq. 2)  isotcode = '28'  
-!     if (isot .eq. 3)  isotcode = '36'  
-!     if (isot .eq. 4)  isotcode = '27'  
-!     if (isot .eq. 5)  isotcode = '62'  
-!     if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5     
-!     @    	.or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then      
-!     write (ibcode1,101) ib           
-!     else       
-!     write (ibcode2,202) ib           
-!     endif      
-!     write (*,'( 30h calculating curtis matrix :  ,2x,         
-!     @    	8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot         
-            
-c integration in angle !!!!!!!!!!!!!!!!!!!!     
-c------- diffusivity approx.        
-      if (iimu.eq.1) then    
-!	  write (*,*)  ' diffusivity approx. beta = ',beta 
-         mu(1) = 1.0d0        
-         amu(1)= 1.0d0        
-c-------data for 8 points integration           
-      elseif (iimu.eq.4) then            
-         write (*,*)' 4 points for the gauss-legendre angle quadrature.'
-         mu(1)=(1.0d0+0.339981043584856)/2.0d0 	      
-         mu(2)=(1.0d0-0.339981043584856)/2.0d0 	      
-         mu(3)=(1.0d0+0.861136311594053)/2.0d0        
-         mu(4)=(1.0d0-0.861136311594053)/2.0d0 	      
-         amu(1)=0.652145154862546 	       
-         amu(2)=amu(1) 	      
-         amu(3)=0.347854845137454 	       
-         amu(4)=amu(3)        
-         beta=1.0d0           
-c-------data for 8 points integration           
-      elseif(iimu.eq.8) then             
-         write (*,*)' 8 points for the gauss-legendre angle quadrature.'
-         mu(1)=(1.0d0+0.183434642495650)/2.0d0        
-         mu(2)=(1.0d0-0.183434642495650)/2.0d0        
-         mu(3)=(1.0d0+0.525532409916329)/2.0d0        
-         mu(4)=(1.0d0-0.525532409916329)/2.0d0        
-         mu(5)=(1.0d0+0.796666477413627)/2.0d0        
-         mu(6)=(1.0d0-0.796666477413627)/2.0d0        
-         mu(7)=(1.0d0+0.960289856497536)/2.0d0        
-         mu(8)=(1.0d0-0.960289856497536)/2.0d0        
-         amu(1)=0.362683783378362         
-         amu(2)=amu(1)        
-         amu(3)=0.313706645877887         
-         amu(4)=amu(3)        
-         amu(5)=0.222381034453374         
-         amu(6)=amu(5)        
-         amu(7)=0.101228536290376         
-         amu(8)=amu(7)        
-         beta=1.0d0           
-      end if     
-c!!!!!!!!!!!!!!!!!!!!!!!            
-            
-ccc         
-ccc  determine abundances included in the absorber amount   
-ccc         
-            
-c first, set up the grid ready for interpolation.           
-      do i=1,nzy              
-         zyd(i) = dble(zy(i))             
-      enddo      
-      do i=1,nl              
-         zld(i) = dble(zl(i))             
-      enddo      
-c vibr. temp of the bending mode :  
-      if (isot.eq.1) call interdp(tvtbs,zyd,nzy, v626t1,zld,nl,1)  
-      if (isot.eq.2) call interdp(tvtbs,zyd,nzy, v628t1,zld,nl,1)  
-      if (isot.eq.3) call interdp(tvtbs,zyd,nzy, v636t1,zld,nl,1)  
-      if (isot.eq.4) call interdp(tvtbs,zyd,nzy, v627t1,zld,nl,1)  
-        !if (isot.eq.5) call interdp ( tvtbs,zxd,nz, vcot1,zld,nl, 1 )  
-        
-c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
-c por similitud a la que se hace en cza.for ; esto solo se hace para CO2    
-	!write (*,*) 'imr(isot) = ', isot, imr(isot)
-      do i=1,nzy              
-         if (isot.eq.5) then  
-	    con(i) = dble( coy(i) * imrco )            
-         else     
-	    con(i) =  dble( co2y(i) * imr(isot) )      
-	    correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )            
-	    con(i) = con(i) * ( 1.d0 - correc )       
-!	    write (*,*) ' iz, correc, co2y(i), con(i) =', 
-!     @            i,correc,co2y(i),con(i) 
-         endif    
-
-	    !-----------------------------------------------------------------
-	    ! mlp & cristina. 17 july 1996    change the calculation of mr.  
-	    ! it is used for calculating partial press
-	    !       alpha = alpha(self,co2)*pp +alpha(n2)*(pt-pp) 
-	    ! for an isotope, if mr is obtained by 
-	    !       co2*imr(iso)/nt 
-	    ! we are considerin collisions with other co2 isotopes 
-	    ! (including the major one, 626) as if they were with n2. 
-	    ! assuming mr as co2/nt, we consider collisions
-	    ! of type 628-626 as of 626-626 instead of as 626-n2.       
-	    !	  mrx(i)=con(i)/ntx(i) ! old malv 
-	    !	  mrx(i)= dble(co2x(i)/ntx(i))	! mlp & crs    
-
-	    ! jan 98:   
-	    ! esta modif de mlp implica anular el correc (deberia revisar esto)
-		      
-         mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98  
-
-	    !-----------------------------------------------------------------
-
-      end do     
-
-! como  beta y 1.d5 son comunes a todas las weighted absorber amounts,  
-! los simplificamos:    
-!	coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )      
-	!write (*,*)  ' con(nz), con(nz-1)  =', con(nz), con(nz-1)        
-      coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )      
-	!write (*,*)  ' coninf =', coninf       
-            
-ccc         
-ccc  temp dependence of the band strength and   
-ccc  nlte correction factor for the absorber amount         
-ccc         
-      call mztf_correccion ( coninf, con, ib, isot, itableout ) 
-ccc         
-ccc reads histogrammed spectral data (strength for lte and vmr=1)       
-ccc         
-	!hfile1 = dirspec//'hi'//dn      !Ya no hacemos distincion d/n en esto
-!	hfile1 = dirspec//'hid'          !(see why in his.for)
-!	hfile1='hid'
-!!	if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his' 
-!	if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his' 
-            
-!	if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5     
-!     @    	.or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then      
-!	   if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat' 
-!	   if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat' 
-!	   if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat' 
-!	   if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat' 
-!	   if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
-!	else       
-!	   if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat' 
-!	   if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat' 
-!	   if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat' 
-!	   if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat' 
-!	   if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat' 
-!	endif      
-      if(ib.eq.1) then
-         if(isot.eq.1) then     !Case 1
-            mm=mm_c1
-            nbox=nbox_c1
-            tmin=tmin_c1
-            tmax=tmax_c1
-            do i=1,nbox_max
-               no(i)=no_c1(i)
-               dist(i)=dist_c1(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c1(j,i)
-                  xls1(j,i)=xls1_c1(j,i)
-                  xln1(j,i)=xln1_c1(j,i)
-                  xld1(j,i)=xld1_c1(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c1(j)
-            enddo
-         else if(isot.eq.2) then !Case 2
-            mm=mm_c2
-            nbox=nbox_c2
-            tmin=tmin_c2
-            tmax=tmax_c2
-            do i=1,nbox_max
-               no(i)=no_c2(i)
-               dist(i)=dist_c2(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c2(j,i)
-                  xls1(j,i)=xls1_c2(j,i)
-                  xln1(j,i)=xln1_c2(j,i)
-                  xld1(j,i)=xld1_c2(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c2(j)
-            enddo
-         else if(isot.eq.3) then !Case 3
-            mm=mm_c3
-            nbox=nbox_c3
-            tmin=tmin_c3
-            tmax=tmax_c3
-            do i=1,nbox_max
-               no(i)=no_c3(i)
-               dist(i)=dist_c3(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c3(j,i)
-                  xls1(j,i)=xls1_c3(j,i)
-                  xln1(j,i)=xln1_c3(j,i)
-                  xld1(j,i)=xld1_c3(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c3(j)
-            enddo
-         else if(isot.eq.4) then !Case 4
-            mm=mm_c4
-            nbox=nbox_c4
-            tmin=tmin_c4
-            tmax=tmax_c4
-            do i=1,nbox_max
-               no(i)=no_c4(i)
-               dist(i)=dist_c4(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c4(j,i)
-                  xls1(j,i)=xls1_c4(j,i)
-                  xln1(j,i)=xln1_c4(j,i)
-                  xld1(j,i)=xld1_c4(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c4(j)
-            enddo
-         else
-            write(*,*)'isot must be 2,3 or 4 for ib=1!!'
-            write(*,*)'stop at mztud/324'
-            stop
-         endif
-      else if (ib.eq.2) then
-         if(isot.eq.1) then	!Case 5
-            mm=mm_c5
-            nbox=nbox_c5
-            tmin=tmin_c5
-            tmax=tmax_c5
-            do i=1,nbox_max
-               no(i)=no_c5(i)
-               dist(i)=dist_c5(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c5(j,i)
-                  xls1(j,i)=xls1_c5(j,i)
-                  xln1(j,i)=xln1_c5(j,i)
-                  xld1(j,i)=xld1_c5(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c5(j)
-            enddo
-         else
-            write(*,*)'isot must be 1 for ib=2!!'
-            write(*,*)'stop at mztud/348'
-            stop
-         endif
-      else if (ib.eq.3) then
-         if(isot.eq.1) then	!Case 6
-            mm=mm_c6
-            nbox=nbox_c6
-            tmin=tmin_c6
-            tmax=tmax_c6
-            do i=1,nbox_max
-               no(i)=no_c6(i)
-               dist(i)=dist_c6(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c6(j,i)
-                  xls1(j,i)=xls1_c6(j,i)
-                  xln1(j,i)=xln1_c6(j,i)
-                  xld1(j,i)=xld1_c6(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c6(j)
-            enddo
-         else
-            write(*,*)'isot must be 1 for ib=3!!'
-            write(*,*)'stop at mztud/372'
-            stop
-         endif
-      else if (ib.eq.4) then
-         if(isot.eq.1) then	!Case 7
-            mm=mm_c7
-            nbox=nbox_c7
-            tmin=tmin_c7
-            tmax=tmax_c7
-            do i=1,nbox_max
-               no(i)=no_c7(i)
-               dist(i)=dist_c7(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c7(j,i)
-                  xls1(j,i)=xls1_c7(j,i)
-                  xln1(j,i)=xln1_c7(j,i)
-                  xld1(j,i)=xld1_c7(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c7(j)
-            enddo
-         else
-            write(*,*)'isot must be 1 for ib=4!!'
-            write(*,*)'stop at mztud/396'
-            stop
-         endif
-      else 
-         write(*,*)'ib must be 1,2,3 or 4!!'
-         write(*,*)'stop at mztud/401'
-      endif
-		 
-	      
-	   
- 
-!	write (*,*) 'hisfile: ', hisfile       
-! the argument to rhist is to make this compatible with mztf_comp.f,    
-! which is a useful modification of mztf.f (to change strengths of bands
-!	call rhist (1.0)       
-      if (isot.ne.5) deltanux = deltanu(isot,ib)     
-      if (isot.eq.5) deltanux = deltanuco            
-      
-c******     
-c****** calculation of tauinf(nl)   
-c******     
-      call initial           
-      ff=1.0e10              
-            
-      do i=nl,1,-1           
-            
-         if(i.eq.nl)then      
-            
-            call intz (zl(i),c2,p2,mr2,t2, con)           
-            do kr=1,nbox          
-               ta(kr)=t2            
-            end do              
-!	write (*,*)  ' i, t2 =', i, t2         
-            call interstrength (st2,t2,ka,ta) 
-            aa = p2 * coninf * mr2 * (st2 * ff)           
-            bb = p2 * coninf * st2            
-            cc = coninf * st2     
-            dd = t2 * coninf * st2            
-            do kr=1,nbox          
-               ccbox(kr) = coninf * ka(kr)          
-               ddbox(kr) = t2 * ccbox(kr)      
-!		  c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5   
-               c2box(kr) = c2 * ka(kr) * dble(deltaz)      
-            end do    
-!		c2 = c2 * st2 * beta * dble(deltaz) * 1.d5   
-            c2 = c2 * st2 * dble(deltaz)      
-            
-         else     
-            call intz (zl(i),c1,p1,mr1,t1, con)           
-            do kr=1,nbox          
-               ta(kr)=t1            
-            end do              
-!	write (*,*)  ' i, t1 =', i, t1         
-            call interstrength (st1,t1,ka,ta) 
-            do kr=1,nbox          
-!		  c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5   
-               c1box(kr) = c1 * ka(kr) * dble(deltaz)      
-            end do    
-!		c1 = c1 * st1 * beta * dble(deltaz) * 1.d5   
-            c1 = c1 * st1 * dble(deltaz)      
-            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0       
-            bb = bb + ( p1*c1 + p2*c2 ) / 2.d0            
-            cc = cc + ( c1 + c2 ) / 2.d0      
-            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0            
-            do kr=1,nbox          
-               ccbox(kr) = ccbox(kr) + 
-     @              ( c1box(kr) + c2box(kr) )/2.d0       
-               ddbox(kr) = ddbox(kr) + 
-     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0 
-            end do    
-            
-            mr2 = mr1             
-            c2=c1     
-            do kr=1,nbox	         
-               c2box(kr) = c1box(kr)           
-            end do    
-            t2=t1     
-            p2=p1     
-         end if   
-
-         pt = bb / cc         
-         pp = aa / (cc*ff)    
-            
-!	  ta=dd/cc            
-!	  tdop = ta           
-         ts = dd/cc           
-         do kr=1,nbox  
-   	    ta(kr) = ddbox(kr) / ccbox(kr)          
-         end do   
-!	write (*,*)  ' i, ts =', i, ts         
-         call interstrength(st,ts,ka,ta)  
-!	  call intershape(alsa,alna,alda,tdop)        
-         call intershape(alsa,alna,alda,ta)           
-*	  ua = cc/st          
-
-c  	next loop calculates the eqw for an especified path uapp,pt,ta     
-            
-         eqwmu = 0.0d0        
-         do im = 1,iimu       
-	    eqw=0.0d0          
-            do  kr=1,nbox           
-               ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
-               if(ua(kr).lt.0.)write(*,*)'mztud/504',ua(kr),ccbox(kr),
-     $              ka(kr),beta,mu(im),kr,im,i,nl
-               call findw(ig,iirw, idummy,c1,p1,Desp,wsL)           
-               if ( i_supersat .eq. 0 ) then     
-	          eqw=eqw+no(kr)*w         
-               else      
-                  eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )           
-               endif     
-	    end do             
-	    eqwmu = eqwmu + eqw * mu(im)*amu(im)       
-         end do   
-	  
-         tauinf(i) = exp( - eqwmu / dble(deltanux) ) 
-            
-      end do               
-!	if ( isot.eq.1 .and. ib.eq.2 ) then           
-!		write (*,*)  ' tauinf(nl) = ', tauinf(nl)         
-!		write (*,*)  ' tauinf(1) = ', tauinf(1)           
-!	endif     
-            
-c******     
-c****** calculation of tau(in,ir) for n<=r      
-c******     
-        
-      do 1 in=1,nl-1         
-         call initial          
-         call intz (zl(in), c1,p1,mr1,t1, con)          
-         do kr=1,nbox           
-            ta(kr) = t1          
-         end do     
-         call interstrength (st1,t1,ka,ta)  
-         do kr=1,nbox           
-!     c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5    
-            c1box(kr) = c1 * ka(kr) * dble(deltaz)       
-         end do     
-!     c1 = c1 * st1 * beta * dble(deltaz) * 1.d5    
-         c1 = c1 * st1 * dble(deltaz)       
-            
-         do 2 ir=in,nl-1        
-            
-            if (ir.eq.in) then     
-               tau(in,ir) = 1.d0    
-               goto 2   
-            end if     
-            
-            call intz (zl(ir), c2,p2,mr2,t2, con)          
-            do kr=1,nbox           
-               ta(kr) = t2          
-            end do     
-            call interstrength (st2,t2,ka,ta)  
-            do kr=1,nbox           
-!	  c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5    
-               c2box(kr) = c2 * ka(kr) * dble(deltaz)       
-            end do     
-!	c2 = c2 * st2 * beta * dble(deltaz) * 1.e5    
-            c2 = c2 * st2 * dble(deltaz)       
-            
-c	aa = aa + ( p1*mr1*c1 + p2*mr2*c2 ) / 2.d0    
-            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0        
-            bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 
-            cc = cc + ( c1 + c2 ) / 2.d0       
-            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 
-            do kr=1,nbox           
-               ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0  
-               ddbox(kr) = ddbox(kr) + 
-     $              ( t1*c1box(kr) + t2*c2box(kr) ) /2.d0        
-            end do     
-            
-            mr1=mr2    
-            t1=t2      
-            c1=c2      
-            p1=p2      
-            do kr=1,nbox	          
-               c1box(kr) = c2box(kr)            
-            end do     
-            
-            pt = bb / cc           
-            pp = aa / (cc * ff)    
-            
-*	ta=dd/cc              
-*	tdop = ta             
-            ts = dd/cc             
-            do kr=1,nbox    
-               ta(kr) = ddbox(kr) / ccbox(kr)          
-            end do     
-            call interstrength(st,ts,ka,ta)    
-            call intershape(alsa,alna,alda,ta) 
-*	ua = cc/st            
-            
-c  	next loop calculates the eqw for an especified path ua,pp,pt,ta     
-            
-            eqwmu = 0.0d0          
-            do im = 1,iimu         
-               eqw=0.0d0            
-               do kr=1,nbox  
-                  ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
-
-                  call findw(ig,iirw, idummy,c1,p1,Desp,wsL)           
-                  if ( i_supersat .eq. 0 ) then     
-                     eqw=eqw+no(kr)*w         
-                  else      
-                     eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )           
-                  endif     
-               end do   
-               eqwmu = eqwmu + eqw * mu(im)*amu(im)         
-            end do     
-
-            tau(in,ir) = exp( - eqwmu / dble(deltanux) )  
-            
- 2       continue             
-            
- 1    continue             
-!	if ( isot.eq.1 .and. ib.eq.2 ) then           
-!		write (*,*)  ' tau(1,*) , *=1,20 '    
-!		write (*,*)  ( sngl(tau(1,k)), k=1,20 )           
-!	endif     
-            
-            
-c**********             
-c**********  calculation of tau(in,ir) for n>r  
-c**********             
-            
-      in=nl      
-            
-      call initial           
-      call intz (zl(in), c1,p1,mr1,t1, con)          
-      do kr=1,nbox           
-         ta(kr) = t1          
-      end do     
-      call interstrength (st1,t1,ka,ta)  
-      do kr=1,nbox           
-!	  c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5    
-         c1box(kr) = c1 * ka(kr) * dble(deltaz)       
-      end do     
-!	c1 = c1 * st1 * beta * dble(deltaz) * 1.d5    
-      c1 = c1 * st1 * dble(deltaz)       
-            
-      do 4 ir=in-1,1,-1      
-            
-         call intz (zl(ir), c2,p2,mr2,t2, con)          
-         do kr=1,nbox           
-            ta(kr) = t2          
-         end do     
-         call interstrength (st2,t2,ka,ta)  
-         do kr=1,nbox           
-!	  c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5    
-            c2box(kr) = c2 * ka(kr) * dble(deltaz)       
-         end do     
-!	c2 = c2 * st2 * beta * dble(deltaz) * 1.d5    
-         c2 = c2 * st2 * dble(deltaz)       
-            
-         aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0        
-         bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 
-         cc = cc + ( c1 + c2 ) / 2.d0       
-         dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 
-         do kr=1,nbox           
-            ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0  
-            ddbox(kr) = ddbox(kr) + 
-     $           ( t1*c1box(kr) + t2*c2box(kr) ) /2.d0        
-         end do     
-         
-         mr1=mr2    
-         c1=c2      
-         t1=t2      
-         p1=p2      
-         do kr=1,nbox           
-            c1box(kr) = c2box(kr)            
-         end do     
-         
-         pt = bb / cc           
-         pp = aa / (cc * ff)    
-         ts = dd / cc           
-         do kr=1,nbox           
-            ta(kr) = ddbox(kr) / ccbox(kr)   
-         end do     
-         call interstrength (st,ts,ka,ta)   
-         call intershape (alsa,alna,alda,ta)            
-            
-*	ua = cc/st            
-            
-c  	next loop calculates the eqw for an especified path ua,pp,pt,ta     
-            
-         eqwmu = 0.0d0          
-         do im = 1,iimu         
-            eqw=0.0d0            
-            do kr=1,nbox  
-               ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)       
-               if(ua(kr).lt.0.)write(*,*)'mztud/691',ua(kr),ccbox(kr),
-     $              ka(kr),beta,mu(im),kr,im,i,nl
-
-               call findw(ig,iirw, idummy,c1,p1,Desp,wsL)           
-               if ( i_supersat .eq. 0 ) then     
-	          eqw=eqw+no(kr)*w         
-               else      
-                  eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )           
-               endif     
-            end do   
-            eqwmu = eqwmu + eqw * mu(im)*amu(im)         
-         end do     
-            
-         tau(in,ir) = exp( - eqwmu / dble(deltanux) )  
-         
- 4    continue             
-            
-c           
-c due to the simmetry of the transmittances     
-c           
-      do in=nl-1,2,-1       
-         do ir=in-1,1,-1      
-            tau(in,ir) = tau(ir,in)           
-         end do   
-      end do     
-
-            
-ccc         
-ccc  writing out transmittances     
-ccc         
-      if (itauout.eq.1) then             
-            
-!	        if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5         
-!     @    	 .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then     
-!  	         open( 1, file=          
-!     @    	   dircurtis//'taul'//isotcode//dn//ibcode1//'.dat',     
-!     @    	   access='sequential', form='unformatted' ) 
-!	        else           
-!  	         open( 1, file=          
-!     @    	   dircurtis//'taul'//isotcode//dn//ibcode2//'.dat',     
-!     @    	   access='sequential', form='unformatted' ) 
-!	        endif          
-            
-!		write(1) dummy        
-!		write(1)' format: (tauinf(n),(tau(n,r),r=1,nl),n=1,nl)'   
-!		do in=1,nl            
-!		    write (1) tauinf(in), ( tau(in,ir), ir=1,nl )         
-!		end do    
-!		close(unit=1)         
-            
-      elseif (itauout.eq.2) then         
-  	         
-!	   if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5  
-!     @    	.or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then	     
-!	     open( 1, file=    
-!     @        dircurtis//'taul'//isotcode//dn//ibcode1//'.dat')     
-!	   else    
-!	     open( 1, file=    
-!     @        dircurtis//'taul'//isotcode//dn//ibcode2//'.dat')     
-!	   endif   
-            
-!		!write(1,*) dummy     
-!		!write(1,*) 'tij for curtis matrix calculations '         
-!		!write(1,*)' cira mars model atmosphere '     
-!		!write(1,*)' beta= ',beta,'deltanu= ',deltanux 
-!		write(1,*) nl
-!		write(1,*)
-!     @             ' format: (tauinf(in),(tau(in,ir),ir=1,nl),in=1,nl)'
-	           
-!		do in=1,nl            
-!		    write (1,*) tauinf(in)        
-!		    write (1,*) (tau(in,ir), ir=1,nl)    
-!		end do    
-!		close(unit=1)         
-            
-!	   if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5  
-!     @    	.or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then      
-!	      write (*,'(1x, 31htransmitances written out in: ,a22)')          
-!     @    	'taul'//isotcode//dn//ibcode1    
-!	   else    
-!	      write (*,'(1x, 31htransmitances written out in: ,a22)')          
-!     @    	'taul'//isotcode//dn//ibcode2    
-!	   endif   
-            
-      end if    
-            
-c cleaning of transmittances        
-!	call elimin_tau(tau,tauinf,nl,nan,itableout,nw,dummy,     
-!     @    					isotcode,dn,ibcode2)        
-            
-c construction of the curtis matrix 
-      
-      call mzcud ( tauinf,tau, cf,cfup,cfdw, vc,taugr,            
-     @     ib,isot,icfout,itableout )            
-            
-c end       
-      return     
-      end
-
-
-
-
-
-c***********************************************************************
-c     mzcud.f  
-c***********************************************************************
-                                                
-      subroutine mzcud( tauinf,tau, c,cup,cdw,vc,taugr,           
-     @     ib,isot,icfout,itableout )           
- 
-c     old times       mlp     first version of mzcf                
-c     a.k.murphy method to avoid extrapolation in the curtis matrix         
-c     feb-89 	        malv 	AKM method to avoid extrapolation in C.M.
-c     25-sept-96  cristina 	dejar las matrices en doble precision 
-c     jan 98		malv	version para mz1d                
-c     oct 01		malv	update version for fluxes for hb and fb
-c     jul 2011        malv+fgg Adapted to LMD-MGCM
-c***********************************************************************
-                                                
-      implicit none                                  
-                 
-      include 'comcstfi.h'
-      include 'nlte_paramdef.h'
-      include 'nlte_commons.h'
-
-c arguments                                     
-      real*8		c(nl,nl), cup(nl,nl), cdw(nl,nl) ! o   
-      real*8 		vc(nl), taugr(nl) ! o       
-      real*8 		tau(nl,nl) ! i                     
-      real*8		tauinf(nl) ! i                      
-      integer		ib      ! i                            
-      integer 	isot		! i                          
-      integer		icfout, itableout ! i               
-                                                
-c external                                      
-      external 	bandid                               
-      character*2 	bandid                            
-                                                
-c local variables                               
-      integer 	i, in, ir, iw, itblout                         
-      real*8		cfup(nl,nl), cfdw(nl,nl)               
-      real*8		a(nl,nl), cf(nl,nl)                    
-      character	isotcode*2, bcode*2                  
-                                                
-c formats                                       
- 101  format(i1)                                 
- 202  format(i2)                                 
- 180  format(a80)                                
- 181  format(a80)                                
-c***********************************************************************
-                                                
-      if (isot.eq.1)  isotcode = '26'               
-      if (isot.eq.2)  isotcode = '28'               
-      if (isot.eq.3)  isotcode = '36'               
-      if (isot.eq.4)  isotcode = '27'               
-      if (isot.eq.5)  isotcode = 'co'               
-      bcode = bandid( ib )                           
-                                                
-!	write (*,*)  ' '                                               
-                                                
-      do in=1,nl                                     
-                                                
-         do ir=1,nl                             
-                                                
-            cf(in,ir) = 0.0d0                     
-            cfup(in,ir) = 0.0d0                   
-            cfdw(in,ir) = 0.0d0                   
-            c(in,ir) = 0.0d0                      
-            cup(in,ir) = 0.0d0                    
-            cdw(in,ir) = 0.0d0                    
-            a(in,ir) = 0.0d0                      
-                                                
-         end do                                 
-                                                
-         vc(in) = 0.0d0                         
-         taugr(in) = 0.0d0                      
-                                                
-      end do                                  
-                                                
-                                                
-c	the next lines are a reduced and equivalent way of calculating        
-c	the c(in,ir) elements for n=2,nl1 and r=1,nl  
-                                                
-                                                
-c	do in=2,nl1                                   
-c	do ir=1,nl                                    
-c	if(ir.eq.1)then                               
-c	c(in,ir)=tau(in-1,1)-tau(in+1,1)              
-c	elseif(ir.eq.nl)then                          
-c	c(in,ir)=tau(in+1,nl1)-tauinf(in+1)-tau(in-1,nl1)+tauinf(in-1)        
-c	else                                          
-c	c(in,ir)=tau(in+1,ir-1)-tau(in+1,ir)-tau(in-1,ir-1)+tau(in-1,ir)      
-c	end if                                        
-c	c(in,ir)=c(in,ir)*pi*deltanu(ib)/(2.*deltaz*1.0e5)	       
-c	end do                                        
-c	end do	                                       
-c	go to 1000                                    
-                                                
-c calculation of the matrix cfup(nl,nl)         
-                                                
-      cfup(1,1) = 1.d0 - tau(1,1)             
-                                                
-      do in=2,nl                              
-         do ir=1,in                              
-                                                
-            if (ir.eq.1) then                       
-               cfup(in,ir) = tau(in,ir) - tau(in,1)        
-            elseif (ir.eq.in) then                  
-               cfup(in,ir) = 1.d0 - tau(in,ir-1)           
-            else                                    
-               cfup(in,ir) = tau(in,ir) - tau(in,ir-1)     
-            end if                                  
-                                                
-         end do                                  
-      end do                                  
-                                                
-! contribution to upwards fluxes from bb at bottom :        
-      do in=1,nl                              
-         taugr(in) =  tau(in,1)                
-      enddo                                   
-                                                
-c calculation of the matrix cfdw(nl,nl)         
-                                                
-      cfdw(nl,nl) = 1.d0 - tauinf(nl)         
-                                                
-      do in=1,nl-1                            
-         do ir=in,nl                             
-                                                
-            if (ir.eq.in) then                      
-               cfdw(in,ir) = 1.d0 - tau(in,ir)             
-            elseif (ir.eq.nl) then                  
-               cfdw(in,ir) = tau(in,ir-1) - tauinf(in)     
-            else                                    
-               cfdw(in,ir) = tau(in,ir-1) - tau(in,ir)     
-            end if                                  
-                                                
-         end do                                  
-      end do                                  
-                                                
-                                                
-c calculation of the matrix cf(nl,nl)           
-                                                
-      do in=1,nl                                     
-         do ir=1,nl                                     
-                                                
-            if (ir.eq.1) then                              
-	    ! version con l_bb(tg)  =  l_bb(t(1))=j(1) (see also vc below)     
-	    !	cf(in,ir) = tau(in,ir)                   
-	    ! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see also vc below)     
-               cf(in,ir) = tau(in,ir) - tau(in,1)            
-            elseif (ir.eq.nl) then                         
-               cf(in,ir) = tauinf(in) - tau(in,ir-1)         
-            else                                           
-               cf(in,ir) = tau(in,ir) - tau(in,ir-1)         
-            end if                                         
-                                                
-         end do                                         
-      end do                                         
-                                                
-                                                
-c  definition of the a(nl,nl) matrix            
-                                                
-      do in=2,nl-1                                   
-         do ir=1,nl                                      
-            if (ir.eq.in+1) a(in,ir) = -1.d0              
-            if (ir.eq.in-1) a(in,ir) = +1.d0              
-            a(in,ir) = a(in,ir) / ( 2.d0*deltaz*1.d5 )         
-         end do                                       
-      end do                                         
-! this is not needed anymore in the akm scheme  
-!	a(1,1) = +3.d0                                
-!	a(1,2) = -4.d0                                
-!	a(1,3) = +1.d0                                
-!	a(nl,nl)   = -3.d0                            
-!	a(nl,nl1) = +4.d0                             
-!	a(nl,nl2) = -1.d0                             
-                                                
-c calculation of the final curtis matrix ("reduced" by murphy's method) 
-                                                
-      if (isot.ne.5) then                            
-         do in=1,nl                                   
-            do ir=1,nl                                  
-               cf(in,ir) = cf(in,ir) * pi*deltanu(isot,ib)            
-               cfup(in,ir) = cfup(in,ir) * pi*deltanu(isot,ib)        
-               cfdw(in,ir) = cfdw(in,ir) * pi*deltanu(isot,ib)        
-            end do                                      
-            taugr(in) = taugr(in) * pi*deltanu(isot,ib) 
-         end do                                       
-      else                                           
-         do in=1,nl                                   
-            do ir=1,nl                                  
-               cf(in,ir) = cf(in,ir) * pi*deltanuco       
-            enddo                                       
-            taugr(in) = taugr(in) * pi*deltanuco        
-         enddo                                        
-      endif                                          
-                                                
-      do in=2,nl-1                                   
-                                                
-         do ir=1,nl                                   
-                                                
-	    do i=1,nl                                  
-	      ! only c contains the matrix a. matrixes cup,cdw dont because
-	      ! these two will be used for flux calculations, not  
-	      ! only for flux divergencies             
-                                                
-               c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir) 
-		! from this matrix we will extract (see below) the        
-		! nl2 x nl2 "core" for the "reduced" final curtis matrix. 
-                                                
-	    end do                                     
-	    cup(in,ir) = cfup(in,ir)                   
-	    cdw(in,ir) = cfdw(in,ir)                   
-                                                
-         end do			                                    
-	  ! version con l_bb(tg)  =  l_bb(t(1))=j(1)  (see cf above)           
-	  !vc(in) = c(in,1)                            
-	  ! version con l_bb(tg) =/= l_bb(t(1))=j(1)  (see cf above)           
-         if (isot.ne.5) then                            
-            vc(in) =  pi*deltanu(isot,ib)/( 2.d0*deltaz*1.d5 ) *     
-     @           ( tau(in-1,1) - tau(in+1,1) )         
-         else
-            vc(in) =  pi*deltanuco/( 2.d0*deltaz*1.d5 ) *     
-     @           ( tau(in-1,1) - tau(in+1,1) )         
-         endif
-                                    
-      end do			                                      
-		                                              
- 5    continue                                     
-                                                
-!	write (*,*)  'mztf/1/ c(2,*) =', (c(2,i), i=1,nl)              
-                                                
-!	call elimin_dibuja(c,nl,itableout)            
-                                                
-c ventana del smoothing de c es nw=3 y de vc es 5 (puesto en lisa):     
-c subroutine elimin_mz4(c,vc,ilayer,nl,nan,iw, nw)          
-                                                
-      iw = nan                                       
-      if (isot.eq.4)  iw = 5    ! eliminates values < 1.d-19
-      if (itableout.eq.30) then 
-         itblout = 0 
-      else 
-         itblout = itableout
-      endif
-      call elimin_mz1d (c,vc,0,iw,itblout,nw)      
-                                                
-! upper boundary condition                      
-!   j'(nl) = j'(nl1) ==> j(nl) = 2j(nl1) - j(nl2) ==>       
-      do in=2,nl-1                                   
-         c(in,nl-2) = c(in,nl-2) - c(in,nl)           
-         c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl)      
-         cup(in,nl-2) = cup(in,nl-2) - cup(in,nl)     
-         cup(in,nl-1) = cup(in,nl-1) + 2.d0*cup(in,nl)            
-         cdw(in,nl-2) = cdw(in,nl-2) - cdw(in,nl)     
-         cdw(in,nl-1) = cdw(in,nl-1) + 2.d0*cdw(in,nl)            
-      end do			                                      
-!   j(nl) = j(nl1) ==>                          
-!	do in=2,nl1                                   
-!	  c(in,nl1) = c(in,nl1) + c(in,nl)            
-!	end do			                                     
-                                                
-! 1000	continue                                 
-                                                
-
-      if (icfout.eq.1) then                          
-                                                
-!	   if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then     
-!		codmatrx = codmatrx_fb                        
-!	   else                                           
-!		codmatrx = codmatrx_hot                       
-!	   end if                                         
-!	   if (ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5    
-!     @        .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then
-!             ibcode2 = '0'//ibcode1
-!           else 
-!             write ( ibcode2, 202) ib
-!           endif
-                                                
-! 	   open ( 1, access='sequential', form='unformatted', file=           
-!     @    dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat')         
-! 	   open ( 2, access='sequential', form='unformatted', file=           
-!     @    dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat')       
-! 	   open ( 3, access='sequential', form='unformatted', file=           
-!     @    dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat')       
-! 	   open ( 4, access='sequential', form='unformatted', file=           
-!     @    dircurtis//'cflgr'//isotcode//dn//ibcode2//codmatrx//'.dat')       
-                                                
-!	    write(1) dummy                             
-!	    write(1) ' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)' 
-!	    do in=2,nl-1                               
-!	     write(1) vc(in), (c(in,ir)  , ir=2,nl-1 )		           
-!!             write (*,*) in, vc(in)
-!	    end do                                     
-                                                
-!	    write(2) dummy                             
-!	    write(2)' format: (cfup(n,r),r=1,nl), n=1,nl)'  
-!	    do in=1,nl                                 
-!	     write(2) ( cup(in,ir)  , ir=1,nl )		      
-!	    end do                                     
-                                                
-!	    write(3) dummy                             
-!	    write(3)' format: (cfdw(n,r),r=1,nl), n=1,nl)'         
-!	    do in=1,nl                                 
-!	     write(3) (cdw(in,ir)  , ir=1,nl )		       
-!	    end do                                     
-                                                
-!	    write(4) dummy   
-!	    write(4)' format: (taugr(n), n=1,nl)'         
-!	    do in=1,nl                                 
-!	     write(4) (taugr(in), ir=1,nl )		       
-!	    end do                 
-!            !write (*,*) ' Last value in file: ', taugr(nl)
-
-!	   write (*,'(1x,30hcurtis matrix written out in: ,a,a,a,a)' )
-!     @     dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat',
-!     @     dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat',
-!     @     dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat',
-!     @     dircurtis//'cflgr'//isotcode//dn//ibcode2//codmatrx//'.dat'
-                                                
-!           close (1)
-!           close (2)
-!           close (3)
-!           close (4)
-
-      else                                           
-	                                               
-	 ! write (*,*)  ' no curtis matrix output file ', char(10)     
-                                                
-      end if                                         
-
-      if (itableout.eq.30) then ! Force output of C.M. in ascii file 
-
-!	   if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then     
-!		codmatrx = codmatrx_fb                        
-!	   else                                           
-!		codmatrx = codmatrx_hot                       
-!	   end if                                         
-!	   if (ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5    
-!     @        .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then
-!             ibcode2 = '0'//ibcode1
-!           else 
-!             write ( ibcode2, 202) ib
-!           endif
-
-! 	   open (10, file=           
-!     &      dircurtis//'table'//isotcode//dn//ibcode2//codmatrx//'.dat')
-!            write(10,*) nl, ' = number of layers '
-!            write(10,*) ' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)'
-!            do in=2,nl-1
-!              write(10,*) vc(in), (c(in,ir)  , ir=2,nl-1 )
-!            enddo
-!           close (10)                              
-      endif
-                               
-c end                                           
-      return                                         
-      end  
-
-
-
-
-
-c***********************************************************************
-c     mztvc
-c***********************************************************************
-
-      subroutine mztvc ( ig,vc, ib,isot,          
-     @     iirw,iimu,itauout,icfout,itableout )   
-
-c     jul 2011 malv+fgg           
-c***********************************************************************
-            
-      implicit none      
-      
-      include 'comcstfi.h'
-      include 'nlte_paramdef.h'
-      include 'nlte_commons.h'
-
-c arguments             
-      integer         ig        ! ADDED FOR TRACEBACK
-      real*8  	cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o    
-      real*8		vc(nl),  taugr(nl) ! o        
-      integer		ib      ! i    
-      integer		isot    ! i  
-      integer		iirw    ! i  
-      integer		iimu    ! i  
-      integer		itauout ! i           
-      integer		icfout  ! i            
-      integer		itableout ! i          
-            
-c local variables and constants     
-      integer 	i, in, ir, im, k ,j         
-      integer 	nmu           
-      parameter 	(nmu = 8)  
-      real*8 		tau(nl,nl)    
-      real*8 		tauinf(nl)    
-      real*8 		con(nzy), coninf           
-      real*8 		c1, c2        
-      real*8 		t1, t2        
-      real*8 		p1, p2        
-      real*8		mr1, mr2       
-      real*8 		st1, st2      
-      real*8 		c1box(70), c2box(70)      
-      real*8		ff      ! to avoid too small numbers      
-      real*8		tvtbs(nzy)      
-      real*8 		st, beta, ts, eqwmu       
-      real*8 		mu(nmu), amu(nmu)         
-      real*8  	zld(nl), zyd(nzy)
-      real*8 		correc        
-      real 		deltanux ! width of vib-rot band (cm-1)    
-      character	isotcode*2 
-      integer         idummy
-      real*8          Desp,wsL
-      
-c     formats   
- 111  format(a1)         
- 112  format(a2)         
- 101  format(i1)         
- 202  format(i2)         
- 180  format(a80)        
- 181  format(a80)        
-c***********************************************************************
-            
-c some needed values    
-!     rl=sqrt(log(2.d0))     
-!     pi2 = 3.14159265358989d0           
-      beta = 1.8d0           
-!     beta = 1.0d0           
-      idummy = 0
-      Desp = 0.0d0
-      wsL = 0.0d0
-      
-                                !write (*,*) ' MZTUD/ iirw = ', iirw
-
-
-c  esto es para que las subroutines de mztfsub calculen we  
-c  de la forma apropiada para mztf, no para fot 
-      icls=icls_mztf         
-            
-c codigos para filenames            
-!	if (isot .eq. 1)  isotcode = '26'  
-!	if (isot .eq. 2)  isotcode = '28'  
-!	if (isot .eq. 3)  isotcode = '36'  
-!	if (isot .eq. 4)  isotcode = '27'  
-!	if (isot .eq. 5)  isotcode = '62'  
-!	if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5     
-!     @    	.or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then      
-!		write (ibcode1,101) ib           
-!	else       
-!		write (ibcode2,202) ib           
-!	endif      
-!	write (*,'( 30h calculating curtis matrix :  ,2x,         
-!     @    	8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot         
-            
-c integration in angle !!!!!!!!!!!!!!!!!!!!     
-            
-c------- diffusivity approx.        
-      if (iimu.eq.1) then    
-!	  write (*,*)  ' diffusivity approx. beta = ',beta 
-         mu(1) = 1.0d0        
-         amu(1)= 1.0d0        
-c-------data for 8 points integration           
-      elseif (iimu.eq.4) then            
-         write (*,*)' 4 points for the gauss-legendre angle quadrature.'
-         mu(1)=(1.0d0+0.339981043584856)/2.0d0 	      
-         mu(2)=(1.0d0-0.339981043584856)/2.0d0 	      
-         mu(3)=(1.0d0+0.861136311594053)/2.0d0        
-         mu(4)=(1.0d0-0.861136311594053)/2.0d0 	      
-         amu(1)=0.652145154862546 	       
-         amu(2)=amu(1) 	      
-         amu(3)=0.347854845137454 	       
-         amu(4)=amu(3)        
-         beta=1.0d0           
-c-------data for 8 points integration           
-      elseif(iimu.eq.8) then             
-         write (*,*)' 8 points for the gauss-legendre angle quadrature.'
-         mu(1)=(1.0d0+0.183434642495650)/2.0d0        
-         mu(2)=(1.0d0-0.183434642495650)/2.0d0        
-         mu(3)=(1.0d0+0.525532409916329)/2.0d0        
-         mu(4)=(1.0d0-0.525532409916329)/2.0d0        
-         mu(5)=(1.0d0+0.796666477413627)/2.0d0        
-         mu(6)=(1.0d0-0.796666477413627)/2.0d0        
-         mu(7)=(1.0d0+0.960289856497536)/2.0d0        
-         mu(8)=(1.0d0-0.960289856497536)/2.0d0        
-         amu(1)=0.362683783378362         
-         amu(2)=amu(1)        
-         amu(3)=0.313706645877887         
-         amu(4)=amu(3)        
-         amu(5)=0.222381034453374         
-         amu(6)=amu(5)        
-         amu(7)=0.101228536290376         
-         amu(8)=amu(7)        
-         beta=1.0d0           
-      end if     
-c!!!!!!!!!!!!!!!!!!!!!!!            
-            
-ccc         
-ccc  determine abundances included in the absorber amount   
-ccc         
-            
-c first, set up the grid ready for interpolation.           
-      do i=1,nzy              
-         zyd(i) = dble(zy(i))             
-      enddo      
-      do i=1,nl              
-         zld(i) = dble(zl(i))             
-      enddo      
-            
-c vibr. temp of the bending mode :  
-      if (isot.eq.1) call interdp(tvtbs,zyd,nzy, v626t1,zld,nl,1)  
-      if (isot.eq.2) call interdp(tvtbs,zyd,nzy, v628t1,zld,nl,1)  
-      if (isot.eq.3) call interdp(tvtbs,zyd,nzy, v636t1,zld,nl,1)  
-      if (isot.eq.4) call interdp(tvtbs,zyd,nzy, v627t1,zld,nl,1)  
-        !if (isot.eq.5) call interdp ( tvtbs,zxd,nz, vcot1,zld,nl, 1 )  
-            
-c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
-c por similitud a la que se hace en cza.for ; esto solo se hace para CO2    
-            
-	!write (*,*) 'imr(isot) = ', isot, imr(isot)
-      do i=1,nzy              
-         if (isot.eq.5) then  
-	    con(i) = dble( coy(i) * imrco )            
-         else     
-	    con(i) =  dble( co2y(i) * imr(isot) )      
-	    correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )            
-	    con(i) = con(i) * ( 1.d0 - correc )       
-!	    write (*,*) ' iz, correc, co2y(i), con(i) =', 
-!     @            i,correc,co2y(i),con(i) 
-         endif    
-
-	    !-----------------------------------------------------------------
-	    ! mlp & cristina. 17 july 1996    change the calculation of mr.  
-	    ! it is used for calculating partial press
-	    !       alpha = alpha(self,co2)*pp +alpha(n2)*(pt-pp) 
-	    ! for an isotope, if mr is obtained by 
-	    !       co2*imr(iso)/nt 
-	    ! we are considerin collisions with other co2 isotopes 
-	    ! (including the major one, 626) as if they were with n2. 
-	    ! assuming mr as co2/nt, we consider collisions
-	    ! of type 628-626 as of 626-626 instead of as 626-n2.       
-	    !	  mrx(i)=con(i)/ntx(i) ! old malv 
-	    !	  mrx(i)= dble(co2x(i)/ntx(i))	! mlp & crs    
-
-	    ! jan 98:   
-	    ! esta modif de mlp implica anular el correc (deberia revisar esto)
-		      
-         mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98  
-
-	    !-----------------------------------------------------------------
-
-      end do     
-            
-! como  beta y 1.d5 son comunes a todas las weighted absorber amounts,  
-! los simplificamos:    
-!	coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )      
-	!write (*,*)  ' con(nz), con(nz-1)  =', con(nz), con(nz-1)        
-      coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )      
-	!write (*,*)  ' coninf =', coninf       
-            
-ccc         
-ccc  temp dependence of the band strength and   
-ccc  nlte correction factor for the absorber amount         
-ccc         
-      call mztf_correccion ( coninf, con, ib, isot, itableout ) 
-            
-ccc         
-ccc reads histogrammed spectral data (strength for lte and vmr=1)       
-ccc         
-	!hfile1 = dirspec//'hi'//dn      !Ya no hacemos distincion d/n en esto
-!!	hfile1 = dirspec//'hid'          !(see why in his.for)
-!	hfile1='hid'
-!!	if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'       
-!	if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
-            
-!	if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5     
-!     @    	.or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then      
-!	   if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat' 
-!	   if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat' 
-!	   if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat' 
-!	   if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat' 
-!	   if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat' 
-!	else       
-!	   if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat' 
-!	   if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat' 
-!	   if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat' 
-!	   if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat' 
-!	   if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat' 
-!	endif      
-!	write (*,*) 'hisfile: ', hisfile       
-            
-! the argument to rhist is to make this compatible with mztf_comp.f,    
-! which is a useful modification of mztf.f (to change strengths of bands
-!	call rhist (1.0)       
-      if(ib.eq.1) then
-         if(isot.eq.1) then     !Case 1
-            mm=mm_c1
-            nbox=nbox_c1
-            tmin=tmin_c1
-            tmax=tmax_c1
-            do i=1,nbox_max
-               no(i)=no_c1(i)
-               dist(i)=dist_c1(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c1(j,i)
-                  xls1(j,i)=xls1_c1(j,i)
-                  xln1(j,i)=xln1_c1(j,i)
-                  xld1(j,i)=xld1_c1(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c1(j)
-            enddo
-         else if(isot.eq.2) then !Case 2
-            mm=mm_c2
-            nbox=nbox_c2
-            tmin=tmin_c2
-            tmax=tmax_c2
-            do i=1,nbox_max
-               no(i)=no_c2(i)
-               dist(i)=dist_c2(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c2(j,i)
-                  xls1(j,i)=xls1_c2(j,i)
-                  xln1(j,i)=xln1_c2(j,i)
-                  xld1(j,i)=xld1_c2(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c2(j)
-            enddo
-         else if(isot.eq.3) then !Case 3
-            mm=mm_c3
-            nbox=nbox_c3
-            tmin=tmin_c3
-            tmax=tmax_c3
-            do i=1,nbox_max
-               no(i)=no_c3(i)
-               dist(i)=dist_c3(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c3(j,i)
-                  xls1(j,i)=xls1_c3(j,i)
-                  xln1(j,i)=xln1_c3(j,i)
-                  xld1(j,i)=xld1_c3(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c3(j)
-            enddo
-         else if(isot.eq.4) then !Case 4
-            mm=mm_c4
-            nbox=nbox_c4
-            tmin=tmin_c4
-            tmax=tmax_c4
-            do i=1,nbox_max
-               no(i)=no_c4(i)
-               dist(i)=dist_c4(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c4(j,i)
-                  xls1(j,i)=xls1_c4(j,i)
-                  xln1(j,i)=xln1_c4(j,i)
-                  xld1(j,i)=xld1_c4(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c4(j)
-            enddo
-         else
-            write(*,*)'isot must be 2,3 or 4 for ib=1!!'
-            write(*,*)'stop at mztvc/310'
-            stop
-         endif
-      else if (ib.eq.2) then
-         if(isot.eq.1) then	!Case 5
-            mm=mm_c5
-            nbox=nbox_c5
-            tmin=tmin_c5
-            tmax=tmax_c5
-            do i=1,nbox_max
-               no(i)=no_c5(i)
-               dist(i)=dist_c5(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c5(j,i)
-                  xls1(j,i)=xls1_c5(j,i)
-                  xln1(j,i)=xln1_c5(j,i)
-                  xld1(j,i)=xld1_c5(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c5(j)
-            enddo
-         else
-            write(*,*)'isot must be 1 for ib=2!!'
-            write(*,*)'stop at mztvc/334'
-            stop
-         endif
-      else if (ib.eq.3) then
-         if(isot.eq.1) then	!Case 6
-            mm=mm_c6
-            nbox=nbox_c6
-            tmin=tmin_c6
-            tmax=tmax_c6
-            do i=1,nbox_max
-               no(i)=no_c6(i)
-               dist(i)=dist_c6(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c6(j,i)
-                  xls1(j,i)=xls1_c6(j,i)
-                  xln1(j,i)=xln1_c6(j,i)
-                  xld1(j,i)=xld1_c6(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c6(j)
-            enddo
-         else
-            write(*,*)'isot must be 1 for ib=3!!'
-            write(*,*)'stop at mztvc/358'
-            stop
-         endif
-      else if (ib.eq.4) then
-         if(isot.eq.1) then	!Case 7
-            mm=mm_c7
-            nbox=nbox_c7
-            tmin=tmin_c7
-            tmax=tmax_c7
-            do i=1,nbox_max
-               no(i)=no_c7(i)
-               dist(i)=dist_c7(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c7(j,i)
-                  xls1(j,i)=xls1_c7(j,i)
-                  xln1(j,i)=xln1_c7(j,i)
-                  xld1(j,i)=xld1_c7(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c7(j)
-            enddo
-         else
-            write(*,*)'isot must be 1 for ib=4!!'
-            write(*,*)'stop at mztvc/382'
-            stop
-         endif
-      else 
-         write(*,*)'ib must be 1,2,3 or 4!!'
-         write(*,*)'stop at mztvc/387'
-      endif
-            
-            
-c******     
-c****** calculation of tau(1,ir) for 1<=r      
-c******     
-      call initial           
-            
-      ff=1.0e10              
-
-      in=1         
-            
-      tau(in,1) = 1.d0 
-
-      call initial          
-      call intz (zl(in), c1,p1,mr1,t1, con)          
-      do kr=1,nbox           
-         ta(kr) = t1          
-      end do     
-      call interstrength (st1,t1,ka,ta)  
-      do kr=1,nbox           
-         c1box(kr) = c1 * ka(kr) * dble(deltaz)       
-      end do     
-      c1 = c1 * st1 * dble(deltaz)       
-      
-      do 2 ir=2,nl        
-            
-         call intz (zl(ir), c2,p2,mr2,t2, con)          
-         do kr=1,nbox           
-            ta(kr) = t2          
-         end do     
-         call interstrength (st2,t2,ka,ta)  
-         do kr=1,nbox           
-            c2box(kr) = c2 * ka(kr) * dble(deltaz)       
-         end do     
-         c2 = c2 * st2 * dble(deltaz)       
-         
-         aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0        
-         bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 
-         cc = cc + ( c1 + c2 ) / 2.d0       
-         dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 
-         do kr=1,nbox           
-            ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0  
-            ddbox(kr) = ddbox(kr) + 
-     $           ( t1*c1box(kr) + t2*c2box(kr) ) /2.d0        
-         end do     
-            
-         mr1=mr2    
-         t1=t2      
-         c1=c2      
-         p1=p2      
-         do kr=1,nbox	          
-            c1box(kr) = c2box(kr)            
-         end do     
-         
-         pt = bb / cc           
-         pp = aa / (cc * ff)    
-         
-         ts = dd/cc             
-         do kr=1,nbox    
-   	    ta(kr) = ddbox(kr) / ccbox(kr)          
-         end do     
-         call interstrength(st,ts,ka,ta)    
-         call intershape(alsa,alna,alda,ta) 
-         
-         
-         eqwmu = 0.0d0          
-         do im = 1,iimu         
-            eqw=0.0d0            
-            do kr=1,nbox  
-               ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
-               call findw(ig,iirw, idummy,c1,p1,Desp,wsL)           
-               if ( i_supersat .eq. 0 ) then     
-	          eqw=eqw+no(kr)*w         
-               else      
-                  eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )           
-               endif     
-            end do   
-            eqwmu = eqwmu + eqw * mu(im)*amu(im)         
-         end do     
-            
-         tau(in,ir) = exp( - eqwmu / dble(deltanu(isot,ib)) )  
-            
- 2    continue             
-            
-            
-           
-c           
-c due to the simmetry of the transmittances     
-c           
-      do in=nl,2,-1  
-         tau(in,1) = tau(1,in)           
-      end do           
-      
-      vc(1) = 0.0d0                         
-      vc(nl) = 0.0d0                         
-      do in=2,nl-1              ! poner aqui nl-1 luego          
-         vc(in) =  pi*deltanu(isot,ib)/( 2.d0*deltaz*1.d5 ) *     
-     @        ( tau(in-1,1) - tau(in+1,1) )         
-      end do			                                      
-      
-            
-c end       
-      return     
-      end 
-
-
-
-
-
-c***********************************************************************
-c     mztvc_626fh.F
-c***********************************************************************
-                                                            
-      subroutine mztvc_626fh(ig)
-
-c     jul 2011 malv+fgg
-c***********************************************************************
-                                                            
-      implicit none                                  
-                                                            
-!!!!!!!!!!!!!!!!!!!!!!!                         
-! common variables & constants                  
-                                                            
-      include 'nlte_paramdef.h'
-      include 'nlte_commons.h'
-
-!!!!!!!!!!!!!!!!!!!!!!!                         
-! arguments                                     
-                 
-      integer   ig              ! ADDED FOR TRACEBACK
-                                                            
-!!!!!!!!!!!!!!!!!!!!!!!                         
-! local variables                               
-                                                            
-      real*4 cdummy(nl,nl), csngl(nl,nl)             
-      
-      real*8 cax1(nl,nl), cax2(nl,nl), cax3(nl,nl)   
-      real*8 v1(nl), v2(nl), v3(nl), cm_factor, vc_factor        
-                                                            
-      integer itauout,icfout,itableout, interpol,ismooth, isngldble          
-      integer i,j,ik,ist,isot,ib,itt                 
-                                                            
-	!character 	bandcode*2
-      character 	isotcode*2
-	!character 	codmatrx_hot*5                      
-                                                            
-!!!!!!!!!!!!!!!!!!!!!!!                         
-! external functions                            
-                                                            
-      external bandid                                
-      character*2 bandid                             
-                                                            
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!               
-! subroutines called:                           
-! 	mz4sub, dmzout, readc_mz4, readcupdw, mztf   
-                                                            
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!               
-! formatos                                      
- 132  format(i2)                                 
-                                                            
-************************************************************************
-************************************************************************
-                                                            
-      isngldble = 1		! =1 --> dble precission        
-                                                            
-      fileroot = 'cfl'                               
-                                                            
-      ist = 1                                        
-      isot = 26                                      
-      write (isotcode,132) isot  
-                             
-      call zerov( vc121, nl )
-      
-      do 11, ik=1,3                                  
-                                                            
-         ib=ik+1                                      
-                                                            
-         call mztvc (ig,v1, ib, 1, irw_mztf, imu, 0,0,0 ) 
-                                                            
-         do i=1,nl                                    
-                                                            
-	    if(ik.eq.1)then                            
-               vc_factor = dble(667.75/618.03)               
-	    elseif(ik.eq.2)then                        
-               vc_factor = 1.d0                              
-	    elseif(ik.eq.3)then                        
-               vc_factor = dble(667.75/720.806)              
-	    end if                                     
-                                                            
-	    vc121(i) = vc121(i) + v1(i) * vc_factor    
-
-         end do          
-                                                            
- 11   continue                                     
-                                                            
-                                                            
-      return                                         
-      end  
-
-
-
-
-
-c***********************************************************************
-c     mztf_correccion
-c***********************************************************************
-                                                
-      subroutine mztf_correccion (coninf, con, ib, isot, icurt_pop)  
-                                                
-c including the dependence of the absort. coeff. on temp., vibr. temp., 
-c function, etc.., when neccessary. imr is already corrected in his.for 
-c we follow pg.39b-43a (l5):                    
-c  tvt1 is the vibr temp of the upper level     
-c  tvt  is the vibr temp of the transition itself           
-c  tvtbs is the vibr temp of the bending mode (used in qv)  
-c  for fundamental bands, they are not used at the moment.  
-c  for the 15 fh and sh bands, only tvt0 is used at the moment.         
-c  for the laser band, all of them are used following pg. 41a -l5- :    
-c    we need s(z) and we can read s(tk) from the histogram (also called 
-c    what we have to calculate now is the factor s(z)/s(tk) or following
-c    l5 notebook notation, s_nlte/s_lte.        
-c           s_nlte/s_lte = xfactor = xlower * xqv * xes     
-                                                
-c  icurt_pop = 30 -> Output of populations of the 0200,0220,1000 states
-c            = otro -> no output of these populations
-
-c     oct 92          malv                    
-c     jan 98		malv		version for mz1d          
-c     jul 2011        malv+fgg        adapted to LMD-MGCM
-c***********************************************************************
-                                                
-      implicit none                                  
-                                                
-      include 'nlte_paramdef.h'
-      include 'nlte_commons.h'
-                                                
-c arguments                                     
-      integer		ib, isot                              
-      integer 	icurt_pop       ! output of Fermi states population
-      real*8		con(nzy), coninf                        
-                                                
-! local variables                               
-      integer 	i                                     
-      real*8	tvt0(nzy),tvt1(nzy),tvtbs(nzy), zld(nl),zyd(nzy) 	       
-      real	xalfa, xbeta, xtv1000, xtv0200, xtv0220, xfactor      
-      real	xqv, xnu_trans, xtv_trans, xes, xlower    
-c***********************************************************************
-                                  
-      xfactor = 1.0
-
+      real*8 		correc
+      real*8  	deltanudbl
+      integer         isot
+      real*8          yy
+
+c     external function
+      external        we_clean
+      real*8          we_clean
+
+
+c     formats
+ 101  format(i1)
+c***********************************************************************
+
+c     some values
+      beta = 1.8d5
+      isot = 1
+      write (ibcode1,101) ib
+      deltanudbl = dble( deltanu(isot,ib) )
+      ff=1.0d10
+      deltazdbl = dble(deltaz)
+
+ccc   
+ccc   
+ccc   
+      do i=1,nl
+         zld(i) = dble(zl(i))
+      enddo
       do i=1,nzy
          zyd(i) = dble(zy(i))
       enddo
-      do i=1,nl                                      
-         zld(i) = dble( zl(i) )                
-      end do                                  
-                                                
-! tvtbs is the bending mode of the molecule. used in xqv.   
-      if (isot.eq.1) call interdp (tvtbs,zyd,nzy,v626t1,zld,nl,1)  
-      if (isot.eq.2) call interdp (tvtbs,zyd,nzy,v628t1,zld,nl,1)  
-      if (isot.eq.3) call interdp (tvtbs,zyd,nzy,v636t1,zld,nl,1)  
-      if (isot.eq.4) call interdp (tvtbs,zyd,nzy,v627t1,zld,nl,1)  
-      if (isot.eq.5) call interdp (tvtbs,zyd,nzy,vcot1,zld,nl,1)   
-      
-! tvt0 is the lower level of the transition. used in xlower.            
-      if (ib.eq.2 .or. ib.eq.3 .or. ib.eq.4 .or. ib.eq.15) then  
-         if (isot.eq.1) call interdp(tvt0,zyd,nzy,v626t1,zld,nl,1) 
-         if (isot.eq.2) call interdp(tvt0,zyd,nzy,v628t1,zld,nl,1) 
-         if (isot.eq.3) call interdp(tvt0,zyd,nzy,v636t1,zld,nl,1) 
-         if (isot.eq.4) call interdp(tvt0,zyd,nzy,v627t1,zld,nl,1) 
-      elseif (ib.eq.6 .or. ib.eq.8 .or. ib.eq.10     
-     @        .or. ib.eq.13 .or. ib.eq.14                  
-     @        .or. ib.eq.17 .or. ib.eq.19 .or. ib.eq.20) then          
-         if (isot.eq.1) call interdp(tvt0,zyd,nzy,v626t2,zld,nl,1)
-         if (isot.eq.2) call interdp(tvt0,zyd,nzy,v628t2,zld,nl,1)
-         if (isot.eq.3) call interdp(tvt0,zyd,nzy,v636t2,zld,nl,1)
-         if (isot.eq.4) then  
-            call interdp ( tvt0,zyd,nzy, v627t2,zld,nl, 1 )
-         endif                                        
-      else                                           
-         do i=1,nzy                                    
-	    tvt0(i) = dble( ty(i) )                    
-         end do                                       
-      end if                                         
-                                                
-c tvt is the vt of the transition. used in xes. 
-c since xes=1.0 except for the laser bands, tvt is only needed for  them
-c but it is actually calculated from the tv of the upper and lower level
-c of the transition. hence, only tvt1 remains to be read for the laser b
-c tvt1 is the upper level of the transition.    
-      if (ib.eq.13 .or. ib.eq.14) then 
-         if (isot.eq.1) call interdp(tvt1,zyd,nzy,v626t4,zld,nl,1)
-         if (isot.eq.2) call interdp(tvt1,zyd,nzy,v628t4,zld,nl,1) 
-         if (isot.eq.3) call interdp(tvt1,zyd,nzy,v636t4,zld,nl,1)
-         if (isot.eq.4) call interdp(tvt1,zyd,nzy,v627t4,zld,nl,1)
-      end if
-      
-c  here we weight the absorber amount by a factor which compensate the l
-c  value of the strength read from hitran. we use that factor in order t
-c  correct the product s*m when we later multiply those two variables.  
-                                                
-!        if ( isot.eq.1 .and. icurt_pop.eq.30 ) then 
-!           open (30, file='020populations.dat')
-!           write (30,*) ' z  tv(020) tv0200 tv0220 tv1000 '
-!        endif
-
-      do i=1,nzy                                      
-                                                
-         if (isot.eq.1) then                  
-
-	    !!! vt of the 3 levels in (020)  (see pag. 36a-sn1 for this)       
-	    xalfa = 1.d0/2.d0*exp(dble(-ee*(nu12_1000-nu(1,2))/ty(i)))      
-	    xbeta = 1.d0/2.d0*exp(dble(-ee*(nu12_0200-nu(1,2))/ty(i)))      
-	    xtv0200 = dble( - ee * nu12_0200 ) /       
-     @           ( log( xbeta/(1.d0+xalfa+xbeta) ) - 
-     @           dble(ee*nu(1,2))/tvt0(i) )    
-	    xtv0220 = dble( - ee * nu(1,2) ) /         
-     @           ( log( 1.d0/(1.d0+xalfa+xbeta) ) - 
-     @           dble(ee*nu(1,2))/tvt0(i) )     
-	    xtv1000 = dble( - ee * nu12_1000 ) /       
-     @           ( log( xalfa/(1.d0+xalfa+xbeta) ) - 
-     @           dble(ee*nu(1,2))/tvt0(i) )
-            !!! correccion 8-Nov-04 (see pag.9b-Marte4-)
-	    xtv0200 = dble( - ee * nu12_0200 /       
-     @           (log(4.*xbeta/(1.d0+xalfa+xbeta))-ee*nu(1,2)/tvt0(i)))    
-	    xtv0220 = dble( - ee * nu(1,2) /         
-     @           ( log(2./(1.d0+xalfa+xbeta)) - ee*nu(1,2)/tvt0(i) ) )     
-	    xtv1000 = dble( - ee * nu12_1000 /       
-     @           (log(4.*xalfa/(1.d0+xalfa+xbeta))-ee*nu(1,2)/tvt0(i)))
-
-!            if ( icurt_pop.eq.30 ) then 
-!               write (30,'( 1x,f7.2, 3x,f8.3, 2x,3(1x,f8.3) )')
-!     @           zx(i), tvt0(i), xtv0200, xtv0220, xtv1000 
-!            endif
-               
-	    !!! xlower and xes for the band            
-	    if (ib.eq.19) then                         
-               xlower = exp( dble(ee*elow(isot,ib)) *   
-     @              ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )        
-               xes = 1.0d0                              
-	    elseif (ib.eq.17) then                     
-               xlower = exp( dble(ee*elow(isot,ib)) *   
-     @              ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )        
-               xes = 1.0d0                              
-	    elseif (ib.eq.20) then                     
-               xlower = exp( dble(ee*elow(isot,ib)) *   
-     @              ( 1.d0/dble(ty(i))-1.d0/xtv0220 ) )        
-               xes = 1.0d0                              
-	    elseif (ib.eq.14) then                     
-               xlower = exp( dble(ee*nu12_1000) *       
-     @              ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )        
-               xnu_trans = dble( nu(1,4)-nu12_1000 )    
-               xtv_trans = xnu_trans / dble(nu(1,4)/tvt1(i)-
-     @              nu12_1000/xtv1000)  
-               xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /  
-     @              (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
-	    elseif (ib.eq.13) then                     
-               xlower = exp( dble(ee*nu12_0200) *       
-     @              ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )        
-               xnu_trans = dble(nu(1,4)-nu12_0200)      
-               xtv_trans = xnu_trans / dble(nu(1,4)/tvt1(i)-
-     @              nu12_0200/xtv0200)  
-               xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /  
-     @              (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
-	    else                                       
-               xlower = exp( dble(ee*elow(isot,ib)) *   
-     @              ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )        
-               xes = 1.0d0                              
-	    end if                                     
-	    xqv = (1.d0-exp( dble(-ee*667.3801/tvtbs(i)) )) /      
-     @           (1.d0-exp( dble(-ee*667.3801/ty(i)) )) 
-	    xfactor = xlower * xqv**2.d0 * xes         
-            
-         elseif (isot.eq.2) then                      
-            
-	    xalfa = 1.d0/2.d0* exp( dble(-ee*(nu22_1000-nu(2,2))/
-     @           ty(i)) )      
-	    xbeta = 1.d0/2.d0* exp( dble(-ee*(nu22_0200-nu(2,2))/
-     @           ty(i)) )      
-	    xtv0200 = dble( - ee * nu22_0200 ) /       
-     @           ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(2,2))/
-     @           tvt0(i) )    
-	    xtv1000 = dble( - ee * nu22_1000 ) /       
-     @           ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(2,2))/
-     @           tvt0(i) )    
-                                                
-	    if (ib.eq.14) then                         
-               xlower = exp( dble(ee*nu22_1000) *       
-     @              ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )        
-               xnu_trans = dble(nu(2,4)-nu22_1000)      
-               xtv_trans = xnu_trans / dble(nu(2,4)/tvt1(i)-nu22_1000/
-     @              xtv1000)  
-               xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /  
-     @              (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
-	    elseif (ib.eq.13) then                     
-               xlower = exp( dble(ee*nu22_0200) *       
-     @              ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )        
-               xnu_trans = dble( nu(2,4)-nu22_0200 )    
-               xtv_trans = xnu_trans / dble(nu(2,4)/tvt1(i)-nu22_0200/
-     @              xtv0200)  
-               xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /  
-     @              (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
-	    else                                       
-               xlower = exp( dble(ee*elow(isot,ib)) *   
-     @              ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )        
-               xes = 1.0d0                              
-	    end if                                     
-	    xqv = (1.d0-exp( dble(-ee*662.3734/tvtbs(i)) )) /      
-     @           (1.d0-exp( dble(-ee*662.3734/ty(i))  ))            
-	    xfactor = xlower * xqv**2.d0 * xes         
-            
-         elseif (isot.eq.3) then                      
-                                                
-	    xalfa = 1.d0/2.d0* exp( dble(-ee*(nu32_1000-nu(3,2))/
-     @           ty(i)) )      
-	    xbeta = 1.d0/2.d0* exp( dble(-ee*(nu32_0200-nu(3,2))/
-     @           ty(i)) )      
-	    xtv0200 = dble( - ee * nu32_0200 ) /       
-     @           ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(3,2))/
-     @           tvt0(i) )    
-	    xtv1000 = dble( - ee * nu32_1000 ) /       
-     @           ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(3,2))/
-     @           tvt0(i) )    
-            
-	    if (ib.eq.14) then                         
-               xlower = exp( dble(ee*nu32_1000) *       
-     @              ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )        
-               xnu_trans = dble( nu(3,4)-nu32_1000 )    
-               xtv_trans = xnu_trans / dble(nu(3,4)/tvt1(i)-nu32_1000/
-     @              xtv1000)  
-               xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /  
-     @              (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
-	    elseif (ib.eq.13) then                     
-               xlower = exp( dble(ee*nu32_0200) *       
-     @              ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )        
-               xnu_trans = dble( nu(3,4)-nu32_0200 )    
-               xtv_trans = xnu_trans / dble(nu(3,4)/tvt1(i)-nu32_0200/
-     @              xtv0200)  
-               xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /  
-     @              (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
-	    else                                       
-               xlower = exp( dble(ee*elow(isot,ib)) *   
-     @              ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )        
-               xes = 1.0d0                              
-	    end if                                     
-	    xqv = (1.d0-exp( dble(-ee*648.4784/tvtbs(i)) )) /      
-     @           (1.d0-exp( dble(-ee*648.4784/ty(i))  ))            
-	    xfactor = xlower * xqv**2.d0 * xes         
-                                                
-         elseif (isot.eq.4) then                      
-            
-	    xalfa = 1.d0/2.d0* exp( dble(-ee*(nu42_1000-nu(4,2))/
-     @           ty(i)) )      
-	    xbeta = 1.d0/2.d0* exp( dble(-ee*(nu42_0200-nu(4,2))/
-     @           ty(i)) )      
-	    xtv0200 = dble( - ee * nu42_0200 ) /       
-     @           ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(4,2))/
-     @           tvt0(i) )    
-	    xtv1000 = dble( - ee * nu42_1000 ) /       
-     @           ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(4,2))/
-     @           tvt0(i) )    
-            
-	    if (ib.eq.14) then                         
-               xlower = exp( dble(ee*nu42_1000) *       
-     @              ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )        
-               xnu_trans = dble( nu(4,4)-nu42_1000 )    
-               xtv_trans = xnu_trans / dble(nu(4,4)/tvt1(i)-nu42_1000/
-     @              xtv1000)  
-               xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /  
-     @              (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
-	    elseif (ib.eq.13) then                     
-               xlower = exp( dble(ee*nu42_0200) *       
-     $              ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )     
-               xnu_trans = dble( nu(4,4)-nu42_0200 )    
-               xtv_trans = xnu_trans / dble(nu(4,4)/tvt1(i)-nu42_0200/
-     @              xtv0200)  
-               xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) /  
-     @              (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
-	    else                                       
-               xlower = exp( dble(ee*elow(isot,ib)) *   
-     @              ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )        
-               xes = 1.0d0                              
-	    end if                                     
-	    xqv = (1.d0-exp( dble(-ee*664.7289/tvtbs(i)) )) /      
-     @           (1.d0-exp( dble(-ee*664.7289/ty(i))  ))            
-	    xfactor = xlower * xqv**2.d0 * xes         
-                                                
-         elseif (isot.eq.5 .and. ib.eq.1) then        
-            
-	    xlower = 1.d0                              
-	    xes = 1.0d0                                
-	    xqv = (1.d0-exp( dble(-ee*nuco_10/tvtbs(i)) )) /          
-     @           (1.d0-exp( dble(-ee*nuco_10/ty(i))  ))    
-	    xfactor = xlower * xqv * xes         
-            
-         end if                                       
-         
-         con(i) = con(i) * xfactor                    
-         if (i.eq.nzy) coninf = coninf * xfactor       
-                                                
-      end do                                         
-                    
-!        if ( isot.eq.1 .and. icurt_pop.eq.30 ) then 
-!           close (30)
-!        endif
-                            
-      return                                         
-      end  
-
-
-
-
-
-c***********************************************************************
-c     mztf.f                                        
-c***********************************************************************
-c                       
-c     program  for calculating atmospheric transmittances       
-c     to be used in the calculation of curtis matrix coefficients           
-            
-      subroutine mztf ( ig,cf,cfup,cfdw,vc,taugr, ib,isot,          
-     @     iirw,iimu,itauout,icfout,itableout )   
-            
-c     i*out = 1	output of data          
-c     i*out = 0	no output   
-c
-c     jul 2011        malv+fgg adapted to LMD-MGCM           
-c     nov 98          mavl    allow for overlaping in the lorentz line
-c     jan 98		malv	version for mz1d. based on curtis/mztf.for   
-c     17-jul-96	mlp&crs	change the calculation of mr.     
-c				evitar: divide por cero. anhadiendo: ff    
-c     oct-92		malv 	correct s(t) dependence for all histogr bands           
-c     june-92		malv	proper lower levels for laser bands         
-c     may-92		malv 	new temperature dependence for laser bands  
-c     @    991 		malv 	boxing for the averaged absorber amount and t            
-c     ?		malv	extension up to 200 km altitude in mars          
-c     13-nov-86	mlp	include the temperature weighted to match         
-c				the eqw in the strong doppler limit.       
-c***********************************************************************
-            
-      implicit none      
-            
+
+      call interhuntdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
+
+      do i=1,nzy
+         con(i) =  dble( co2y(i) * imr(isot) )
+         correc = 2.d0 * exp( -ee*dble(elow(isot,2))/tvtbs(i) )
+         con(i) = con(i) * ( 1.d0 - correc )
+         mr(i) = dble( co2y(i) / nty(i) )
+      end do
+
+      coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
+      call mztf_correccion ( coninf, con, ib )
+
+ccc   
+      call gethist_03 ( ib )
+
+
+c     
+c     tauinf(nl)
+c     
+      call initial
+
+      iaquiZ = nzy - 2
+      iaquiHIST = nhist / 2
+
+      do i=nl,1,-1
+
+         if(i.eq.nl)then
+
+            call intzhunt ( iaquiZ, zl(i),c2,p2,mr2,t2, con)
+            do kr=1,nbox
+               ta(kr)=t2
+            end do
+            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
+            aa = p2 * coninf * mr2 * (st2 * ff)
+            cc = coninf * st2
+            dd = t2 * coninf * st2
+            do kr=1,nbox
+               ccbox(kr) = coninf * ka(kr)
+               ddbox(kr) = t2 * ccbox(kr)
+               c2box(kr) = c2 * ka(kr) * deltazdbl
+            end do
+            c2 = c2 * st2 * deltazdbl
+
+         else
+            call intzhunt ( iaquiZ, zl(i),c1,p1,mr1,t1, con)
+            do kr=1,nbox
+               ta(kr)=t1
+            end do
+            call interstrhunt (iaquiHIST, st1,t1,ka,ta)
+            do kr=1,nbox
+               c1box(kr) = c1 * ka(kr) * deltazdbl
+            end do
+            c1 = c1 * st1 * deltazdbl
+            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
+            cc = cc + ( c1 + c2 ) / 2.d0
+            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
+            do kr=1,nbox
+               ccbox(kr) = ccbox(kr) +
+     @              ( c1box(kr) + c2box(kr) )/2.d0
+               ddbox(kr) = ddbox(kr) +
+     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0
+            end do
+
+            mr2 = mr1
+            c2=c1
+            do kr=1,nbox
+               c2box(kr) = c1box(kr)
+            end do
+            t2=t1
+            p2=p1
+         end if
+
+         pp = aa / (cc*ff)
+
+         ts = dd/cc
+         do kr=1,nbox
+   	    ta(kr) = ddbox(kr) / ccbox(kr)
+         end do
+         call interstrhunt(iaquiHIST, st,ts,ka,ta)
+         call intershphunt(iaquiHIST, alsa,alda,ta)
+
+c     
+
+         eqw = 0.0d0
+         do  kr=1,nbox
+            yy = ccbox(kr) * beta
+            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
+            eqw = eqw + no(kr)*w
+         end do
+
+         argumento = eqw / deltanudbl
+         tauinf(i) = dexp( - argumento )
+
+
+      end do                    ! i continue
+
+
+c     
+c     tau(in,ir) for n<=r
+c     
+
+      iaquiZ = 2
+      do 1 in=1,nl-1
+
+         call initial
+         call intzhunt ( iaquiZ, zl(in), c1,p1,mr1,t1, con)
+         do kr=1,nbox
+            ta(kr) = t1
+         end do
+         call interstrhunt (iaquiHIST, st1,t1,ka,ta)
+         do kr=1,nbox
+            c1box(kr) = c1 * ka(kr) * deltazdbl
+         end do
+         c1 = c1 * st1 * deltazdbl
+
+         do 2 ir=in,nl-1
+
+            if (ir.eq.in) then
+               tau(in,ir) = 1.d0
+               goto 2
+            end if
+
+            call intzhunt ( iaquiZ, zl(ir), c2,p2,mr2,t2, con)
+            do kr=1,nbox
+               ta(kr) = t2
+            end do
+            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
+            do kr=1,nbox
+               c2box(kr) = c2 * ka(kr) * deltazdbl
+            end do
+            c2 = c2 * st2 * deltazdbl
+
+            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
+            cc = cc + ( c1 + c2 ) / 2.d0
+            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
+            do kr=1,nbox
+               ccbox(kr) = ccbox(kr) + 
+     $              ( c1box(kr) + c2box(kr) ) / 2.d0
+               ddbox(kr) = ddbox(kr) + 
+     $              ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
+            end do
+
+            mr1=mr2
+            t1=t2
+            c1=c2
+            p1=p2
+            do kr=1,nbox
+               c1box(kr) = c2box(kr)
+            end do
+
+            pp = aa / (cc * ff)
+
+            ts = dd/cc
+            do kr=1,nbox
+               ta(kr) = ddbox(kr) / ccbox(kr)
+            end do
+            call interstrhunt(iaquiHIST, st,ts,ka,ta)
+            call intershphunt(iaquiHIST, alsa,alda,ta)
+
+c     
+
+            eqw = 0.0d0
+            do kr=1,nbox
+               yy = ccbox(kr) * beta
+               w = we_clean ( yy, pp, alsa(kr),alda(kr) )
+               eqw = eqw + no(kr)*w
+            end do
+
+            argumento = eqw / deltanudbl
+            tau(in,ir) = dexp( - argumento )
+
+ 2       continue
+
+ 1    continue
+
+c     
+c     tau(in,ir) for n>r
+c     
+
+      in=nl
+
+      call initial
+      iaquiZ = nzy - 2
+      call intzhunt ( iaquiZ, zl(in), c1,p1,mr1,t1, con)
+      do kr=1,nbox
+         ta(kr) = t1
+      end do
+      call interstrhunt (iaquiHIST, st1,t1,ka,ta)
+      do kr=1,nbox
+         c1box(kr) = c1 * ka(kr) * deltazdbl
+      end do
+      c1 = c1 * st1 * deltazdbl
+
+      do 4 ir=in-1,1,-1
+
+         call intzhunt ( iaquiZ, zl(ir), c2,p2,mr2,t2, con)
+         do kr=1,nbox
+            ta(kr) = t2
+         end do
+         call interstrhunt (iaquiHIST, st2,t2,ka,ta)
+         do kr=1,nbox
+            c2box(kr) = c2 * ka(kr) * deltazdbl
+         end do
+         c2 = c2 * st2 * deltazdbl
+
+         aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
+         cc = cc + ( c1 + c2 ) / 2.d0
+         dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
+         do kr=1,nbox
+            ccbox(kr) = ccbox(kr) + 
+     $           ( c1box(kr) + c2box(kr) ) / 2.d0
+            ddbox(kr) = ddbox(kr) + 
+     $           ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
+         end do
+
+         mr1=mr2
+         c1=c2
+         t1=t2
+         p1=p2
+         do kr=1,nbox
+            c1box(kr) = c2box(kr)
+         end do
+
+         pp = aa / (cc * ff)
+         ts = dd / cc
+         do kr=1,nbox
+            ta(kr) = ddbox(kr) / ccbox(kr)
+         end do
+         call interstrhunt (iaquiHIST, st,ts,ka,ta)
+         call intershphunt (iaquiHIST, alsa,alda,ta)
+
+c     
+         eqw=0.0d0
+         do kr=1,nbox
+            yy = ccbox(kr) * beta
+            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
+            eqw = eqw + no(kr)*w
+         end do
+
+         argumento = eqw / deltanudbl
+         tau(in,ir) = dexp( - argumento )
+
+ 4    continue
+
+c     
+c     
+c     
+      do in=nl-1,2,-1
+         do ir=in-1,1,-1
+            tau(in,ir) = tau(ir,in)
+         end do
+      end do
+
+c     
+      call MZCUD121 ( tauinf,tau, cf, vc, ib )
+
+
+c     end
+      return
+      end
+
+
+
+c     *** Old MZCUD121_dlvr11.f ***
+
+c***********************************************************************
+
+      subroutine MZCUD121 ( tauinf,tau, c,vc, ib )
+
+c***********************************************************************
+
+      implicit none
+
       include 'nlte_paramdef.h'
       include 'nlte_commons.h'
-            
-            
-c arguments             
-      integer         ig        !ADDED FOR TRACEBACK
-      real*8  	cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o.          
-      real*8		vc(nl),  taugr(nl) ! o        
-      integer		ib      ! i    
-      integer		isot    ! i  
-      integer		iirw    ! i  
-      integer		iimu    ! i  
-      integer		itauout ! i           
-      integer		icfout  ! i            
-      integer		itableout ! i          
-            
-c local variables and constants     
-      integer 	i, in, ir, im, k ,j         
-      integer 	nmu           
-      parameter 	(nmu = 8)  
-      real*8 		tau(nl,nl)    
-      real*8 		tauinf(nl)    
-      real*8 		con(nzy), coninf           
-      real*8 		c1, c2        
-      real*8 		t1, t2        
-      real*8 		p1, p2        
-      real*8		mr1, mr2       
-      real*8 		st1, st2      
-      real*8 		c1box(70), c2box(70)      
-      real*8		ff      ! to avoid too small numbers      
-      real*8		tvtbs(nzy)      
-      real*8 		st, beta, ts, eqwmu       
-      real*8 		mu(nmu), amu(nmu)         
-      real*8  	zld(nl), zyd(nzy)      
-      real*8 		correc        
-      real 		deltanux ! width of vib-rot band (cm-1) 
-!	character	isotcode*2
-      integer         idummy
-      real*8          Desp,wsL
-       
-c formats   
-! 111	format(a1)         
-! 112	format(a2)         
- 101  format(i1)         
- 202  format(i2)         
-! 180	format(a80)        
-! 181	format(a80)        
-c***********************************************************************
-            
-c some needed values    
-!	rl=sqrt(log(2.d0))     
-!	pi2 = 3.14159265358989d0           
-      beta = 1.8d0           
-      idummy = 0
-      Desp = 0.d0
-      wsL = 0.d0
-
-c  esto es para que las subroutines de mztfsub calculen we  
-c  de la forma apropiada para mztf, no para fot 
-      icls=icls_mztf         
-            
-c codigos para filenames            
-!	if (isot .eq. 1)  isotcode = '26'  
-!	if (isot .eq. 2)  isotcode = '28'  
-!	if (isot .eq. 3)  isotcode = '36'  
-!	if (isot .eq. 4)  isotcode = '27'  
-!	if (isot .eq. 5)  isotcode = '62'  
-!	if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5     
-!     @    	.or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then      
-!		write (ibcode1,101) ib           
-!	else       
-!		write (ibcode2,202) ib           
-!	endif      
-!	write (*,'( 30h calculating curtis matrix :  ,2x,         
-!     @    	8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot         
-            
-c integration in angle !!!!!!!!!!!!!!!!!!!!     
-            
-c------- diffusivity approx.        
-      if (iimu.eq.1) then    
-!	  write (*,*)  ' diffusivity approx. beta = ',beta 
-         mu(1) = 1.0d0        
-         amu(1)= 1.0d0        
-c-------data for 8 points integration           
-      elseif (iimu.eq.4) then            
-         write (*,*)' 4 points for the gauss-legendre angle quadrature.'
-         mu(1)=(1.0d0+0.339981043584856)/2.0d0 	      
-         mu(2)=(1.0d0-0.339981043584856)/2.0d0 	      
-         mu(3)=(1.0d0+0.861136311594053)/2.0d0        
-         mu(4)=(1.0d0-0.861136311594053)/2.0d0 	      
-         amu(1)=0.652145154862546 	       
-         amu(2)=amu(1) 	      
-         amu(3)=0.347854845137454 	       
-         amu(4)=amu(3)        
-         beta=1.0d0           
-c-------data for 8 points integration           
-      elseif(iimu.eq.8) then             
-         write (*,*)' 8 points for the gauss-legendre angle quadrature.'
-         mu(1)=(1.0d0+0.183434642495650)/2.0d0        
-         mu(2)=(1.0d0-0.183434642495650)/2.0d0        
-         mu(3)=(1.0d0+0.525532409916329)/2.0d0        
-         mu(4)=(1.0d0-0.525532409916329)/2.0d0        
-         mu(5)=(1.0d0+0.796666477413627)/2.0d0        
-         mu(6)=(1.0d0-0.796666477413627)/2.0d0        
-         mu(7)=(1.0d0+0.960289856497536)/2.0d0        
-         mu(8)=(1.0d0-0.960289856497536)/2.0d0        
-         amu(1)=0.362683783378362         
-         amu(2)=amu(1)        
-         amu(3)=0.313706645877887         
-         amu(4)=amu(3)        
-         amu(5)=0.222381034453374         
-         amu(6)=amu(5)        
-         amu(7)=0.101228536290376         
-         amu(8)=amu(7)        
-         beta=1.0d0           
-      end if     
-c!!!!!!!!!!!!!!!!!!!!!!!            
-            
-ccc         
-ccc  determine abundances included in the absorber amount   
-ccc         
-            
-c first, set up the grid ready for interpolation.           
-      do i=1,nzy              
-         zyd(i) = dble(zy(i))             
-      enddo      
-      do i=1,nl              
-         zld(i) = dble(zl(i))             
-      enddo      
-            
-            
-c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
-c por similitud a la que se hace en cza.for     
-            
-      do i=1,nzy              
-         if (isot.eq.5) then  
-	    con(i) = dble( coy(i) * imrco )            
-         else     
-	    con(i) =  dble( co2y(i) * imr(isot) )      
-c vibr. temp of the bending mode :  
-            if(isot.eq.1) call interdp(tvtbs,zyd,nzy,v626t1,zld,nl,1)  
-            if(isot.eq.2) call interdp(tvtbs,zyd,nzy,v628t1,zld,nl,1)  
-            if(isot.eq.3) call interdp(tvtbs,zyd,nzy,v636t1,zld,nl,1)  
-            if(isot.eq.4) call interdp(tvtbs,zyd,nzy,v627t1,zld,nl,1)  
-	    correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )            
-	    con(i) = con(i) * ( 1.d0 - correc )        
-         endif    
-c-----------------------------------------------------------------------
-c mlp & cristina. 17 july 1996      
-c change the calculation of mr. it is used for calculating partial press
-c alpha = alpha(self,co2)*pp +alpha(n2)*(pt-pp) 
-c for an isotope, if mr is obtained by co2*imr(iso)/nt we are considerin
-c collisions with other co2 isotopes (including the major one, 626)     
-c as if they were with n2. assuming mr as co2/nt, we consider collisions
-c of type 628-626 as of 626-626 instead of as 626-n2.       
-c	  mrx(i)=con(i)/ntx(i) ! old malv 
-            
-!	  mrx(i)= dble(co2x(i)/ntx(i))	! mlp & crs    
-            
-c jan 98:   
-c esta modif de mlp implica anular el correc (deberia revisar esto)     
-         mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98  
-            
-c-----------------------------------------------------------------------
-            
-      end do     
-            
-! como  beta y 1.d5 son comunes a todas las weighted absorber amounts,  
-! los simplificamos:    
-!	coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )      
-      coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )      
-            
-!	write (*,*)  ' coninf =', coninf       
-            
-ccc         
-ccc  temp dependence of the band strength and   
-ccc  nlte correction factor for the absorber amount         
-ccc         
-      call mztf_correccion ( coninf, con, ib, isot, itableout ) 
-            
-ccc         
-ccc reads histogrammed spectral data (strength for lte and vmr=1)       
-ccc         
-	!hfile1 = dirspec//'hi'//dn   ! ya no distinguimos entre d/n     
-!!	hfile1 = dirspec//'hid'       ! (see why in his.for)
-!        hfile='hid'
-!!	if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
-!        if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
-!            
-!	if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5     
-!     @    	.or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then      
-!	   if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat' 
-!	   if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat' 
-!	   if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat' 
-!	   if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat' 
-!	   if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat' 
-!	else       
-!	   if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat' 
-!	   if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat' 
-!	   if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat' 
-!	   if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat' 
-!	   if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat' 
-!	endif      
-!	write (*,*) 'hisfile: ', hisfile       
-            
-! the argument to rhist is to make this compatible with mztf_comp.f,    
-! which is a useful modification of mztf.f (to change strengths of bands
-!	call rhist (1.0)       
-      if(ib.eq.1) then
-         if(isot.eq.1) then     !Case 1
-            mm=mm_c1
-            nbox=nbox_c1
-            tmin=tmin_c1
-            tmax=tmax_c1
-            do i=1,nbox_max
-               no(i)=no_c1(i)
-               dist(i)=dist_c1(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c1(j,i)
-                  xls1(j,i)=xls1_c1(j,i)
-                  xln1(j,i)=xln1_c1(j,i)
-                  xld1(j,i)=xld1_c1(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c1(j)
-            enddo
-         else if(isot.eq.2) then !Case 2
-            mm=mm_c2
-            nbox=nbox_c2
-            tmin=tmin_c2
-            tmax=tmax_c2
-            do i=1,nbox_max
-               no(i)=no_c2(i)
-               dist(i)=dist_c2(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c2(j,i)
-                  xls1(j,i)=xls1_c2(j,i)
-                  xln1(j,i)=xln1_c2(j,i)
-                  xld1(j,i)=xld1_c2(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c2(j)
-            enddo
-         else if(isot.eq.3) then !Case 3
-            mm=mm_c3
-            nbox=nbox_c3
-            tmin=tmin_c3
-            tmax=tmax_c3
-            do i=1,nbox_max
-               no(i)=no_c3(i)
-               dist(i)=dist_c3(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c3(j,i)
-                  xls1(j,i)=xls1_c3(j,i)
-                  xln1(j,i)=xln1_c3(j,i)
-                  xld1(j,i)=xld1_c3(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c3(j)
-            enddo
-         else if(isot.eq.4) then !Case 4
-            mm=mm_c4
-            nbox=nbox_c4
-            tmin=tmin_c4
-            tmax=tmax_c4
-            do i=1,nbox_max
-               no(i)=no_c4(i)
-               dist(i)=dist_c4(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c4(j,i)
-                  xls1(j,i)=xls1_c4(j,i)
-                  xln1(j,i)=xln1_c4(j,i)
-                  xld1(j,i)=xld1_c4(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c4(j)
-            enddo
-         else
-            write(*,*)'isot must be 2,3 or 4 for ib=1!!'
-            write(*,*)'stop at mztf_overlap/317'
-            stop
-         endif
-      else if (ib.eq.2) then
-         if(isot.eq.1) then	!Case 5
-            mm=mm_c5
-            nbox=nbox_c5
-            tmin=tmin_c5
-            tmax=tmax_c5
-            do i=1,nbox_max
-               no(i)=no_c5(i)
-               dist(i)=dist_c5(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c5(j,i)
-                  xls1(j,i)=xls1_c5(j,i)
-                  xln1(j,i)=xln1_c5(j,i)
-                  xld1(j,i)=xld1_c5(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c5(j)
-            enddo
-         else
-            write(*,*)'isot must be 1 for ib=2!!'
-            write(*,*)'stop at mztf_overlap/341'
-            stop
-         endif
-      else if (ib.eq.3) then
-         if(isot.eq.1) then	!Case 6
-            mm=mm_c6
-            nbox=nbox_c6
-            tmin=tmin_c6
-            tmax=tmax_c6
-            do i=1,nbox_max
-               no(i)=no_c6(i)
-               dist(i)=dist_c6(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c6(j,i)
-                  xls1(j,i)=xls1_c6(j,i)
-                  xln1(j,i)=xln1_c6(j,i)
-                  xld1(j,i)=xld1_c6(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c6(j)
-            enddo
-         else
-            write(*,*)'isot must be 1 for ib=3!!'
-            write(*,*)'stop at mztf_overlap/365'
-            stop
-         endif
-      else if (ib.eq.4) then
-         if(isot.eq.1) then	!Case 7
-            mm=mm_c7
-            nbox=nbox_c7
-            tmin=tmin_c7
-            tmax=tmax_c7
-            do i=1,nbox_max
-               no(i)=no_c7(i)
-               dist(i)=dist_c7(i)
-               do j=1,nhist
-                  sk1(j,i)=sk1_c7(j,i)
-                  xls1(j,i)=xls1_c7(j,i)
-                  xln1(j,i)=xln1_c7(j,i)
-                  xld1(j,i)=xld1_c7(j,i)
-               enddo
-            enddo
-            do j=1,nhist 
-               thist(j)=thist_c7(j)
-            enddo
-         else
-            write(*,*)'isot must be 1 for ib=4!!'
-            write(*,*)'stop at mztf_overlap/389'
-            stop
-         endif
-      else 
-         write(*,*)'ib must be 1,2,3 or 4!!'
-         write(*,*)'stop at mztf_overlap/394'
-      endif
-      
-      if (isot.ne.5) deltanux = deltanu(isot,ib)     
-      if (isot.eq.5) deltanux = deltanuco            
-            
-c******     
-c****** calculation of tauinf(nl)   
-c******     
-      call initial           
-            
-      ff=1.0e10              
-            
-      do i=nl,1,-1           
-            
-         if(i.eq.nl)then      
-            
-            call intz (zl(i),c2,p2,mr2,t2, con)           
-            do kr=1,nbox          
-               ta(kr)=t2            
-            end do              
-!     write (*,*)  ' i, t2 =', i, t2         
-            call interstrength (st2,t2,ka,ta) 
-            aa = p2 * coninf * mr2 * (st2 * ff)           
-            bb = p2 * coninf * st2            
-            cc = coninf * st2     
-            dd = t2 * coninf * st2            
-            do kr=1,nbox          
-               ccbox(kr) = coninf * ka(kr)          
-               ddbox(kr) = t2 * ccbox(kr)      
-!		  c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5   
-               c2box(kr) = c2 * ka(kr) * dble(deltaz)      
-            end do    
-!		c2 = c2 * st2 * beta * dble(deltaz) * 1.d5   
-            c2 = c2 * st2 * dble(deltaz)      
-            
-         else     
-            call intz (zl(i),c1,p1,mr1,t1, con)           
-            do kr=1,nbox          
-               ta(kr)=t1            
-            end do              
-!	write (*,*)  ' i, t1 =', i, t1         
-            call interstrength (st1,t1,ka,ta) 
-            do kr=1,nbox          
-!     c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5   
-               c1box(kr) = c1 * ka(kr) * dble(deltaz)      
-            end do    
-!		c1 = c1 * st1 * beta * dble(deltaz) * 1.d5   
-            c1 = c1 * st1 * dble(deltaz)      
-            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0       
-            bb = bb + ( p1*c1 + p2*c2 ) / 2.d0            
-            cc = cc + ( c1 + c2 ) / 2.d0      
-            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0            
-            do kr=1,nbox          
-               ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) )/2.d0       
-               ddbox(kr) = ddbox(kr) + (t1*c1box(kr)+t2*c2box(kr))/2.d0
-            end do    
-            
-            mr2 = mr1             
-            c2=c1     
-            do kr=1,nbox	         
-               c2box(kr) = c1box(kr)           
-            end do    
-            t2=t1     
-            p2=p1     
-         end if   
-         
-         pt = bb / cc         
-         pp = aa / (cc*ff)    
-         
-!	  ta=dd/cc            
-!	  tdop = ta           
-         ts = dd/cc           
-         do kr=1,nbox  
-   	    ta(kr) = ddbox(kr) / ccbox(kr)          
-         end do   
-!	write (*,*)  ' i, ts =', i, ts         
-         call interstrength(st,ts,ka,ta)  
-!	  call intershape(alsa,alna,alda,tdop)        
-         call intershape(alsa,alna,alda,ta)           
-            
-*	  ua = cc/st          
-            
-c  	next loop calculates the eqw for an especified path ua,pp,pt,ta     
-            
-         eqwmu = 0.0d0        
-         do im = 1,iimu       
-	    eqw=0.0d0          
-            do  kr=1,nbox           
-               ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)  
-               if(ua(kr).lt.0.)write(*,*)'mztf_overlap/483',ua(kr),
-     $              ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
-               
-               call findw(ig,iirw, idummy,c1,p1,Desp,wsL)           
-               if ( i_supersat .eq. 0 ) then     
-	          eqw=eqw+no(kr)*w         
-               else      
-                  eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )           
-               endif     
-	    end do             
-	    eqwmu = eqwmu + eqw * mu(im)*amu(im)       
-         end do   
-            
-         tauinf(i) = exp( - eqwmu / dble(deltanux) ) 
-            
-      end do                    ! i continue   
-            
-!	if ( isot.eq.1 .and. ib.eq.2 ) then           
-!		write (*,*)  ' tauinf(nl) = ', tauinf(nl)         
-!		write (*,*)  ' tauinf(1) = ', tauinf(1)           
-!	endif     
-            
-c******     
-c****** calculation of tau(in,ir) for n<=r      
-c******     
-            
-      do 1 in=1,nl-1         
-            
-         call initial          
-         call intz (zl(in), c1,p1,mr1,t1, con)          
-         do kr=1,nbox           
-            ta(kr) = t1          
-         end do     
-         call interstrength (st1,t1,ka,ta)  
-         do kr=1,nbox           
-!	  c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5    
-            c1box(kr) = c1 * ka(kr) * dble(deltaz)       
-         end do     
-!     c1 = c1 * st1 * beta * dble(deltaz) * 1.d5    
-         c1 = c1 * st1 * dble(deltaz)       
-            
-         do 2 ir=in,nl-1        
-            
-            if (ir.eq.in) then     
-               tau(in,ir) = 1.d0    
-               goto 2   
-            end if     
-            
-            call intz (zl(ir), c2,p2,mr2,t2, con)          
-            do kr=1,nbox           
-               ta(kr) = t2          
-            end do     
-            call interstrength (st2,t2,ka,ta)  
-            do kr=1,nbox           
-!	  c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5    
-               c2box(kr) = c2 * ka(kr) * dble(deltaz)       
-            end do     
-!	c2 = c2 * st2 * beta * dble(deltaz) * 1.e5    
-            c2 = c2 * st2 * dble(deltaz)       
-            
-c	aa = aa + ( p1*mr1*c1 + p2*mr2*c2 ) / 2.d0    
-            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0        
-            bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 
-            cc = cc + ( c1 + c2 ) / 2.d0       
-            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 
-            do kr=1,nbox           
-               ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0  
-               ddbox(kr) = ddbox(kr) + 
-     $              ( t1*c1box(kr) + t2*c2box(kr) ) /2.d0        
-            end do     
-            
-            mr1=mr2    
-            t1=t2      
-            c1=c2      
-            p1=p2      
-            do kr=1,nbox	          
-               c1box(kr) = c2box(kr)            
-            end do     
-            
-            pt = bb / cc           
-            pp = aa / (cc * ff)    
-            
-*	ta=dd/cc              
-*	tdop = ta             
-            ts = dd/cc             
-            do kr=1,nbox    
-               ta(kr) = ddbox(kr) / ccbox(kr)          
-            end do     
-            call interstrength(st,ts,ka,ta)    
-            call intershape(alsa,alna,alda,ta) 
-*     ua = cc/st            
-            
-c  	next loop calculates the eqw for an especified path ua,pp,pt,ta     
-            
-            eqwmu = 0.0d0          
-            do im = 1,iimu         
-               eqw=0.0d0            
-               do kr=1,nbox  
-                  ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)       
-                  if(ua(kr).lt.0.)write(*,*)'mztf_overlap/581',ua(kr),
-     $                 ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
-
-                  call findw(ig,iirw, idummy,c1,p1,Desp,wsL)           
-                  if ( i_supersat .eq. 0 ) then     
-                     eqw=eqw+no(kr)*w         
-                  else      
-                     eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )           
-                  endif     
-               end do   
-               eqwmu = eqwmu + eqw * mu(im)*amu(im)         
-            end do     
-            
-            tau(in,ir) = exp( - eqwmu / dble(deltanux) )  
-            
- 2       continue             
-            
- 1    continue             
-            
-!	if ( isot.eq.1 .and. ib.eq.2 ) then           
-!		write (*,*)  ' tau(1,*) , *=1,20 '    
-!		write (*,*)  ( sngl(tau(1,k)), k=1,20 )           
-!	endif     
-            
-            
-c**********             
-c**********  calculation of tau(in,ir) for n>r  
-c**********             
-            
-      in=nl      
-            
-      call initial           
-      call intz (zl(in), c1,p1,mr1,t1, con)          
-      do kr=1,nbox           
-         ta(kr) = t1          
-      end do     
-      call interstrength (st1,t1,ka,ta)  
-      do kr=1,nbox           
-!     c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5    
-         c1box(kr) = c1 * ka(kr) * dble(deltaz)       
-      end do     
-!     c1 = c1 * st1 * beta * dble(deltaz) * 1.d5    
-      c1 = c1 * st1 * dble(deltaz)       
-            
-      do 4 ir=in-1,1,-1      
-            
-         call intz (zl(ir), c2,p2,mr2,t2, con)          
-         do kr=1,nbox           
-            ta(kr) = t2          
-         end do     
-         call interstrength (st2,t2,ka,ta)  
-         do kr=1,nbox           
-!     c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5    
-            c2box(kr) = c2 * ka(kr) * dble(deltaz)       
-         end do     
-!	c2 = c2 * st2 * beta * dble(deltaz) * 1.d5    
-         c2 = c2 * st2 * dble(deltaz)       
-            
-         aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0        
-         bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 
-         cc = cc + ( c1 + c2 ) / 2.d0       
-         dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 
-         do kr=1,nbox           
-            ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0  
-            ddbox(kr) = ddbox(kr) + ( t1*c1box(kr) + t2*c2box(kr) )/2.d0        
-         end do     
-            
-         mr1=mr2    
-         c1=c2      
-         t1=t2      
-         p1=p2      
-         do kr=1,nbox           
-            c1box(kr) = c2box(kr)            
-         end do     
-            
-         pt = bb / cc           
-         pp = aa / (cc * ff)    
-         ts = dd / cc           
-         do kr=1,nbox           
-            ta(kr) = ddbox(kr) / ccbox(kr)   
-         end do     
-         call interstrength (st,ts,ka,ta)   
-         call intershape (alsa,alna,alda,ta)            
-            
-*	ua = cc/st            
-            
-c  	next loop calculates the eqw for an especified path ua,pp,pt,ta     
-            
-         eqwmu = 0.0d0          
-         do im = 1,iimu         
-            eqw=0.0d0            
-            do kr=1,nbox  
-               ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
-               if(ua(kr).lt.0.)write(*,*)'mztf_overlap/674',ua(kr),
-     $              ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl
-               
-               call findw(ig,iirw, idummy,c1,p1,Desp,wsL)           
-               if ( i_supersat .eq. 0 ) then     
-	          eqw=eqw+no(kr)*w         
-               else      
-                  eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )           
-               endif     
-            end do   
-            eqwmu = eqwmu + eqw * mu(im)*amu(im)         
-         end do     
-            
-         tau(in,ir) = exp( - eqwmu / dble(deltanux) )  
-            
- 4    continue             
-            
-c           
-c due to the simmetry of the transmittances     
-c           
-      do in=nl-1,2,-1       
-         do ir=in-1,1,-1      
-            tau(in,ir) = tau(ir,in)           
-         end do   
-      end do     
-            
-            
-ccc         
-ccc  writing out transmittances     
-ccc         
-      if (itauout.eq.1) then             
-            
-!	        if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5         
-!     @    	 .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then     
-!  	         open( 1, file=          
-!     @    	   dircurtis//'taul'//isotcode//dn//ibcode1//'.dat',     
-!     @    	   access='sequential', form='unformatted' ) 
-!	        else           
-!  	         open( 1, file=          
-!     @    	   dircurtis//'taul'//isotcode//dn//ibcode2//'.dat',     
-!     @    	   access='sequential', form='unformatted' ) 
-!	        endif          
-            
-!		write(1) dummy        
-!		write(1)' format: (tauinf(n),(tau(n,r),r=1,nl),n=1,nl)'   
-!		do in=1,nl            
-!		    write (1) tauinf(in), ( tau(in,ir), ir=1,nl )         
-!		end do    
-!		close(unit=1)         
-            
-      elseif (itauout.eq.2) then         
-  	         
-!	   if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5  
-!     @    	.or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then	     
-!	     open( 1, file=    
-!     @        dircurtis//'taul'//isotcode//dn//ibcode1//'.dat')     
-!	   else    
-!	     open( 1, file=    
-!     @        dircurtis//'taul'//isotcode//dn//ibcode2//'.dat')     
-!	   endif   
-            
-!		!write(1,*) dummy     
-!		!write(1,*) 'tij for curtis matrix calculations '         
-!		!write(1,*)' cira mars model atmosphere '     
-!		write(1,*)' beta= ',beta,'deltanu= ',deltanux 
-!		write(1,*)' number of elements (in,ir)= ',nl,nl           
-!		write(1,*)' format: (tauinf(in),(tau(in,ir),ir=1,nl),in=1,nl)'
-	           
-!		do in=1,nl            
-!		    write (1,*) tauinf(in)        
-!		    do ir=1,nl        
-!			write(1,*) tau(in,ir)            
-!		    end do            
-!		end do    
-!		close(unit=1)         
-            
-!	   if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5  
-!     @    	.or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then      
-!	      write (*,'(1x, 31htransmitances written out in: ,a22)')          
-!     @    	'taul'//isotcode//dn//ibcode1    
-!	   else    
-!	      write (*,'(1x, 31htransmitances written out in: ,a22)')          
-!     @    	'taul'//isotcode//dn//ibcode2    
-!	   endif   
-            
-      end if     
-            
-c cleaning of transmittances        
-!	call elimin_tau(tau,tauinf,nl,nan,itableout,nw,dummy,     
-!     @    					isotcode,dn,ibcode2)        
-            
-c construction of the curtis matrix 
-            
-      call mzcf ( tauinf,tau, cf,cfup,cfdw, vc,taugr,            
-     @     ib,isot,icfout,itableout )            
-            
-            
-c end       
-      return     
-      end    
-
-
-
-
-c***********************************************************************
-c      mzcf
-c***********************************************************************
-                                                
-      subroutine mzcf( tauinf,tau, c,cup,cdw,vc,taugr,           
-     @     ib,isot,icfout,itableout )            
-                                                
-c     a.k.murphy method to avoid extrapolation in the curtis matrix         
-c     feb-89 	    m. angel 	granada                 
-c     25-sept-96  cristina 	dejar las matrices en doble precision           
-c     jan 98		malv	version para mz1d                
-c     jul 2011 malv+fgg       adapted to LMD-MGCM
-c***********************************************************************
-                                                
-      implicit none                                  
-
-      include 'comcstfi.h'
+
+
+c     arguments
+      real*8		c(nl,nl) ! o
+      real*8 		vc(nl)  ! o
+      real*8 		tau(nl,nl) ! i
+      real*8		tauinf(nl) ! i
+      integer		ib      ! i
+
+c     local variables
+      integer 	i, in, ir, isot
+      real*8		a(nl,nl), cf(nl,nl), pideltanu, deltazdbl,pi
+
+c***********************************************************************
+
+      pi=3.141592
+      isot = 1
+      pideltanu = pi*dble(deltanu(isot,ib))
+      deltazdbl = dble(deltaz)
+c     
+      do in=1,nl
+
+         do ir=1,nl
+
+            cf(in,ir) = 0.0d0
+            c(in,ir) = 0.0d0
+            a(in,ir) = 0.0d0
+
+         end do
+
+         vc(in) = 0.0d0
+
+      end do
+
+
+c     
+      do in=1,nl
+         do ir=1,nl
+
+            if (ir.eq.1) then
+               cf(in,ir) = tau(in,ir) - tau(in,1)
+            elseif (ir.eq.nl) then
+               cf(in,ir) = tauinf(in) - tau(in,ir-1)
+            else
+               cf(in,ir) = tau(in,ir) - tau(in,ir-1)
+            end if
+
+         end do
+      end do
+
+
+c     
+      do in=2,nl-1
+         do ir=1,nl
+            if (ir.eq.in+1) a(in,ir) = -1.d0
+            if (ir.eq.in-1) a(in,ir) = +1.d0
+            a(in,ir) = a(in,ir) / ( 2.d5*deltazdbl )
+         end do
+      end do
+
+c     
+      do in=1,nl
+         do ir=1,nl
+	    cf(in,ir) = cf(in,ir) * pideltanu
+         end do
+      end do
+
+
+      do in=2,nl-1
+         do ir=1,nl
+	    do i=1,nl
+               c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir)
+	    end do
+         end do
+         vc(in) =  pideltanu /( 2.d5*deltazdbl ) *
+     @        ( tau(in-1,1) - tau(in+1,1) )
+      end do
+
+c     
+      do in=2,nl-1
+         c(in,nl-2) = c(in,nl-2) - c(in,nl)
+         c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl)
+      end do
+
+
+c     end
+      return
+      end
+
+
+
+c     *** Old MZESC121_dlvr11_03.f ***
+
+c***********************************************************************
+      subroutine MZESC121
+c***********************************************************************
+
+      implicit none
+
       include 'nlte_paramdef.h'
       include 'nlte_commons.h'
-                                                
-c arguments                                     
-      real*8		c(nl,nl), cup(nl,nl), cdw(nl,nl) ! o   
-      real*8 		vc(nl), taugr(nl) ! o       
-      real*8 		tau(nl,nl) ! i                     
-      real*8		tauinf(nl) ! i                      
-      integer		ib      ! i                            
-      integer 	isot		! i                          
-      integer		icfout, itableout ! i               
-                                                
-c external                                      
-      external 	bandid                               
-      character*2 	bandid                            
-                                                
-c local variables                               
-      integer 	i, in, ir, iw                         
-      real*8		cfup(nl,nl), cfdw(nl,nl)               
-      real*8		a(nl,nl), cf(nl,nl)                    
-      character	isotcode*2, bcode*2                  
-                                                
-c formats                                       
- 101  format(i1)                                 
- 202  format(i2)                                 
- 180  format(a80)                                
- 181  format(a80)                                
-c***********************************************************************
-                                                
-      if (isot.eq.1)  isotcode = '26'               
-      if (isot.eq.2)  isotcode = '28'               
-      if (isot.eq.3)  isotcode = '36'               
-      if (isot.eq.4)  isotcode = '27'               
-      if (isot.eq.5)  isotcode = 'co'               
-      bcode = bandid( ib )                           
-                                                
-!	write (*,*)  ' '                                               
-                                                
-      do in=1,nl                                     
-                                                
-         do ir=1,nl                             
-                                                
-            cf(in,ir) = 0.0d0                     
-            cfup(in,ir) = 0.0d0                   
-            cfdw(in,ir) = 0.0d0                   
-            c(in,ir) = 0.0d0                      
-            cup(in,ir) = 0.0d0                    
-            cdw(in,ir) = 0.0d0                    
-            a(in,ir) = 0.0d0                      
-                                                
-         end do                                 
-                                                
-         vc(in) = 0.0d0                         
-         taugr(in) = 0.0d0                      
-                                                
-      end do                                  
-                                                
-                                                
-c	the next lines are a reduced and equivalent way of calculating        
-c	the c(in,ir) elements for n=2,nl1 and r=1,nl  
-                                                
-                                                
-c	do in=2,nl1                                   
-c	do ir=1,nl                                    
-c	if(ir.eq.1)then                               
-c	c(in,ir)=tau(in-1,1)-tau(in+1,1)              
-c	elseif(ir.eq.nl)then                          
-c	c(in,ir)=tau(in+1,nl1)-tauinf(in+1)-tau(in-1,nl1)+tauinf(in-1)        
-c	else                                          
-c	c(in,ir)=tau(in+1,ir-1)-tau(in+1,ir)-tau(in-1,ir-1)+tau(in-1,ir)      
-c	end if                                        
-c	c(in,ir)=c(in,ir)*pi*deltanu(ib)/(2.*deltaz*1.0e5)	       
-c	end do                                        
-c	end do	                                       
-c	go to 1000                                    
-                                                
-c calculation of the matrix cfup(nl,nl)         
-                                                
-      cfup(1,1) = 1.d0 - tau(1,1)             
-                                                
-      do in=2,nl                              
-         do ir=1,in                              
-                                                
-            if (ir.eq.1) then                       
-               cfup(in,ir) = tau(in,ir) - tau(in,1)        
-            elseif (ir.eq.in) then                  
-               cfup(in,ir) = 1.d0 - tau(in,ir-1)           
-            else                                    
-               cfup(in,ir) = tau(in,ir) - tau(in,ir-1)     
-            end if                                  
-            
-         end do                                  
-      end do                                  
-                                                
-! contribution to upwards fluxes from bb at bottom :        
-      do in=1,nl                              
-         taugr(in) =  tau(in,1)                
-      enddo                                   
-                                                
-c calculation of the matrix cfdw(nl,nl)         
-                                                
-      cfdw(nl,nl) = 1.d0 - tauinf(nl)         
-                                                
-      do in=1,nl-1                            
-         do ir=in,nl                             
-                                                
-            if (ir.eq.in) then                      
-               cfdw(in,ir) = 1.d0 - tau(in,ir)             
-            elseif (ir.eq.nl) then                  
-               cfdw(in,ir) = tau(in,ir-1) - tauinf(in)     
-            else                                    
-               cfdw(in,ir) = tau(in,ir-1) - tau(in,ir)     
-            end if                                  
-                                                
-         end do                                  
-      end do                                  
-                                                
-                                                
-c calculation of the matrix cf(nl,nl)           
-                                                
-      do in=1,nl                                     
-         do ir=1,nl                                     
-                                                
-            if (ir.eq.1) then                              
-	    ! version con l_bb(tg)  =  l_bb(t(1))=j(1) (see also vc below)     
-	    !	cf(in,ir) = tau(in,ir)                   
-	    ! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see also vc below)     
-               cf(in,ir) = tau(in,ir) - tau(in,1)            
-            elseif (ir.eq.nl) then                         
-               cf(in,ir) = tauinf(in) - tau(in,ir-1)         
-            else                                           
-               cf(in,ir) = tau(in,ir) - tau(in,ir-1)         
-            end if                                         
-                                                
-         end do                                         
-      end do                                         
-                                                
-                                                
-c  definition of the a(nl,nl) matrix            
-                                                
-      do in=2,nl-1                                   
-         do ir=1,nl                                      
-            if (ir.eq.in+1) a(in,ir) = -1.d0              
-            if (ir.eq.in-1) a(in,ir) = +1.d0              
-            a(in,ir) = a(in,ir) / ( 2.d0*deltaz*1.d5 )         
-         end do                                       
-      end do                                         
-! this is not needed anymore in the akm scheme  
-!	a(1,1) = +3.d0                                
-!	a(1,2) = -4.d0                                
-!	a(1,3) = +1.d0                                
-!	a(nl,nl)   = -3.d0                            
-!	a(nl,nl1) = +4.d0                             
-!	a(nl,nl2) = -1.d0                             
-                                                
-c calculation of the final curtis matrix ("reduced" by murphy's method) 
-                                                
-      if (isot.ne.5) then                            
-         do in=1,nl                                   
-            do ir=1,nl                                  
-               cf(in,ir) = cf(in,ir) * pi*deltanu(isot,ib)            
-               cfup(in,ir) = cfup(in,ir) * pi*deltanu(isot,ib)        
-               cfdw(in,ir) = cfdw(in,ir) * pi*deltanu(isot,ib)        
-            end do                                      
-            taugr(in) = taugr(in) * pi*deltanu(isot,ib) 
-         end do                                       
-      else                                           
-         do in=1,nl                                   
-            do ir=1,nl                                  
-               cf(in,ir) = cf(in,ir) * pi*deltanuco       
-            enddo                                       
-            taugr(in) = taugr(in) * pi*deltanuco        
-         enddo                                        
-      endif                                          
-                                                
-      do in=2,nl-1                                   
-                                                
-         do ir=1,nl                                   
-                                                
-	    do i=1,nl                                  
-	      ! only c contains the matrix a. matrixes cup,cdw dont because
-	      ! these two will be used for flux calculations, not  
-	      ! only for flux divergencies             
-                                                
-               c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir) 
-		! from this matrix we will extract (see below) the        
-		! nl2 x nl2 "core" for the "reduced" final curtis matrix. 
-                                                
-	    end do                                     
-	    cup(in,ir) = cfup(in,ir)                   
-	    cdw(in,ir) = cfdw(in,ir)                   
-                                                
-         end do			                                    
-	  ! version con l_bb(tg)  =  l_bb(t(1))=j(1)  (see cf above)           
-	  !vc(in) = c(in,1)                            
-	  ! version con l_bb(tg) =/= l_bb(t(1))=j(1)  (see cf above)           
-         vc(in) =  pi*deltanu(isot,ib)/( 2.d0*deltaz*1.d5 ) *     
-     @        ( tau(in-1,1) - tau(in+1,1) )         
-                                                
-      end do			                                      
-		                                              
- 5    continue                                     
-                                                
-!	write (*,*)  'mztf/1/ c(2,*) =', (c(2,i), i=1,nl)              
-                                                
-!	call elimin_dibuja(c,nl,itableout)            
-                                                
-c ventana del smoothing de c es nw=3 y de vc es 5 (puesto en lisa):     
-c subroutine elimin_mz4(c,vc,ilayer,nl,nan,iw, nw)          
-                                                
-      iw = nan                                       
-      if (isot.eq.4)  iw = 5                         
-      call elimin_mz1d (c,vc,0,iw,itableout,nw)      
-                                                
-! upper boundary condition                      
-!   j'(nl) = j'(nl1) ==> j(nl) = 2j(nl1) - j(nl2) ==>       
-      do in=2,nl-1                                   
-         c(in,nl-2) = c(in,nl-2) - c(in,nl)           
-         c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl)      
-         cup(in,nl-2) = cup(in,nl-2) - cup(in,nl)     
-         cup(in,nl-1) = cup(in,nl-1) + 2.d0*cup(in,nl)            
-         cdw(in,nl-2) = cdw(in,nl-2) - cdw(in,nl)     
-         cdw(in,nl-1) = cdw(in,nl-1) + 2.d0*cdw(in,nl)            
-      end do			                                      
-!   j(nl) = j(nl1) ==>                          
-!	do in=2,nl1                                   
-!	  c(in,nl1) = c(in,nl1) + c(in,nl)            
-!	end do			                                     
-                                                
-! 1000	continue                                 
-        
-      if (icfout.eq.1) then                          
-                                                
-! 	 if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then  
-!		codmatrx = codmatrx_fb                        
-!	 else                                           
-!		codmatrx = codmatrx_hot                       
-!	 end if                                         
-                                                
-!	 if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5    
-!     @        .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then      
-                                                
-! 	   open ( 1, access='sequential', form='unformatted', file=           
-!     @    dircurtis//'cfl'//isotcode//dn//ibcode1//codmatrx//'.dat')         
-! 	   open ( 2, access='sequential', form='unformatted', file=           
-!     @    dircurtis//'cflup'//isotcode//dn//ibcode1//codmatrx//'.dat')       
-! 	   open ( 3, access='sequential', form='unformatted', file=           
-!     @    dircurtis//'cfldw'//isotcode//dn//ibcode1//codmatrx//'.dat')       
-                                                
-!	 else                                          
-                                                
-! 	   open ( 1, access='sequential', form='unformatted', file=           
-!     @    dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat')         
-! 	   open ( 2, access='sequential', form='unformatted', file=           
-!     @    dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat')       
-! 	   open ( 3, access='sequential', form='unformatted', file=           
-!     @    dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat')       
-                                                
-!	 endif                                         
-                                                
-!	    write(1) dummy                             
-!	    write(1)' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)' 
-!	    do in=2,nl-1                               
-!	     write(1) vc(in), (c(in,ir)  , ir=2,nl-1 )		           
-	! es mas importante la precision que ocupar mucho espacio asi que      
-	! escribiremos las matrices en doble precision y por tanto en          
-	! [lib]readc_mz4.for no hay que reconvertirlas a doble precision.      
-		! ch is stored in single prec. to save storage space.     
-!	    end do                                     
-                                                
-!	    write(2) dummy                             
-!	    write(2)' format: (cfup(n,r),r=1,nl), n=1,nl)'         
-!	    do in=1,nl                                 
-!	     write(2) ( cup(in,ir)  , ir=1,nl )		      
-!	    end do                                     
-                                                
-!	    write(3) dummy                             
-!	    write(3)' format: (cfdw(n,r),r=1,nl), n=1,nl)'         
-!	    do in=1,nl                                 
-!	     write(3) (cdw(in,ir)  , ir=1,nl )		       
-!	    end do                                     
-                                                
-!	   if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5  
-!     @        .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9)  then      
-!	     write (*,'(1x,30hcurtis matrix written out in: ,a50)' )           
-!     @     dircurtis//'cfl'//isotcode//dn//ibcode1//codmatrx//'.dat'         
-!	   else                                        
-!	     write (*,'(1x,30hcurtis matrix written out in: ,a50)' )           
-!     @     dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat'         
-!	   endif                                       
-                                                
-      else                                           
-           
-	 ! write (*,*)  ' no curtis matrix output file ', char(10)     
-                                                
-      end if                                         
-                                                
-                                                
-c end                                           
-      return                                         
+
+
+c     local variables
+      integer 	i
+      real*8          factor0200, factor0220, factor1000
+      real*8          aux_0200(nl), aux2_0200(nl)
+      real*8          aux_0220(nl), aux2_0220(nl)
+      real*8          aux_1000(nl), aux2_1000(nl)
+
+c***********************************************************************
+
+!      call zerov (taustar12, nl)
+      taustar12(1:nl)=0.d0
+      call zero2v(aux_0200,aux2_0200, nl)
+      call zero2v(aux_0220,aux2_0220, nl)
+      call zero2v(aux_1000,aux2_1000, nl)
+
+      call MZESC121sub (aux_0200,aux2_0200, 2 )
+      call MZESC121sub (aux_0220,aux2_0220, 3 )
+      call MZESC121sub (aux_1000,aux2_1000, 4 )
+
+      factor0220 = 1.d0
+      factor0200 = dble( (nu(1,2)-nu(1,1)) / (nu12_0200-nu(1,1)) )
+      factor1000 = dble( (nu(1,2)-nu(1,1)) / (nu12_1000-nu(1,1)) )
+      do i=1,nl
+         taustar12(i) = taustar12(i)
+     @        + aux_0200(i) * factor0200
+     @        + aux_0220(i) * factor0220
+     @        + aux_1000(i) * factor1000
+      enddo
+
+      call mzescape_normaliz ( taustar12, 2 )
+
+c     end
+      return
       end
 
 
-
-
-
-c***********************************************************************
-c     cm15um_hb_simple
-c***********************************************************************
-                                                            
-      subroutine cm15um_hb_simple (ig,icurt)            
-                                                            
-c     computing the curtix matrixes for the 15 um hot bands   
-c     (las de las bandas fudnamentales las calcula cm15um_fb) 
-                                                            
-c     jan 98 		malv		version de mod3/cm_15um.f para mz1d 
-c     jul 2011 malv+fgg               adapted to LMD-MGCM
-c***********************************************************************
-                                                            
-      implicit none                                  
-                                                            
-!!!!!!!!!!!!!!!!!!!!!!!                         
-! common variables & constants                  
-                                                            
+c     *** Old MZESC121sub_dlvr11_03.f ***
+
+c***********************************************************************
+
+      subroutine MZESC121sub (taustar,tauinf, ib )
+
+c***********************************************************************
+
+      implicit none
+
+      include 'datafile.h'
       include 'nlte_paramdef.h'
       include 'nlte_commons.h'
-                                                            
-!!!!!!!!!!!!!!!!!!!!!!!                         
-! arguments                                     
-                   
-      integer ig                ! ADDED FOR TRACEBACK
-      integer icurt             ! icurt=0,1,2                  
-					! new calculations? (see caa.f heads)
-                                                            
-!!!!!!!!!!!!!!!!!!!!!!!                         
-! local variables                               
-                                                            
-      real*4 cdummy(nl,nl), csngl(nl,nl)             
-                                                            
-      real*8 cax1(nl,nl), cax2(nl,nl), cax3(nl,nl)   
-      real*8 v1(nl), v2(nl), v3(nl), cm_factor, vc_factor        
-                                                            
-      integer itauout,icfout,itableout, interpol,ismooth, isngldble          
-      integer i,j,ik,ist,isot,ib,itt                 
-                                                            
-	!character 	bandcode*2
-      character 	isotcode*2
-      character 	codmatrx_hot*5                      
-                                                            
-!!!!!!!!!!!!!!!!!!!!!!!                         
-! external functions                            
-                                                            
-      external bandid                                
-      character*2 bandid                             
-                                                            
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!               
-! subroutines called:                           
-! 	mz4sub, dmzout, readc_mz4, readcupdw, mztf   
-                                                            
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!               
-! formatos                                      
- 132  format(i2)                                 
-                                                            
+
+
+c     arguments
+      real*8 		taustar(nl) ! o
+      real*8 		tauinf(nl) ! o
+      integer		ib      ! i
+
+
+c     local variables and constants
+      integer 	i, iaquiHIST, iaquiZ, isot
+      real*8 		con(nzy), coninf
+      real*8 		c1, c2, ccc
+      real*8 		t1, t2
+      real*8 		p1, p2
+      real*8		mr1, mr2
+      real*8 		st1, st2
+      real*8 		c1box(70), c2box(70)
+      real*8		ff      ! to avoid too small numbers
+      real*8		tvtbs(nzy)
+      real*8 		st, beta, ts
+      real*8  	zld(nl), zyd(nzy)
+      real*8 		correc
+      real*8 		deltanudbl, deltazdbl
+      real*8          yy
+
+c     external function
+      external        we_clean
+      real*8  we_clean
+
+c     formats
+ 101  format(i1)
+
+c***********************************************************************
+
+c     
+      beta = 1.8d5
+      isot = 1
+      write ( ibcode1, 101) ib
+      deltanudbl = dble( deltanu(isot,ib) )
+      ff=1.0d10
+      deltazdbl = dble(deltaz)
+
+c     
+      do i=1,nzy
+         zyd(i) = dble(zy(i))
+      enddo
+      do i=1,nl
+         zld(i) = dble(zl(i))
+      enddo
+
+      call interhuntdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
+
+      do i=1,nzy
+         con(i) =  dble( co2y(i) * imr(isot) )
+         correc = 2.d0 * dexp( -ee*dble(elow(isot,2))/tvtbs(i) )
+         con(i) = con(i) * ( 1.d0 - correc )
+         mr(i) = dble(co2y(i)/nty(i))
+      end do
+
+      coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
+      call mztf_correccion ( coninf, con, ib )
+
+c     
+      call gethist_03 ( ib )
+
+c     
+c     tauinf
+c     
+      call initial
+
+      iaquiHIST = nhist/2
+      iaquiZ = nzy - 2
+
+      do i=nl,1,-1
+
+         if(i.eq.nl)then
+
+            call intzhunt (iaquiZ, zl(i),c2,p2,mr2,t2, con)
+            do kr=1,nbox
+               ta(kr)=t2
+	    end do
+            call interstrhunt (iaquiHIST, st2,t2,ka,ta)
+            aa = p2 * coninf * mr2 * (st2 * ff)
+            cc = coninf * st2
+            dd = t2 * coninf * st2
+            do kr=1,nbox
+               ccbox(kr) = coninf * ka(kr)
+               ddbox(kr) = t2 * ccbox(kr)
+               c2box(kr) = c2 * ka(kr) * deltazdbl
+            end do
+            c2 = c2 * st2 * deltazdbl
+
+         else
+            call intzhunt (iaquiZ, zl(i),c1,p1,mr1,t1, con)
+            do kr=1,nbox
+               ta(kr)=t1
+            end do
+            call interstrhunt (iaquiHIST,st1,t1,ka,ta)
+            do kr=1,nbox
+               c1box(kr) = c1 * ka(kr) * deltazdbl
+            end do
+            c1 = c1 * st1 * deltazdbl
+            aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
+            cc = cc + ( c1 + c2 ) / 2.d0
+            ccc = ( c1 + c2 ) / 2.d0
+            dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
+            do kr=1,nbox
+               ccbox(kr) = ccbox(kr) +
+     @              ( c1box(kr) + c2box(kr) )/2.d0
+               ddbox(kr) = ddbox(kr) +
+     @              ( t1*c1box(kr)+t2*c2box(kr) )/2.d0
+            end do
+
+            mr2 = mr1
+            c2=c1
+            do kr=1,nbox
+               c2box(kr) = c1box(kr)
+            end do
+            t2=t1
+            p2=p1
+         end if
+
+         pp = aa / (cc*ff)
+
+         ts = dd/cc
+         do kr=1,nbox
+   	    ta(kr) = ddbox(kr) / ccbox(kr)
+         end do
+         call interstrhunt(iaquiHIST,st,ts,ka,ta)
+         call intershphunt(iaquiHIST,alsa,alda,ta)
+
+c     
+         eqw=0.0d0
+         do  kr=1,nbox
+            yy = ccbox(kr) * beta
+            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
+            eqw = eqw + no(kr)*w
+         end do
+         tauinf(i) = dexp( - eqw / deltanudbl )
+         if (tauinf(i).lt.0.d0) tauinf(i) = 0.0d0
+
+         if (i.eq.nl) then
+            taustar(i) = 0.0d0
+         else
+            taustar(i) = deltanudbl * (tauinf(i+1)-tauinf(i))
+     @           / ( beta * ccc  )
+         endif
+
+      end do
+
+
+
+c     end
+      return
+      end
+
+
+c     *** Old MZTVC121_dlvr11.f *** 
+
+c***********************************************************************
+
+      subroutine MZTVC121
+
+c***********************************************************************
+
+      implicit none
+
+!!!!!!!!!!!!!!!!!!!!!!!
+!     common variables & constants
+
+      include 'nlte_paramdef.h'
+      include 'nlte_commons.h'
+
+
+!     arguments
+      integer ierr
+      real*8 varerr
+
+
+!     local variables
+
+      real*8 v1(nl), vc_factor
+      integer i,ik,ib
+
 ************************************************************************
-************************************************************************
-                                                            
-      call zerom (c121,nl)                           
-      
-      call zerov (vc121,nl)                          
-      
-      call zerom (cup121,nl)                         
-      call zerom (cdw121,nl)                         
-      
-      call zerov (taugr121,nl)                       
-                                                            
-                                                            
-      itauout = 0		! =1 --> with output of tau       
-      icfout = 0		! =1 --> with output of cf         
-      itableout = 0		! =1 --> with output of table of taus       
-      isngldble = 1		! =1 --> dble precission        
-                                                            
-      codmatrx_hot='     '
-      if (icurt.eq.2) then                           
-         icfout=1                                      
-      elseif (icurt.eq.0) then                       
-         write (*,'(a,a$)')                            
-     @        ' hot bands. code for old matrixes (5 chars): '     
-         read (*,'(a)')  codmatrx_hot                  
-      endif                                          
-                                                            
-      fileroot = 'cfl'                               
-                                                            
-! ====================== curtis matrixes for fh bands ==================
-                                                            
-                                                            
-! una piedra en el camino ...                   
-!	write (*,*)  ' cm15um_hb/1 '                                   
-                                                            
-ccc                                             
-      if ( input_cza.ge.1 ) then                     
-ccc                                             
-                                                            
-         if (icurt.eq.2) then                           
-            write (*,'(a,a$)')                           
-     @           '  new calculation of curt. mat. for fh bands.',         
-     @           '    code for new matrixes : '               
-            read (*,'(a)') codmatrx_hot                  
-         elseif (icurt.eq.0) then                       
-            write (*,'(a,a$)')                           
-     @           '  reading in curt. mat. for fh bands.',     
-     @           '    code for old matrixes : '               
-            read (*,'(a)') codmatrx_hot                  
-         else                                           
-!     write (*,'(a)')                             
-!     @        '  new calculation of curt. mat. for fh bands.'         
-         end if                                         
-                                                            
-!       fh bands for the 626 isotope ================================-  
-                                                            
-         ist = 1                                        
-         isot = 26                                      
-! 	encode (2,132,isotcode) isot                  
-         write (isotcode,132) isot  
-                             
-         do 11, ik=1,3                                  
-            
-            ib=ik+1                                      
-                                                            
-            if (icurt.gt.0) then                         
-               call zero3m (cax1,cax2,cax3,nl)            
-! una piedra en el camino ...                   
-	    !write (*,*)  ' cm15um_hb/11 '                                  
-	    !write (*,*)  ' ib, ist, irw, imu =', ib, ist, irw_mztf, imu     
-               call mztf(ig,cax1,cax2,cax3,v1,v2,ib,ist,irw_mztf,imu,
-     @              itauout,icfout,itableout)                    
-!	  else                                         
-!	    bandcode = bandid(ib)                      
-!	    filend=isotcode//dn//bandcode//codmatrx_hot            
-!!	    write (*,*)  char(9), fileroot//filend                     
-!	    call zero3m (cax1,cax2,cax3,nl)            
-!	    call readcud_mz1d ( cax1,cax2,cax3, v1, v2,            
-!     @        	fileroot,filend, csngl, nl,nan,0,isngldble) 
-            end if                                       
-                                                            
-c 	  calculating the total c121(n,r) matrix for the first hot band      
-            do i=1,nl                                    
-                                                            
-               if ( ib .eq. 4 ) then                       
-!	    write (*,*)  ' '                                           
-!	    write (*,*)  i, ' ib,ist, altura :', ib,ist, zl(i)         
-               endif                                       
-                                                            
-!	    if ( v1(i) .le. 1.d-99 ) v1(i) = 0.0d0    
-!	    if ( v2(i) .le. 1.d-99 ) v2(i) = 0.0d0    
-                                                            
-                                                            
-               if(ik.eq.1)then                            
-                  cm_factor = (dble(618.03/667.75))**2.d0*      
-     @                 exp( dble(ee*(667.75-618.03)/t(i)) )     
-                  vc_factor = dble(667.75/618.03)               
-               elseif(ik.eq.2)then                        
-                  cm_factor = 1.d0                              
-                  vc_factor = 1.d0                              
-               elseif(ik.eq.3)then                        
-                  cm_factor = ( dble(720.806/667.75) )**2.d0*   
-     @                 exp( dble(ee*(667.75-720.806)/t(i)) )    
-                  vc_factor = dble(667.75/720.806)              
-               end if                                     
-               do j=1,nl                                  
-!	       if (cax1(i,j) .le. 1.d-99 ) cax1(i,j) = 0.0d0      
-!	       if (cax2(i,j) .le. 1.d-99 ) cax2(i,j) = 0.0d0      
-!	       if (cax3(i,j) .le. 1.d-99 ) cax3(i,j) = 0.0d0      
-                  c121(i,j) = c121(i,j) + cax1(i,j) * cm_factor        
-                  cup121(i,j) = cup121(i,j) + cax2(i,j) * cm_factor    
-                  cdw121(i,j) = cdw121(i,j) + cax3(i,j) * cm_factor    
-               end do                                     
-                                                            
-!	    write (*,*)  ' i =', i                                     
-!	    write (*,*)  ' vc_factor =', vc_factor                     
-!	    write (*,*)  ' v1 =', v1(i)                                
-!	    write (*,*)  ' v2 =', v2(i)                                
-!	    write (*,*)  vc121(i), taugr121(i)                         
-!	    write (*,*)  v1(i) * vc_factor                             
-!	    write (*,*)  vc121(i) + v1(i) * vc_factor                  
-                                                            
-               vc121(i) = vc121(i) + v1(i) * vc_factor    
-            
-                                      
-!	    write (*,*)  v2(i) * vc_factor                             
-!	    write (*,*)  taugr121(i) + v2(i) * vc_factor               
-                                                            
-               taugr121(i) = taugr121(i) + v2(i) * vc_factor          
-                                                            
-            end do                                       
- 11      continue                                     
-                                                            
-ccc                                             
-      end if                                         
-ccc                                             
-                                                            
-                                                            
-      return                                         
-      end 
-
-
-
-
-
+
+!      call zerov( vc121, nl )
+      vc121(1:nl)=0.d0
+
+      do 11, ik=1,3
+
+         ib=ik+1
+
+         call MZTVC121sub (v1, ib, ierr,varerr )
+
+         do i=1,nl
+
+	    if(ik.eq.1)then
+               vc_factor =
+     @              dble( (nu(1,2)-nu(1,1)) / (nu12_0200-nu(1,1)) )
+	    elseif(ik.eq.2)then
+               vc_factor = 1.d0
+	    elseif(ik.eq.3)then
+               vc_factor =
+     @              dble( (nu(1,2)-nu(1,1)) / (nu12_1000-nu(1,1)) )
+	    end if
+
+	    vc121(i) = vc121(i) + v1(i) * vc_factor
+
+         end do
+
+ 11   continue
+
+
+      return
+      end
+
+
+c     *** Old MZTVC121sub_dlvr11_03.f ***
+
+c***********************************************************************
+c     mztf.f
+c***********************************************************************
+
+      subroutine MZTVC121sub  ( vc, ib,  ierr, varerr )
+
+c***********************************************************************
+
+      implicit none
+
+      include 'datafile.h'
+      include 'nlte_paramdef.h'
+      include 'nlte_commons.h'
+
+
+c     arguments
+      real*8		vc(nl)  ! o
+      integer		ib      ! i
+      integer ierr              ! o
+      real*8  varerr            ! o
+
+c     local variables and constants
+      integer 	i, in, ir, iaquiHIST , iaquiZ, isot
+      integer 	nmu
+      parameter 	(nmu = 8)
+      real*8 		tau(nl,nl), argumento
+      real*8 		con(nzy), coninf
+      real*8 		c1, c2
+      real*8 		t1, t2
+      real*8 		p1, p2
+      real*8		mr1, mr2
+      real*8 		st1, st2
+      real*8 		c1box(70), c2box(70)
+      real*8		ff      ! to avoid too small numbers
+      real*8		tvtbs(nzy)
+      real*8 		st, beta, ts
+      real*8  	zld(nl), zyd(nzy), deltazdbl
+      real*8 		correc
+      real*8 		deltanudbl, pideltanu,pi
+      real*8          yy
+      real*8          minvc, maxtau
+
+c     external function
+      external        we_clean
+      real*8          we_clean
+
+c     formats
+ 101  format(i1)
+
+c***********************************************************************
+
+c     
+      pi=3.141592
+      isot = 1
+      beta = 1.8d5
+      write (ibcode1,101) ib
+      deltanudbl = dble( deltanu(isot,ib) )
+      pideltanu = pi*deltanudbl
+      ff=1.0d10
+      deltazdbl = dble(deltaz)
+c     
+c     
+c     
+
+      do i=1,nzy
+         zyd(i) = dble(zy(i))
+      enddo
+      do i=1,nl
+         zld(i) = dble(zl(i))
+      enddo
+
+      call interhuntdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
+
+      do i=1,nzy
+         con(i) =  dble( co2y(i) * imr(isot) )
+         correc = 2.d0 * dexp( -ee*dble(elow(isot,2))/tvtbs(i) )
+         con(i) = con(i) * ( 1.d0 - correc )
+         mr(i) = dble(co2y(i)/nty(i))
+      end do
+
+      coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
+      call mztf_correccion ( coninf, con, ib )
+
+ccc   
+      call gethist_03 ( ib )
+
+c     
+c     tau(1,ir)
+c     
+      call initial
+
+      iaquiHIST = nhist/2
+
+      in=1
+
+      tau(in,1) = 1.d0
+
+      call initial
+      iaquiZ = 2
+      call intzhunt ( iaquiZ, zl(in), c1,p1,mr1,t1, con)
+      do kr=1,nbox
+         ta(kr) = t1
+      end do
+      call interstrhunt (iaquiHIST, st1,t1,ka,ta)
+      do kr=1,nbox
+         c1box(kr) = c1 * ka(kr) * deltazdbl
+      end do
+      c1 = c1 * st1 * deltazdbl
+                                ! Check interpolation errors :
+      if (c1.le.0.0d0) then
+         ierr=15
+         varerr=c1
+         return
+      elseif (p1.le.0.0d0) then
+         ierr=16
+         varerr=p1
+         return
+      elseif (mr1.le.0.0d0) then
+         ierr=17
+         varerr=mr1
+         return
+      elseif (t1.le.0.0d0) then
+         ierr=18
+         varerr=t1
+         return
+      elseif (st1.le.0.0d0) then
+         ierr=19
+         varerr=st1
+         return
+      endif
+                                !
+
+      do 2 ir=2,nl
+
+         call intzhunt (iaquiZ, zl(ir), c2,p2,mr2,t2, con)
+         do kr=1,nbox
+            ta(kr) = t2
+         end do
+         call interstrhunt (iaquiHIST, st2,t2,ka,ta)
+         do kr=1,nbox
+            c2box(kr) = c2 * ka(kr) * deltazdbl
+         end do
+         c2 = c2 * st2 * deltazdbl
+
+         aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0
+         cc = cc + ( c1 + c2 ) / 2.d0
+         dd = dd + ( t1*c1 + t2*c2 ) / 2.d0
+         do kr=1,nbox
+            ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
+            ddbox(kr) = ddbox(kr) + 
+     $           ( t1*c1box(kr) + t2*c2box(kr) ) / 2.d0
+         end do
+
+         mr1=mr2
+         t1=t2
+         c1=c2
+         p1=p2
+         do kr=1,nbox
+            c1box(kr) = c2box(kr)
+         end do
+
+         pp = aa / (cc * ff)
+
+         ts = dd/cc
+         do kr=1,nbox
+   	    ta(kr) = ddbox(kr) / ccbox(kr)
+         end do
+         call interstrhunt(iaquiHIST, st,ts,ka,ta)
+         call intershphunt(iaquiHIST, alsa,alda,ta)
+
+         eqw=0.0d0
+         do kr=1,nbox
+            yy = ccbox(kr) * beta
+            w = we_clean ( yy, pp, alsa(kr),alda(kr) )
+            eqw = eqw + no(kr)*w
+         end do
+
+         argumento = eqw / deltanudbl
+         tau(in,ir) = dexp( - argumento )
+
+ 2    continue
+
+
+c     
+c     
+c     
+      do in=nl,2,-1
+         tau(in,1) = tau(1,in)
+      end do
+
+c     
+      vc(1) = 0.0d0
+      vc(nl) = 0.0d0
+      do in=2,nl-1
+         vc(in) =  pideltanu /( 2.d5*deltazdbl ) *
+     @        ( tau(in-1,1) - tau(in+1,1) )
+         if (vc(in) .lt. 0.0d0) vc(in) = vc(in-1)
+      end do
+
+c     
+c     Tracking potential numerical errors
+c     
+      minvc = 1.d6
+      maxtau = tau(nl,1)
+      do in=2,nl-1
+         minvc = min( minvc, vc(in) )
+         maxtau = max( maxtau, tau(in,1) )
+      end do
+      if (maxtau .gt. 1.0d0) then
+         ierr = 13
+         varerr = maxtau
+         return
+      else if (minvc .lt. 0.0d0) then
+         ierr = 14
+         varerr = minvc
+         return
+      endif
+
+c     end
+      return
+      end
+
+
+
+
+
+
+
+
+
