!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Fast scheme for NLTE cooling rates at 15um by CO2 in a Martian GCM !
!                 Version dlvr11_03. 2012.                           !
! Software written and provided by IAA/CSIC, Granada, Spain,         !
! under ESA contract "Mars Climate Database and Physical Models"     !
! Person of contact: Miguel Angel Lopez Valverde  valverde@iaa.es    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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,ierr,varerr)
c***********************************************************************

      implicit none

#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, iaquiHIST , iaquiZ
      integer   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***********************************************************************
      ierr = 0
      varerr = 0.d0

      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
      if ( con(nzy_cts_real) .le. 0.0d0 ) then
         ierr = 33
         varerr = con(nzy_cts_real)
         return
      elseif ( con(nzy_cts_real-1) .le. con(nzy_cts_real) ) then
         write (*,*) ' WARNING in MZESC110 '
         write (*,*) '    [CO2] growing with altitude at TOA.'
         write (*,*) '    [CO2] @ TOA = ', con(nzy_cts_real)
         coninf = dble( con(nzy_cts_real) )
      else
         coninf = dble( con(nzy_cts_real) /
     @        log( con(nzy_cts_real-1) / con(nzy_cts_real) ) )
      endif
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 "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***********************************************************************

      ierr = 0
      varerr = 0.d0

      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
      if ( con(nzy) .le. 0.0d0 ) then
         ierr = 43
         varerr = con(nzy)
         return
      elseif ( con(nzy-1) .le. con(nzy) ) then
         write (*,*) ' WARNING in MZTUD110 '
         write (*,*) '    [CO2] grows with height at CurtisMatrix top.'
         write (*,*) '    [CO2] @ top = ', con(nzy)
         coninf = dble( con(nzy) )
      else
         coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
      endif
      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=45
               varerr=c2
               return
            elseif (p2.le.0.0d0) then
               ierr=46
               varerr=p2
               return
            elseif (mr2.le.0.0d0) then
               ierr=47
               varerr=mr2
               return
            elseif (t2.le.0.0d0) then
               ierr=48
               varerr=t2
               return
            elseif (st2.le.0.0d0) then
               ierr=49
               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
            ! Check interpolation errors :
            if (c1.le.0.0d0) then
               ierr=75
               varerr=c1
               return
            elseif (p1.le.0.0d0) then
               ierr=76
               varerr=p1
               return
            elseif (mr1.le.0.0d0) then
               ierr=77
               varerr=mr1
               return
            elseif (t1.le.0.0d0) then
               ierr=78
               varerr=t1
               return
            elseif (st1.le.0.0d0) then
               ierr=79
               varerr=st1
               return
            endif
            !
            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 = 42
         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
      
#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
      real*8  varerr

      integer i,j,ik,ib
      integer ierr      

************************************************************************

      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, ierr, varerr )
         if (ierr .gt. 0) call ERRORS (ierr,varerr)

         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, ierr, varerr )
c***********************************************************************

      implicit none

#include "nlte_paramdef.h"
#include "nlte_commons.h"
      

c     arguments
      real*8  	cf(nl,nl) ! o.
      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
      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
      real*8  	zld(nl), zyd(nzy)
      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***********************************************************************

      ierr = 0
      varerr = 0.d0

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

      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

      if ( con(nzy) .le. 0.0d0 ) then
         ierr = 83
         varerr = con(nzy)
         return
      elseif ( con(nzy-1) .le. con(nzy) ) then
         write (*,*) ' WARNING in MZTUD121 '
         write (*,*) '    [CO2] grows with height at CurtisMatrix top.'
         write (*,*) '    [CO2] @ top = ', con(nzy)
         coninf = dble( con(nzy) )
      else
         coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
      endif
      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
            ! Check interpolation errors :
            if (c1.le.0.0d0) then
               ierr=85
               varerr=c1
               return
            elseif (p1.le.0.0d0) then
               ierr=86
               varerr=p1
               return
            elseif (mr1.le.0.0d0) then
               ierr=87
               varerr=mr1
               return
            elseif (t1.le.0.0d0) then
               ierr=88
               varerr=t1
               return
            elseif (st1.le.0.0d0) then
               ierr=89
               varerr=st1
               return
            endif
	    !
            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
      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     local variables
      integer 	      i,ierr
      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)
      real*8          varerr 

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 , ierr, varerr)
      if (ierr .gt. 0) call ERRORS (ierr,varerr)
      call MZESC121sub (aux_0220,aux2_0220, 3 , ierr, varerr)
      if (ierr .gt. 0) call ERRORS (ierr,varerr)
      call MZESC121sub (aux_1000,aux2_1000, 4 , ierr, varerr)
      if (ierr .gt. 0) call ERRORS (ierr,varerr)

      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     *** Old MZESC121sub_dlvr11_03.f ***

c***********************************************************************

      subroutine MZESC121sub (taustar,tauinf, ib, ierr, varerr )

c***********************************************************************

      implicit none

#include "nlte_paramdef.h"
#include "nlte_commons.h"


c     arguments
      real*8 	      taustar(nl) ! o
      real*8 	      tauinf(nl) ! o
      integer	      ib      ! i
      integer         ierr        ! o
      real*8          varerr      ! o


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***********************************************************************

      ierr = 0
      varerr = 0.d0

      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
      if ( con(nzy) .le. 0.0d0 ) then
         ierr = 63
         varerr = con(nzy)
         return
      elseif ( con(nzy-1) .le. con(nzy) ) then
         write (*,*) ' WARNING in MZESC121sub '
         write (*,*) '    [CO2] grows with height at CurtisMatrix top.'
         write (*,*) '    [CO2] @ top = ', con(nzy)
         coninf = dble( con(nzy) )
      else
         coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
      endif
      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)
            ! Check interpolation errors :
            if (c2.le.0.0d0) then
               ierr=65
               varerr=c2
               return
            elseif (p2.le.0.0d0) then
               ierr=66
               varerr=p2
               return
            elseif (mr2.le.0.0d0) then
               ierr=67
               varerr=mr2
               return
            elseif (t2.le.0.0d0) then
               ierr=68
               varerr=t2
               return
            elseif (st2.le.0.0d0) then
               ierr=69
               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
            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 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***********************************************************************

      subroutine MZTVC121sub  ( vc, ib,  ierr, varerr )

c***********************************************************************

      implicit none

#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***********************************************************************

      ierr = 0
      varerr = 0.d0
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

      if ( con(nzy) .le. 0.0d0 ) then
         ierr = 53
         varerr = con(nzy)
         return
      elseif ( con(nzy-1) .le. con(nzy) ) then
         write (*,*) ' WARNING in MZTVC121sub '
         write (*,*) '    [CO2] grows with height at CurtisMatrix top.'
         write (*,*) '    [CO2] @ top = ', con(nzy)
         coninf = dble( con(nzy) )
      else
         coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
      endif
      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=55
         varerr=c1
         return
      elseif (p1.le.0.0d0) then
         ierr=56
         varerr=p1
         return
      elseif (mr1.le.0.0d0) then
         ierr=57
         varerr=mr1
         return
      elseif (t1.le.0.0d0) then
         ierr=58
         varerr=t1
         return
      elseif (st1.le.0.0d0) then
         ierr=59
         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 = 52
         varerr = maxtau
         return
      else if (minvc .lt. 0.0d0) then
         ierr = 51
         varerr = minvc
         return
      endif

c     end
      return
      end









