c $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/traceupr.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
	subroutine traceupr(z1,z2,a1,re,z,rfndx,nz,th2,a2,
     &losdist,rpl,mupl,nplpts,x,layindx)
c trace the ray from altitude z1 to z2 in the direction t1 at z1
c given by the sine of t1 (st1).
c
c input:
c z1 - initial alt.
c z2 - final alt.
c a1 - direction at z1
c
c output:
c th2 - angular distance traveled [radians]
c a2 - direction at z2
c
c losdist - segment distance
c rpl  -  radial distance from the center of the earth of the pts along the path
c mupl -  cosine of the local angle at pts along the path
c nplpts - # of path points
c x    -  angle between local zenith at A and local zenith at pts along the path
c         (radians)
c layindx - index of the layer for each segment
c
c
      implicit none
c the rcs revision #
      character*20 rcsrev
      parameter (rcsrev='$Revision: 2.22 $')
      include 'consts.cmn'
c passed
      integer nz,nplpts,layindx(2*nz)
      real*8 losdist(2*nz),rpl(2*nz+1),mupl(2*nz+1)
      real*8 x(nz*2+1)
      real*8 z1,z2,a1,re,rfndx(nz),th2,z(nz),a2
c local
      integer klo1,khi1,klo2,khi2,i,k
      real*8 repz1,c,cc,nz1,Hrf1,n01,Hrf2,n02,dtht,coem(3,3),xp(3)
      real*8 st1,st2,ds,nz2,dsN
      logical up,debug,lusedsN
c external
      real*8 rfintp
      external rfintp
c
      data debug/.FALSE./
      data lusedsN/.false./
      save debug,lusedsN
c
c check to see that z1 and z2 differ by amount larger than smallx
      if((z2-z1).gt.smallx)then
c
c in the upward direction
c
c compute the sine of the inital direction
c
        st1=sin(a1)
c
c find the index of refraction at the 2 altitudes
c find the layer of first index of refraction
c
c get n(r1)
        repz1=re+z1
c
c find the 2 pts that straddle z1
c
        call findex(z,nz,z1,klo1,khi1)
c check to make sure that z(klo1) is not equal to z1 (then we won't need this segment)
        if(klo1.ne.1)then
          if(abs(z1-z(klo1)).lt.smallx)then
            klo1=klo1-1
            khi1=khi1-1
          endif
        endif
        if(debug)write(6,*)'traceup. klo1 khi1',klo1,khi1
        if(debug)write(6,*)'traceup. zklo1 zkhi1',z(klo1),z(khi1)
        k=1
        mupl(k)=cos(a1)
        rpl(k)=re+z1
        x(k)=0.0d0
c
c calc the constant of the motion
c
        if(z1.gt.z(1))then
          nz1=1.0d0
          Hrf1=0.0d0
          n01=0.0d0
        else
          call rfscal(rfndx(klo1),rfndx(khi1),z(klo1),z(khi1),Hrf1,n01)
          nz1=1.0d0+rfintp(Hrf1,n01,z1)
        endif
c
c constant of the motion
        c=st1*nz1*repz1
        cc=c*c
c
c find the 2 pts that straddle z2
c
        call findex(z,nz,z2,klo2,khi2)
        if(z2.gt.z(1))then
          nz2=1.0d0
          Hrf2=0.0d0
          n02=0.0d0
        else
          call rfscal(rfndx(klo2),rfndx(khi2),z(klo2),z(khi2),Hrf2,n02)
          nz2=1.0d0+rfintp(Hrf2,n01,z2)
        endif
c
c check to see if z1 and z2 are within the same layer
c
        if(klo1.eq.klo2)then
c
c check to see if z1&z2 are below the top of the atmosphere
c
          if(((z1-z(1)).lt.smallx).and.((z2-z(1)).lt.smallx))then
            call soltht(z1,z2,re,Hrf1,n01,c,th2,xp,coem,st2,ds,dsN)
            if(lusedsN)then
              losdist(k)=dsN
            else
              losdist(k)=ds
            endif
            layindx(k)=klo1
          else if((z1-z(1)).lt.smallx)then
c from z1 to the top of atmo
            call soltht(z1,z(1),re,Hrf1,n01,c,th2,xp,coem,st2,ds,dsN)
            if(lusedsN)then
              losdist(k)=dsN
            else
              losdist(k)=ds
            endif
            layindx(k)=klo1
            k=k+1
            if(st2.lt.1.0d0)then
              mupl(k)=cos(asin(st2))
            else
              mupl(k)=0.0d0
            endif
            rpl(k)=re+z(1)
            x(k)=th2
c and then from top to z2
            call soltht(z(1),z2,re,Hrf1,n01,c,th2,xp,coem,st2,ds,dsN)
            if(lusedsN)then
              losdist(k)=dsN
            else
              losdist(k)=ds
            endif
            layindx(k)=klo1
          else
c from z1 to z2
            write(6,*)'h n',hrf1,n01,z1,z2
            call soltht(z1,z2,re,Hrf1,n01,c,th2,xp,coem,st2,ds,dsN)
            write(6,*)z1,z2,re,Hrf1,n01,c,th2,st2
            write(6,*)'traceup.3 st2 ',st2
            if(lusedsN)then
              losdist(k)=dsN
            else
              losdist(k)=ds
            endif
            layindx(k)=klo1
          endif
          if(debug)write(6,*)' ',z1,z2,th2/dr
        else
c
c start at layer klo1 and go up
c
c first segment from point to the top of the layer
          call soltht(z1,z(klo1),re,Hrf1,n01,c,th2,xp,coem,st2,ds,dsN)
          if(lusedsN)then
            losdist(k)=dsN
          else
            losdist(k)=ds
          endif
          layindx(k)=klo1
          k=k+1
          if(st2.lt.1.0d0)then
            mupl(k)=cos(asin(st2))
          else
            mupl(k)=0.0d0
          endif
          rpl(k)=re+z(klo1)
          x(k)=th2
          if(debug)write(6,*)' ',z1,z(klo1),th2/dr
c now from the top of the first layer to the top of the atmosphere
          do i=klo1-1,khi2,-1
             call rfscal(rfndx(i+1),rfndx(i),z(i+1),z(i),Hrf1,n01)
             call soltht(z(i+1),z(i),re,Hrf1,n01,c,dtht,xp,coem,st2,ds,
     &dsN)
             th2=th2+dtht
             if(debug)write(6,*)' ',z(i+1),z(i),dtht/dr
             layindx(k)=i
             if(lusedsN)then
               losdist(k)=dsN
             else
               losdist(k)=ds
             endif
             k=k+1
             st2=(c/(re+z(i))/(1.0d0+rfintp(Hrf1,n01,z(i))))
             if(st2.lt.1.0d0)then
               mupl(k)=cos(asin(st2))
             else
               mupl(k)=0.0d0
             endif
             rpl(k)=re+z(i)
             x(k)=th2
          enddo
c now the last segment from z(khi2) to z2
          call soltht(z(khi2),z2,re,Hrf2,n02,c,dtht,xp,coem,st2,ds,dsN)
          th2=th2+dtht
          layindx(k)=klo2
          if(lusedsN)then
            losdist(k)=dsN
          else
            losdist(k)=ds
          endif
          if(debug)write(6,*)' ',z(khi2),z2,dtht/dr
        endif
c compute sine(direction at z2)
        k=k+1
        if(st2.lt.1.0d0)then
          a2=asin(st2)
        else
          a2=pi_2
        endif
        mupl(k)=cos(a2)
        rpl(k)=re+z2
        x(k)=th2
      else
        th2=0.0d0
        a2=a1
        k=0
      endif
      nplpts=k
      if(debug)stop
      return
      end