! Biliear interpolation using binary search
! Avoid masked points (e.g. Land ) during interpolation
! Author: Biju Thomas  on 2013/06/13 

MODULE hbin2dinpolmd
IMPLICIT NONE
  TYPE pair
     INTEGER   :: lo, up
  END TYPE pair
  CONTAINS

  SUBROUTINE hbin2dinpol(maskin,maskout,datain,dataout,   &
                              xo,yo,xn,yn,io,jo,in,jn)

! Doing  interpolation 2dim fields
! Author: Biju Thomas  on 2007/10/31 ....

  IMPLICIT NONE
    REAL, INTENT(IN)                     :: maskin,maskout
    INTEGER, INTENT(IN)                  :: io,jo,in,jn
    REAL, INTENT(IN), DIMENSION(io)      :: xo
    REAL, INTENT(IN), DIMENSION(jo)      :: yo
    REAL, INTENT(IN), DIMENSION(in)      :: xn
    REAL, INTENT(IN), DIMENSION(jn)      :: yn
    REAL, INTENT(IN), DIMENSION(io, jo)  :: datain
    REAL, INTENT(OUT),DIMENSION(in, jn)  :: dataout

    TYPE(pair),       DIMENSION(in)      :: ig
    TYPE(pair),       DIMENSION(jn)      :: jg 
      
    INTEGER                              :: i,j, i1,i2,j1,j2,ik,jk
    REAL                                 :: x, y, x1,x2,y1,y2
    REAL                                 :: sum, dxy
    REAL                                 :: flag(4)

    IF(yo(1) > yo(jo) .OR. xo(1) > xo(io)) THEN
      PRINT*, 'Input X, Y are not acending order'
      PRINT*,yo(1),yo(io),xo(1),xo(io)
      RETURN
    END IF

!   IF ( /= PRESENT(ig) ) THEN 
       CALL search_bin(xo,xn,io,in,ig)
!   ENDIF 
!   IF ( /= PRESENT(jg) ) THEN
       CALL search_bin(yo,yn,jo,jn,jg)
!   ENDIF 

      
       
    jloop: DO j = 1, jn
    iloop: DO i = 1, in
 
      i1 = ig(i)%lo
      i2 = ig(i)%up
      j1 = jg(j)%lo
      j2 = jg(j)%up
 
! j=1 is south end point and j=jn is north end point 
      IF (i1 == -99 .OR. i2 == -99 .OR.                  &
          j1 == -99 .OR. j2 == -99 ) THEN 
          dataout(i,j) = maskout
          CYCLE
      ENDIF
      x1 = xo(i1)
      x2 = xo(i2)
      y1 = yo(j1)
      y2 = yo(j2)

      x = xn(i)
      y = yn(j)

      IF ( x1 == x ) THEN
        x2 = x1
        i2 = i1
      END IF
      IF ( x1 /= x .AND. x2 == x ) THEN
        x1 = x2
        i1 = i2
      END IF
      IF ( y1 == y ) THEN
        y2 = y1
        j2 = j1
      END IF
      IF (y1 /= y .AND. y2 == y) THEN
        y1 = y2
        j1 = j2
      END IF

      IF(datain(i1,j1) == maskin) THEN
          flag(1)=0.
      ELSE
          flag(1)=1.
      ENDIF
      IF(datain(i2,j1) == maskin) THEN
          flag(2)=0.
      ELSE
          flag(2)=1.
      ENDIF
      IF(datain(i1,j2) == maskin) THEN
          flag(3)=0.
      ELSE
          flag(3)=1.
      ENDIF
      IF(datain(i2,j2) == maskin) THEN
          flag(4)=0.
      ELSE
          flag(4)=1.
      ENDIF

      IF( x2 /= x1 .AND. y2 /= y1) THEN
         sum=flag(1) + flag(2) + flag(3) + flag(4)
         IF(sum == 0.) THEN
           dataout(i,j)=maskin
         ELSE
           dxy = flag(1)*(x2-x)*(y2-y) + flag(2)*(x-x1)*(y2-y)       & 
               + flag(3)*(x2-x)*(y-y1) + flag(4)*(x-x1)*(y-y1)

           dataout(i,j) = flag(1)*datain(i1,j1)*(x2-x)*(y2-y)/dxy    &
                        + flag(2)*datain(i2,j1)*(x-x1)*(y2-y)/dxy    &     
                        + flag(3)*datain(i1,j2)*(x2-x)*(y-y1)/dxy    &
                        + flag(4)*datain(i2,j2)*(x-x1)*(y-y1)/dxy
         ENDIF
      ELSEIF( x2 /= x1 .AND. y2 == y1) THEN
         sum = flag(1) + flag(2) + flag(3) + flag(4)
         IF(sum == 0.) THEN
            dataout(i,j)=maskin
         ELSE
            dxy = flag(1)*(x2-x) + flag(2)*(x-x1) + flag(3)*(x2-x)  &
                                                  + flag(4)*(x-x1)
            dataout(i,j) = flag(1)*datain(i1,j1)*(x2-x)/dxy         &
                         + flag(2)*datain(i2,j1)*(x-x1)/dxy         &
                         + flag(3)*datain(i1,j2)*(x2-x)/dxy         &
                         + flag(4)*datain(i2,j2)*(x-x1)/dxy
         ENDIF
      ELSEIF( x2 == x1 .AND. y2 /= y1) THEN
         sum = flag(1) + flag(2) + flag(3) + flag(4)
         IF(sum == 0.) THEN
            dataout(i,j)=maskin
         ELSE
            dxy = flag(1)*(y2-y) + flag(2)*(y2-y) + flag(3)*(y-y1)  &
                                                  + flag(4)*(y-y1)
            dataout(i,j) = flag(1)*datain(i1,j1)*(y2-y)/dxy         &
                         + flag(2)*datain(i2,j1)*(y2-y)/dxy         &
                         + flag(3)*datain(i1,j2)*(y-y1)/dxy         &
                         + flag(4)*datain(i2,j2)*(y-y1)/dxy
         ENDIF
      ELSE
         sum = flag(1)
         IF(sum == 0.) THEN
            dataout(i,j)=maskin
         ELSE
           dataout(i,j) = datain(i1,j1)
         ENDIF
      ENDIF
      
    END DO iloop
    END DO jloop 

    RETURN
    END SUBROUTINE hbin2dinpol

  SUBROUTINE search_bin(yo,yn,jo,jn,jg)

! Doing binary searching  neighboring points.....
! Author: Biju Thomas  on 2013/05/23 ....

  IMPLICIT NONE
    INTEGER,   INTENT(IN)                 :: jo,jn
    TYPE(pair), DIMENSION(jn), INTENT(OUT):: jg
    REAL, DIMENSION(jo), INTENT(IN)       :: yo
    REAL, DIMENSION(jn), INTENT(IN)       :: yn
    INTEGER                               :: j, jk

    IF(yo(1) > yo(jo) .OR. yo(1) > yo(2) ) THEN
      PRINT*, 'Input X, Y are not acending order'
      PRINT*,yo(1),yo(jo)
      RETURN
    END IF
      
    jloop1: DO j=1,jn
       
       IF (yn(j) == yo(1)) THEN
           jg(j)%lo=1
           jg(j)%up=1
       ELSEIF (yn(j) == yo(jo)) THEN
           jg(j)%lo=jo
           jg(j)%up=jo
       ELSEIF (yn(j) < yo(1)) THEN
           jg(j)%lo=-99
           jg(j)%up=-99
       ELSEIF (yn(j) > yo(jo)) THEN
           jg(j)%lo=-99
           jg(j)%up=-99
       ELSE
           jg(j)%lo=1
           jg(j)%up=jo
           DO WHILE( jg(j)%up - jg(j)%lo > 1 )
              CALL binsrch( yo, jo, yn(j), jg(j) )
           ENDDO
       END IF
      
    END DO jloop1
      
  RETURN
  END  SUBROUTINE search_bin

  SUBROUTINE binsrch(yo, jo, y1, jg)

! Doing binary searching  neighboring points.....
! Author: Biju Thomas  on 2013/05/23 ....

  IMPLICIT NONE
    INTEGER,   INTENT(IN)          :: jo
    TYPE(pair), INTENT(INOUT)      :: jg
    REAL                           :: y1   
    REAL, DIMENSION(jo), INTENT(IN):: yo
    INTEGER                        :: j, jk, jm

    IF ( y1 > yo(jo) .OR. y1 < yo(1) ) THEN
        PRINT*, ' xinp is out of range ' 
        RETURN
    END IF

    IF ( jg%lo <= jg%up ) THEN
      jm = INT( 0.5 *( jg%lo + jg%up ) )
      IF ( yo(jm) >= y1 ) THEN
           jg%up = jm
      ELSE
           jg%lo = jm
      END IF
    ELSE
      PRINT*, ' data is not ordered '
    END IF
        
  RETURN
  END  SUBROUTINE binsrch
      
END MODULE hbin2dinpolmd