!//////////////////////////////////////////////////////////////
!/  ABEL INVERSION of bending angle to determin refraction
!/  profile. Refraction profile is used to derive density and
!/  density is used to determine Pressure via hydrostatics,
!/  Temperature is determined from the ideal gas law.
!/
!/  B.T. Marshall March, 2008
!//////////////////////////////////////////////////////////////    
CTM This routine performs a geometric inversion on bending angle
CTM to retrieve refractive index and density. It then retrieves
CTM T(P) from an assumed known upper boundary (TI,PI).
CTM The required bending angle derivative "Deltabend" is
CTM determined from Deltabend(i)=bend(i)/Scale(i) where ScaleH(i)
CTM was pre-determined using an exponential function fit to the
CTM bending angle data (bend=bend(i)*exp(deltaZ/ScaleH)).
CTM The routine "bend_fit" is used for the fitting.   
      subroutine abel_inversion_refraction(nl,nsr,ner,za,xlat,Rearth,
     &bend,scaleH,refr_indx,refractivity,density,pi,ti,zt,pt,tt,
     &i_top,z_top,i_bot,z_mid)
      implicit none
      real*4 P0,T0,MWair,Runiversal
      Parameter(P0=1013.25,T0=273.15)
      Parameter(MWair=28.964,Runiversal=8314.0)
      integer nl,nsr,ner,i,j,jj,jj1,jj2,num,num3,i_top,i_bot
      real*4 Rearth,ZA(nl),BEND(nl),ScaleH(nl),refractivity(nl),
     &density(nl),zt(nl),pt(nl),tt(nl),Factor,xlat,g,Rg,ALTITUDE,
     &dd,pi,ti,z_top,z_mid
C      real*4 denfit(nl),Scaleden(nl)
      REAL*8 pi_inverse,R,Rold,REFR_INDX(nl),Deltah,psum,dense,ScaleD,
     &SUM,DeltaZ,ZZZ,Deltabend,H,H1,W
      Data pi_inverse/0.3183099d0/
c      Data Factor/78.425/
CTM Factor*P0/T0 gives density in atm
C      Factor=Factor*P0/T0
CTM Factor*Runiversal/MWair gives density in 100Kg/M^3
CTM the 100 comes from pascals/mb this is so we can use p in mb
C
c       print*, za(nsr), za(ner),pi,ti
c      do i=1,nl
c        print*,'# ',za(i), bend(i), scaleH(i)
c      enddo
c        stop
c        ti=ti*2
c        bend(1) = bend(1) * 0.8

      num3=3
      Factor=78.425*Runiversal/MWair
      DeltaZ=1.0d-3
C      write(*,*) xlat,Rearth
      do i=1,ner                   ! geometric refr_indx inversion.
        R=DBLE(Rearth+ZA(i))
        SUM=0.0d0
        do j=1,i                   ! geometric refr_indx inversion.     
          ZZZ=ZA(J)+0.5*DeltaZ
          if(j.eq.1) THEN
CTM Top layer is 5 scale heights
            NUM=INT(ABS(5.0*ScaleH(j)/DeltaZ))
            H=DBLE(ScaleH(j))
            DO jj=1,num
              Rold=DBLE(Rearth+ZZZ)
              Deltabend=DeltaZ*bend(j)*EXP(-(ZZZ-ZA(j))/H)/H
              sum=sum+log(Rold/R+SQRT((Rold/R)**2-1.0d0))*Deltabend
              ZZZ=ZZZ+DeltaZ
            ENDDO
          ELSE
            NUM=INT(ABS((ZA(J-1)-ZA(J))/DeltaZ))
C            H=DBLE(0.5d0*(ScaleH(j)+ScaleH(j-1)))
            H=DBLE(ScaleH(J))
            H1=DBLE(ScaleH(J-1))
            DO jj=1,num
              W=1.0d0-((ZZZ-ZA(J))/(ZA(J-1)-ZA(J)))
              Rold=DBLE(Rearth+ZZZ)
              Deltabend=DeltaZ*((W*bend(j)*EXP(-(ZZZ-ZA(j))/H)/H)
     &                 +((1.0d0-W)*bend(j-1)*EXP(-(ZZZ-ZA(j-1))/H1)/H1))
c         if(i .eq. 510) print*, i, j, deltabend, bend(j), bend(j-1) 
              sum=sum+log(Rold/R+SQRT((Rold/R)**2-1.0d0))*Deltabend
              ZZZ=ZZZ+DeltaZ
            ENDDO
          ENDIF
c          print*, i, sum
C          DO jj=1,num
C            Rold=DBLE(Rearth+ZZZ)
C            Deltabend=DeltaZ*bend(j)*EXP(-(ZZZ-ZA(j))/H)/H
C            sum=sum+log(Rold/R+SQRT((Rold/R)**2-1.0d0))*Deltabend
C            ZZZ=ZZZ+DeltaZ
C          ENDDO
        enddo
c         stop
        refr_indx(i)=exp(sum*pi_inverse)
        rold=R/refr_indx(i)
c       print*, i, refr_indx(i)
        zt(i)=real(rold-Rearth)
        refractivity(i)=1.0e6*real(refr_indx(i)-1.0d0)
        density(i)=refractivity(i)/Factor
      enddo
c       stop
C      call bend_fit(nl,za,density,num3,num3,za(nsr),za(ner),denfit,
C     &scaleden)
      dd=MWair*PI/TI/Runiversal
      pt(nsr)=PI*density(nsr)/dd
      tt(nsr)=TI
      psum=pt(nsr)
CTM Calculate P and T 
      do i=nsr+1,ner
CTM Altitude is at layer middle and Density is for the layer not tangent point
CTM Use density scale height ScaleH for integration of density to get pressure.
CTM Doing this gives each layer contributing:
CTM gave*density(i)*ScaleH(i)*(exp((z(i-1)-z(i))/ScaleH(i))-1.0)
        IF(i.eq.1) then
          Deltah=DBLE(zt(i)-zt(i+1))
          ALTITUDE=zt(i)+0.5*(zt(i)-zt(i+1))
          ScaleD=DBLE(-ScaleH(i))
        else
          Deltah=DBLE(zt(i-1)-zt(i)) 
          ALTITUDE=zt(i)+0.5*(zt(i-1)-zt(i))
          ScaleD=DBLE(-0.5d0*(ScaleH(i)+ScaleH(i-1)))
        endif
CTM Covert Deltah to meters
        dense=1000.0d0*dble(density(i))
        call radgrav(ALTITUDE,xlat,g,Rg)
        psum=psum+DBLE(g)*dense*ScaleD*(DEXP(DeltaH/ScaleD)-1.0d0)
        pt(i)=real(psum)
CTM Calculate T from ideal gas law.
        tt(i)=MWair*pt(i)/Runiversal/density(i)
      enddo
CTM convert density to kg/m^3, needs factor of 100
      do i=1,nl
        density(i)=density(i)*100.0
c      write(*,101) za(i),zt(i),bend(i),refr_indx(i),refractivity(i),
c     & density(i),pt(i),tt(i)
      enddo
      z_top = zt(i_top)
      z_mid = 0.5*(zt(i_top)+zt(i_bot))
      write(*,*) i_top,i_bot,z_top,z_mid
c        print*, '# refrac '
c       do i = nsr,ner
c          print*, tt(i),zt(i)
c       enddo
c       print*,'&'
c      stop
      return
 101  Format(F9.4,F9.4,E14.6,F16.12,F12.6,2E14.6,F9.4)
      end