!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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     Includes the following old 1-D model files/subroutines

c     -MZTCRSUB_dlvr11.f
c     *dinterconnection
c     *planckd
c     *leetvt
c     -MZTFSUB_dlvr11_02.f
c     *initial
c     *intershphunt
c     *interstrhunt
c     *intzhunt
c     *intzhunt_cts
c     *rhist
c     *we_clean
c     *mztf_correccion
c     *mzescape_normaliz
c     *mzescape_normaliz_02
c     -interdpESCTVCISO_dlvr11.f
c     -hunt_cts.f
c     -huntdp.f
c     -hunt.f
c     -interdp_limits.f
c     -interhunt2veces.f
c     -interhunt5veces.f
c     -interhuntdp3veces.f
c     -interhuntdp4veces.f
c     -interhuntdp.f
c     -interhunt.f
c     -interhuntlimits2veces.f
c     -interhuntlimits5veces.f
c     -interhuntlimits.f
c     -lubksb_dp.f
c     -ludcmp_dp.f
c     -LUdec.f
c     -mat_oper.f
c     *unit
c     *diago
c     *invdiag
c     *samem
c     *mulmv
c     *trucodiag
c     *trucommvv
c     *sypvmv
c     *mulmm
c     *resmm
c     *sumvv
c     *sypvvv
c     *zerom
c     *zero4m
c     *zero3m
c     *zero2m
c     *zerov
c     *zero4v
c     *zero3v
c     *zero2v
c     -suaviza.f

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


c     *** Old MZTCRSUB_dlvr11.f ***

!************************************************************************

!      subroutine dinterconnection ( v, vt )


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

!      implicit none
!#include "nlte_paramdef.h"

c     argumentos
!      real*8 vt(nl), v(nl)

c     local variables
!      integer 	i

c     *************
!
!      do i=1,nl
!         v(i) = vt(i)
!      end do

!      return
!      end

c***********************************************************************
      function planckdp(tp,xnu)
c***********************************************************************

      implicit none

#include "nlte_paramdef.h"

      real*8 planckdp
      real*8 xnu
      real tp

      planckdp = gamma*xnu**3.0d0 / exp( ee*xnu/dble(tp) )
                                !erg cm-2.sr-1/cm-1.

c     end
      return
      end

c***********************************************************************
      subroutine leetvt

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

      implicit none

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

c     local variables
      integer i
      real*8	zld(nl), zyd(nzy)
      real*8  xvt11(nzy), xvt21(nzy), xvt31(nzy), xvt41(nzy) 

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

      do i=1,nzy
         zyd(i) = dble(zy(i))
         xvt11(i)= dble( ty(i) )
         xvt21(i)= dble( ty(i) )
         xvt31(i)= dble( ty(i) )
         xvt41(i)= dble( ty(i) )
      end do

      do i=1,nl
         zld(i) = dble( zl(i) )
      enddo
      call interhuntdp4veces ( v626t1,v628t1,v636t1,v627t1, zld,nl,
     $     xvt11, xvt21, xvt31, xvt41, zyd,nzy, 1 )


c     end
      return
      end


c     *** MZTFSUB_dlvr11_02.f ***


c     ****************************************************************
      subroutine initial

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

      implicit none

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

c     local variables
      integer 	i

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

      eqw = 0.0d00
      aa = 0.0d00
      cc = 0.0d00
      dd = 0.0d00

      do i=1,nbox
         ccbox(i) = 0.0d0
         ddbox(i) = 0.0d0
      end do

      return
      end

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

      subroutine intershphunt (i, alsx,adx,xtemp)

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

      implicit none

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

c     arguments
      real*8 alsx(nbox_max),adx(nbox_max) ! Output
      real*8 xtemp(nbox_max)    ! Input
      integer    i              ! I , O

c     local variables
      integer 	k
      real*8          factor
      real*8    temperatura     ! para evitar valores ligeramnt out of limits
      character(len=100) text

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

      do 1, k=1,nbox
         temperatura = xtemp(k)
         if (abs(xtemp(k)-thist(1)).le.0.01d0) then
            temperatura=thist(1)
         elseif (abs(xtemp(k)-thist(nhist)).le.0.01d0) then
            temperatura=thist(nhist)
         elseif (xtemp(k).lt.thist(1)) then
            temperatura=thist(1)
c            write (*,*) ' WARNING intershphunt/ Too low atmosph Tk:'
c            write (*,*) ' WARNING      k,xtemp = ', k,xtemp(k)
c            write (*,*) ' Minimum Tk in histogram used : ', thist(1)
         elseif (xtemp(k).gt.thist(nhist)) then
            temperatura=thist(nhist)
c            write (*,*) ' WARNING intershphunt/ Very high atmosph Tk:'
c            write (*,*) ' WARNING      k,xtemp = ', k,xtemp(k)
c            write (*,*) ' Max Tk in histogram used : ', thist(nhist)
         endif
         call huntdp ( thist,nhist, temperatura, i )
         if ( i.eq.0 .or. i.eq.nhist ) then
c	    write (*,*) ' HUNT/ Limits input grid:',
c     @           thist(1),thist(nhist)
c	    write (*,*) ' HUNT/ location in grid:', xtemp(k)
            text='INTERSHP/ Interpolation error. T out of Histogram.'
            CALL abort_physic("nlte_aux",text, 1)
         endif
         factor = 1.d0 /  (thist(i+1)-thist(i))
         alsx(k) = (( xls1(i,k)*(thist(i+1)-xtemp(k)) +
     @        xls1(i+1,k)*(xtemp(k)-thist(i)) )) * factor
         adx(k)  = (( xld1(i,k)*(thist(i+1)-xtemp(k)) +
     @        xld1(i+1,k)*(xtemp(k)-thist(i)) )) * factor
 1    continue

      return
      end

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

      subroutine interstrhunt (i, stx, ts, sx, xtemp )

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

      implicit none

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

c     arguments
      real*8 		stx     ! output, total band strength
      real*8		ts      ! input, temp for stx
      real*8		sx(nbox_max) ! output, strength for each box
      real*8		xtemp(nbox_max) ! input, temp for sx
      integer 	i

c     local variables
      integer 	k
      real*8          factor
      real*8    temperatura
      character(len=100) text

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

      do 1, k=1,nbox
         temperatura = xtemp(k)
         if (abs(xtemp(k)-thist(1)).le.0.01d0) then
            temperatura=thist(1)
         elseif (abs(xtemp(k)-thist(nhist)).le.0.01d0) then
            temperatura=thist(nhist)
         elseif (xtemp(k).lt.thist(1)) then
            temperatura=thist(1)
c            write (*,*) ' WARNING interstrhunt/ Too low atmosph Tk:'
c	    write (*,*) ' WARNING     k,xtemp(k) = ', k,xtemp(k)
c	    write (*,*) ' Minimum Tk in histogram used : ', thist(1)
         elseif (xtemp(k).gt.thist(nhist)) then
            temperatura=thist(nhist)
c            write (*,*) ' WARNING interstrhunt/ Very high atmosph Tk:'
c	    write (*,*) ' WARNING     k,xtemp(k) = ', k,xtemp(k)
c	    write (*,*) ' Max Tk in histogram used : ', thist(nhist)
         endif
         call huntdp ( thist,nhist, temperatura, i )
         if ( i.eq.0 .or. i.eq.nhist ) then
c	    write(*,*)'HUNT/ Limits input grid:',
c     $           thist(1),thist(nhist)
c	    write(*,*)'HUNT/ location in grid:',xtemp(k)
            text='INTERSTR/1/ Interpolation error. T out of Histogram.'
            CALL abort_physic("nlte_aux",text, 1)
         endif
         factor = 1.d0 /  (thist(i+1)-thist(i))
         sx(k) = ( sk1(i,k)   * (thist(i+1)-xtemp(k))
     @        + sk1(i+1,k) * (xtemp(k)-thist(i))  ) * factor
 1    continue


      temperatura = ts
      if (abs(ts-thist(1)).le.0.01d0) then
         temperatura=thist(1)
      elseif (abs(ts-thist(nhist)).le.0.01d0) then
         temperatura=thist(nhist)
      elseif (ts.lt.thist(1)) then
         temperatura=thist(1)
c         write (*,*) ' WARNING interstrhunt/ Too low atmosph Tk:'
c         write (*,*) ' WARNING            ts = ', temperatura
c         write (*,*) ' Minimum Tk in histogram used : ', thist(1)
      elseif (ts.gt.thist(nhist)) then
         temperatura=thist(nhist)
c         write (*,*) ' WARNING interstrhunt/ Very high atmosph Tk:'
c         write (*,*) ' WARNING            ts = ', temperatura
c         write (*,*) ' Max Tk in histogram used : ', thist(nhist)
      endif
      call huntdp ( thist,nhist, temperatura, i )
      if ( i.eq.0 .or. i.eq.nhist ) then
c         write (*,*) ' HUNT/ Limits input grid:',
c     @        thist(1),thist(nhist)
c         write (*,*) ' HUNT/ location in grid:', ts
         text='INTERSTR/2/ Interpolat error. T out of Histogram.'
         CALL abort_physic("nlte_aux",text, 1)
      endif
      factor = 1.d0 /  (thist(i+1)-thist(i))
      stx = 0.d0
      do k=1,nbox
         stx = stx + no(k) * ( sk1(i,k)*(thist(i+1)-ts) +
     @        sk1(i+1,k)*(ts-thist(i)) ) * factor
      end do


      return
      end

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

      subroutine intzhunt (k, h, aco2,ap,amr,at, con)

c     k lleva la posicion de la ultima llamada a intz , necesario para
c     que esto represente una aceleracion real.
c     **********************************************************************

      implicit none
#include "nlte_paramdef.h"
#include "nlte_commons.h"

c     arguments
      real		h       ! i
      real*8		con(nzy) ! i
      real*8		aco2, ap, at, amr ! o
      integer		k       ! i

c     local variables
      real          factor

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

      call hunt ( zy,nzy, h, k )
      factor =  (h-zy(k)) /  (zy(k+1)-zy(k))
      ap = dble( exp( log(py(k)) + log(py(k+1)/py(k)) * factor ) )
      aco2 = dlog(con(k)) + dlog( con(k+1)/con(k) ) * dble(factor)
      aco2 = exp( aco2 )
      at = dble( ty(k) + (ty(k+1)-ty(k)) * factor )
      amr = dble( mr(k) + (mr(k+1)-mr(k)) * factor )


      return
      end

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

      subroutine intzhunt_cts (k, h, nzy_cts_real, 
     @     aco2,ap,amr,at, con)

c     k lleva la posicion de la ultima llamada a intz , necesario para
c     que esto represente una aceleracion real.
c     **********************************************************************

      implicit none
#include "nlte_paramdef.h"
#include "nlte_commons.h"

c     arguments
      real		h       ! i
      real*8		con(nzy_cts) ! i
      real*8		aco2, ap, at, amr ! o
      integer		k       ! i
      integer         nzy_cts_real ! i

c     local variables
      real          factor

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

      call hunt_cts ( zy_cts,nzy_cts, nzy_cts_real, h, k )
      factor =  (h-zy_cts(k)) /  (zy_cts(k+1)-zy_cts(k))
      ap = dble( exp( log(py_cts(k)) +
     @     log(py_cts(k+1)/py_cts(k)) * factor ) )
      aco2 = dlog(con(k)) + dlog( con(k+1)/con(k) ) * dble(factor)
      aco2 = exp( aco2 )
      at = dble( ty_cts(k) + (ty_cts(k+1)-ty_cts(k)) * factor )
      amr = dble( mr_cts(k) + (mr_cts(k+1)-mr_cts(k)) * factor )


      return
      end


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

      real*8 function we_clean  ( y,pl, xalsa, xalda )

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

      implicit none

#include "nlte_paramdef.h"

c     arguments
      real*8		y       ! I. path's absorber amount * strength
      real*8          pl        ! I. path's partial pressure of CO2
      real*8          xalsa     ! I.  Self lorentz linewidth for 1 isot & 1 box
      real*8          xalda     ! I.  Doppler linewidth        "           "

c     local variables
      integer 	i
      real*8 		x,wl,wd,wvoigt
      real*8		cn(0:7),dn(0:7)
      real*8          factor, denom
      real*8          pi, pi2, sqrtpi

c     data blocks
      data cn/9.99998291698d-1,-3.53508187098d-1,9.60267807976d-2,
     @     -2.04969011013d-2,3.43927368627d-3,-4.27593051557d-4,
     @     3.42209457833d-5,-1.28380804108d-6/
      data dn/1.99999898289,5.774919878d-1,-5.05367549898d-1,
     @     8.21896973657d-1,-2.5222672453,6.1007027481,
     @     -8.51001627836,4.6535116765/

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

      pi = 3.141592
      pi2= 6.28318531
      sqrtpi = 1.77245385

      x=y / ( pi2 * xalsa*pl )


c     Lorentz
      wl=y/sqrt(1.0d0+pi*x/2.0d0)

c     Doppler
      x = y / (xalda*sqrtpi)
      if (x .lt. 5.0d0) then
         wd = cn(0)
         factor = 1.d0
         do i=1,7
            factor = factor * x
	    wd = wd + cn(i) * factor
         end do
         wd = xalda * x * sqrtpi * wd
      else
         wd = dn(0)
         factor = 1.d0 / log(x)
         denom = 1.d0
         do i=1,7
            denom = denom * factor
	    wd = wd + dn(i) * denom
         end do
         wd = xalda * sqrt(log(x)) * wd
      end if

c     Voigt
      wvoigt = wl*wl + wd*wd - (wd*wl/y)*(wd*wl/y)

      if ( wvoigt .lt. 0.0d0 ) then
       write (*,*) ' Subroutine WE/ Error in Voift EQS calculation'
       write (*,*) '  WL, WD, X, Y = ', wl, wd, x, y
       stop '  ERROR : Imaginary EQW. Revise spectral data. '
      endif

      we_clean = sqrt( wvoigt )


      return
      end


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

      subroutine mztf_correccion (coninf, con, ib )

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

      implicit none

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

c     arguments
      integer		ib
      real*8		con(nzy), coninf

!     local variables
      integer 	i, isot
      real*8	tvt0(nzy), tvtbs(nzy), zld(nl),zyd(nzy)
      real*8  xqv, xes, xlower, xfactor

c     *********

      isot = 1
      nu11 = dble( nu(1,1) )

      do i=1,nzy
         zyd(i) = dble(zy(i))
      enddo
      do i=1,nl
         zld(i) = dble( zl(i) )
      end do

!     tvtbs
      call interhuntdp (tvtbs,zyd,nzy, v626t1,zld,nl, 1 )

!     tvt0
      if (ib.eq.2 .or. ib.eq.3 .or. ib.eq.4) then
         call interhuntdp (tvt0,zyd,nzy, v626t1,zld,nl, 1 )
      else
         do i=1,nzy
            tvt0(i) = dble( ty(i) )
         end do
      end if

c     factor
      do i=1,nzy

         xlower = exp( ee*dble(elow(isot,ib)) *
     @        ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )
         xes = 1.0d0
         xqv = ( 1.d0-exp( -ee*nu11/tvtbs(i) ) ) /
     @        (1.d0-exp( -ee*nu11/dble(ty(i)) ))
         xfactor = xlower * xqv**2.d0 * xes

         con(i) = con(i) * xfactor
         if (i.eq.nzy) coninf = coninf * xfactor

      end do


      return
      end


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

      subroutine mzescape_normaliz ( taustar, istyle )

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

      implicit none
#include "nlte_paramdef.h"

c     arguments
      real*8 		taustar(nl) ! o
      integer         istyle    ! i

c     local variables and constants
      integer 	i, imaximum
      real*8          maximum

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

      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

      do i=1,nl
         taustar(i) = taustar(i) / maximum
      enddo


c     end
      return
      end

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

      subroutine mzescape_normaliz_02 ( taustar, nn, istyle )

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

      implicit none

c     arguments
      real*8 		taustar(nn) ! i,o
      integer         istyle    ! i
      integer         nn        ! i

c     local variables and constants
      integer 	i, imaximum
      real*8          maximum

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

      taustar(nn) = taustar(nn-1)

      if ( istyle .eq. 1 ) then
         imaximum = nn
         maximum = taustar(nn)
         do i=1,nn-1
	    if (taustar(i).gt.maximum) taustar(i) = taustar(nn)
         enddo
      elseif ( istyle .eq. 2 ) then
         imaximum = nn
         maximum = taustar(nn)
         do i=nn-1,1,-1
	    if (taustar(i).gt.maximum) then
	       maximum = taustar(i)
	       imaximum = i
	    endif
         enddo
         do i=imaximum,nn
	    if (taustar(i).lt.maximum) taustar(i) = maximum
         enddo
      endif

      do i=1,nn
         taustar(i) = taustar(i) / maximum
      enddo


c     end
      return
      end


c     *** interdp_ESCTVCISO_dlvr11.f ***

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

      subroutine interdp_ESCTVCISO

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

      implicit none

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

c     local variables
      integer    i
      real*8     lnpnb(nl)


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

c     Use pressure in the NLTE grid but in log and in nb
      do i=1,nl
         lnpnb(i) = log( dble( pl(i) * 1013.25 * 1.e6) )
      enddo

c     Interpolations

      call interhuntdp3veces
     @     ( taustar21,taustar31,taustar41,    lnpnb, nl,
     @     tstar21tab,tstar31tab,tstar41tab, lnpnbtab, nztabul,
     @     1 )

      call interhuntdp3veces ( vc210,vc310,vc410, lnpnb, nl,
     @     vc210tab,vc310tab,vc410tab, lnpnbtab, nztabul, 2 )

c     end
      return
      end


c     *** hunt_cts.f ***

cccc  
      SUBROUTINE hunt_cts(xx,n,n_cts,x,jlo)  
c     
c     La dif con hunt es el uso de un indice superior (n_cts) mas bajito que (n)
c     
c     Arguments
      INTEGER jlo               ! O
      INTEGER n                 ! I
      INTEGER n_cts             ! I
      REAL  xx(n)               ! I
      REAL  x                   ! I

c     Local variables
      INTEGER inc,jhi,jm  
      LOGICAL ascnd  
c     
cccc  
c     
      ascnd=xx(n_cts).ge.xx(1)  
      if(jlo.le.0.or.jlo.gt.n_cts)then  
         jlo=0  
         jhi=n_cts+1  
         goto 3  
      endif  
      inc=1  
      if(x.ge.xx(jlo).eqv.ascnd)then  
 1       jhi=jlo+inc  
!     write (*,*) jlo
         if(jhi.gt.n_cts)then  
            jhi=n_cts+1  
!     write (*,*) jhi-1
         else if(x.ge.xx(jhi).eqv.ascnd)then  
            jlo=jhi  
            inc=inc+inc  
!     write (*,*) jlo
            goto 1  
         endif  
      else  
         jhi=jlo  
 2       jlo=jhi-inc  
!     write (*,*) jlo
         if(jlo.lt.1)then  
            jlo=0  
         else if(x.lt.xx(jlo).eqv.ascnd)then  
            jhi=jlo  
            inc=inc+inc  
            goto 2  
         endif  
      endif  
 3    if(jhi-jlo.eq.1)then  
         if(x.eq.xx(n_cts))jlo=n_cts-1  
         if(x.eq.xx(1))jlo=1  
!     write (*,*) jlo
         return  
      endif  
      jm=(jhi+jlo)/2  
      if(x.ge.xx(jm).eqv.ascnd)then  
         jlo=jm  
      else  
         jhi=jm  
      endif  
!     write (*,*) jhi-1
      goto 3  
c     
      END  

      
c     *** huntdp.f ***

cccc  
      SUBROUTINE huntdp(xx,n,x,jlo)  
c     
c     Arguments
      INTEGER jlo               ! O
      INTEGER n                 ! I
      REAL*8  xx(n)             ! I
      REAL*8  x                 ! I

c     Local variables
      INTEGER inc,jhi,jm  
      LOGICAL ascnd  
c     
cccc  
c     
      ascnd=xx(n).ge.xx(1)  
      if(jlo.le.0.or.jlo.gt.n)then  
         jlo=0  
         jhi=n+1  
         goto 3  
      endif  
      inc=1  
      if(x.ge.xx(jlo).eqv.ascnd)then  
 1       jhi=jlo+inc  
         if(jhi.gt.n)then  
            jhi=n+1  
         else if(x.ge.xx(jhi).eqv.ascnd)then  
            jlo=jhi  
            inc=inc+inc  
            goto 1  
         endif  
      else  
         jhi=jlo  
 2       jlo=jhi-inc  
         if(jlo.lt.1)then  
            jlo=0  
         else if(x.lt.xx(jlo).eqv.ascnd)then  
            jhi=jlo  
            inc=inc+inc  
            goto 2  
         endif  
      endif  
 3    if(jhi-jlo.eq.1)then  
         if(x.eq.xx(n))jlo=n-1  
         if(x.eq.xx(1))jlo=1  
         return  
      endif  
      jm=(jhi+jlo)/2  
      if(x.ge.xx(jm).eqv.ascnd)then  
         jlo=jm  
      else  
         jhi=jm  
      endif  
      goto 3  
c     
      END  

      
c     *** hunt.f ***

cccc  
      SUBROUTINE hunt(xx,n,x,jlo)  
c     
c     Arguments
      INTEGER jlo               ! O
      INTEGER n                 ! I
      REAL  xx(n)               ! I
      REAL  x                   ! I

c     Local variables
      INTEGER inc,jhi,jm  
      LOGICAL ascnd  
c     
cccc  
c     
      ascnd=xx(n).ge.xx(1)  
      if(jlo.le.0.or.jlo.gt.n)then  
         jlo=0  
         jhi=n+1  
         goto 3  
      endif  
      inc=1  
      if(x.ge.xx(jlo).eqv.ascnd)then  
 1       jhi=jlo+inc  
!     write (*,*) jlo
         if(jhi.gt.n)then  
            jhi=n+1  
!     write (*,*) jhi-1
         else if(x.ge.xx(jhi).eqv.ascnd)then  
            jlo=jhi  
            inc=inc+inc  
!     write (*,*) jlo
            goto 1  
         endif  
      else  
         jhi=jlo  
 2       jlo=jhi-inc  
!     write (*,*) jlo
         if(jlo.lt.1)then  
            jlo=0  
         else if(x.lt.xx(jlo).eqv.ascnd)then  
            jhi=jlo  
            inc=inc+inc  
            goto 2  
         endif  
      endif  
 3    if(jhi-jlo.eq.1)then  
         if(x.eq.xx(n))jlo=n-1  
         if(x.eq.xx(1))jlo=1  
!     write (*,*) jlo
         return  
      endif  
      jm=(jhi+jlo)/2  
      if(x.ge.xx(jm).eqv.ascnd)then  
         jlo=jm  
      else  
         jhi=jm  
      endif  
!     write (*,*) jhi-1
      goto 3  
c     
      END  

      
c     *** interdp_limits.f ***

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

      subroutine interdp_limits ( yy, zz, m,   i1,i2, 
     @     y,  z, n,   j1,j2,  opt)

c     Interpolation soubroutine. 
c     Returns values between indexes i1 & i2, donde  1 =< i1 =< i2 =< m
c     Solo usan los indices de los inputs entre j1,j2, 1 =< j1 =< j2 =< n    
c     Input values: y(n) , z(n)  (solo se usarann los valores entre j1,j2)
c     zz(m) (solo se necesita entre i1,i2)
c     Output values: yy(m) (solo se calculan entre i1,i2)
c     Options:    opt=1 -> lineal ,,  opt=2 -> logarithmic
c     Difference with interdp:  
c     here interpolation proceeds between indexes i1,i2 only 
c     if i1=1 & i2=m, both subroutines are exactly the same
c     thus previous calls to interdp or interdp2 could be easily replaced

c     JAN 98 	MALV 		Version for mz1d
c     ***********************************************************************

      implicit none

!     Arguments 
      integer 	n,m             ! I. Dimensions
      integer 	i1, i2, j1, j2, opt ! I
      real*8 		zz(m)   ! I
      real*8 		yy(m)   ! O
      real*8		z(n),y(n) ! I

!     Local variables
      integer 	i,j
      real*8 		zmin,zzmin,zmax,zzmax

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

!     write (*,*) ' d interpolating '
!     call mindp_limits (z,n,zmin, j1,j2)
!     call mindp_limits (zz,m,zzmin, i1,i2)
!     call maxdp_limits (z,n,zmax, j1,j2)
!     call maxdp_limits (zz,m,zzmax, i1,i2)
      zmin=minval(z(j1:j2))
      zzmin=minval(zz(i1:i2))
      zmax=maxval(z(j1:j2))
      zzmax=maxval(zz(i1:i2))

      if(zzmin.lt.zmin)then
         write (*,*) 'from d interp: new variable out of limits'
         write (*,*) zzmin,'must be .ge. ',zmin
         stop
!     elseif(zzmax.gt.zmax)then
!     type *,'from interp: new variable out of limits'
!     type *,zzmax, 'must be .le. ',zmax
!     stop
      end if

      do 1,i=i1,i2

         do 2,j=j1,j2-1
            if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3
 2       continue
c     in this case (zz(i2).eq.z(j2)) and j leaves the loop with j=j2-1+1=j2
         if(opt.eq.1)then
            yy(i)=y(j2-1)+(y(j2)-y(j2-1))*(zz(i)-z(j2-1))/
     $           (z(j2)-z(j2-1))
         elseif(opt.eq.2)then
            if(y(j2).eq.0.0d0.or.y(j2-1).eq.0.0d0)then
               yy(i)=0.0d0
            else
               yy(i)=exp(log(y(j2-1))+log(y(j2)/y(j2-1))*
     @              (zz(i)-z(j2-1))/(z(j2)-z(j2-1)))
            end if
         else
            write (*,*) ' d interp : opt must be 1 or 2, opt= ',opt
         end if
         goto 1
 3       continue
         if(opt.eq.1)then
            yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j))
!     type *, ' '
!     type *, ' z(j),z(j+1) =', z(j),z(j+1)
!     type *, ' t(j),t(j+1) =', y(j),y(j+1)
!     type *, ' zz, tt =  ', zz(i), yy(i)
         elseif(opt.eq.2)then
            if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then
               yy(i)=0.0d0
            else
               yy(i)=exp(log(y(j))+log(y(j+1)/y(j))*
     @              (zz(i)-z(j))/(z(j+1)-z(j)))
            end if
         else
            write (*,*) ' interp : opt must be 1 or 2, opt= ',opt
         end if
 1    continue
      return
      end



c     *** interhunt2veces.f ***

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

      subroutine interhunt2veces ( y1,y2,  zz,m, 
     @     x1,x2,  z,n,  opt)

c     interpolation soubroutine basada en Numerical Recipes HUNT.FOR
c     input values: y(n) at z(n) 
c     output values: yy(m) at zz(m)
c     options: 1 -> lineal
c     2 -> logarithmic
c     ***********************************************************************

      implicit none

!     Arguments
      integer	n,m,opt         ! I
      real	zz(m),z(n)      ! I
      real    y1(m),y2(m)       ! O
      real    x1(n),x2(n)       ! I


!     Local variables 
      integer i, j 
      real    factor 
      real    zaux
      character(len=100) text

!!!!  

      j = 1                     ! initial first guess (=n/2 is anothr pssblty)

      do 1,i=1,m                ! 

                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
         zaux = zz(i)
         if (abs(zaux-z(1)).le.0.01) then 
            zaux=z(1)
         elseif (abs(zaux-z(n)).le.0.01) then
            zaux=z(n)
         endif
         call hunt ( z,n, zaux, j )
         if ( j.eq.0 .or. j.eq.n ) then 
	    write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
	    write (*,*) ' HUNT/ location in new grid:', zz(i)
            text='interhunt2/ Interpolat error. zz out of limits.'
            CALL abort_physic("nlte_aux",text, 1)
         endif

                                ! Perform interpolation
         factor = (zz(i)-z(j))/(z(j+1)-z(j))
         if (opt.eq.1) then
	    y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
	    y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
         else
	    y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	    y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
         end if

 1    continue

      return
      end


c     *** interhunt5veces.f ***

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

      subroutine interhunt5veces ( y1,y2,y3,y4,y5,  zz,m, 
     @     x1,x2,x3,x4,x5,  z,n,  opt)

c     interpolation soubroutine basada en Numerical Recipes HUNT.FOR
c     input values: y(n) at z(n) 
c     output values: yy(m) at zz(m)
c     options: 1 -> lineal
c     2 -> logarithmic
c     ***********************************************************************

      implicit none
!     Arguments
      integer	n,m,opt         ! I
      real	zz(m),z(n)      ! I
      real    y1(m),y2(m),y3(m),y4(m),y5(m) ! O
      real    x1(n),x2(n),x3(n),x4(n),x5(n) ! I


!     Local variables 
      integer i, j 
      real    factor 
      real    zaux
      character(len=100) text

!!!!  

      j = 1                     ! initial first guess (=n/2 is anothr pssblty)

      do 1,i=1,m                ! 

                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
         zaux = zz(i)
         if (abs(zaux-z(1)).le.0.01) then 
            zaux=z(1)
         elseif (abs(zaux-z(n)).le.0.01) then
            zaux=z(n)
         endif
         call hunt ( z,n, zaux, j )
         if ( j.eq.0 .or. j.eq.n ) then 
	    write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
	    write (*,*) ' HUNT/ location in new grid:', zz(i)
            text='interhunt5/ Interpolat error. zz out of limits.'
            CALL abort_physic("nlte_aux",text, 1)
         endif

                                ! Perform interpolation
         factor = (zz(i)-z(j))/(z(j+1)-z(j))
         if (opt.eq.1) then
	    y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
	    y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
	    y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor
	    y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor
	    y5(i) = x5(j) + (x5(j+1)-x5(j)) * factor
         else
	    y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	    y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
	    y3(i) = exp( log(x3(j)) + log(x3(j+1)/x3(j)) * factor )
	    y4(i) = exp( log(x4(j)) + log(x4(j+1)/x4(j)) * factor )
	    y5(i) = exp( log(x5(j)) + log(x5(j+1)/x5(j)) * factor )
         end if

 1    continue

      return
      end



c     *** interhuntdp3veces.f ***

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

      subroutine interhuntdp3veces ( y1,y2,y3, zz,m, 
     @     x1,x2,x3,  z,n,  opt)

c     interpolation soubroutine basada en Numerical Recipes HUNT.FOR
c     input values: x(n) at z(n) 
c     output values: y(m) at zz(m)
c     options: opt = 1 -> lineal
c     opt=/=1 -> logarithmic
c     ***********************************************************************
!     Arguments
      integer	n,m,opt         ! I
      real*8	zz(m),z(n)      ! I
      real*8    y1(m),y2(m),y3(m) ! O
      real*8    x1(n),x2(n),x3(n) ! I


!     Local variables 
      integer i, j 
      real*8    factor 
      real*8    zaux
      character(len=100) text

!!!!  

      j = 1                     ! initial first guess (=n/2 is anothr pssblty)

      do 1,i=1,m                ! 

                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
         zaux = zz(i)
         if (abs(zaux-z(1)).le.0.01d0) then 
            zaux=z(1)
         elseif (abs(zaux-z(n)).le.0.01d0) then
            zaux=z(n)
         endif
         call huntdp ( z,n, zaux, j )
         if ( j.eq.0 .or. j.eq.n ) then 
	    write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
	    write (*,*) ' HUNT/ location in new grid:', zz(i)
            text='INTERHUNTDP3/ Interpolat error. zz out of limits.'
            CALL abort_physic("nlte_aux",text, 1)
         endif

                                ! Perform interpolation
         factor = (zz(i)-z(j))/(z(j+1)-z(j))
         if (opt.eq.1) then
	    y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
	    y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
	    y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor
         else
	    y1(i) = dexp( dlog(x1(j)) + dlog(x1(j+1)/x1(j)) * factor )
	    y2(i) = dexp( dlog(x2(j)) + dlog(x2(j+1)/x2(j)) * factor )
	    y3(i) = dexp( dlog(x3(j)) + dlog(x3(j+1)/x3(j)) * factor )
         end if

 1    continue

      return
      end


c     *** interhuntdp4veces.f ***

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

      subroutine interhuntdp4veces ( y1,y2,y3,y4, zz,m, 
     @     x1,x2,x3,x4,  z,n,  opt)

c     interpolation soubroutine basada en Numerical Recipes HUNT.FOR
c     input values: x1(n),x2(n),x3(n),x4(n) at z(n) 
c     output values: y1(m),y2(m),y3(m),y4(m) at zz(m)
c     options: 1 -> lineal
c     2 -> logarithmic
c     ***********************************************************************

      implicit none

!     Arguments
      integer	n,m,opt         ! I
      real*8	zz(m),z(n)      ! I
      real*8    y1(m),y2(m),y3(m),y4(m) ! O
      real*8    x1(n),x2(n),x3(n),x4(n) ! I


!     Local variables 
      integer i, j 
      real*8    factor 
      real*8    zaux
      character(len=100) text

!!!!  

      j = 1                     ! initial first guess (=n/2 is anothr pssblty)

      do 1,i=1,m                ! 

                                ! Caza del indice j donde ocurre que zz(i) esta entre [z(j),z(j+1)]
         zaux = zz(i)
         if (abs(zaux-z(1)).le.0.01d0) then 
            zaux=z(1)
         elseif (abs(zaux-z(n)).le.0.01d0) then
            zaux=z(n)
         endif
         call huntdp ( z,n, zaux, j )
         if ( j.eq.0 .or. j.eq.n ) then 
	    write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
	    write (*,*) ' HUNT/ location in new grid:', zz(i)
            text='INTERHUNTDP4/ Interpolat error. zz out of limits.'
            CALL abort_physic("nlte_aux",text, 1)
         endif

                                ! Perform interpolation
         factor = (zz(i)-z(j))/(z(j+1)-z(j))
         if (opt.eq.1) then
	    y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
	    y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
	    y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor
	    y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor
         else
	    y1(i) = dexp( dlog(x1(j)) + dlog(x1(j+1)/x1(j)) * factor )
	    y2(i) = dexp( dlog(x2(j)) + dlog(x2(j+1)/x2(j)) * factor )
	    y3(i) = dexp( dlog(x3(j)) + dlog(x3(j+1)/x3(j)) * factor )
	    y4(i) = dexp( dlog(x4(j)) + dlog(x4(j+1)/x4(j)) * factor )
         end if

 1    continue

      return
      end


c     *** interhuntdp.f ***

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

      subroutine interhuntdp ( y1, zz,m, 
     @     x1,  z,n,  opt)

c     interpolation soubroutine basada en Numerical Recipes HUNT.FOR
c     input values: x1(n) at z(n) 
c     output values: y1(m) at zz(m)
c     options: 1 -> lineal
c     2 -> logarithmic
c     ***********************************************************************

      implicit none 

!     Arguments
      integer	n,m,opt         ! I
      real*8	zz(m),z(n)      ! I
      real*8    y1(m)           ! O
      real*8    x1(n)           ! I


!     Local variables 
      integer i, j 
      real*8    factor 
      real*8    zaux
      character(len=100) text

!!!!  

      j = 1                     ! initial first guess (=n/2 is anothr pssblty)

      do 1,i=1,m                ! 

                                ! Caza del indice j donde ocurre que zz(i) esta entre [z(j),z(j+1)]
         zaux = zz(i)
         if (abs(zaux-z(1)).le.0.01d0) then 
            zaux=z(1)
         elseif (abs(zaux-z(n)).le.0.01d0) then
            zaux=z(n)
         endif
         call huntdp ( z,n, zaux, j )
         if ( j.eq.0 .or. j.eq.n ) then 
	    write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
	    write (*,*) ' HUNT/ location in new grid:', zz(i)
            text='INTERHUNT/ Interpolat error. zz out of limits.'
            CALL abort_physic("nlte_aux",text, 1)
         endif

                                ! Perform interpolation
         factor = (zz(i)-z(j))/(z(j+1)-z(j))
         if (opt.eq.1) then
	    y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
         else
	    y1(i) = dexp( dlog(x1(j)) + dlog(x1(j+1)/x1(j)) * factor )
         end if

 1    continue

      return
      end


c     *** interhunt.f ***

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

      subroutine interhunt ( y1, zz,m, 
     @     x1,  z,n,  opt)

c     interpolation soubroutine basada en Numerical Recipes HUNT.FOR
c     input values: x1(n) at z(n) 
c     output values: y1(m) at zz(m)
c     options: 1 -> lineal
c     2 -> logarithmic
c***********************************************************************

      implicit none

!     Arguments
      integer	n,m,opt         ! I
      real	zz(m),z(n)      ! I
      real    y1(m)             ! O
      real    x1(n)             ! I


!     Local variables 
      integer i, j 
      real    factor 
      real    zaux
      character(len=100) text

!!!!  

      j = 1                     ! initial first guess (=n/2 is anothr pssblty)

      do 1,i=1,m                ! 

                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
         zaux = zz(i)
         if (abs(zaux-z(1)).le.0.01) then 
            zaux=z(1)
         elseif (abs(zaux-z(n)).le.0.01) then
            zaux=z(n)
         endif
         call hunt ( z,n, zaux, j )
         if ( j.eq.0 .or. j.eq.n ) then 
	    write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
	    write (*,*) ' HUNT/ location in new grid:', zz(i)
            text='interhunt/ Interpolat error. z out of limits.'
            CALL abort_physic("nlte_aux",text, 1)
         endif

                                ! Perform interpolation
         factor = (zz(i)-z(j))/(z(j+1)-z(j))
         if (opt.eq.1) then
	    y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
         else
	    y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
         end if


 1    continue

      return
      end


c     *** interhuntlimits2veces.f ***

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

      subroutine interhuntlimits2veces 
     @     ( y1,y2, zz,m, limite1,limite2, 
     @     x1,x2,  z,n, opt)

c     Interpolation soubroutine basada en Numerical Recipes HUNT.FOR
c     Input values:  x1,x2(n) at z(n) 
c     Output values: 
c     y1,y2(m) at zz(m)   pero solo entre los indices de zz
c     siguientes: [limite1,limite2]
c     Options: 1 -> linear in z and linear in x
c     2 -> linear in z and logarithmic in x
c     3 -> logarithmic in z and linear in x
c     4 -> logarithmic in z and logaritmic in x
c     
c     NOTAS: Esta subrutina extiende y generaliza la usual  
c     "interhunt5veces" en 2 direcciones: 
c     - la condicion en los limites es que zz(limite1:limite2)
c     esté dentro de los limites de z (pero quizas no todo zz)
c     - se han añadido 3 opciones mas al caso de interpolacion
c     logaritmica, ahora se hace en log de z, de x o de ambos.
c     Notese que esta subrutina engloba a la interhunt5veces 
c     ( esta es reproducible haciendo  limite1=1  y  limite2=m 
c     y usando una de las 2 primeras opciones  opt=1,2 )
c***********************************************************************

      implicit none

!     Arguments
      integer	n,m,opt, limite1,limite2 ! I
      real	zz(m),z(n)      ! I
      real    y1(m),y2(m)       ! O
      real    x1(n),x2(n)       ! I


!     Local variables 
      integer i, j 
      real    factor 
      real    zaux
      character(len=100) text

!!!!  

      j = 1                     ! initial first guess (=n/2 is anothr pssblty)

      do 1,i=limite1,limite2              

                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
         zaux = zz(i)
         if (abs(zaux-z(1)).le.0.01) then 
            zaux=z(1)
         elseif (abs(zaux-z(n)).le.0.01) then
            zaux=z(n)
         endif
         call hunt ( z,n, zaux, j )
         if ( j.eq.0 .or. j.eq.n ) then 
	    write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
	    write (*,*) ' HUNT/ location in new grid:', zz(i)
            text='interhuntlimits/ Interpolat error. z out of limits.'
            CALL abort_physic("nlte_aux",text, 1)
         endif

                                ! Perform interpolation
         if (opt.eq.1) then
	    factor = (zz(i)-z(j))/(z(j+1)-z(j))
	    y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
	    y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor 

         elseif (opt.eq.2) then 
	    factor = (zz(i)-z(j))/(z(j+1)-z(j))
	    y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	    y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
         elseif (opt.eq.3) then 
	    factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
	    y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
	    y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor 
         elseif (opt.eq.4) then 
	    factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
	    y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	    y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
         end if


 1    continue

      return
      end


c     *** interhuntlimits5veces.f ***

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

      subroutine interhuntlimits5veces 
     @     ( y1,y2,y3,y4,y5, zz,m, limite1,limite2, 
     @     x1,x2,x3,x4,x5,  z,n, opt)

c     Interpolation soubroutine basada en Numerical Recipes HUNT.FOR
c     Input values:  x1,x2,..,x5(n) at z(n) 
c     Output values: 
c     y1,y2,...,y5(m) at zz(m)   pero solo entre los indices de zz
c     siguientes: [limite1,limite2]
c     Options: 1 -> linear in z and linear in x
c     2 -> linear in z and logarithmic in x
c     3 -> logarithmic in z and linear in x
c     4 -> logarithmic in z and logaritmic in x
c     
c     NOTAS: Esta subrutina extiende y generaliza la usual  
c     "interhunt5veces" en 2 direcciones: 
c     - la condicion en los limites es que zz(limite1:limite2)
c     esté dentro de los limites de z (pero quizas no todo zz)
c     - se han añadido 3 opciones mas al caso de interpolacion
c     logaritmica, ahora se hace en log de z, de x o de ambos.
c     Notese que esta subrutina engloba a la interhunt5veces 
c     ( esta es reproducible haciendo  limite1=1  y  limite2=m 
c     y usando una de las 2 primeras opciones  opt=1,2 )
c***********************************************************************

      implicit none

!     Arguments
      integer	n,m,opt, limite1,limite2 ! I
      real	zz(m),z(n)      ! I
      real    y1(m),y2(m),y3(m),y4(m),y5(m) ! O
      real    x1(n),x2(n),x3(n),x4(n),x5(n) ! I


!     Local variables 
      integer i, j 
      real    factor 
      real    zaux
      character(len=100) text

!!!!  

      j = 1                     ! initial first guess (=n/2 is anothr pssblty)

      do 1,i=limite1,limite2              

                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
         zaux = zz(i)
         if (abs(zaux-z(1)).le.0.01) then 
            zaux=z(1)
         elseif (abs(zaux-z(n)).le.0.01) then
            zaux=z(n)
         endif
         call hunt ( z,n, zaux, j )
         if ( j.eq.0 .or. j.eq.n ) then 
	    write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
	    write (*,*) ' HUNT/ location in new grid:', zz(i)
            text='interhuntlimits/ Interpolat error. z out of limits.'
            CALL abort_physic("nlte_aux",text, 1)
         endif

                                ! Perform interpolation
         if (opt.eq.1) then
	    factor = (zz(i)-z(j))/(z(j+1)-z(j))
	    y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
	    y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor 
	    y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor 
	    y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor 
	    y5(i) = x5(j) + (x5(j+1)-x5(j)) * factor 
         elseif (opt.eq.2) then 
	    factor = (zz(i)-z(j))/(z(j+1)-z(j))
	    y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	    y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
	    y3(i) = exp( log(x3(j)) + log(x3(j+1)/x3(j)) * factor )
	    y4(i) = exp( log(x4(j)) + log(x4(j+1)/x4(j)) * factor )
	    y5(i) = exp( log(x5(j)) + log(x5(j+1)/x5(j)) * factor )
         elseif (opt.eq.3) then 
	    factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
	    y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
	    y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor 
	    y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor 
	    y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor 
	    y5(i) = x5(j) + (x5(j+1)-x5(j)) * factor 
         elseif (opt.eq.4) then 
	    factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
	    y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
	    y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
	    y3(i) = exp( log(x3(j)) + log(x3(j+1)/x3(j)) * factor )
	    y4(i) = exp( log(x4(j)) + log(x4(j+1)/x4(j)) * factor )
	    y5(i) = exp( log(x5(j)) + log(x5(j+1)/x5(j)) * factor )
         end if


 1    continue

      return
      end



c     *** interhuntlimits.f ***

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

      subroutine interhuntlimits ( y, zz,m, limite1,limite2, 
     @     x,  z,n, opt)

c     Interpolation soubroutine basada en Numerical Recipes HUNT.FOR
c     Input values:  x(n) at z(n) 
c     Output values: y(m) at zz(m)   pero solo entre los indices de zz
c     siguientes: [limite1,limite2]
c     Options: 1 -> linear in z and linear in x
c     2 -> linear in z and logarithmic in x
c     3 -> logarithmic in z and linear in x
c     4 -> logarithmic in z and logaritmic in x
c     
c     NOTAS: Esta subrutina extiende y generaliza la usual  "interhunt"
c     en 2 direcciones: 
c     - la condicion en los limites es que zz(limite1:limite2)
c     esté dentro de los limites de z (pero quizas no todo zz)
c     - se han añadido 3 opciones mas al caso de interpolacion
c     logaritmica, ahora se hace en log de z, de x o de ambos.
c     Notese que esta subrutina engloba a la usual interhunt 
c     ( esta es reproducible haciendo  limite1=1  y  limite2=m 
c     y usando una de las 2 primeras opciones  opt=1,2 )
c***********************************************************************

      implicit none

!     Arguments
      integer	n,m,opt, limite1,limite2 ! I
      real	zz(m),z(n)      ! I
      real    y(m)              ! O
      real    x(n)              ! I


!     Local variables 
      integer i, j 
      real    factor 
      real    zaux
      character(len=100) text

!!!!  

      j = 1                     ! initial first guess (=n/2 is anothr pssblty)

      do 1,i=limite1,limite2              

                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
         zaux = zz(i)
         if (abs(zaux-z(1)).le.0.01) then 
            zaux=z(1)
         elseif (abs(zaux-z(n)).le.0.01) then
            zaux=z(n)
         endif
         call hunt ( z,n, zaux, j )
         if ( j.eq.0 .or. j.eq.n ) then 
	    write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
	    write (*,*) ' HUNT/ location in new grid:', zz(i)
            text='interhuntlimits/ Interpolat error. z out of limits.'
            CALL abort_physic("nlte_aux",text, 1)
         endif

                                ! Perform interpolation
         if (opt.eq.1) then
	    factor = (zz(i)-z(j))/(z(j+1)-z(j))
	    y(i) = x(j) + (x(j+1)-x(j)) * factor 
         elseif (opt.eq.2) then 
	    factor = (zz(i)-z(j))/(z(j+1)-z(j))
	    y(i) = exp( log(x(j)) + log(x(j+1)/x(j)) * factor )
         elseif (opt.eq.3) then 
	    factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
	    y(i) = x(j) + (x(j+1)-x(j)) * factor 
         elseif (opt.eq.4) then 
	    factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
	    y(i) = exp( log(x(j)) + log(x(j+1)/x(j)) * factor )
         end if


 1    continue

      return
      end


c     *** lubksb_dp.f ***

      subroutine lubksb_dp(a,n,np,indx,b)     

      implicit none

      integer,intent(in) :: n,np
      real*8,intent(in) :: a(np,np) 
      integer,intent(in) :: indx(n)
      real*8,intent(out) :: b(n) 

      real*8 sum
      integer ii, ll, i, j

      ii=0            
      do 12 i=1,n           
         ll=indx(i)          
         sum=b(ll)           
         b(ll)=b(i)          
         if (ii.ne.0)then    
            do 11 j=ii,i-1        
               sum=sum-a(i,j)*b(j)     
 11         continue                  
         else if (sum.ne.0.0) then  
            ii=i                      
         endif                       
         b(i)=sum                    
 12   continue                      
      do 14 i=n,1,-1                
         sum=b(i)                    
         if(i.lt.n)then              
            do 13 j=i+1,n             
               sum=sum-a(i,j)*b(j)     
 13         continue                  
         endif                       
         b(i)=sum/a(i,i)             
 14   continue                      
      return   
      end      


c     *** ludcmp_dp.f ***

      subroutine ludcmp_dp(a,n,np,indx,d)

      implicit none

      integer,intent(in) :: n, np
      real*8,intent(inout) :: a(np,np)
      real*8,intent(out) :: d
      integer,intent(out) :: indx(n)
      
      integer nmax, i, j, k, imax
      real*8 tiny
      parameter (nmax=100,tiny=1.0d-20)   
      real*8 vv(nmax), aamax, sum, dum


      d=1.0d0
      do 12 i=1,n             
         aamax=0.0d0
         do 11 j=1,n           
            if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))             
 11      continue              
         if (aamax.eq.0.0) then
            write(*,*) 'ludcmp_dp: singular matrix!'
            stop
         endif            
         vv(i)=1.0d0/aamax  
 12   continue                
      do 19 j=1,n             
         if (j.gt.1) then      
            do 14 i=1,j-1       
               sum=a(i,j)        
               if (i.gt.1)then               
                  do 13 k=1,i-1               
                     sum=sum-a(i,k)*a(k,j)     
 13               continue    
                  a(i,j)=sum      
               endif             
 14         continue            
         endif                 
         aamax=0.0d0           
         do 16 i=j,n           
            sum=a(i,j)          
            if (j.gt.1)then     
               do 15 k=1,j-1                     
                  sum=sum-a(i,k)*a(k,j)                     
 15            continue              
               a(i,j)=sum            
            endif                   
            dum=vv(i)*abs(sum)      
            if (dum.ge.aamax) then  
               imax=i                
               aamax=dum             
            endif                   
 16      continue                  
         if (j.ne.imax)then        
            do 17 k=1,n             
               dum=a(imax,k)         
               a(imax,k)=a(j,k)      
               a(j,k)=dum            
 17         continue                
            d=-d                    
            vv(imax)=vv(j)          
         endif                     
         indx(j)=imax              
         if(j.ne.n)then            
            if(a(j,j).eq.0.0)a(j,j)=tiny               
            dum=1.0d0/a(j,j)      
            do 18 i=j+1,n           
               a(i,j)=a(i,j)*dum     
 18         continue                
         endif                     
 19   continue                    
      if(a(n,n).eq.0.0)a(n,n)=tiny               
      return                      
      end   


c     *** LUdec.f ***

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     
c     Solution of linear equation without inverting matrix
c     using LU decomposition: 
c     AA * xx = bb         AA, bb: known
c     xx: to be found
c     AA and bb are not modified in this subroutine
c     
c     MALV , Sep 2007
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine LUdec(xx,aa,bb,m,n)

      implicit none

!     Arguments 
      integer,intent(in) ::     m, n
      real*8,intent(in) ::      aa(m,m), bb(m)
      real*8,intent(out) ::     xx(m)


!     Local variables
      real*8      a(n,n), b(n), x(n), d
      integer    i, j, indx(n)      


!     Subrutinas utilizadas
!     ludcmp_dp, lubksb_dp

!!!!!!!!!!!!!!!Comienza el programa !!!!!!!!!!!!!!
      
      do i=1,n
         b(i) = bb(i+1)
         do j=1,n
            a(i,j) = aa(i+1,j+1)
         enddo
      enddo

                                ! Descomposicion de auxm1
      call ludcmp_dp ( a, n, n, indx, d)

                                ! Sustituciones foward y backwards para hallar la solucion
      do i=1,n
         x(i) = b(i)
      enddo
      call lubksb_dp( a, n, n, indx, x )

      do i=1,n
         xx(i+1) = x(i)
      enddo

      return
      end


c     *** mat_oper.f ***

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c     ***********************************************************************
      subroutine unit(a,n)
c     store the unit value in the diagonal of a 
c     ***********************************************************************
      implicit none
      real*8 a(n,n)
      integer n,i,j,k
      do 1,i=2,n-1
         do 2,j=2,n-1
	    if(i.eq.j) then
               a(i,j) = 1.d0
	    else
               a(i,j)=0.0d0
	    end if
 2       continue
 1    continue
      do k=1,n
         a(n,k) = 0.0d0
         a(1,k) = 0.0d0
         a(k,1) = 0.0d0
         a(k,n) = 0.0d0
      end do
      return
      end

c     ***********************************************************************
      subroutine diago(a,v,n)
c     store the vector v in the diagonal elements of the square matrix a
c     ***********************************************************************
      implicit none

      integer n,i,j,k
      real*8 a(n,n),v(n)

      do 1,i=2,n-1
         do 2,j=2,n-1
	    if(i.eq.j) then
               a(i,j) = v(i)
	    else
               a(i,j)=0.0d0
	    end if
 2       continue
 1    continue
      do k=1,n
         a(n,k) = 0.0d0
         a(1,k) = 0.0d0
         a(k,1) = 0.0d0
         a(k,n) = 0.0d0
      end do
      return
      end  

c     ***********************************************************************
      subroutine invdiag(a,b,n)
c     inverse of a diagonal matrix 
c     ***********************************************************************
      implicit none

      integer n,i,j,k
      real*8 a(n,n),b(n,n)

      do 1,i=2,n-1
         do 2,j=2,n-1
	    if (i.eq.j) then
               a(i,j) = 1.d0/b(i,i)
	    else
               a(i,j)=0.0d0
	    end if
 2       continue
 1    continue
      do k=1,n
         a(n,k) = 0.0d0
         a(1,k) = 0.0d0
         a(k,1) = 0.0d0
         a(k,n) = 0.0d0
      end do
      return
      end


c     ***********************************************************************
      subroutine samem (a,m,n)
c     store the matrix m in the matrix a 
c     ***********************************************************************
      implicit none
      real*8 a(n,n),m(n,n)
      integer n,i,j,k
      do 1,i=2,n-1
         do 2,j=2,n-1
            a(i,j) = m(i,j)  
 2       continue
 1    continue 	
      do k=1,n
         a(n,k) = 0.0d0
         a(1,k) = 0.0d0
         a(k,1) = 0.0d0
         a(k,n) = 0.0d0
      end do
      return
      end 


c     ***********************************************************************
      subroutine mulmv(a,b,c,n)
c     do a(i)=b(i,j)*c(j). a, b, and c must be distint
c     ***********************************************************************
      implicit none

      integer n,i,j
      real*8 a(n),b(n,n),c(n),sum

      do 1,i=2,n-1
         sum=0.0d0
         do 2,j=2,n-1
	    sum = sum + b(i,j) * c(j)
 2       continue
         a(i)=sum
 1    continue
      a(1) = 0.0d0
      a(n) = 0.0d0
      return
      end


c     ***********************************************************************
      subroutine trucodiag(a,b,c,d,e,n)
c     inputs: matrices b,c,d,e
c     output: matriz diagonal a
c     Operacion a realizar:  a = b * c^(-1) * d + e
c     La matriz c va a ser invertida
c     Todas las matrices de entrada son diagonales excepto b
c     Aprovechamos esa condicion para invertir c, acelerar el calculo, y 
c     ademas, para forzar que a sea diagonal
c     ***********************************************************************
      implicit none
      real*8 a(n,n),b(n,n),c(n,n),d(n,n),e(n,n), sum
      integer n,i,j,k
      do 1,i=2,n-1
         sum=0.0d0
         do 2,j=2,n-1
	    sum=sum+ b(i,j) * d(j,j)/c(j,j)
 2       continue
         a(i,i) = sum + e(i,i)
 1    continue
      do k=1,n
         a(n,k) = 0.0d0
         a(1,k) = 0.0d0
         a(k,1) = 0.0d0
         a(k,n) = 0.0d0
      end do
      return
      end


c     ***********************************************************************
      subroutine trucommvv(v,b,c,u,w,n)
c     inputs: matrices b,c , vectores u,w
c     output: vector v 
c     Operacion a realizar:  v = b * c^(-1) * u + w
c     La matriz c va a ser invertida 
c     c es diagonal, b no
c     Aprovechamos esa condicion para invertir c, y acelerar el calculo
c     ***********************************************************************
      implicit none
      real*8 v(n),b(n,n),c(n,n),u(n),w(n), sum
      integer n,i,j
      do 1,i=2,n-1
         sum=0.0d0
         do 2,j=2,n-1
	    sum=sum+ b(i,j) * u(j)/c(j,j)
 2       continue
         v(i) = sum + w(i)
 1    continue
      v(1) = 0.d0
      v(n) = 0.d0
      return
      end


c     ***********************************************************************
      subroutine sypvmv(v,u,c,w,n)
c     inputs: matriz diagonal c , vectores u,w
c     output: vector v 
c     Operacion a realizar:  v = u + c * w 
c     ***********************************************************************
      implicit none
      real*8 v(n),u(n),c(n,n),w(n)
      integer n,i
      do 1,i=2,n-1
         v(i)= u(i) + c(i,i) * w(i)
 1    continue
      v(1) = 0.0d0
      v(n) = 0.0d0
      return
      end


c     ***********************************************************************
      subroutine sumvv(a,b,c,n)
c     a(i)=b(i)+c(i)
c     ***********************************************************************
      implicit none

      integer n,i
      real*8 a(n),b(n),c(n)

      do 1,i=2,n-1
         a(i)= b(i) + c(i)
 1    continue
      a(1) = 0.0d0
      a(n) = 0.0d0
      return
      end


c     ***********************************************************************
      subroutine sypvvv(a,b,c,d,n)
c     a(i)=b(i)+c(i)*d(i)
c     ***********************************************************************
      implicit none
      real*8 a(n),b(n),c(n),d(n)
      integer n,i
      do 1,i=2,n-1
         a(i)= b(i) + c(i) * d(i)
 1    continue
      a(1) = 0.0d0
      a(n) = 0.0d0
      return
      end


c     ***********************************************************************
!      subroutine zerom(a,n)
c     a(i,j)= 0.0
c     ***********************************************************************
!      implicit none
!      integer n,i,j
!      real*8 a(n,n)

!      do 1,i=1,n
!         do 2,j=1,n
!	    a(i,j) = 0.0d0
! 2       continue
! 1    continue
!      return
!      end


c     ***********************************************************************
      subroutine zero4m(a,b,c,d,n)
c     a(i,j) = b(i,j) = c(i,j) = d(i,j) = 0.0 
c     ***********************************************************************
      implicit none
      real*8 a(n,n), b(n,n), c(n,n), d(n,n)
      integer n
      a(1:n,1:n)=0.d0
      b(1:n,1:n)=0.d0
      c(1:n,1:n)=0.d0
      d(1:n,1:n)=0.d0
!      do 1,i=1,n
!         do 2,j=1,n
!	    a(i,j) = 0.0d0
!	    b(i,j) = 0.0d0
!	    c(i,j) = 0.0d0
!	    d(i,j) = 0.0d0
! 2       continue
! 1    continue
      return
      end


c     ***********************************************************************
      subroutine zero3m(a,b,c,n)
c     a(i,j) = b(i,j) = c(i,j) = 0.0 
c     **********************************************************************
      implicit none
      real*8 a(n,n), b(n,n), c(n,n)
      integer n
      a(1:n,1:n)=0.d0
      b(1:n,1:n)=0.d0
      c(1:n,1:n)=0.d0
!      do 1,i=1,n
!         do 2,j=1,n
!	    a(i,j) = 0.0d0
!	    b(i,j) = 0.0d0
!	    c(i,j) = 0.0d0
! 2       continue
! 1    continue
      return
      end


c     ***********************************************************************
      subroutine zero2m(a,b,n)
c     a(i,j) = b(i,j) = 0.0 
c     ***********************************************************************
      implicit none
      real*8 a(n,n), b(n,n)
      integer n
      a(1:n,1:n)=0.d0
      b(1:n,1:n)=0.d0
!      do 1,i=1,n
!         do 2,j=1,n
!	    a(i,j) = 0.0d0
!	    b(i,j) = 0.0d0
! 2       continue
! 1    continue
      return
      end


c     ***********************************************************************
!      subroutine zerov(a,n)
c     a(i)= 0.0
c     ***********************************************************************
!      implicit none
!      real*8 a(n)
!      integer n,i
!      do 1,i=1,n
!         a(i) = 0.0d0
! 1    continue
!      return
!      end


c     ***********************************************************************
      subroutine zero4v(a,b,c,d,n)
c     a(i) = b(i) = c(i) = d(i,j) = 0.0
c     ***********************************************************************
      implicit none
      real*8 a(n), b(n), c(n), d(n)
      integer n
      a(1:n)=0.d0
      b(1:n)=0.d0
      c(1:n)=0.d0
      d(1:n)=0.d0
!      do 1,i=1,n
!         a(i) = 0.0d0
!         b(i) = 0.0d0
!         c(i) = 0.0d0
!         d(i) = 0.0d0
! 1    continue
      return
      end


c     ***********************************************************************
      subroutine zero3v(a,b,c,n)
c     a(i) = b(i) = c(i) = 0.0
c     ***********************************************************************
      implicit none
      real*8 a(n), b(n), c(n) 
      integer n
      a(1:n)=0.d0
      b(1:n)=0.d0
      c(1:n)=0.d0
!      do 1,i=1,n
!         a(i) = 0.0d0
!         b(i) = 0.0d0
!         c(i) = 0.0d0
! 1    continue
      return
      end


c     ***********************************************************************
      subroutine zero2v(a,b,n)
c     a(i) = b(i) = 0.0
c     ***********************************************************************
      implicit none
      real*8 a(n), b(n) 
      integer n
      a(1:n)=0.d0
      b(1:n)=0.d0
!      do 1,i=1,n
!         a(i) = 0.0d0
!         b(i) = 0.0d0
! 1    continue
      return
      end

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


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

c     *** suaviza.f ***

c*****************************************************************************
c     
      subroutine suaviza ( x, n, ismooth, y )
c     
c     x - input and return values 
c     y - auxiliary vector
c     ismooth = 0  --> no smoothing is performed
c     ismooth = 1  --> weak smoothing (5 points, centred weighted)
c     ismooth = 2  --> normal smoothing (3 points, evenly weighted)
c     ismooth = 3  --> strong smoothing (5 points, evenly weighted)


c     august 1991
c*****************************************************************************

      implicit none

      integer	n, imax, imin, i, ismooth
      real*8	x(n), y(n)
c*****************************************************************************

      imin=1
      imax=n

      if (ismooth.eq.0) then

         return 

      elseif (ismooth.eq.1) then ! 5 points, with central weighting 

         do i=imin,imax 
	    if(i.eq.imin)then
               y(i)=x(imin)
	    elseif(i.eq.imax)then
               y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0
	    elseif(i.gt.(imin+1) .and. i.lt.(imax-1) )then
               y(i) = ( x(i+2)/4.d0 + x(i+1)/2.d0 + 2.d0*x(i)/3.d0 + 
     @              x(i-1)/2.d0 + x(i-2)/4.d0 )* 6.d0/13.d0
	    else
               y(i)=(x(i+1)/2.d0+x(i)+x(i-1)/2.d0)/2.d0
	    end if
         end do	

      elseif (ismooth.eq.2) then ! 3 points, evenly spaced

         do i=imin,imax 
	    if(i.eq.imin)then
               y(i)=x(imin)
	    elseif(i.eq.imax)then
               y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0
	    else
               y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0
	    end if
         end do	
         
      elseif (ismooth.eq.3) then ! 5 points, evenly spaced

         do i=imin,imax 
	    if(i.eq.imin)then
               y(i) = x(imin)
	    elseif(i.eq.(imin+1) .or. i.eq.(imax-1))then
               y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0 
	    elseif(i.eq.imax)then
               y(i) = ( x(imax-1) + x(imax-1) + x(imax-2) ) / 3.d0
	    else
               y(i) = ( x(i+2)+x(i+1)+x(i)+x(i-1)+x(i-2) )/5.d0
	    end if
         end do	

      else

         write (*,*) ' Error in suaviza.f   Wrong ismooth value.'
         stop

      endif

c     rehago el cambio, para devolver x(i)
      do i=imin,imax 
         x(i)=y(i)
      end do

      return
      end


c     ***********************************************************************
      subroutine mulmmf90(a,b,c,n)
c     ***********************************************************************
      implicit none
      real*8 a(n,n), b(n,n), c(n,n)
      integer n

      a=matmul(b,c)
      a(1,:)=0.d0
      a(:,1)=0.d0
      a(n,:)=0.d0
      a(:,n)=0.d0

      return
      end


c     ***********************************************************************
      subroutine resmmf90(a,b,c,n)
c     ***********************************************************************
      implicit none
      real*8 a(n,n), b(n,n), c(n,n)
      integer n

      a=b-c
      a(1,:)=0.d0
      a(:,1)=0.d0
      a(n,:)=0.d0
      a(:,n)=0.d0

      return
      end


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

      subroutine gethist_03 (ihist) 

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

      implicit none

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


c     arguments 
      integer         ihist

c     local variables
      integer 	      j, r

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

      nbox = nbox_stored(ihist)
      do j=1,mm_stored(ihist)
         thist(j) = thist_stored(ihist,j) 
         do r=1,nbox_stored(ihist)
	    no(r) = no_stored(ihist,r) 
            sk1(j,r) = sk1_stored(ihist,j,r) 
            xls1(j,r) = xls1_stored(ihist,j,r)
            xld1(j,r) = xld1_stored(ihist,j,r)
         enddo
      enddo


      return
      end


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

      subroutine rhist_03 (ihist) 

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

      implicit none

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


c     arguments 
      integer         ihist

c     local variables
      integer 	      j, r
      real*8          xx

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

      open(unit=3,file=hisfile,status='old')

      read(3,*)
      read(3,*)
      read(3,*) mm_stored(ihist)
      read(3,*)
      read(3,*) nbox_stored(ihist)
      read(3,*)

      if ( nbox_stored(ihist) .gt. nbox_max ) then
         write (*,*) ' nbox too large in input file ', hisfile
         stop ' Check maximum number nbox_max in mz1d.par '
      endif

      do j=mm_stored(ihist),1,-1
         read(3,*) thist_stored(ihist,j) 
         do r=1,nbox_stored(ihist)
	    read(3,*) no_stored(ihist,r), 
     &           sk1_stored(ihist,j,r), 
     &           xls1_stored(ihist,j,r),
     &           xx,
     &           xld1_stored(ihist,j,r)
         enddo

      enddo

      close(unit=3)


      return
      end
