*-----------------------------------------------------------------------
*     G2IBL2
*-----------------------------------------------------------------------
*     Copyright (C) 2000-2016 GFD Dennou Club. All rights reserved.
*-----------------------------------------------------------------------
*     G2IBL2  Inverse of the bilinear interpolation of 2 params
*             (see G2FBL2)
*-----------------------------------------------------------------------
      SUBROUTINE G2IBL2(X,Y,X00,X10,X01,X11,Y00,Y10,Y01,Y11,
     &                  P,Q)

      EXTERNAL LREQ1, LRLE1, LRGE1
      LOGICAL  LREQ1, LRLE1, LRGE1

      A = X00 - X10 - X01 + X11
      E = Y00 - Y10 - Y01 + Y11

      IF ( LREQ1(A,0.0) .AND. LREQ1(E,0.0) ) THEN
        B = (-X00+X10)
        C = (-X00+X01)
        D = (X00-X)
        F = (-Y00+Y10)
        G = (-Y00+Y01)
        H = (Y00-Y)
        DET = B*G-C*F
        IF (LREQ1(DET,0.0)) THEN
          CALL MSGDMP('W','G2IBLI','SOLUTION NON EXSISTENT')
          P = 0
          Q = 0
        ELSE
          P = (-G*D+C*H)/DET
          Q = (F*D-B*H)/DET
        ENDIF
      ELSE IF( LREQ1(A,0.0) ) THEN
        B = (-X00+X10)
        C = (-X00+X01)
        D = (X00-X)
        F = (-Y00+Y10)/E
        G = (-Y00+Y01)/E
        H = (Y00-Y)/E
        IF ( LREQ1(B,0.0) .AND. LREQ1(C,0.0) ) THEN
          CALL MSGDMP('W','G2IBLI','SOLUTION NON EXSISTENT')
          P = 0
          Q = 0
        ELSE IF ( LREQ1(B,0.0) ) THEN
          Q = -D/C
          P = -(G*Q+H)/(Q+F)
        ELSE IF ( LREQ1(C,0.0) ) THEN
          P = -D/B
          Q = -(F*P+H)/(P+G)
        ELSE
          SQ = (C*F-B*G-D)**2 + 4*B*(C*H-D*G)
          IF ( LRLE1(SQ,0.0) ) SQ=0.0
          P1 = ( (C*F-B*G-D) + SQRT( SQ ) ) / (2*B)
          P2 = ( (C*F-B*G-D) - SQRT( SQ ) ) / (2*B)
          IF ( LRGE1(P1,0.0) .AND. LRLE1(P1,1.0) ) THEN
            P =P1
          ELSE IF ( LRGE1(P2,0.0) .AND. LRLE1(P2,1.0) ) THEN
            P = P2
          ELSE
            CALL MSGDMP('W','G2IBLI','NO SOLUTION IN THE DOMAIN')
            IF (ABS(P1).LT.ABS(P2)) THEN
              P = P1
            ELSE
              P = P2
            ENDIF
          ENDIF
          Q = -(B*P+D)/C
        ENDIF
      ELSE IF(LREQ1(E,0.0)) THEN
        B = (-X00+X10)/A
        C = (-X00+X01)/A
        D = (X00-X)/A
        F = (-Y00+Y10)
        G = (-Y00+Y01)
        H = (Y00-Y)
        IF (LREQ1(F,0.0) .AND. LREQ1(G,0.0)) THEN
          CALL MSGDMP('W','G2IBLI','SOLUTION NON EXSISTENT')
          P = 0
          Q = 0
        ELSE IF ( LREQ1(F,0.0) ) THEN
          Q = -H/G
          P = -(C*Q+D)/(Q+B)
        ELSE IF ( LREQ1(G,0.0) ) THEN
          P = -H/F
          Q = -(B*P+D)/(P+C)
        ELSE
          SQ = (G*B-F*C-H)**2 + 4*F*(G*D-H*C)
          IF ( LRLE1(SQ,0.0) ) SQ=0.0
          P1 = ( (G*B-F*C-H) + SQRT( SQ ) ) / (2*F)
          P2 = ( (G*B-F*C-H) - SQRT( SQ ) ) / (2*F)
          IF ( LRGE1(P1,0.0) .AND. LRLE1(P1,1.0) ) THEN
            P =P1
          ELSE IF ( LRGE1(P2,0.0) .AND. LRLE1(P2,1.0) ) THEN
            P = P2
          ELSE
            CALL MSGDMP('W','G2IBLI','NO SOLUTION IN THE DOMAIN')
            IF (ABS(P1).LT.ABS(P2)) THEN
              P = P1
            ELSE
              P = P2
            ENDIF
          ENDIF
          Q = -(F*P+H)/G
        ENDIF
      ELSE
        B = (-X00+X10)/A
        C = (-X00+X01)/A
        D = (X00-X)/A
        F = (-Y00+Y10)/E
        G = (-Y00+Y01)/E
        H = (Y00-Y)/E
        IF ( LREQ1(B-F,0.0) .AND. LREQ1(C-G,0.0) ) THEN
          CALL MSGDMP('W','G2IBLI','SOLUTION NON EXSISTENT')
          P = 0
          Q = 0
        ELSE IF ( LREQ1(B-F,0.0) ) THEN
          Q = (-D+H)/(C-G)
          P = -(C*Q+D)/(Q+B)
        ELSE IF( LREQ1(C-G,0.0) ) THEN
          P = -(D-H)/(B-F)
          Q = -(B*P+D)/(P+C)
        ELSE
          SQ = (C*F-B*G-D+H)**2 + 4*(B-F)*(C*H-D*G)
          IF ( LRLE1(SQ,0.0) ) SQ=0.0
          P1 = ( (C*F-B*G-D+H) + SQRT( SQ ) ) / (2*(B-F))
          P2 = ( (C*F-B*G-D+H) - SQRT( SQ ) ) / (2*(B-F))
          IF ( LRGE1(P1,0.0) .AND. LRLE1(P1,1.0) ) THEN
            P =P1
          ELSE IF ( LRGE1(P2,0.0) .AND. LRLE1(P2,1.0) ) THEN
            P = P2
          ELSE
            CALL MSGDMP('W','G2IBLI','NO SOLUTION IN THE DOMAIN')
            IF (ABS(P1).LT.ABS(P2)) THEN
              P = P1
            ELSE
              P = P2
            ENDIF
          ENDIF
          Q = -((B-F)*P+(D-H))/(C-G)
        ENDIF
      ENDIF

      END
