! Inverse Distance Square interpolation
! Avoid masked points (e.g. Land ) during interpolation
! Author: Biju Thomas  on 2016/01/12 

MODULE invdistinpmd
IMPLICIT NONE
  INTEGER,  PARAMETER :: flagint = -999999
  INTEGER :: istep = 0
  REAL, PARAMETER :: flagreal = 1.0e29, noinpval = 0.0
  REAL, PARAMETER :: pi = 3.141592653589
  REAL, PARAMETER :: d2rad = pi/180.0
  REAL, PARAMETER :: pow = 2.0
!  REAL, PARAMETER :: eps = 8.0e-02
  REAL, PARAMETER :: eps = 1.0e-01

  CONTAINS

  SUBROUTINE invdistinp(maskin,maskout,fin,fout,   &
                       i1, j1, i2, j2,w11, w12, w21, w22, nx,ny)

! Doing  interpolation 2dim fields
! Author: Biju Thomas  on 2015/01/12 ....

  IMPLICIT NONE
    REAL, INTENT(IN) :: maskin,maskout
    INTEGER, INTENT(IN) :: nx, ny, i1, j1, i2, j2
    REAL, INTENT(IN) :: w11, w12, w21, w22
    REAL, INTENT(IN), DIMENSION(nx,ny) :: fin
    REAL, INTENT(OUT) :: fout
    REAL :: wt11, wt12, wt21, wt22
    INTEGER :: i, j
    REAL :: sum
      
    IF ( w11 == noinpval ) THEN
      fout = fin(i1,j1)
      RETURN
    ELSEIF ( w12 == noinpval ) THEN
      fout = fin(i1,j2)
      RETURN
    ELSEIF ( w21 == noinpval ) THEN
      fout = fin(i2,j1)
      RETURN
    ELSEIF ( w22 == noinpval ) THEN
      fout = fin(i2,j2)
      RETURN
    ENDIF

    IF ( w11 == flagreal ) THEN
      wt11 = 0.0
    ELSE
      wt11 = w11
    ENDIF
    IF ( w12 == flagreal ) THEN
      wt12 = 0.0
    ELSE
      wt12 = w12
    ENDIF
    IF ( w21 == flagreal ) THEN
      wt21 = 0.0
    ELSE
      wt21 = w21
    ENDIF
    IF ( w22 == flagreal ) THEN
       wt22 = 0.0
    ELSE
       wt22 = w22 
    ENDIF

    sum = wt11 + wt12 + wt21 + wt22
    IF ( sum /= 0.0 ) THEN
       fout = ( wt11*fin(i1,j1) + wt12*fin(i1,j2) +                    &
                wt21*fin(i2,j1) + wt22*fin(i2,j2) ) / sum
    ELSE
       fout = maskout
    ENDIF

    
    RETURN
    END SUBROUTINE invdistinp

 
  FUNCTION monotonic (nn, g1d) RESULT(check)
! Check grid points are monotonically increaing
! Author: Biju Thomas  on 2016/01/12 ....

  IMPLICIT NONE
    INTEGER, INTENT(IN) :: nn
    REAL, DIMENSION(nn), INTENT(IN) :: g1d
    LOGICAL :: check
    INTEGER :: k

    check = .TRUE.
    DO k = 1, nn-1
      IF ( g1d(k) > g1d(k+1) ) THEN
        check = .FALSE.
        RETURN
      ENDIF
    ENDDO 
  END FUNCTION

  SUBROUTINE nearpts(nx, ny, g2dx, g2dy, xi, yi, is, ie,            &
                     js, je, fin, ilo, jlo, iup, jup, maskin )
! Doing  searching  4  neighboring points.....
! Author: Biju Thomas  on 2016/01/12 ....
! g2dx, g2dy need to monotonically increaing

  IMPLICIT NONE
    INTEGER, INTENT(IN) :: nx, ny, is, ie, js, je
    INTEGER, INTENT(OUT) :: ilo, jlo, iup, jup
    REAL, INTENT(IN) :: xi, yi, maskin
    REAL, DIMENSION(nx, ny), INTENT(IN) :: fin, g2dx, g2dy
    REAL :: d11, d12, d21, d22, mindist 
    REAL :: nbrpts(4)
    INTEGER :: i, j, m11, m12, m21, m22, indx
    INTEGER :: ilo1, jlo1, iup1, jup1
    LOGICAL :: loopexit
    LOGICAL :: asearch = .FALSE.
    
    IF ( istep == 0 ) THEN
      DO j = 1, ny
        IF ( .NOT. monotonic (nx, g2dx(:,j)) )  THEN
          PRINT*, 'INVALID INPUT GRID(lono: invdistinp) EXIT', j
          RETURN
        ENDIF
      ENDDO
      DO i = 1, nx
        IF ( .NOT. monotonic (ny, g2dy(i,:)) )  THEN
          PRINT*, 'INVALID INPUT GRID(lato: invdistinp) EXIT', i
          RETURN
        ENDIF
      ENDDO
    ENDIF
    istep =  1
    IF ( is > ie .OR. js > je) THEN
      PRINT*, 'INVALID is, ie, js, je (invdistinp) EXIT', is, ie, js, je
      RETURN
    ENDIF
    loopexit = .FALSE.
    ilo = flagint 
    jlo = flagint
    iup = flagint
    jup = flagint
    DO j = js, je-1
    DO i = is, ie-1
      IF ( xi >= g2dx(i,j) .AND. xi <= g2dx(i+1,j) .AND.       &
           yi >= g2dy(i,j) .AND. yi <= g2dy(i,j+1) ) THEN
        ilo = i
        jlo = j
        iup = i+1
        jup = j+1
        loopexit = .TRUE.
        EXIT
      ENDIF
    ENDDO
    IF ( loopexit ) EXIT
    ENDDO
    IF ( ilo == flagint .OR. jlo == flagint .OR. iup == flagint .OR. &
         jup == flagint ) THEN
      DO j = js, je
      DO i = is, ie
        IF ( xi >= g2dx(i,j)-eps .AND. xi <= g2dx(i+1,j)+eps .AND.       &
             yi >= g2dy(i,j)-eps .AND. yi <= g2dy(i,j+1)+eps ) THEN
          ilo = i
          jlo = j
          iup = i+1
          jup = j+1
          loopexit = .TRUE.
          EXIT
        ENDIF
      ENDDO
      IF ( loopexit ) EXIT
      ENDDO
    ENDIF
    IF ( ilo == flagint .OR. jlo == flagint .OR. iup == flagint .OR. jup == flagint ) THEN
      PRINT*,'FAILED TO FIND NEAR POINTS', ilo, iup, jlo, jup, xi, yi
      PRINT*,'EXIT(SUB: nearpts <= invdistinp.f90)'
      RETURN
    ENDIF
    IF ( .NOT. asearch ) RETURN
    m11 = 1
    m12 = 1
    m21 = 1
    m22 = 1
    IF ( fin(ilo,jlo) == maskin ) m11 = 0
    IF ( fin(ilo,jup) == maskin ) m12 = 0
    IF ( fin(iup,jlo) == maskin ) m21 = 0
    IF ( fin(iup,jup) == maskin ) m22 = 0

    IF ( m11 +  m12 + m21 + m22 == 1 ) THEN 
      loopexit = .FALSE.
      ilo1 = flagint
      jlo1 = flagint
      iup1 = flagint
      jup1 = flagint
      DO j = MAX(2,js), MAX(je-2,1)
      DO i = MAX(2,is), MAX(ie-2,1)
        IF ( xi >= g2dx(i-1,j) .AND. xi <= g2dx(i+2,j) .AND.       &
             yi >= g2dy(i,j-1) .AND. yi <= g2dy(i,j+2) ) THEN
          ilo1 = i-1
          jlo1 = j-1
          iup1 = i+2
          jup1 = j+2
          loopexit = .TRUE.
          EXIT
        ENDIF
      ENDDO
      IF ( loopexit ) EXIT
      ENDDO
      d11 = maskin
      d12 = maskin
      d21 = maskin
      d22 = maskin
      IF ( ilo1 /= flagint .AND. jlo1 /= flagint ) THEN
        IF ( fin(ilo1,jlo1) /= maskin )                             &
             d11 = gcirdist(g2dx(ilo1,jlo1), g2dy(ilo1,jlo1), xi, yi)
      ENDIF 
      IF ( ilo1 /= flagint .AND. jup1 /= flagint ) THEN
        IF ( fin(ilo1,jup1) /= maskin )                             &
             d12 = gcirdist(g2dx(ilo1,jup1), g2dy(ilo1,jup1), xi, yi)
      ENDIF
      IF ( iup1 /= flagint .AND. jlo1 /= flagint ) THEN
        IF ( fin(iup1,jlo1) /= maskin )                             &
            d21 = gcirdist(g2dx(iup1,jlo1), g2dy(iup1,jlo1), xi, yi)
      ENDIF
      IF ( iup1 /= flagint .AND. jup1 /= flagint ) THEN
         IF ( fin(iup1,jup1) /= maskin )                            &
            d22 = gcirdist(g2dx(iup1,jup1), g2dy(iup1,jup1), xi, yi)
         ENDIF
      IF ( m11 == 1 ) THEN
         nbrpts = (/d11, d12, d21, d22/)
         CALL indxminv(4, nbrpts, indx, mindist)
         IF ( mindist > 0 .AND. mindist /= maskin ) THEN
            IF ( indx == 1 ) THEN
               iup = ilo1
               jup = jlo1
            ELSE IF ( indx == 2 ) THEN
               iup = ilo1
               jup = jup1
            ELSE IF ( indx == 3 ) THEN
               iup = iup1
               jup = jlo1
            ELSE IF ( indx == 4 ) THEN
               iup = iup1
               jup = jup1
            END IF
            IF ( fin(iup,jup) == maskin ) PRINT*, 'ERROR IN REMAP1 (SUBR: nearpts)'
         END IF
      ELSE IF ( m12 == 1 ) THEN
         nbrpts = (/d11, d12, d21, d22/)
         CALL indxminv(4, nbrpts, indx, mindist)
         IF ( mindist > 0 .AND. mindist /= maskin ) THEN
            IF ( indx == 1 ) THEN
               iup = ilo1
               jup = jlo1
            ELSE IF ( indx == 2 ) THEN
               iup = ilo1
               jup = jup1
            ELSE IF ( indx == 3 ) THEN
               iup = iup1
               jup = jlo1
            ELSE IF ( indx == 4 ) THEN
               iup = iup1
               jup = jup1
            END IF
            IF ( fin(iup,jup) == maskin ) PRINT*, 'ERROR IN REMAP2 (SUBR: nearpts)'
         END IF
      ELSE IF ( m21 == 1 ) THEN
         nbrpts = (/d11, d12, d21, d22/)
         CALL indxminv(4, nbrpts, indx, mindist)
         IF ( mindist > 0 .AND. mindist /= maskin ) THEN
            IF ( indx == 1 ) THEN
               iup = ilo1
               jup = jlo1
            ELSE IF ( indx == 2 ) THEN
               iup = ilo1
               jup = jup1
            ELSE IF ( indx == 3 ) THEN
               iup = iup1
               jup = jlo1
            ELSE IF ( indx == 4 ) THEN
               iup = iup1
               jup = jup1
            END IF
            IF ( fin(iup,jup) == maskin ) PRINT*, 'ERROR IN REMAP3 (SUBR: nearpts)'
         END IF
      ELSE IF ( m22 == 1 ) THEN
         nbrpts = (/d11, d12, d21, d22/)
         CALL indxminv(4, nbrpts, indx, mindist)
         IF ( mindist > 0 .AND. mindist /= maskin ) THEN
            IF ( indx == 1 ) THEN
               ilo = ilo1
               jlo = jlo1
            ELSE IF ( indx == 2 ) THEN
               ilo = ilo1
               jlo = jup1
            ELSE IF ( indx == 3 ) THEN
               ilo = iup1
               jlo = jlo1
            ELSE IF ( indx == 4 ) THEN
               ilo = iup1
               jlo = jup1
            END IF
            IF ( fin(ilo,jlo) == maskin ) PRINT*, 'ERROR IN REMAP4 (SUBR: nearpts)'
         END IF
      END IF
    END IF
  END SUBROUTINE  nearpts
  
  SUBROUTINE weights(ii, jj, nx, ny, g2dx, g2dy, xi, yi, fin,     &
            ilo, jlo, iup, jup, w11, w12, w21, w22, maskin)
! Calculating weights.....
! Author: Biju Thomas  on 2016/01/12 ....
! g2dx, g2dy need to monotonically increaing

  IMPLICIT NONE 
    INTEGER, INTENT(IN) :: ii, jj, nx, ny
    INTEGER, INTENT(INOUT) :: ilo, jlo, iup, jup
    REAL, INTENT(IN) :: xi, yi, maskin 
    REAL, INTENT(OUT) :: w11, w12, w21, w22
    REAL, DIMENSION(nx, ny), INTENT(IN) :: fin, g2dx, g2dy
    INTEGER, PARAMETER :: mnipts = 1
    INTEGER :: mlo, mup, nlo, nup
    REAL :: nbrpts(3)
    REAL :: d1, d2, d3, mindist, msum
    REAL :: ds(12)
    INTEGER :: i, j, iix, jjy, m11, m12, m21, m22, indx
    LOGICAL :: asearch = .FALSE.

    IF ( ilo == flagint .OR. jlo == flagint .OR. iup == flagint .OR. jup == flagint ) THEN
     PRINT*,'FAILED TO FIND NEAR POINTS', ilo, iup, jlo, jup, xi, yi
     PRINT*,'EXIT(invdistinp.f90)'
     RETURN
    ENDIF
    
    IF ( fin(ilo,jlo) == maskin ) THEN
       w11 =  flagreal
       m11 = 0
    ELSE
       w11 = invdiswt(g2dx(ilo,jlo), g2dy(ilo,jlo), xi, yi, pow)
       m11 = 1       
    ENDIF
    IF ( fin(ilo,jup) == maskin ) THEN
       w12 =  flagreal
       m12 = 0
    ELSE
       w12 = invdiswt(g2dx(ilo,jup), g2dy(ilo,jup), xi, yi, pow)
       m12 = 1
    ENDIF
    IF ( fin(iup,jlo) == maskin ) THEN
       w21 =  flagreal
       m21 = 0
    ELSE
       w21 = invdiswt(g2dx(iup,jlo), g2dy(iup,jlo), xi, yi, pow)
       m21 = 1
    ENDIF
    IF ( fin(iup,jup) == maskin ) THEN
       w22 =  flagreal
       m22 = 0
    ELSE
       w22 = invdiswt(g2dx(iup,jup), g2dy(iup,jup), xi, yi, pow)
       m22 = 1
    ENDIF

    msum = m11 +  m12 + m21 + m22
    IF ( msum /= mnipts ) RETURN 
!    IF ( m11 +  m12 + m21 + m22 < 1 .OR. m11 +  m12 + m21 + m22 > 2 )  RETURN
    mlo = MAX(ilo -1, 1)
    mup = MIN(iup+1, nx)
    nlo = MAX(jlo -1, 1)
    nup = MIN(jup+1, ny)

    IF ( fin(mlo,nlo) == maskin ) THEN
      ds(1) = maskin
    ELSE
      ds(1) = gcirdist(g2dx(mlo,nlo), g2dy(mlo,nlo), xi, yi)
    ENDIF
   
    IF ( fin(ilo,nlo) == maskin ) THEN
      ds(2) = maskin
    ELSE
      ds(2) = gcirdist(g2dx(ilo,nlo), g2dy(ilo,nlo), xi, yi)
    ENDIF

    IF ( fin(iup,nlo) == maskin ) THEN
      ds(3) = maskin
    ELSE
      ds(3) = gcirdist(g2dx(iup,nlo), g2dy(iup,nlo), xi, yi)
    ENDIF

    IF ( fin(mup,nlo) == maskin ) THEN
      ds(4) = maskin
    ELSE
      ds(4) = gcirdist(g2dx(mup,nlo), g2dy(mup,nlo), xi, yi)
    ENDIF

    IF ( fin(mup,jlo) == maskin ) THEN
      ds(5) = maskin
    ELSE
      ds(5) = gcirdist(g2dx(mup,jlo), g2dy(mup,nlo), xi, yi)
    ENDIF

    IF ( fin(mup,jup) == maskin ) THEN
      ds(6) = maskin
    ELSE
      ds(6) = gcirdist(g2dx(mup,jup), g2dy(mup,jup), xi, yi)
    ENDIF

    IF ( fin(mup,nup) == maskin ) THEN
      ds(7) = maskin
    ELSE
      ds(7) = gcirdist(g2dx(mup,nup), g2dy(mup,nup), xi, yi)
    ENDIF

    IF ( fin(iup,nup) == maskin ) THEN
      ds(8) = maskin
    ELSE
      ds(8) = gcirdist(g2dx(iup,nup), g2dy(iup,nup), xi, yi)
    ENDIF

    IF ( fin(ilo,nup) == maskin ) THEN
      ds(9) = maskin
    ELSE
      ds(9) = gcirdist(g2dx(ilo,nup), g2dy(ilo,nup), xi, yi)
    ENDIF

    IF ( fin(mlo,nup) == maskin ) THEN
      ds(10) = maskin
    ELSE
      ds(10) = gcirdist(g2dx(mlo,nup), g2dy(mlo,nup), xi, yi)
    ENDIF

    IF ( fin(mlo,jup) == maskin ) THEN
      ds(11) = maskin
    ELSE
      ds(11) = gcirdist(g2dx(mlo,jup), g2dy(mlo,jup), xi, yi)
    ENDIF

    IF ( fin(mlo,jlo) == maskin ) THEN
      ds(12) = maskin
    ELSE
      ds(12) = gcirdist(g2dx(mlo,jlo), g2dy(mlo,jlo), xi, yi)
    ENDIF
    
    CALL indxminv(12, ds, indx, mindist)
    IF ( mindist > 0 .AND. mindist /= maskin ) THEN
      SELECT CASE(indx)
        CASE(1)
          iix = mlo
          jjy = nlo
        CASE(2)
          iix = ilo
          jjy = nlo
        CASE(3)
          iix = iup
          jjy = nlo
        CASE(4)
          iix = mup
          jjy = nlo
        CASE(5)
          iix = mup
          jjy = jlo
        CASE(6)
          iix = mup
          jjy = jup
        CASE(7)
          iix = mup
          jjy = nup
        CASE(8)
          iix = iup
          jjy = nup
        CASE(9)
          iix = ilo
          jjy = nup
        CASE(10)
          iix = mlo
          jjy = nup
        CASE(11)
          iix = mlo
          jjy = jup
        CASE(12)
          iix = mlo
          jjy = jlo
        CASE DEFAULT
          PRINT*,'INVALID **indx** weights(invdistinp)'
      END SELECT
      IF ( m11 == 1 ) THEN
        iup = iix
        jup = jjy
        IF ( fin(ilo,jup) /= maskin ) THEN
           w12 = invdiswt(g2dx(ilo,jup), g2dy(ilo,jup), xi, yi, pow)
           m12 = 1
        ENDIF
        IF ( fin(iup,jlo) /= maskin ) THEN
           w21 = invdiswt(g2dx(iup,jlo), g2dy(iup,jlo), xi, yi, pow) 
           m21 = 1
        ENDIF
        IF ( fin(iup,jup) /= maskin ) THEN
           w22 = invdiswt(g2dx(iup,jup), g2dy(iup,jup), xi, yi, pow)
           m22 = 1
        ENDIF
      ELSE IF ( m12 == 1 ) THEN
        iup = iix
        jlo = jjy
        IF ( fin(ilo,jlo) /= maskin )  THEN
           w11 = invdiswt(g2dx(ilo,jlo), g2dy(ilo,jlo), xi, yi, pow)
           m11 = 1
        ENDIF
        IF ( fin(iup,jlo) /= maskin )  THEN
           w21 = invdiswt(g2dx(iup,jlo), g2dy(iup,jlo), xi, yi, pow)
           m21 = 1
        ENDIF
        IF ( fin(iup,jup) /= maskin )  THEN
           w22 = invdiswt(g2dx(iup,jup), g2dy(iup,jup), xi, yi, pow)
           m22 = 1
        ENDIF
      ELSE IF ( m21 == 1 ) THEN
        ilo = iix
        jup = jjy
        IF ( fin(ilo,jlo) /= maskin )  THEN
           w11 = invdiswt(g2dx(ilo,jlo), g2dy(ilo,jlo), xi, yi, pow)
           m11 = 1
        ENDIF
        IF ( fin(ilo,jup) /= maskin )  THEN
           w12 = invdiswt(g2dx(ilo,jup), g2dy(ilo,jup), xi, yi, pow)
           m12 = 1
        ENDIF
        IF ( fin(iup,jup) /= maskin )  THEN
           w22 = invdiswt(g2dx(iup,jup), g2dy(iup,jup), xi, yi, pow)
           m22 = 1
        ENDIF
      ELSE IF ( m22 == 1 ) THEN
        ilo = iix
        jlo = jjy
        IF ( fin(ilo,jlo) /= maskin )  THEN
           w11 = invdiswt(g2dx(ilo,jlo), g2dy(ilo,jlo), xi, yi, pow)
           m11 = 1
        ENDIF
        IF ( fin(ilo,jup) /= maskin )  THEN
           w12 = invdiswt(g2dx(ilo,jup), g2dy(ilo,jup), xi, yi, pow)
           m12 = 1
        ENDIF
        IF ( fin(iup,jlo) /= maskin )  THEN
           w21 = invdiswt(g2dx(iup,jlo), g2dy(iup,jlo), xi, yi, pow)
           m21 = 1
        ENDIF
      ENDIF
    ENDIF
    msum = m11 + m12 + m21 + m22
    IF ( msum == mnipts ) THEN
       IF ( m11 == 1 ) THEN  
         w11 = flagreal
       ELSEIF ( m12 == 1 ) THEN
         w12 = flagreal
       ELSEIF ( m21 == 1 ) THEN
         w21 = flagreal
       ELSEIF ( m22 == 1 ) THEN
         w22 = flagreal
       ENDIF
    ENDIF
    IF ( msum >= 0 ) RETURN
    IF ( m11 == 0 ) THEN
       IF ( fin(mlo,nlo) == maskin ) THEN
          d1 = maskin
       ELSE
          d1 = gcirdist(g2dx(mlo,nlo), g2dy(mlo,nlo), xi, yi)
       ENDIF
       IF ( fin(mlo,jlo) == maskin ) THEN
          d2 = maskin
       ELSE
          d2 = gcirdist(g2dx(mlo,jlo), g2dy(mlo,jlo), xi, yi)
       ENDIF
       IF ( fin(ilo,nlo) == maskin ) THEN
          d3 = maskin
       ELSE
          d3 = gcirdist(g2dx(ilo,nlo), g2dy(ilo,nlo), xi, yi)
       ENDIF
       IF ( d1 == maskin .AND. d2 == maskin .AND. d3 == maskin ) THEN
          w11 =  flagreal
       ELSE
          nbrpts = (/d1, d2, d3/)
          CALL indxminv(3, nbrpts, indx, mindist)
          IF ( mindist > 0 .AND. mindist /= maskin ) THEN
             IF ( indx == 1 ) THEN
               ilo = mlo 
               jlo = nlo
               w11 = invdiswt(g2dx(ilo,jlo), g2dy(ilo,jlo), xi, yi, pow)
             ELSE IF ( indx == 2 ) THEN
               ilo = mlo   
               w11 = invdiswt(g2dx(ilo,jlo), g2dy(ilo,jlo), xi, yi, pow)
             ELSE IF ( indx == 3 ) THEN
               jlo = nlo
               w11 = invdiswt(g2dx(ilo,jlo), g2dy(ilo,jlo), xi, yi, pow)
             ENDIF
          ELSE
             w11 =  flagreal
          ENDIF
       ENDIF
    ENDIF
    IF ( m12 == 0 ) THEN
       IF ( fin(mlo,nup) == maskin ) THEN
          d1 = maskin
       ELSE
          d1 = gcirdist(g2dx(mlo,nup), g2dy(mlo,nup), xi, yi)
       ENDIF
       IF ( fin(mlo,jup) == maskin ) THEN
          d2 = maskin
       ELSE
          d2 = gcirdist(g2dx(mlo,jup), g2dy(mlo,jup), xi, yi)
       ENDIF
       IF ( fin(ilo,nup) == maskin ) THEN
          d3 = maskin
       ELSE
          d3 = gcirdist(g2dx(ilo,nup), g2dy(ilo,nup), xi, yi)
       ENDIF
       IF ( d1 == maskin .AND. d2 == maskin .AND. d3 == maskin ) THEN
          w12 =  flagreal
       ELSE
          nbrpts = (/d1, d2, d3/)
          CALL indxminv(3, nbrpts, indx, mindist)
          IF ( mindist > 0 .AND. mindist /= maskin ) THEN
             IF ( indx == 1 ) THEN
               ilo = mlo
               jup = nup
               w12 = invdiswt(g2dx(ilo,jup), g2dy(ilo,jup), xi, yi, pow)
             ELSE IF ( indx == 2 ) THEN
               ilo = mlo   
               w12 = invdiswt(g2dx(ilo,jup), g2dy(ilo,jup), xi, yi, pow)
             ELSE IF ( indx == 3 ) THEN
               jup = nup
               w12 = invdiswt(g2dx(ilo,jup), g2dy(ilo,jup), xi, yi, pow)
             ENDIF
          ELSE
             w12 =  flagreal
          ENDIF
       ENDIF
    ENDIF
    IF ( m21 == 0 ) THEN
       IF ( fin(mup,nlo) == maskin ) THEN
          d1 = maskin
       ELSE
          d1 = gcirdist(g2dx(mup,nlo), g2dy(mup,nlo), xi, yi)
       ENDIF
       IF ( fin(mup,jlo) == maskin ) THEN
          d2 = maskin
       ELSE
          d2 = gcirdist(g2dx(mup,jlo), g2dy(mup,jlo), xi, yi)
       ENDIF
       IF ( fin(iup,nlo) == maskin ) THEN
          d3 = maskin
       ELSE
          d3 = gcirdist(g2dx(iup,nlo), g2dy(iup,nlo), xi, yi)
       ENDIF
       IF ( d1 == maskin .AND. d2 == maskin .AND. d3 == maskin ) THEN
          w21 =  flagreal
       ELSE
          nbrpts = (/d1, d2, d3/)
          CALL indxminv(3, nbrpts, indx, mindist)
          IF ( mindist > 0 .AND. mindist /= maskin ) THEN
             IF ( indx == 1 ) THEN
               iup = mup
               jlo = nlo
               w21 = invdiswt(g2dx(iup,jlo), g2dy(iup,jlo), xi, yi, pow)
             ELSE IF ( indx == 2 ) THEN
               iup = mup
               w21 = invdiswt(g2dx(iup,jlo), g2dy(iup,jlo), xi, yi, pow)
             ELSE IF ( indx == 3 ) THEN
               jlo = nlo
               w21 = invdiswt(g2dx(iup,jlo), g2dy(iup,jlo), xi, yi, pow)
             ENDIF
          ELSE
             w21 =  flagreal
          ENDIF
       ENDIF
    ENDIF
    IF ( m22 == 0 ) THEN
       IF ( fin(mup,nup) == maskin ) THEN
          d1 = maskin
       ELSE
          d1 = gcirdist(g2dx(mup,nup), g2dy(mup,nup), xi, yi)
       ENDIF
       IF ( fin(mup,jup) == maskin ) THEN
          d2 = maskin
       ELSE
          d2 = gcirdist(g2dx(mup,jup), g2dy(mup,jup), xi, yi)
       ENDIF
       IF ( fin(iup,nup) == maskin ) THEN
          d3 = maskin
       ELSE
          d3 = gcirdist(g2dx(iup,nup), g2dy(iup,nup), xi, yi)
       ENDIF
       IF ( d1 == maskin .AND. d2 == maskin .AND. d3 == maskin ) THEN
          w22 =  flagreal
       ELSE
          nbrpts = (/d1, d2, d3/)
          CALL indxminv(3, nbrpts, indx, mindist)
          IF ( mindist > 0 .AND. mindist /= maskin ) THEN
             IF ( indx == 1 ) THEN
               iup = mup
               jup = nup
               w22 = invdiswt(g2dx(iup,jup), g2dy(iup,jup), xi, yi, pow)
             ELSE IF ( indx == 2 ) THEN
               iup = mup
               w22 = invdiswt(g2dx(iup,jup), g2dy(iup,jup), xi, yi, pow)
             ELSE IF ( indx == 3 ) THEN
               jup = nup
               w22 = invdiswt(g2dx(iup,jup), g2dy(iup,jup), xi, yi, pow)
             ENDIF
          ELSE
             w22 =  flagreal
          ENDIF
       ENDIF
    ENDIF

  END SUBROUTINE weights

  FUNCTION gcirdist(lon1, lat1, lon2, lat2) RESULT(dist)
  IMPLICIT NONE
     REAL, INTENT(IN) :: lon1, lat1, lon2, lat2
     REAL :: rlat1, rlat2
     REAL :: dlonr, dist
     IF ( lon1 == lon2 .AND. lat1 == lat2 ) THEN
        dist = 0.0
        RETURN
     ENDIF    
     dlonr = ABS(lon1 - lon2) * d2rad
     rlat1 = lat1 * d2rad
     rlat2 = lat2 * d2rad
     dist = ACOS(MAX(MIN(SIN(rlat1)*SIN(rlat2)+COS(rlat1)*COS(rlat2)*  &
            COS(dlonr),1.0), -1.0))/d2rad
  END FUNCTION 

  FUNCTION invdiswt(lon1, lat1, lon2, lat2, pow) RESULT(wt)
  IMPLICIT NONE
     REAL, INTENT(IN) :: lon1, lat1, lon2, lat2, pow
     REAL :: dist, distpow, wt
     dist = gcirdist(lon1, lat1, lon2, lat2)
     distpow = dist**pow
     IF (distpow <= 0.0) THEN 
        wt = noinpval
     ELSE
        wt = 1.0 / distpow
     ENDIF 
  END FUNCTION

  SUBROUTINE check_weights(maskin, f11, f12, f21, f22, w11, w12, w21, w22 )
  IMPLICIT NONE    
     REAL, INTENT(IN) :: maskin
     REAL, INTENT(INOUT) :: w11, w12, w21, w22
     REAL, INTENT(IN) :: f11, f12, f21, f22
       IF (f11 == maskin ) THEN
          w11 =  flagreal
       ENDIF
       IF (f12 == maskin ) THEN
          w12 =  flagreal
       ENDIF
       IF (f21 == maskin ) THEN
          w21 =  flagreal
       ENDIF
       IF (f22 == maskin ) THEN
          w22 =  flagreal
       ENDIF
  END SUBROUTINE check_weights

  SUBROUTINE indxminv(np, array, indx, mnval)
    INTEGER, INTENT(IN) :: np
    REAL, DIMENSION(np), INTENT(IN) :: array
    INTEGER, INTENT(OUT) :: indx
    REAL, INTENT(OUT) :: mnval
    INTEGER :: n
    mnval = MINVAL(array)
    indx = -999
    DO n = 1, np
      IF ( mnval == array(n) ) THEN
        indx = n
        EXIT
      END IF
    ENDDO
  END SUBROUTINE indxminv
END MODULE invdistinpmd