      SUBROUTINE XYIND(XIND,YIND,XGI,XGO,YGI,YGO,II,IO,JI,JO,IFX,IFY)
      DIMENSION XIND(IO),YIND(JO),XGI(II),XGO(IO),YGI(JI),YGO(JO)

      IF (IFX.GE.0) THEN
         IF (XGO(1).LT.XGI(1).OR.XGO(IO).GT.XGI(II)) THEN
            PRINT 3
    3       FORMAT(' IMPROPER X GRID DOMAIN IN XYIND')
            PRINT 4,XGI
    4       FORMAT(' XGI=',1P10E12.5)
            PRINT 5,XGO
    5       FORMAT(' XGO=',1P10E12.5)
         ENDIF
      ENDIF
      IF (IFY.GE.0) THEN
         IF (YGO(1).LT.YGI(1).OR.YGO(JO).GT.YGI(JI)) THEN
            PRINT 6
    6       FORMAT(' IMPROPER Y GRID DOMAIN IN XYIND')
            PRINT 7,YGI
    7       FORMAT(' YGI=',1P10E12.5)
            PRINT 8,YGO
    8       FORMAT(' YGO=',1P10E12.5)
         ENDIF
      ENDIF

      IF (IFX.LT.0) THEN
         RATX=(FLOAT(II)-1.)/(FLOAT(IO)-1.)
         XIND(1)=1.
         DO 10 I=2,IO
         XIND(I)=FLOAT(I-1)*RATX+1.
   10    CONTINUE
      ELSE
         I1=1
         DO 40 I=1,IO
   20    CONTINUE
         IF (XGI(I1+1).GE.XGO(I)) GO TO 30
         I1=I1+1
         IF (I1.GT.II) THEN
            PRINT 22,I1,II
   22       FORMAT(' XIND SEARCH OUT OF BOUNDS IN XYIND,I1,II=',2I4)
            PRINT 4,XGI
            PRINT 5,XGO
         ENDIF
         GO TO 20
   30    CONTINUE
         XIND(I)=(XGO(I)-XGI(I1))/(XGI(I1+1)-XGI(I1))+I1
	 print *,'xgo(',I,')= ',xgo(i),', xgi(',I,')= ',xgi(i)
   40    CONTINUE
      ENDIF

      IF (IFY.LT.0) THEN
         RATY=(FLOAT(JI)-1.)/(FLOAT(JO)-1.)
         YIND(1)=1.
         DO 50 J=2,JO
         YIND(J)=FLOAT(J-1)*RATY+1.
   50    CONTINUE
      ELSE
         J1=1
         IF (J1.GT.JI) THEN
            PRINT 52,J1,JI
   52       FORMAT(' YIND SEARCH OUT OF BOUNDS IN XYIND,I1,II=',2I4)
            PRINT 7,YGI
            PRINT 8,YGO
         ENDIF
         DO 80 J=1,JO
   60    CONTINUE
         IF (YGI(J1+1).GE.YGO(J)) GO TO 70
         J1=J1+1
         GO TO 60
   70    CONTINUE
         YIND(J)=(YGO(J)-YGI(J1))/(YGI(J1+1)-YGI(J1))+J1
   80    CONTINUE
      ENDIF
      RETURN
      END

      SUBROUTINE LININT(Z,NZ,F,NF,XIND,YIND,II,IO,JI,JO)
      DIMENSION Z(NZ,JI),F(NF,JO),XIND(IO),YIND(JO)

      DO J=1,JO
        JZ=INT(YIND(J))
        JZP1=MIN0(JZ+1,JI)
        Q=YIND(J)-JZ
        DO I=1,IO
          IZ=INT(XIND(I))
          IZP1=MIN0(IZ+1,II)
          P=XIND(I)-IZ
          F(I,J)=(1.-Q)*((1.-P)*Z(IZ,JZ)+P*Z(IZP1,JZ))+
     &       Q*(P*Z(IZP1,JZP1)+(1.-P)*Z(IZ,JZP1))
	enddo
      enddo

      RETURN
      END
