C MARCH 5, 2001; MARCH 24, 2005
C FUNCTION DECLARATION FOR THE APPROXIMATE INVERSE OF THE DERIVATIVE
C FUNCTION SPECIFIED IN  DC(I,F).
C IN PARTICULAR, FOR A GIVEN ARC I AND T, IT COMPUTES AN F SATISFYING
C            -EPS3 <= DC(I,F)-T <= EPS3,
C WHERE EPS3 < EPS2 IS SPECIFIED IN THE CALLING PROGRAM
C AND PASSED VIA COMMON BLOCK. 
C THE FUNCTION USES F(I) AND COST(I) AND SURPI FOR REFERENCE
      
      FUNCTION DCINV(I,T)
      IMPLICIT NONE
      INTEGER*4 MAXARCS,MAXNDAT,I,N,NA
      PARAMETER (MAXARCS=20000, MAXNDAT=5)
      REAL*8 A(MAXARCS),CDAT(MAXARCS,MAXNDAT),F(MAXARCS),
     $       U(MAXARCS),COST(MAXARCS)
      REAL*8 T,DC,DCINV,LARGE,ZERO,EPS3
      REAL*8 DERIV,DU,D0,DA,DB,DE,FA,FB,FE,FD,C1,C2
      COMMON /SCALARS/N,NA,LARGE,ZERO
      COMMON /LCOS/A
      COMMON /QCO1/CDAT
      COMMON /UBOUND/U
      COMMON /COST/COST
      COMMON /FLOW/F
      COMMON /EPS3/EPS3

C USING F(I) AS REFERENCE, FIND FA & FB THAT BRACKET DESIRED F AND 
C WHERE DERIVATIVES ARE FINITE.  THE LATTER IS NEEDED FOR INTERPOLATING.

c	print*,'dcinv0',cost(i),t,eps3,f(i),u(i)

      DERIV=COST(I)
      
c      print*,'2',deriv-t
      
      IF (DERIV-T.LT.-EPS3) THEN 
        DU=DC(I,U(I))
        IF (DU-T.GT.EPS3) THEN
          FA=F(I)
          DA=DERIV
          FB=U(I)
          IF (DU.EQ.LARGE) THEN
            FB=0.5*F(I)+0.5*U(I)
            DB=DC(I,FB)
1           IF (DB-T.LE.EPS3) THEN
              FB=0.5*FB+0.5*U(I)
              DB=DC(I,FB)
              GO TO 1
            END IF
          ELSE
            DB=DU
          END IF
          
c            print*,'1: a,b,da-t,db-t=',fa,fb,da-t,db-t
            
        ELSE
          DCINV=U(I)
          RETURN
        END IF
      ELSE 

c      print*,'3',deriv-t
      
        IF (DERIV-T.GT.EPS3) THEN
          D0=DC(I,0.D0)
          
c          print*,d0
          
          IF (D0-T.LT.-EPS3) THEN
            FB=F(I)
            DB=DERIV
            FA=0.0
            IF (D0.EQ.-LARGE) THEN 
              FA=0.5*F(I)
              DA=DC(I,FA)
2             IF (DA-T.GE.-EPS3) THEN
                FA=0.5*FA
                DA=DC(I,FA)
                GO TO 2
              END IF
            ELSE
              DA=D0
            END IF
          
c            print*,'2: a,b,da-t,db-t=',fa,fb,da-t,db-t
            
          ELSE
            DCINV=0.0
            RETURN
          END IF
        ELSE
          DCINV=F(I)
          RETURN
        END IF
      END IF

C LINEAR INTERPOLATE BETWEEN FA & FB TO OBTAIN FE  
    
      FE=FA-(DA-T)*(FB-FA)/(DB-DA)
      DE=DC(I,FE)
      
c      print*,'dcinv0.5=',fe,DE-T,eps3
            
3     IF ((DE-T.LT.-EPS3).OR.(DE-T.GE.EPS3)) THEN

c      print*,'dcinv1=',fA,fe,fB,DA-t,DE-t,DB-t


C QUADRATIC INTERPOLATE BETWEEN FA & FB & FE TO OBTAIN FD 

        C2=( (DB-DA)/(FB-FA) - (DE-DA)/(FE-FA) )/(FB-FE)
        C1=(DB-DA)/(FB-FA) - (FB-FA)*C2
        
C SOLVE QUADRATIC  Q(F) = DA + C1*(F-FA) + C2*(F-FA)**2 

        IF (C2.EQ.0.0) THEN
        
C CAN THIS CASE EVER OCCUR?      
        
          FD=FA-(DA-T)*(FB-FA)/(DB-DA)
        ELSE
          FD=FA + (DSQRT(C1**2-4*C2*(DA-T))-C1)/(2*C2)
          
c        print*,'D,c1,c2=',FD,c1,c2
c        print*,FA - 2*(DA-T)/(DSQRT(C1**2-4*C2*(DA-T))+C1)
        
          
        END IF

        IF (FD.LT.FE) THEN
          FB=FE
          DB=DE
        ELSE
          FA=FE
          DA=DE
        END IF
        FE=FD
        DE=DC(I,FE)
        
c        print*,'de-t=',de-t
        
        GO TO 3
      END IF
      DCINV=FE
      
c      print*,'type 1 to continue'
c      read(*,*)fe
      

      RETURN
      END