C
C  FK
C  FK IS A FUNCTION THAT USES LINEAR INTERPOLATION TO FIND  FK,
C  GIVEN  ZK  AND THE ARRAYS  Z  AND  F. FOR EXAMPLE, IN ROUTINE
C  DERIVS,  BVF = FK(ALT,BVFZ,BVFF)  WILL BE USED TO GET  BVF  AT  ALT
C  FROM THE ARRAYS  BVFZ  AND  BVFF.
      REAL FUNCTION FK(ZK,Z,F,NZ,IZTOP)
      REAL Z(*), F(*)
C
      COMMON /IFLAGB/IFLAG
C
C  MAKE SURE Z(IZTOP) .GE. ZK
    5 CONTINUE
      IF(ZK.GT.Z(IZTOP)) THEN
         IZTOP = IZTOP + 1
         IF(IZTOP.GT.NZ) THEN
            WRITE(2,9) ZK,Z(NZ)
    9       FORMAT(' ZK TOO HIGH. ZK,Z(NZ) = ',2(1PE15.5))
            IZTOP = NZ
            FK = F(NZ)
            IFLAG = 1
            RETURN
CCC           CLOSE(2)
CCC           CLOSE(3)
CCC           STOP
         ENDIF
         GO TO 5
      ENDIF
C
      IF(ZK.EQ.Z(IZTOP)) THEN
         FK = F(IZTOP)
         RETURN
      ENDIF
C
C  FIND IZK SUCH THAT  Z(IZK) .GT. ZK .GE. Z(IZK-1)
      DO 100 IZ=IZTOP,2,-1
      IZM1 = IZ - 1
      IF(Z(IZ).GT.ZK .AND. ZK.GE.Z(IZM1)) THEN
         ZRATIO = (ZK - Z(IZM1)) / (Z(IZ) - Z(IZM1))
         FK = F(IZM1) + ZRATIO * (F(IZ) - F(IZM1))
         IZTOP = IZM1 + 1
         RETURN
      ENDIF
  100 CONTINUE
C
C  ZK TOO LOW
!     WRITE(2,109) ZK,Z(1)
c      WRITE(6,109) ZK,Z(1)
  109 FORMAT(' ZK TOO LOW. ZK,Z(1) = ',2(1PE15.5))
      IZTOP = 2
      FK = F(1)
      IFLAG = 1
      RETURN
CCC      CLOSE(2)
CCC      CLOSE(3)
CCC      STOP
C
      END
C
CCC      FUNCTION RTSAFE(EPS,FUNCD,X1,X2,XACC)
CCC      PARAMETER(MAXIT=100)
CCC      CALL FUNCD(EPS,X1,FL,DF)
CCC      CALL FUNCD(EPS,X2,FH,DF)
CCC      IF(FL*FH .GE. 0.) PAUSE 'root not bracketed'
CCC      IF(FL .LT. 0.) THEN
CCC         XL = X1
CCC         XH = X2
CCC      ELSE
CCC         XH = X1
CCC         XL = X2
CCC      ENDIF
CCC      RTSAFE = .5*(X1+X2)
CCC      DXOLD = ABS(X2-X1)
CCC      DX = DXOLD
CCC      CALL FUNCD(EPS,RTSAFE,F,DF)
CCC      DO 10 J=1,MAXIT
CCC      IF(((RTSAFE-XH)*DF-F)*((RTSAFE-XL)*DF-F) .GE. 0.
CCC X        .OR. ABS(2.*F) .GT. ABS(DXOLD*DF)) THEN
CCC            DXOLD = DX
CCC            DX = 0.5*(XH-XL)
CCC            RTSAFE = XL + DX
CCC            IF(XL.EQ.RTSAFE) RETURN
CCC         ELSE
CCC            DXOLD = DX
CCC            DX = F/DF
CCC            TEMP = RTSAFE
CCC            RTSAFE = RTSAFE - DX
CCC            IF(TEMP.EQ.RTSAFE) RETURN
CCC         ENDIF
CCC         IF(ABS(DX).LT.XACC) RETURN
CCC         CALL FUNCD(EPS,RTSAFE,F,DF)
CCC         IF(F .LT. 0.) THEN
CCC            XL = RTSAFE
CCC         ELSE
CCC            XH = RTSAFE
CCC         ENDIF
CCC   10    CONTINUE
CCC         PAUSE 'RTSAFE exceeding maximum iterations'
CCC         RETURN
CCC         END
c
CCC      subroutine funcd(eps,rtsafe,f,df)
CCC      f = exp(-.7*rtsafe)*rtsafe**0.25 - eps
CCC      df = (f + eps) * (-.7 +.25/rtsafe)
CCC      return
CCC      end