mod_interp.f90

Go to the documentation of this file.

!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

MODULE mod_interp
  USE mod_utils
  USE mod_prec
  USE mod_clock
  USE mod_spherical
  IMPLICIT NONE  



  !/===========================================================================/
  !
  !  LAND MASK IN BILINEAR INTERP!
  !
  !  MASK == 0  OCEAN
  !  MASK == 1  LAND
  !
  !/===========================================================================/

  ! SEARCH RADIUS
  REAL(sp) :: search=80000.0_sp

  ! SLOPE TOLERANCE FOR CURVALINEAR INTERP PROGRAM
  REAL(sp), PARAMETER :: tb_tol=1000.0_sp
  REAL(sp), PARAMETER :: lr_tol=1000.0_sp

  REAL(sp), PARAMETER :: small= 1e-6

  ! NUMBER OF POINTS PER QUADRENT
  !  INTEGER, parameter :: nppq=2
  !  INTEGER, parameter :: quad=4

  TYPE r2pts
     REAL(sp) :: weights(4,2)
     INTEGER :: quad_num(4,2)
  END TYPE r2pts

  TYPE interp_weights
     INTEGER :: nin
     INTEGER :: nout
     INTEGER, ALLOCATABLE :: index(:)
     TYPE(r2pts), ALLOCATABLE:: ptw(:)

     ! HANDLE LAND MASK TO FIND NEAREST NEIGHBOR IN MESH WITH DATA
     ! THIS IS A HACK FOR GETTING VALUES ON AN INLAND RIVER...
     INTEGER :: n_type
     INTEGER, POINTER :: n_index(:,:)
     INTEGER, POINTER :: n_cnt(:)
     INTEGER, POINTER :: n_order(:)
     INTEGER, POINTER :: n_found(:)


  END TYPE interp_weights

  INTEGER, POINTER :: n_found(:)
  
  INTEGER, PARAMETER :: type_node = 1
  INTEGER, PARAMETER :: type_element = 2



  TYPE(watch) :: min_w

  TYPE(watch) :: box_w

  TYPE(watch) :: cond_w

  TYPE(watch) :: case_w

  TYPE(watch) :: tot_w


CONTAINS

  SUBROUTINE kill_weights(MYW)
    IMPLICIT NONE
    TYPE(interp_weights),POINTER :: MYW


    IF(.not.ASSOCIATED(myw)) RETURN

    IF(ALLOCATED(myw%INDEX)) DEALLOCATE(myw%INDEX)

    IF(ALLOCATED(myw%PTW)) DEALLOCATE(myw%PTW)
    
    IF(ASSOCIATED(myw%N_INDEX)) DEALLOCATE(myw%N_INDEX)
    IF(ASSOCIATED(myw%N_CNT)) DEALLOCATE(myw%N_CNT)
    IF(ASSOCIATED(myw%N_ORDER)) DEALLOCATE(myw%N_ORDER)
    IF(ASSOCIATED(myw%N_FOUND)) DEALLOCATE(myw%N_FOUND)

    DEALLOCATE(myw)

  END SUBROUTINE kill_weights


  SUBROUTINE print_weights(MYW,STR)
    IMPLICIT NONE
    TYPE(interp_weights) :: MYW
    CHARACTER(LEN=*) :: STR
    INTEGER :: I

    WRITE(ipt,*) "===== PRINTING WEIGHTS:"//trim(str)//"========"
    WRITE(ipt,*)
    WRITE(ipt,*) "Nin=",myw%Nin
    WRITE(ipt,*) "Nout=",myw%Nout
    IF(ALLOCATED(myw%index)) THEN
       WRITE(ipt,*) "Size(index)=",size(myw%index)
    ELSE
       WRITE(ipt,*) "Index not allocated!"
    END IF

    IF(ALLOCATED(myw%PTW)) THEN
       DO i = 1,myw%Nout
          CALL print_ptw(myw%PTW(i),i)
       END DO
    ELSE
       WRITE(ipt,*) "PTW not ALLOCATED!"       
    END IF

    WRITE(ipt,*) "NType=",myw%N_TYPE

    IF(Associated(myw%N_INDEX)) THEN
       WRITE(ipt,*) "Size(N_INDEX)=",size(myw%N_INDEX)
    ELSE
       WRITE(ipt,*) "N_INDEX not allocated!"
    END IF

    IF(Associated(myw%N_CNT)) THEN
       WRITE(ipt,*) "Size(N_CNT)=",size(myw%N_CNT)
    ELSE
       WRITE(ipt,*) "N_CNT not allocated!"
    END IF

    IF(Associated(myw%N_ORDER)) THEN
       WRITE(ipt,*) "Size(N_ORDER)=",size(myw%N_ORDER)
    ELSE
       WRITE(ipt,*) "N_ORDER not allocated!"
    END IF

    IF(Associated(myw%N_FOUND)) THEN
       WRITE(ipt,*) "Size(N_FOUND)=",size(myw%N_FOUND)
    ELSE
       WRITE(ipt,*) "N_FOUND not allocated!"
    END IF
    WRITE(ipt,*) ""
    WRITE(ipt,*)"========== END WEIGHTS ========"


  END SUBROUTINE print_weights

  SUBROUTINE print_ptw(PTW,CNT)
    IMPLICIT NONE
    TYPE(r2pts) :: PTW
    INTEGER :: CNT
    WRITE(ipt,*) "===== PRINTING PTW#",cnt
    WRITE(ipt,*) "WEIGHTS="
    WRITE(ipt,'(4F12.6)')ptw%WEIGHTS(:,1)
    WRITE(ipt,'(4F12.6)')ptw%WEIGHTS(:,2)
    WRITE(ipt,*) "QUADNUM="
    WRITE(ipt,'(4I8)')ptw%QUAD_NUM(:,1)
    WRITE(ipt,'(4I8)')ptw%QUAD_NUM(:,2)


  END SUBROUTINE print_ptw

  ! THE BILINEAR METHOD USED HERE IS DESIGNED TO FUNCTION ON
  ! CURVALINEAR GRIDDED DATA. ANY 



  ! GET THE LOCATION OF THE NEAREST POINT AND RETURN AN INDEX TO IT
  ! ALLONG WITH ITS NEIGHBORS IN THE MESH FOR THE CELL THAT CONTAINS
  ! IT. THE CONDITION INDICATES HOW TO INTERPOLATE IN THIS CELL TO
  ! HANDLE A LAND MASK OR A VALUE OUTSIDE THE DOMAIN BOUNDARY
  SUBROUTINE get_loc(X_bnd,Y_bnd,msze,nsze,RSQ,MASK,X,Y,BOX,CONDITION)
    IMPLICIT NONE 
    REAL(SP), POINTER, INTENT(IN) :: X_bnd(:,:)
    REAL(SP), POINTER, INTENT(IN) :: Y_bnd(:,:)
    REAL(SP), INTENT(IN) :: X,Y
    INTEGER, INTENT(IN) :: MSZE,NSZE
    REAL(SP), POINTER :: Rsq(:,:)
    INTEGER, POINTER,INTENT(IN) :: MASK(:,:)


    INTEGER, INTENT(OUT) :: BOX(4,2)
    INTEGER, INTENT(OUT) :: CONDITION

    INTEGER :: MISSING(4)

    INTEGER, DIMENSION(2) :: NN, NNT
    INTEGER :: err,I,j,CNT
    REAL(SP) :: Xn,Yn

    REAL(SP) :: PTEST

    REAL(SP),POINTER::XIN(:,:),YIN(:,:)

    xin => x_bnd(1:msze,1:nsze)
    yin => y_bnd(1:msze,1:nsze)

    CALL watch_start_lap(min_w)
    ! Find Nearest point

!!$!---------------------------------------------------------------
!!$# if defined(SPHERICAL)
!!$!---------------------------------------------------------------
!!$
!!$
!!$    IF(ASSOCIATED(MASK))THEN;
!!$       
!!$       DO I=1,msze
!!$          DO J=1,nsze
!!$             IF(MASK(I,J)==0) CALL ARC(Xin(i,j),Yin(i,j),X,Y,RSQ(i,j))
!!$          END DO
!!$       END DO
!!$       
!!$       NN = minloc(Rsq,MASK==0)
!!$
!!$       ! AT THE NORTH POLE, ONLY THE LONGITUDE MATTERS...
!!$       IF( YIN(NN(1),NN(2)) == 90.0_SP .or. Y == 90.0_SP) THEN
!!$          
!!$          NN = minloc(abs(Xin - X),Yin == 90.0_SP)
!!$          
!!$       ELSEIF ( YIN(NN(1),NN(2)) == -90.0_SP .or. Y == -90.0_SP) THEN
!!$          NN = minloc(abs(Xin - X),Yin == -90.0_SP)
!!$          
!!$       END IF
!!$       
!!$    ELSE
!!$       DO I=1,msze
!!$          DO J=1,nsze
!!$             CALL ARC(Xin(i,j),Yin(i,j),X,Y,RSQ(i,j))
!!$          END DO
!!$       END DO
!!$       
!!$       NN = minloc(Rsq)
!!$ 
!!$       ! AT THE NORTH POLE, ONLY THE LONGITUDE MATTERS...
!!$       IF( YIN(NN(1),NN(2)) == 90.0_SP .or. Y == 90.0_SP) THEN
!!$          
!!$          NN = minloc(abs(Xin - X),Yin == 90.0_SP)
!!$          
!!$       ELSEIF ( YIN(NN(1),NN(2)) == -90.0_SP .or. Y == -90.0_SP) THEN
!!$          NN = minloc(abs(Xin - X),Yin == -90.0_SP)
!!$          
!!$       END IF
!!$ 
!!$
!!$   END IF
!!$# else
    IF(ASSOCIATED(mask))then;
       
       WHERE(mask==0)
          rsq = (xin -x)**2 + (yin -y)**2
       END WHERE
       
       nn = minloc(rsq,rsq>-1.0_sp)
       
    ELSE
       
       rsq = (xin -x)**2 + (yin -y)**2
       
       nn = minloc(rsq,rsq>-1.0_sp)
       
    END IF
!!$!---------------------------------------------------------------
!!$# endif
!!$!---------------------------------------------------------------
    
 
    CALL watch_stop_lap(min_w)

    !    PRINT*, "NN=",NN


    CALL watch_start_lap(box_w)
    box = get_box(x,y,nn,x_bnd,y_bnd,err)


    ! IN SOME STRANGE GEOMETRY, THE CLOSEST POINT IS NOT IN THE SAME CELL!
    nnt = nn
    cnt = 0
    DO WHILE(err/=0)
       ptest = rsq(nnt(1),nnt(2))
       nnt = minloc(rsq,rsq>ptest)

       IF(nnt(1) == 0) EXIT
!       IF(NNT(1) == 0) THEN
!          WRITE(IPT,*) "SEARCHED WHOLE GRID - NO VALID CELLS CONTAIN THE POINT"
!          exit
!       END IF
       box = get_box(x,y,nnt,x_bnd,y_bnd,err)
       
       IF(cnt > 12) EXIT
       cnt = cnt +1

    END DO
    
    ! SET THE FINALL VALUE FOR NN OR HANDLE ERRROR
    IF( err == 0) THEN
       nn = nnt
    ELSE
       
       IF(ASSOCIATED(mask)) THEN
          ! IF THE MASK BLANKS OUT ALL THE NODES IN THE CELL
          ! CONTAINING THE POINT, JSUT USE THE VALUE AT THE NEAREST POINT


          ! THIS IS NOT A GOOD SOLUTION - NEED TO HAVE A BACK UP
          ! USING NEAREST NEIGHBOR WITH DATA...
          condition = -7
          box(1,:) = nn
          return
       ELSE
          ! GET THE LOCATION
          xn = x_bnd(nn(1),nn(2))
          yn = y_bnd(nn(1),nn(2))

          ! THIS SHOULD NOT HAPPEN IF THRE IS NO MASK
          write(ipt,*) "Error searching for Bilinear Interpolation to Point:"
          write(ipt,*) "X=",x,"; Y=",y
          write(ipt,*) "Found Nearest Grid index:",nn
          write(ipt,*) "Bounds(",msze,nsze,")"
          write(ipt,*) "Grid location:",xn,yn
          
!          WRITE(IPT,*) RSQ(:,1)
!          WRITE(IPT,*) RSQ(:,2)

          CALL fatal_error("Could not find the triangle of the four nearest nodes",&
               & "containing the point x,y ?")

       END IF
    END IF


    CALL watch_stop_lap(box_w)


    ! ASSESS CONDITION - IS THE POINT OUTSIDE THE DOMAIN?

    ! ASSUME CONDITION 0
    condition = 0

    ! ANY VALUES IN THE BOUNDARY SET TO ZERO
    box(:,1) = mod(box(:,1),msze+1)

    box(:,2) = mod(box(:,2),nsze+1)


    IF(ASSOCIATED(mask)) THEN
       missing = 0
       DO i = 1,4

          IF(box(i,1)*box(i,2)==0) THEN
             box(i,:)=0
             missing(i) = 1
             cycle
          END IF


!         WRITE(IPT,*) "I",I,"; MASK(BOX(I)):",MASK(BOX(I,1),BOX(I,2))

          IF(mask(box(i,1),box(i,2))==1) THEN
             box(i,:)=0
             missing(i) = 1
          END IF
!          WRITE(IPT,*) "I",I,"; missing;",missing

       END DO

       select case(sum(missing))
       CASE(1)
          IF(all(missing ==(/0,0,0,1/)) ) THEN
             condition=-1
          ELSEIF(all(missing ==(/0,0,1,0/)) ) THEN
             condition=-2
          ELSEIF(all(missing ==(/0,1,0,0/)) ) THEN
             condition=-3
          ELSEIF(all(missing ==(/1,0,0,0/)) ) THEN
             condition=-4
          END IF
          RETURN
       CASE(2)

          IF(all(missing == (/1,0,1,0/)) ) THEN
             condition=-5
          ELSEIF(all(missing == (/0,1,0,1/)) ) THEN
             condition=-6
          ELSEIF(all(missing ==(/1,1,0,0/)) ) THEN
             condition=6
          ELSEIF(all(missing ==(/0,1,1,0/)) ) THEN
             condition=4
          ELSEIF(all(missing ==(/0,0,1,1/)) ) THEN
             condition=2
          ELSEIF(all(missing ==(/1,0,0,1/)) ) THEN
             condition=8
          END IF
          RETURN
       CASE(3)

          IF(all(missing ==(/1,1,1,0/)) ) THEN
             condition=5
          ELSEIF(all(missing ==(/1,1,0,1/)) ) THEN
             condition=7
          ELSEIF(all(missing ==(/1,0,1,1/)) ) THEN
             condition=1
          ELSEIF(all(missing ==(/0,1,1,1/)) ) THEN
             condition=3
          END IF
          RETURN
       CASE(4)
          CALL fatal_error&
               ("THE LAND MASK HAS LEFT NO VALID POINTS - THIS SHOULD BE IMPOSSIBLE!")
       END select

    END IF

    CALL watch_start_lap(cond_w)

    ! IF ANY LESS THAN ONE, VALUES ARE OUTSIDE!
    IF(any(box==0)) THEN
       ! OUTSIDE BOUNDARY OR ON THE LAND MASK!!!


       ! TEST FOR 8,6,4,2 while one value will correctly determine value
       IF(count(box==0)==2) THEN

          ! TEST FOR ZERO: TOP OR BOTTOM, LEFT OR RIGHT
          IF(box(1,1) == 0) condition = 8
          IF(box(1,2) == 0) condition = 6
          IF(box(3,1) == 0) condition = 4
          IF(box(3,2) == 0) condition = 2
       END IF


       ! SET BOTH INDICIES TO ZERO IF EITHER IS
       DO i = 1,4
          IF(product(box(i,:))==0) box(i,:)=0
       END DO


       IF(count(box==0)==6)THEN
          ! TEST THE NON ZERO CORNERS
          IF(box(1,1) /= 0) condition = 3
          IF(box(2,1) /= 0) condition = 1
          IF(box(3,1) /= 0) condition = 7
          IF(box(4,1) /= 0) condition = 5


       END IF

    END IF

    CALL watch_stop_lap(cond_w)

  END SUBROUTINE get_loc

  ! GET THE CELL THAT CONTAINS THE POINT, IT COULD BE ONE OF FOUR
  ! THAT SHARE A COMMON NODE....
  FUNCTION get_box(X,Y,NN,X_bnd,Y_bnd,err) RESULT(BOX)
    IMPLICIT NONE
    REAL(sp), POINTER, INTENT(IN) :: x_bnd(:,:)
    REAL(sp), POINTER, INTENT(IN) :: y_bnd(:,:)
    REAL(sp), INTENT(IN) :: x,y
    INTEGER, INTENT(IN):: nn(2)
    INTEGER, INTENT(OUT) :: err

    INTEGER:: box(4,2)

    REAL(sp) :: xn,yn,xt(3),yt(3)
    INTEGER :: i,j,tri


    INTEGER:: tn(4,2)
    INTEGER:: tno(4,2) ! TRIANGLE OPOSITE
    INTEGER:: bn(4,2)
    INTEGER:: dn(4,2)


    ! INDEX OF POINTS TO TRINAGLE
    tn(1,:) = (/ 1 , 0 /)
    tn(2,:) = (/ 0 , 1 /)
    tn(3,:) = (/-1 , 0 /)
    tn(4,:) = (/ 0 ,-1 /)

    ! INDEX OF POINTS TO OPOSITE CORNER
    tno(1,:) = (/ 1 , 1 /)
    tno(2,:) = (/-1 , 1 /)
    tno(3,:) = (/-1 ,-1 /)
    tno(4,:) = (/ 1 ,-1 /)

    ! INDEX OF POINTS TO BOX
    bn(1,:) = (/ 0 , 1 /)
    bn(2,:) = (/ 1 , 1 /)
    bn(3,:) = (/ 1 , 0 /)
    bn(4,:) = (/ 0 , 0 /)

    ! DELTA BOX DEPENDING ON WHICH TRIANGLE
    dn(1,:) = (/ 0 , 0 /)
    dn(2,:) = (/-1 , 0 /)
    dn(3,:) = (/-1 ,-1 /)
    dn(4,:) = (/ 0 ,-1 /)


    xt(1) = x_bnd(nn(1),nn(2))

    yt(1) = y_bnd(nn(1),nn(2))

    err = -1
    box=-1


    ! TEST EACH OF THE FOUR TRIANGLE IN CELLS SURROUNDING THE NEAREST NODE
    DO i=1,4

       j = mod(i,4)+1

       xt(2) = x_bnd(nn(1)+tn(i,1),nn(2)+tn(i,2))
       yt(2) = y_bnd(nn(1)+tn(i,1),nn(2)+tn(i,2))

       xt(3) = x_bnd(nn(1)+tn(j,1),nn(2)+tn(j,2))
       yt(3) = y_bnd(nn(1)+tn(j,1),nn(2)+tn(j,2))

!       WRITE(IPT,*) "======FOUND TRI==========="
!       write(ipt,*) "XT",XT
!       write(ipt,*) "YT",YT
!       write(ipt,*) "X,Y",X,Y


       IF(isintri(x,y,xt,yt)) THEN
          !          print*, "TRI#",I


          err = 0
          bn(:,1) = bn(:,1) +  dn(i,1)
          bn(:,2) = bn(:,2) +  dn(i,2)

          box(:,1) = bn(:,1) + nn(1) 
          box(:,2) = bn(:,2) + nn(2) 


          RETURN
       END IF

    END DO

    ! TEST EACH OF THE FOUR TRIANGLES IN THE OPOSITE CORNER OF CELL
    ! WITH THE NEAREST NODE
    DO i=1,4

       j = mod(i,4)+1

       !OPOSITE CORNER
       xt(1) = x_bnd(nn(1)+tno(i,1),nn(2)+tno(i,2))
       yt(1) = y_bnd(nn(1)+tno(i,1),nn(2)+tno(i,2))


       ! TWO OTHERS
       xt(2) = x_bnd(nn(1)+tn(i,1),nn(2)+tn(i,2))
       yt(2) = y_bnd(nn(1)+tn(i,1),nn(2)+tn(i,2))

       xt(3) = x_bnd(nn(1)+tn(j,1),nn(2)+tn(j,2))
       yt(3) = y_bnd(nn(1)+tn(j,1),nn(2)+tn(j,2))

!          WRITE(IPT,*) "======FOUND TRI==========="
!          write(ipt,*) "XT",XT
!          write(ipt,*) "YT",YT
!          write(ipt,*) "X,Y",X,Y


       IF(isintri(x,y,xt,yt)) THEN
          !          print*, "TRI OPOSITE! #",I


          err = 0
          bn(:,1) = bn(:,1) +  dn(i,1)
          bn(:,2) = bn(:,2) +  dn(i,2)

          box(:,1) = bn(:,1) + nn(1) 
          box(:,2) = bn(:,2) + nn(2) 


          RETURN
       END IF

    END DO




  END FUNCTION get_box

  SUBROUTINE setup_interp_bilinear_a(Xin,Yin,Xout,Yout,WEIGHTS, land_mask)
    IMPLICIT NONE 
    REAL(SP), ALLOCATABLE, TARGET, INTENT(IN) :: Xin(:,:)
    REAL(SP), ALLOCATABLE, TARGET, INTENT(IN) :: Yin(:,:)
    REAL(SP), ALLOCATABLE, TARGET, INTENT(IN) :: Xout(:)
    REAL(SP), ALLOCATABLE, TARGET, INTENT(IN) :: Yout(:)

    TYPE(interp_weights), INTENT(OUT) :: WEIGHTS

    INTEGER, OPTIONAL, ALLOCATABLE, TARGET, INTENT(IN) :: LAND_MASK(:,:)

    REAL(SP), POINTER :: XinP(:,:)
    REAL(SP), POINTER :: YinP(:,:)
    REAL(SP), POINTER :: XoutP(:)
    REAL(SP), POINTER :: YoutP(:)

    INTEGER, POINTER :: LAND_MASKP(:,:)

    NULLIFY(xinp,yinp,xoutp,youtp,land_maskp)


    IF(ALLOCATED(xin)) xinp => xin

    IF(ALLOCATED(yin)) yinp => yin

    IF(ALLOCATED(xout)) xoutp => xout

    IF(ALLOCATED(yout)) youtp => yout

    IF(PRESENT(land_mask)) THEN

       IF(ALLOCATED(land_mask)) THEN
          land_maskp => land_mask
       ELSE
          CALL fatal_error("PASSED LAND_MASK BUT IT IS NOT ALLOCATED") 
       END IF

       CALL setup_interp_bilinear_p(xinp,yinp,xoutp,youtp,weights,land_maskp)
    ELSE
       CALL setup_interp_bilinear_p(xinp,yinp,xoutp,youtp,weights)
    END IF
  END SUBROUTINE setup_interp_bilinear_a


  SUBROUTINE setup_interp_bilinear_p(Xin,Yin,Xout,Yout,WEIGHTS,land_mask)
    USE all_vars, only :nbsn,m,n,mt,nt,par

    IMPLICIT NONE 
    REAL(SP), POINTER, INTENT(IN) :: Xin(:,:)
    REAL(SP), POINTER, INTENT(IN) :: Yin(:,:)
    REAL(SP), POINTER, INTENT(IN) :: Xout(:)
    REAL(SP), POINTER, INTENT(IN) :: Yout(:)

    TYPE(interp_weights), INTENT(OUT) :: WEIGHTS
    INTEGER, OPTIONAL, POINTER, INTENT(IN) :: LAND_MASK(:,:)

    INTEGER, POINTER :: MASK(:,:)

    INTEGER, POINTER :: N_INDEX(:,:)
    INTEGER, POINTER :: N_FOUND(:)
    INTEGER, POINTER :: N_CNT(:)
    INTEGER, POINTER :: N_ORDER(:)
    INTEGER :: mk_cnt


    INTEGER :: N_TYPE

    INTEGER :: I, J, status, lb, ub, msze,nsze

    REAL(SP) :: Xbox(4)
    REAL(SP) :: Ybox(4)
    REAL(SP) :: wghts(4)


    REAL(SP) :: DRA,DRB,DRT

    ! DELTA X/Y for each edge
    REAL(SP) :: DXT, DXB
    REAL(SP) :: DYT, DYB

    REAL(SP) :: ARCX1, ARCX2,ARCX3

    REAL(SP),POINTER::X_bnd(:,:),Y_bnd(:,:),RSQ(:,:)

    LOGICAL :: BOUNDS_FLAG = .false.

    REAL(SP) :: X,Y
    INTEGER :: BOX(4,2)
    INTEGER :: CONDITION

    ! CHECK INPUT ALLOCATION STATUS
    if(.not.associated(xin)) CALL fatal_error("SETUP_INTERP: INPUT ARGU&
         &MENTS MUST BE ALLOCATED!")
    if(.not.associated(yin)) CALL fatal_error("SETUP_INTERP: INPUT ARGU&
         &MENTS MUST BE ALLOCATED!")
    if(.not.associated(xout)) CALL fatal_error("SETUP_INTERP: INPUT ARGU&
         &MENTS MUST BE ALLOCATED!")
    if(.not.associated(yout)) CALL fatal_error("SETUP_INTERP: INPUT ARGU&
         &MENTS MUST BE ALLOCATED!")

    bounds_flag = .false.


    nullify(n_index)
    nullify(n_cnt)
    nullify(n_order)
    nullify(n_found)

    ! GET THE DIMENSIONS
    weights%Nout= ubound(xout,1)
    msze=size(xin,1)
    nsze=size(xin,2)

    weights%Nin  = msze*nsze


    ! ALLOCATE THE SPACE FOR THE WEIGHTS OUTPUT AND INDEX
    ALLOCATE(weights%PTW(weights%Nout), stat=status)
    IF(status /= 0) CALL fatal_error("SETUP_INTERP: COULD NOT ALLOCATE SPACE")


    ALLOCATE(rsq(msze,nsze)); rsq = 0.0_sp
    ALLOCATE(x_bnd(0:msze+1,0:nsze+1)); x_bnd=0.0_sp
    ALLOCATE(y_bnd(0:msze+1,0:nsze+1)); y_bnd=0.0_sp


    ! WE NEED A LOCATION MATRIX WHICH CAN HANDLE VALUES OUTSIDE THE
    ! DOMAIN SO USE THE LIMITS AT EACH EDGE...

!---------------------------------------------------------------
    ! Set center values
    x_bnd(1:msze,1:nsze)=xin

    ! Set edge values
    x_bnd(0,:) = 2*x_bnd(1,:) - x_bnd(msze,:) ! LEFT
    x_bnd(msze+1,:) = 2*x_bnd(msze,:) - x_bnd(1,:) ! RIGHT

    x_bnd(:,0) = x_bnd(:,1) !TOP
    x_bnd(:,nsze+1) = x_bnd(:,nsze) !BOTTOM

    ! Set center values
    y_bnd(1:msze,1:nsze)=yin

    ! Set edge values
    y_bnd(0,:) = y_bnd(1,:) ! LEFT
    y_bnd(msze+1,:) = y_bnd(msze,:) ! RIGHT

    y_bnd(:,0) = 2*y_bnd(:,1) - y_bnd(:,nsze) !TOP
    y_bnd(:,nsze+1) = 2*y_bnd(:,nsze) - y_bnd(:,1) !BOTTOM



    ! CORNERS
    y_bnd(0,0) = y_bnd(1,0)
    y_bnd(msze+1,0) = y_bnd(msze,0)

    y_bnd(0,nsze+1) = y_bnd(1,nsze+1)
    y_bnd(msze+1,nsze+1) = y_bnd(msze,nsze+1)

    x_bnd(0,0) = x_bnd(0,1)
    x_bnd(msze+1,0) = x_bnd(msze+1,1)

    x_bnd(0,nsze+1) = x_bnd(0,nsze)
    x_bnd(msze+1,nsze+1) = x_bnd(msze+1,nsze)

!---------------------------------------------------------------
!---------------------------------------------------------------


! TO PRINT THE X_BND AND Y_BND ARRAY
! Uncommented by Jadon Ge
!    DO I=0,UBOUND(X_bnd,2)
!       write(ipt,*) "X",X_bnd(:,I)
!    END DO

!    DO I=0,UBOUND(X_bnd,2)
!       write(ipt,*) "Y",Y_bnd(:,I)
!    END DO

    !===================================================
    ! TIMING STUFF TO TEST OPTIMIZATION
    !===================================================

    CALL watch_init(min_w)
    CALL watch_init(box_w)
    CALL watch_init(cond_w)
    CALL watch_init(case_w)

    CALL watch_init(tot_w)

    !===================================================

    n_type = 0
    NULLIFY(mask)
    IF(PRESENT(land_mask))THEN
       mask => land_mask

       IF(all(mask==1)) CALL fatal_error&
            &("BILINEAR INTERP CAN NOT WORK WITH EVERY POINT MASKED")

       IF(any(mask<0 .OR. mask>1)) CALL fatal_error&
            &("BILINEAR INTERP MASK VALUES ARE INCORRECT:",&
            & "SET 1 FOR MASKED (LAND)",&
            & "SET 0 FOR UNMASKED (VALID/OCEAN)")

       IF(all(mask==0)) NULLIFY(mask)


       IF(weights%Nout == mt) THEN

          IF (.not. allocated(nbsn)) CALL fatal_error&
               &("CAN NOD USE MASK FOR A GRID VARIABLE WITHOUT RUNNING TGE FIRST")

          n_type = type_node
          ! allocate :: M, max number of nodes surrounding a node
          ALLOCATE(n_index(ubound(nbsn,1),ubound(nbsn,2))); n_index=0
          ALLOCATE(n_found(0:mt)); n_found=0
          ALLOCATE(n_cnt(m)); n_cnt=0

          IF(par) call warning("Possible interpolation using land mask",&
               &"This procedure works in parallel but may not report"&
               &,"Domain boundary errors which stop the code",&
               &"If your code stops with no message try running single processor mode")


       ELSEIF(weights%Nout == nt) THEN

          IF (.not. allocated(nbsn)) CALL fatal_error&
               &("CAN NOD USE MASK FOR A GRID VARIABLE WITHOUT RUNNING TGE FIRST")

          n_type = type_element
          ALLOCATE(n_index(n,3)); n_index=0
          ALLOCATE(n_found(0:nt)); n_found=0
          ALLOCATE(n_cnt(n)); n_cnt=0

          IF(par) call warning("Possible interpolation using land mask",&
               &"This procedure works in parallel but may not report"&
               &,"Domain boundary errors which stop the code",&
               &"If your code stops with no message try running single processor mode")



       END IF


    END IF




    do i=1,weights%Nout

       x = xout(i) ! THE LOCATION WE ARE SOLVING FOR!
       y = yout(i) 

       rsq = -1.0_sp
       CALL get_loc(x_bnd,y_bnd,msze,nsze,rsq,mask,x,y,box,condition)


       ! GET_LOC RETURNS THE INDICIES TO THE GRID CELL CONTAINING X,Y
       ! IN A CLOCKWISE ORDER STARTING FROM UPPER LEFT


       !               |                  |
       !      VII      |       VI         |           V
       !               |                  |
       !---------------------------------------------------------
       !               |                  |
       !               |                  |
       !      VIII     |       box        |           IV        
       !               |                  |
       !---------------------------------------------------------
       !               |                  |
       !      I        |       II         |           III


       ! NEGATIVE CONDITION VALUES OCCUR WHEN A LAND MASK IS APPLIED
       !
       ! o - ocean
       ! x - land
       ! 
       ! BILINEAR INTERP IN TRIANGLE, WEIGHTED AVERAGE IN MISSING QUADRANT
       !           o o
       !   I   -   x o 
       ! 
       !           o o
       !  II   -   o x
       ! 
       !           o x
       ! III   -   o o
       ! 
       !           x o
       !  IV   -   o o
       !
       ! WEIGHTED AVERAGE BETEEN POINTS
       !           o x 
       !   V   -   x o
       ! 
       !           x o 
       !  VI   -   o x
       ! OTHER CASES OF TWO MISSING CORRISPOND TO CASES 2/4/6/8
       !
       ! ALL CASES OF THREE MISSING CORRISPOND TO CASES 1/3/5/7
       !
       ! VII  - IF THE MASK BLANKS OUT ALL NODES IN THE CELL
       ! CONTAINING THE POINT JUST USE THE NEAREST VALUE



!      WRITE(IPT,*) "CONDITION=",CONDITION

       weights%PTW(i)%QUAD_NUM= box
       weights%PTW(i)%WEIGHTS = 0.0_sp

       CALL watch_start_lap(case_w)

       IF(condition/=0) bounds_flag = .true.

       !   Normal situation
       SELECT CASE(condition)

       CASE(0) ! INSIDE BOUNDARIES

!---------------------------------------------------------------
          ! DELTA X
          dxt = (xin(box(2,1),box(2,2)) -  xin(box(1,1),box(1,2)))

          dxb = (xin(box(4,1),box(4,2)) -  xin(box(3,1),box(3,2)))

          ! DELTA Y
          dyt = (yin(box(2,1),box(2,2)) -  yin(box(1,1),box(1,2)))

          dyb = (yin(box(4,1),box(4,2)) -  yin(box(3,1),box(3,2)))

!---------------------------------------------------------------
!---------------------------------------------------------------


          ! CHECK FOR REASONABLE SLOPES!
          IF(abs(dyt) < tb_tol/2._sp *abs(dxt) .and. abs(dyb) < tb_tol/2._sp *abs(dxb)) THEN
             
             
             xbox(1) = xin(box(1,1),box(1,2))
             xbox(2) = xin(box(2,1),box(2,2))
             xbox(3) = xin(box(3,1),box(3,2))
             xbox(4) = xin(box(4,1),box(4,2))
             
             ybox(1) = yin(box(1,1),box(1,2))
             ybox(2) = yin(box(2,1),box(2,2))
             ybox(3) = yin(box(3,1),box(3,2))
             ybox(4) = yin(box(4,1),box(4,2))
             
             
             
             CALL interp_0(xbox,ybox,x,y,wghts)
             
             weights%PTW(i)%WEIGHTS(1,1) = wghts(1)
             weights%PTW(i)%WEIGHTS(2,1) = wghts(2)
             weights%PTW(i)%WEIGHTS(3,1) = wghts(3)
             weights%PTW(i)%WEIGHTS(4,1) = wghts(4)
             
          ELSE ! TRY ROTATION!

              ! CHECK FOR REASONABLE SLOPES!
             IF(abs(dyt) > tb_tol *abs(dxt)) THEN
                CALL fatal_error&
                     &("THIS CURVALINEAR MESH HAS CELLS WHICH ARE TOO EXTREME FOR THIS METHOD",&
                     & "PLEASE EXAMINE TOP_DIR/testing/interp",&
                     & "IF YOU WOULD LIKE TO ADD SUPPORT FOR EXTREME CURVALINEAR INTERPOLATION")
             END IF


             ! ROTATION METHOD DOES NOT APPEAR TO WORK PROPERLY -
             ! LEAVE AS BASIS FOR FUTURE WORK
             
!!$             write(ipt,*) "BOX-pre",BOX
!!$             BOX = cshift(BOX,1)
!!$             write(ipt,*) "BOX-post",BOX
!!$
!!$
!!$             ! DELTA X
!!$             DXT = (XIN(BOX(2,1),BOX(2,2)) -  XIN(BOX(1,1),BOX(1,2)))
!!$             
!!$             DXB = (XIN(BOX(4,1),BOX(4,2)) -  XIN(BOX(3,1),BOX(3,2)))
!!$             
!!$             ! DELTA Y
!!$             DYT = (YIN(BOX(2,1),BOX(2,2)) -  YIN(BOX(1,1),BOX(1,2)))
!!$             
!!$             DYB = (YIN(BOX(4,1),BOX(4,2)) -  YIN(BOX(3,1),BOX(3,2)))
!!$             
!!$             
!!$             
!!$             
!!$             ! CHECK FOR REASONABLE SLOPES!
!!$             IF(ABS(DYT) > TB_TOL *ABS(DXT)) THEN
!!$                CALL FATAL_ERROR("CAN NOT FIND A GOOD ROTATION TO INTE&
!!$                     &RPOLATE THIS CELL!")
!!$             END IF
!!$             
!!$             IF(ABS(DYB) > TB_TOL *ABS(DXB)) THEN
!!$                CALL FATAL_ERROR("CAN NOT FIND A GOOD ROTATION TO INTE&
!!$                     &RPOLATE THIS CELL!")
!!$             END IF
!!$             
!!$
!!$             XBOX(1) = XIN(BOX(1,1),BOX(1,2))
!!$             XBOX(2) = XIN(BOX(2,1),BOX(2,2))
!!$             XBOX(3) = XIN(BOX(3,1),BOX(3,2))
!!$             XBOX(4) = XIN(BOX(4,1),BOX(4,2))
!!$             
!!$             YBOX(1) = YIN(BOX(1,1),BOX(1,2))
!!$             YBOX(2) = YIN(BOX(2,1),BOX(2,2))
!!$             YBOX(3) = YIN(BOX(3,1),BOX(3,2))
!!$             YBOX(4) = YIN(BOX(4,1),BOX(4,2))
!!$             
!!$             
!!$             
!!$             CALL INTERP_0(Xbox,Ybox,X,Y,wghts)
!!$             
!!$             wghts = cshift(wghts,-1)
!!$             
!!$
!!$             WEIGHTS%PTW(i)%WEIGHTS(1,1) = wghts(1)
!!$             WEIGHTS%PTW(i)%WEIGHTS(2,1) = wghts(2)
!!$             WEIGHTS%PTW(i)%WEIGHTS(3,1) = wghts(3)
!!$             WEIGHTS%PTW(i)%WEIGHTS(4,1) = wghts(4)
!!$             
          END IF

          


       CASE(1)

          weights%PTW(i)%WEIGHTS(2,1) = 1.0_sp

       CASE(2)

!---------------------------------------------------------------
          weights%PTW(i)%WEIGHTS(1,1)= (xin(box(2,1),box(2,2)) - x)/(xin(box(2,1),box(2,2)) - xin(box(1,1),box(1,2)))

          weights%PTW(i)%WEIGHTS(2,1)= (x - xin(box(1,1),box(1,2)))/(xin(box(2,1),box(2,2)) - xin(box(1,1),box(1,2)))
!--------------------------------------------------------------------
!--------------------------------------------------------------------

       CASE(3)

          weights%PTW(i)%WEIGHTS(1,1) = 1.0_sp

       CASE(4)

          ! RATIO OF THE DIFFERENCE IN LATITUDE IS THE SAME IN SPHERICAL
          weights%PTW(i)%WEIGHTS(1,1)= (yin(box(4,1),box(4,2)) - y)/(yin(box(4,1),box(4,2)) - yin(box(1,1),box(1,2)))

          weights%PTW(i)%WEIGHTS(4,1)= (y - yin(box(1,1),box(1,2)))/(yin(box(4,1),box(4,2)) - yin(box(1,1),box(1,2)))


       CASE(5)

          weights%PTW(i)%WEIGHTS(4,1) = 1.0_sp

       CASE(6)
!---------------------------------------------------------------
          weights%PTW(i)%WEIGHTS(3,1)= (xin(box(4,1),box(4,2)) - x)/(xin(box(4,1),box(4,2)) - xin(box(3,1),box(3,2)))

          weights%PTW(i)%WEIGHTS(4,1)= (x - xin(box(3,1),box(3,2)))/(xin(box(4,1),box(4,2)) - xin(box(3,1),box(3,2)))
!--------------------------------------------------------------------
!--------------------------------------------------------------------
       CASE(7)

          weights%PTW(i)%WEIGHTS(3,1) = 1.0_sp

       CASE(8)

          ! RATIO OF THE DIFFERENCE IN LATITUDE IS THE SAME IN SPHERICAL
          weights%PTW(i)%WEIGHTS(2,1)= (y - yin(box(3,1),box(3,2)))/(yin(box(2,1),box(2,2)) - yin(box(3,1),box(3,2)))

          weights%PTW(i)%WEIGHTS(3,1)= (yin(box(2,1),box(2,2)) - y)/(yin(box(2,1),box(2,2)) - yin(box(3,1),box(3,2)))

       CASE(-1)

          ! BOTTOM LEFT IS MISSING
          ! FOR THE TWO EDGES WITH DATA
 
          xbox(1) = x_bnd(box(1,1),box(1,2))
          xbox(2) = x_bnd(box(2,1),box(2,2))
          xbox(3) = x_bnd(box(3,1),box(3,2))
          xbox(4) = x_bnd(box(4,1),box(4,2))
          
          ybox(1) = y_bnd(box(1,1),box(1,2))
          ybox(2) = y_bnd(box(2,1),box(2,2))
          ybox(3) = y_bnd(box(3,1),box(3,2))
          ybox(4) = y_bnd(box(4,1),box(4,2))
          
          CALL interp_neg(xbox,ybox,x,y,wghts)
          
          weights%PTW(i)%WEIGHTS(1,1) = wghts(1)
          weights%PTW(i)%WEIGHTS(2,1) = wghts(2)
          weights%PTW(i)%WEIGHTS(3,1) = wghts(3)
          weights%PTW(i)%WEIGHTS(4,1) = wghts(4)

       CASE(-2)

          ! BOTTOM RIGHT IS MISSING
          box = cshift(box,-1)

          xbox(1) = x_bnd(box(1,1),box(1,2))
          xbox(2) = x_bnd(box(2,1),box(2,2))
          xbox(3) = x_bnd(box(3,1),box(3,2))
          xbox(4) = x_bnd(box(4,1),box(4,2))
          
          ybox(1) = y_bnd(box(1,1),box(1,2))
          ybox(2) = y_bnd(box(2,1),box(2,2))
          ybox(3) = y_bnd(box(3,1),box(3,2))
          ybox(4) = y_bnd(box(4,1),box(4,2))
          
          CALL interp_neg(xbox,ybox,x,y,wghts)
          
          wghts = cshift(wghts,1)

          weights%PTW(i)%WEIGHTS(1,1) = wghts(1)
          weights%PTW(i)%WEIGHTS(2,1) = wghts(2)
          weights%PTW(i)%WEIGHTS(3,1) = wghts(3)
          weights%PTW(i)%WEIGHTS(4,1) = wghts(4)
          
 
      CASE(-3)

          ! TOP RIGHT IS MISSING

          box = cshift(box,-2)

          xbox(1) = x_bnd(box(1,1),box(1,2))
          xbox(2) = x_bnd(box(2,1),box(2,2))
          xbox(3) = x_bnd(box(3,1),box(3,2))
          xbox(4) = x_bnd(box(4,1),box(4,2))
          
          ybox(1) = y_bnd(box(1,1),box(1,2))
          ybox(2) = y_bnd(box(2,1),box(2,2))
          ybox(3) = y_bnd(box(3,1),box(3,2))
          ybox(4) = y_bnd(box(4,1),box(4,2))
          
          CALL interp_neg(xbox,ybox,x,y,wghts)
          
          wghts = cshift(wghts,2)

          weights%PTW(i)%WEIGHTS(1,1) = wghts(1)
          weights%PTW(i)%WEIGHTS(2,1) = wghts(2)
          weights%PTW(i)%WEIGHTS(3,1) = wghts(3)
          weights%PTW(i)%WEIGHTS(4,1) = wghts(4)
          

      CASE(-4)

         ! TOP LEFT IS MISSING
 
         box = cshift(box,-3)

 
          xbox(1) = x_bnd(box(1,1),box(1,2))
          xbox(2) = x_bnd(box(2,1),box(2,2))
          xbox(3) = x_bnd(box(3,1),box(3,2))
          xbox(4) = x_bnd(box(4,1),box(4,2))
          
          ybox(1) = y_bnd(box(1,1),box(1,2))
          ybox(2) = y_bnd(box(2,1),box(2,2))
          ybox(3) = y_bnd(box(3,1),box(3,2))
          ybox(4) = y_bnd(box(4,1),box(4,2))
          
          CALL interp_neg(xbox,ybox,x,y,wghts)
          
          wghts = cshift(wghts,3)

          weights%PTW(i)%WEIGHTS(1,1) = wghts(1)
          weights%PTW(i)%WEIGHTS(2,1) = wghts(2)
          weights%PTW(i)%WEIGHTS(3,1) = wghts(3)
          weights%PTW(i)%WEIGHTS(4,1) = wghts(4)
          

       CASE(-5) 

!---------------------------------------------------------------
          dra = sqrt( (x_bnd(box(2,1),box(2,2)) - x)**2 + (y_bnd(box(2,1),box(2,2)) - y)**2)
          drb = sqrt( (x_bnd(box(4,1),box(4,2)) - x)**2 + (y_bnd(box(4,1),box(4,2)) - y)**2)
!---------------------------------------------------------------
!---------------------------------------------------------------

          drt = (dra +drb)
          ! NORMALIZE TO ONE
          weights%PTW(i)%WEIGHTS(2,1) = drb/drt

          weights%PTW(i)%WEIGHTS(4,1) = dra/drt

       CASE(-6) 

!---------------------------------------------------------------
          dra = sqrt( (x_bnd(box(1,1),box(1,2)) - x)**2 + (y_bnd(box(1,1),box(1,2)) - y)**2)
          drb = sqrt( (x_bnd(box(3,1),box(3,2)) - x)**2 + (y_bnd(box(3,1),box(3,2)) - y)**2)

!---------------------------------------------------------------
!---------------------------------------------------------------

          drt = (dra +drb)
          ! NORMALIZE TO ONE
          weights%PTW(i)%WEIGHTS(1,1) = drb/drt

          weights%PTW(i)%WEIGHTS(3,1) = dra/drt


       CASE(-7) 

          ! NOTE: THIS METHOD IS DESIGNED TO SEARCH FOR THE NEAREST
          ! MESH NEIGHBOR THAT HAS AN UNMASKED VALUES AND TAKE THE
          ! AVERAGE OF THESE VALUES, PROPAGATING SEQUENTIALLY INTO
          ! THE MASKED REGION. THIS CAN BE SUCCESSFUL IN SMALL
          ! REGIONS OF MASKED DATA LIKE A RIVER. IT CAN NOT MAKE UP
          ! PHYSICS OR DATA AND SHOULD GENERALLY NOT BE TRUSTED!

          IF(n_type==type_node) THEN
             n_found(i) = -1
             weights%PTW(i)%QUAD_NUM = 0
             IF(i>m) CALL fatal_error&
                  &("INTERPOLATION WITH MASK: ERORR!",&
                  & "NEARST MESH NEIGHBOR SEARCH CAN NOT WORK WHEN THE MASK COVERS A DOMAIN BOUNDARY!")
             
          ELSEIF( n_type==type_element) THEN
             n_found(i) = -1
             weights%PTW(i)%QUAD_NUM = 0
             IF(i>n) CALL fatal_error&
                  &("INTERPOLATION WITH MASK: ERORR!",&
                  & "NEARST MESH NEIGHBOR SEARCH CAN NOT WORK WHEN THE MASK COVERS A DOMAIN BOUNDARY!")

          ELSE
             ! IF NOT USING FVCOM GRID USE NEAREST VALUE BY DISTANCE
             weights%PTW(i)%WEIGHTS(1,1) = 1.0_sp
             weights%PTW(i)%QUAD_NUM(2:4,:) = 0
             
          END IF



       CASE DEFAULT

          CALL fatal_error("SETUP_INTERP_BILINEAR: INVALID CONDITION" )
       END SELECT

       CALL watch_stop_lap(case_w)

    END do

    
    ! TAKE CARE OF ANY MARKED AS MASKED
    IF(n_type==type_node) THEN
       IF(any(n_found < 0)) THEN
          
          CALL warning&
               &("FVCOM ATTEMPTING TO GUESS VALUES FOR A REGION WHERE ALL DATA ARE MASKED",&
               & "PLEASE CHECK RESULTS CAREFULLY AND GET BETTER DATA IF POSSIBLE")
          
          mk_cnt = abs(sum(n_found))


          ALLOCATE(n_order(mk_cnt))

          ! IMPORTANT INTIALIZATION TO PREVENT COUNTING BOUNDARY AS
          ! FOUND CELL!
          n_found(0) = -2

          CALL grid_neighbor_index(n_found,n_index,n_cnt,n_order)

          weights%N_INDEX => n_index
          weights%N_CNT => n_cnt
          weights%N_ORDER => n_order
          weights%N_FOUND => n_found

          nullify(n_index)
          nullify(n_cnt)
          nullify(n_order)
          nullify(n_found)

       ELSE

          DEALLOCATE(n_found)
          DEALLOCATE(n_index)
          DEALLOCATE(n_cnt)

       END IF

    ELSEIF( n_type==type_element) THEN
       IF(any(n_found < 0)) THEN
          CALL warning&
            &("FVCOM ATTEMPTING TO GUESS VALUES FOR A REGION WHERE ALL DATA ARE MASKED",&
            & "PLEASE CHECK RESULTS CAREFULLY AND GET BETTER DATA IF POSSIBLE")
       
          mk_cnt = abs(sum(n_found))
          

          ALLOCATE(n_order(mk_cnt))
          
          ! IMPORTANT INTIALIZATION TO PREVENT COUNTING BOUNDARY AS
          ! FOUND CELL!
          n_found(0) = -2

          CALL grid_neighbor_index(n_found,n_index,n_cnt,n_order)

          weights%N_INDEX => n_index
          weights%N_CNT => n_cnt
          weights%N_ORDER => n_order
          weights%N_FOUND => n_found

          nullify(n_index)
          nullify(n_cnt)
          nullify(n_order)
          nullify(n_found)

       ELSE

          DEALLOCATE(n_found)
          DEALLOCATE(n_index)
          DEALLOCATE(n_cnt)


       END IF

    END IF


    IF(dbg_set(dbg_sbr)) THEN
       WRITE(ipt,*) "///////////////////////////////////////////////////"
       WRITE(ipt,*) "//// TIMING REPORT FROM BILINEAR INTERP ///////////"
       WRITE(ipt,*) "///////////////////////////////////////////////////"
       CALL watch_report(min_w,ipt,"MIN LOCATION")

       CALL watch_report(box_w,ipt,"GET BOX")

       CALL watch_report(cond_w,ipt,"CONDITION")

       CALL watch_report(case_w,ipt,"WEIGHTS")

       CALL watch_lap(tot_w,ipt,"TOTAL")
       WRITE(ipt,*) "///////////////////////////////////////////////////"
    END IF

    DEALLOCATE(x_bnd)
    DEALLOCATE(y_bnd)
    DEALLOCATE(rsq)


    IF(bounds_flag) CALL warning&
         &("FVCOM GRID IS OUT OF BOUNDS FOR FORCING",&
         & "USING CONSTANT VALUE OUTSIDE OF FORCING REGION")

  END SUBROUTINE setup_interp_bilinear_p

  ! CALCULATE WEIGHTS FOR CURVILINEAR INTERP
  SUBROUTINE interp_0(Xbox,Ybox,X,Y,wghts)
    IMPLICIT NONE 

    REAL(SP), INTENT(IN) :: Xbox(4)
    REAL(SP), INTENT(IN) :: Ybox(4)
    REAL(SP), INTENT(IN) :: X
    REAL(SP), INTENT(IN) :: Y
    REAL(SP), INTENT(OUT) :: wghts(4)


    ! FOR GENERAL CASE
    REAL(SP) :: AB,BB,AT,BT ! SLOPE AND INTERCEPT, TOP and BOTTOM
    REAL(SP) :: AL,BL,AR,BR ! SLOPE AND INTERCEPT, TOP and BOTTOM

    ! DELTA X/Y for each edge
    REAL(SP) :: DXL, DXR, DXT, DXB
    REAL(SP) :: DYL, DYR, DYT, DYB

    ! X/Y INTERCEPT ON TOP AND BOTTOM EDGES
    REAL(sp) :: XT,YT,XB,YB

    ! SLOPE, INTERPCEPT and DELTA's FOR AN 'AVERAGE' LINE BETWEEN TOP
    ! AND BOTTOM 
    REAL(SP) :: XLA,XRA
    REAL(SP) :: DYM, DXM
    REAL(SP) :: AM,BM
    REAL(SP) :: RADL,RADR,RADM

    ! DISTANCES
    REAL(SP) :: DR1, DR2,DR3,DR4,DR5,DR6
    REAL(SP) :: DRM,DRT,DRB
    
    ! FOR SPHERICAL MODEL, DO INTERCEPT CALCULATION IN LATITUDE AND LONGITUDE!

    ! DELTA X
    dxt = xbox(2) -  xbox(1)
    
    dxb = xbox(4) -  xbox(3)
    
    dxl = xbox(1) -  xbox(4)
    
    dxr = xbox(2) -  xbox(3)
    
    ! DELTA Y
    dyt = ybox(2) -  ybox(1)
    
    dyb = ybox(4) -  ybox(3)
    
    dyl = ybox(1) -  ybox(4)
    
    dyr = ybox(2) -  ybox(3)
    

    !SLOPE AND INTERCEPT FOR TOP AND BOTTOM OF QUAD
    at = dyt/dxt
    bt = ybox(2) - at*xbox(2)

    !       write(ipt,*) "AT=",AT,"; BT=",BT


    ab = dyb/dxb
    bb = ybox(4) - ab*xbox(4)

    !       write(ipt,*) "AB=",AB,"; BB=",BB


    !GUESS AT X INTERCEPT ON TOP AND BOTTOM LINES
    xt = x

    xb = x


    ! IF THE SIDES ARE NOT NEARLY VERTICAL GET A NEW ESTIMATE FOR XT, XB
    IF(abs(dyl) < lr_tol * abs(dxl) .OR. abs(dyr) < lr_tol*abs(dxr) ) THEN

       !XLA = (XBOX(BOX(1,1),BOX(1,2)) + XBOX(BOX(4,1),BOX(4,2)) ) / 2.0_SP

       IF(abs(dyl) > lr_tol * abs(dxl)) THEN ! IF THE LEFT SIDE IS HORIZONTAL
          xla = max(xbox(1),xbox(4))
       ELSE
          ! USE INTERCEPT AT Y!
          al = dyl/dxl
          bl = ybox(4) - al*xbox(4)

          xla = (y -bl)/al
       END IF


       IF(abs(dyr) > lr_tol * abs(dxr)) THEN ! IF THE LEFT SIDE IS HORIZONTAL
          xra = min(xbox(2),xbox(3))
       ELSE
          ! USE INTERCEPT AT Y!
          ar = dyr/dxr
          br = ybox(2) - ar*xbox(2)

          xra = (y -br)/ar
       END IF


!       WRITE(IPT,*) "XLA=",XLA,"; XRA=",XRA

       !             write(ipt,*) "Left slope=",(DYL)/(DXL)
       !             write(ipt,*) "Right slope=",(DYR)/(DXR)


       radl = atan2(dyl,dxl)
       radr = atan2(dyr,dxr)

       !             write(ipt,*) "Left RAD=",RADL
       !             write(ipt,*) "Right RAD=",RADR


       ! GET AN AVERAGE DX and DY
       radm  = (xla -xt)*(radr)/(xla-xra) + (xt - xra)*radl/(xla-xra)

       !             write(ipt,*) "RADM=",RADM

       am = tan(radm)

       ! IF THE RESULTING SLOPE IS NOT VERTICAL
       IF(am < lr_tol )THEN
          bm = y - am *x

          xt = (bt-bm)/(am - at)

          xb = (bb-bm)/(am - ab)

          !                write(ipt,*) "Middle! AM=",AM


       END IF

    END IF

    !GET YT, YB FROM XB,XT
    yb = ab*xb+bb 

    yt = at*xt+bt 


!    WRITE(IPT,*) "XT=",XT,"; YT=",YT

!    WRITE(IPT,*) "XB=",XB,"; YB=",YB



!---------------------------------------------------------------

    ! TOP AND BOTTOM LENGTH
    drt = sqrt(dxt**2+ dyt**2)
    drb = sqrt(dxb**2+ dyb**2)

    ! TOP SEGMENTS
    dr1 = sqrt( (xbox(1) - xt)**2 + (ybox(1) - yt)**2)
    dr2 = drt - dr1

    ! BOTTOM SEGMENTS
    dr3 = sqrt( (xbox(3) - xb)**2 + (ybox(3) - yb)**2)
    dr4 = drb - dr3

    ! MIDDLE LINE LENGTH
    drm = sqrt((yb-yt)**2 + (xb-xt)**2)

    ! MIDDLE SEGMENTS
    dr5 = sqrt((y-yt)**2 + (x-xt)**2)
    dr6 = drm-dr5

!---------------------------------------------------------------
!---------------------------------------------------------------

!    WRITE(IPT,*) "DRT",DRT
!    WRITE(IPT,*) "DRB",DRB
!    WRITE(IPT,*) "DRM",DRM
!    WRITE(IPT,*) "DR1",DR1
!    WRITE(IPT,*) "DR2",DR2
!    WRITE(IPT,*) "DR3",DR3
!    WRITE(IPT,*) "DR4",DR4
!    WRITE(IPT,*) "DR5",DR5
!    WRITE(IPT,*) "DR6",DR6




    wghts(1) =(dr6*dr2) / (drt*drm)
    wghts(2) =(dr6*dr1) / (drt*drm)
    wghts(3) =(dr5*dr4) / (drb*drm)
    wghts(4) =(dr5*dr3) / (drb*drm)

  END SUBROUTINE interp_0

  SUBROUTINE interp_neg(Xbox,Ybox,X,Y,wghts)
    IMPLICIT NONE 

    REAL(SP), INTENT(IN) :: Xbox(4)
    REAL(SP), INTENT(IN) :: Ybox(4)
    REAL(SP), INTENT(IN) :: X
    REAL(SP), INTENT(IN) :: Y
    REAL(SP), INTENT(OUT) :: wghts(4)


    ! FOR GENERAL CASE
    REAL(SP) :: A1,B1,A2,B2,A3,B3 ! SLOPE AND INTERCEPT, TOP and BOTTOM

    ! DELTA X/Y for each edge
    REAL(SP) :: DX1, DX2, DX3
    REAL(SP) :: DY1, DY2, DY3

    ! X/Y INTERCEPT ON EDGES
    REAL(sp) :: X1,Y1,X2,Y2

    ! SLOPE, INTERPCEPT and DELTA's FOR AN 'AVERAGE' LINE BETWEEN TOP
    ! AND BOTTOM 

    ! DISTANCES
    REAL(SP) :: DRA, DRB,DRC,DRD,DRE,DRF
    REAL(SP) :: DR1,DR2,DR3


    REAL(SP) :: Xt(3),Yt(3)


    xt = xbox(1:3)
    yt = ybox(1:3)

    IF(isintri(x,y,xt,yt))THEN
       ! BOTTOM LEFT IS MISSING
       ! FOR THE TWO EDGES WITH DATA
       ! DELTA X
       dx1 = xbox(2) -  xbox(1) !SIDE ONE
       dx2 = xbox(2) -  xbox(3) !SIDE TWO

       ! DELTA Y 
       dy1 = ybox(2) -  ybox(1)
       dy2 = ybox(2) -  ybox(3)


       ! ACCROSS THE CENTER - THE TWO VALID POINTS!
       dx3 = xbox(3) -  xbox(1) 
       dy3 = ybox(3) -  ybox(1)


       ! USING THE SAME SLOPE AS THE MIDSECTION, GET THE INTERCEPT OF
       ! THE POINT X/Y ON THE VALID EDGES

       ! IF THE MID SECTION IS NEARLY VERTICAL - GREAT
       IF(abs(dy3) > tb_tol *abs(dx3)) THEN
          x1 = x
          x2 = x

          ! assume the other sides are not also vertical!
          a1 = dx1/dx1
          b1 = ybox(2) - a1*xbox(2)

          ! GET YT 
          y1 = a1*x1+b1 


          a2 = dy2/dx2
          b2 = ybox(2) - a2*xbox(2)
          ! GET YB
          y2 = a2*x2+b2 

       ELSE
          ! GET THE SLOPE OF THE MIDDLE LINE, AND THE EQ FOR THE LINE
          ! THROUGH THE POINT X/Y WITH THAT SLOPE
          a3 = dy3 / dx3
          b3 = y - a3 *x

!          WRITE(IPT,*) "A3=",A3,"; BM=",B3

          ! SIDE ONE
          IF(abs(dy1) > tb_tol *abs(dx1)) THEN
             x1 = xbox(2)
             y1 = a3*x1+b3
          ELSE
             ! ASSUME THE LINES ARE NOT NEARLY PARALLEL
             a1 = dy1/dx1
             b1 = ybox(2) - a1*xbox(2)

             x1 = (b1-b3)/(a3 - a1)

             y1 = a1*x1+b1 
          END IF


          IF(abs(dy2) > tb_tol *abs(dx2)) THEN
             x2 = xbox(2)
             y2 = a3*x2+b3
          ELSE
             ! ASSUME THE LINES ARE NOT NEARLY PARALLEL

             a2 = dy2/dx2
             b2 = ybox(2) - a2*xbox(2)
             x2 = (b2-b3)/(a3 - a2)

             y2 = a2*x2+b2 
          END IF

       END IF

!---------------------------------------------------------------
       ! LENGTH OF SIDES
       dr1 = sqrt(dx1**2+ dy1**2) !SIDE ONE
       dr2 = sqrt(dx2**2+ dy2**2) !SIDE TWO

       ! SIDE ONE SEGMENTS
       dra = sqrt( (xbox(1) - x1)**2 + (ybox(1) - y1)**2)
       drb = dr1 - dra

       ! SIDE TWO SEGMENTS
       drc = sqrt( (xbox(2) - x2)**2 + (ybox(2) - y2)**2)
       drd = dr2 - drc


       dr3 = sqrt((y2-y1)**2 + (x2-x1)**2)
       dre = sqrt((y-y1)**2 + (x-x1)**2)
       drf = dr3-dre

!       write(ipt,*) "Middle! A3=",A3

!       WRITE(IPT,*) "X1=",X1,"; YT=",Y1


!       WRITE(IPT,*) "XB=",X2,"; YB=",Y2


!       write(ipt,*) "DR1=",DR1

!       write(ipt,*) "DR2=",DR2
!       write(ipt,*) "DR3=",DR3

       wghts(1) = (drf*drb) / (dr1*dr3)
       wghts(2) = (drf*dra) / (dr1*dr3) + (dre*drd) / (dr2*dr3)
       wghts(3) = (dre*drc) / (dr2*dr3)
       wghts(4) = 0.0_sp

    ELSE

!---------------------------------------------------------------
       dr3 = sqrt( (xbox(1) - x)**2 + (ybox(1) - y)**2)

       dr1 = sqrt( (xbox(3) - x)**2 + (ybox(3) - y)**2)

       dra = (dr1 +dr3)
       ! NORMALIZE TO ONE

       wghts(1) = dr1/dra
       wghts(2) = 0.0_sp
       wghts(3) = dr3/dra
       wghts(4) = 0.0_sp


    END IF


  END SUBROUTINE interp_neg
    

  subroutine interp_bilinear_a(Zin,WEIGHTS,zout)
    implicit none
    TYPE(interp_weights), INTENT(IN) :: WEIGHTS
    real(SP), ALLOCATABLE ,TARGET, INTENT(IN):: zin(:,:)
    real(SP), ALLOCATABLE,TARGET, INTENT(INOUT) :: Zout(:)


    REAL(SP), POINTER :: ZinP(:,:),ZoutP(:)

    IF (allocated(zin)) THEN
       zinp => zin
    ELSE
       CALL fatal_error &
            &("INTERP: ZIN IS NOT ALLOCATED")
    END IF

    IF (allocated(zout)) THEN
       zoutp => zout   !(lbound(zout):ubound(zout))
    ELSE
       CALL fatal_error &
            &("INTERP: ZOUT IS NOT ALLOCATED")
    END IF

    CALL interp_bilinear_p(zinp,weights,zoutp)


  end subroutine interp_bilinear_a

  subroutine interp_bilinear_p(Zin,WEIGHTS,zout)
    implicit none
    TYPE(interp_weights), INTENT(IN) :: WEIGHTS
    real(SP), POINTER , INTENT(IN):: zin(:,:)
    real(SP), POINTER, INTENT(INOUT) :: Zout(:)
    TYPE(r2pts) :: PTW
    integer  :: i,j,x,y, kk
    real(SP) :: tt

    ! WEIGHTS ARE NORMALIZED TO ONE ALREADY: DO NOT DIVIDE BY THE SUM!

    DO j = 1,ubound(zout,1)

       ptw=weights%PTW(j)

       tt=0.0_sp
       do i=1,4
          x=ptw%quad_num(i,1)
          y=ptw%quad_num(i,2)
          IF(x*y == 0) cycle
          tt=zin(x,y)*ptw%weights(i,1)+tt
       end do

       zout(j)=tt

    END DO


    IF(ASSOCIATED(weights%N_INDEX)) THEN
       ! THERE IS MASKED DATA - USE THE NEIGHBOR INDEX!
       DO i = 1,ubound(weights%N_ORDER,1)
          
          j = weights%N_ORDER(i)
          tt=0.0_sp
          
          DO kk = 1, weights%N_CNT(j)
             tt = tt + zout(weights%N_INDEX(j,kk))
          END DO

          zout(j) = tt / real(weights%N_CNT(j),sp)
       END DO
    END IF


  end subroutine interp_bilinear_p









 
  SUBROUTINE setup_interp_quad_a(Xin,Yin,Xout,Yout,WEIGHTS, rzero)
    IMPLICIT NONE 
    REAL(SP), ALLOCATABLE, TARGET, INTENT(IN) :: Xin(:)
    REAL(SP), ALLOCATABLE, TARGET, INTENT(IN) :: Yin(:)
    REAL(SP), ALLOCATABLE, TARGET, INTENT(IN) :: Xout(:)
    REAL(SP), ALLOCATABLE, TARGET, INTENT(IN) :: Yout(:)
    REAL(SP), OPTIONAL :: rzero

    TYPE(interp_weights), INTENT(OUT) :: WEIGHTS

    REAL(SP), POINTER :: XinP(:)
    REAL(SP), POINTER :: YinP(:)
    REAL(SP), POINTER :: XoutP(:)
    REAL(SP), POINTER :: YoutP(:)

    NULLIFY(xinp,yinp,xoutp,youtp)


    IF(ALLOCATED(xin)) xinp => xin

    IF(ALLOCATED(yin)) yinp => yin

    IF(ALLOCATED(xout)) xoutp => xout

    IF(ALLOCATED(yout)) youtp => yout

    IF(PRESENT(rzero)) THEN
       CALL setup_interp_quad_p(xinp,yinp,xoutp,youtp,weights,rzero)
    ELSE
       CALL setup_interp_quad_p(xinp,yinp,xoutp,youtp,weights)
    END IF
  END SUBROUTINE setup_interp_quad_a



  SUBROUTINE setup_interp_quad_p(Xin,Yin,Xout,Yout,WEIGHTS,rzero)
    IMPLICIT NONE 
    REAL(SP), POINTER, INTENT(IN) :: Xin(:)
    REAL(SP), POINTER, INTENT(IN) :: Yin(:)
    REAL(SP), POINTER, INTENT(IN) :: Xout(:)
    REAL(SP), POINTER, INTENT(IN) :: Yout(:)
    REAL(SP), OPTIONAL :: rzero

    REAL(SP), ALLOCATABLE :: Xvec(:),Yvec(:)

    TYPE(interp_weights), INTENT(OUT) :: WEIGHTS

    INTEGER :: I, status, lb, ub

!---------------------------------------------------------------

    ! CHECK INPUT ALLOCATION STATUS
    if(.not.associated(xin)) CALL fatal_error&
         &("SETUP_INTERP: INPUT ARGUMENTS MUST BE ALLOCATED!")
    if(.not.associated(yin)) CALL fatal_error&
         &("SETUP_INTERP: INPUT ARGUMENTS MUST BE ALLOCATED!")
    if(.not.associated(xout)) CALL fatal_error&
         &("SETUP_INTERP: INPUT ARGUMENTS MUST BE ALLOCATED!")
    if(.not.associated(yout)) CALL fatal_error&
         &("SETUP_INTERP: INPUT ARGUMENTS MUST BE ALLOCATED!")
    

    ! GET THE DIMENSIONS
 
    weights%Nout = ubound(xout,1)
    weights%Nin  = ubound(xin,1)

    ALLOCATE(weights%PTW(weights%Nout), stat=status)
    IF(status /= 0) CALL fatal_error("SETUP_INTERP: COULD NOT ALLOCATE SPACE")

    ALLOCATE(weights%INDEX(weights%Nin), stat=status)
    IF(status /= 0) CALL fatal_error("SETUP_INTERP: COULD NOT ALLOCATE SPACE")


    CALL sortrx(weights%Nin,xin,weights%INDEX)


    ALLOCATE(xvec(weights%Nin), stat=status)
    IF(status /= 0) CALL fatal_error&
         & ("SETUP_INTERP: COULD NOT ALLOCATE SPACE")


    ALLOCATE(yvec(weights%Nin), stat=status)
    IF(status /= 0) CALL fatal_error&
         & ("SETUP_INTERP: COULD NOT ALLOCATE SPACE")
    

    DO i = 1,weights%Nin
       xvec(i) = xin(weights%INDEX(i))
       yvec(i) = yin(weights%INDEX(i))
    END DO
    

    IF(PRESENT(rzero)) THEN
       DO i = 1, weights%Nout
          CALL gen_wts(xvec,yvec,xout(i),yout(i),weights%INDEX,weights&
               &%PTW(i),rzero)
       END DO
    ELSE
       DO i = 1, weights%Nout
          CALL gen_wts(xvec,yvec,xout(i),yout(i),weights%INDEX,weights&
               &%PTW(i))
       END DO
    END IF
       

  END SUBROUTINE setup_interp_quad_p


  ! THIS CODE IS SETUP TO WORK FOR REGULAR GRIDDED DATA
  ! TO USE IT ON OBSERVATION DATA IT WOULD NEED TO BE MODIFIED TO
  ! HANDLE STRANGE CASES WHERE THERE ARE CLUSTERS OF DATA POINTS
  SUBROUTINE gen_wts(Xvec,Yvec,Xout,Yout,INDX,PTW,rzero)
    IMPLICIT NONE
    real(SP), INTENT(in) :: Xout,Yout
    real(SP), ALLOCATABLE, INTENT(IN) :: Xvec(:),Yvec(:)
    INTEGER, ALLOCATABLE, INTENT(IN) :: INDX(:)
    TYPE(r2pts), INTENT(OUT) :: PTW
    REAL(SP), OPTIONAL :: rzero

    INTEGER :: I,Q,K,NB,J
    
    Real(SP) :: DELX, DELY, D, SUMWGHT

    REAL(SP) :: D_MAX(4)
    REAL(SP) :: DIST(4,2)
    INTEGER :: NDM

    nb = ubound(xvec,1)

    ptw%WEIGHTS=0
    ptw%QUAD_NUM=0

    d_max = huge(d)
    dist  = huge(d)

    ndm=1

    k = 0
    DO i=1,nb

       delx = xvec(i) - xout

       ! SINCE XVEC IS IN INCREASING ORDER, ONCE DELX IS GREATER THAN
       ! SEARCH WE ARE DONE
       If (delx.GT.search) EXIT

       If (abs(delx).LE.search) Then
          dely = yvec(i) - yout
          d = sqrt(delx**2+dely**2)
          
          ! WE FOUND A POINT IN THE SEARCH RADIUS
          If (d.LE.search) Then
             
             k = k +1
             ! DEAL WITH D == 0.0 here?
             q=quadrant(delx,dely)

             ! REPLACE THAT VALUE IF IT IS CLOSER THAN THE LAST ONE
             IF(d .LT. d_max(q))THEN
                ! RECORD THE DISTANCE AND UPDATE THE MAX FOR THAT QUADRANT
                ndm= maxloc(dist(q,:),1)
                dist(q,ndm) = d
                d_max(q) = maxval(dist(q,:))
                
                ! RECORD THE INDEX IN QUAD_NUM
                ptw%QUAD_NUM(q,ndm) = i
             END IF
          End If
       End If
    END DO

    if (k == 0) call warning("FOUND NO Values in interpolation serach radius!")

    IF (PRESENT(rzero)) THEN
       dist = dist/rzero
    END IF


    DO i = 1,4
       DO k =1 , 2
          
          IF (dist(i,k) < huge(d)) THEN
             ptw%WEIGHTS(i,k) = 1/(.5_sp+dist(i,k)**2)
          END IF
       END DO
    END DO
    
    ! NORMALIZE THE TOTAL WEIGHT TO ONE
    sumwght = sum(ptw%WEIGHTS)
    ptw%WEIGHTS = ptw%WEIGHTS / sumwght

  END SUBROUTINE gen_wts


  INTEGER FUNCTION quadrant(DX,DY)
    IMPLICIT NONE
    REAL(sp), INTENT(IN):: dx,dy
    INTEGER :: ibit1, ibit2
    
    ibit1 = (int(sign(1.,dx))+1) / 2
    ibit2 = (int(sign(1.,dy))+1) / 2
    !    quadrant number:
    quadrant = 2 * ibit1 + ibit2 + 1
  END FUNCTION quadrant
  

!*****************************************************************
!C Here's a hybrid QuickSort I wrote a number of years ago.  It's
!C based on suggestions in Knuth, Volume 3, and performs much better
!C than a pure QuickSort on short or partially ordered input arrays.  

  SUBROUTINE sortrx(N,DATA,INDX)
!C===================================================================
!C
!C     SORTRX -- SORT, Real input, indx output
!C
!C
!C     Input:  N     INTEGER
!C             DATA  REAL
!C
!C     Output: INDX INTEGER (DIMENSION N)
!C
!C This routine performs an in-memory sort of the first N elements of
!C array DATA, returning into array INDX the indices of elements of
!C DATA arranged in ascending order.  Thus,
!C
!C    DATA(INDX(1)) will be the smallest number in array DATA;
!C    DATA(INDX(N)) will be the largest number in DATA.
!C
!C The original data is not physically rearranged.  The original order
!C of equal input values is not necessarily preserved.
!C
!C===================================================================
!C
!C SORTRX uses a hybrid QuickSort algorithm, based on several
!C suggestions in Knuth, Volume 3, Section 5.2.2.  In particular, the
!C "pivot key" [my term] for dividing each subsequence is chosen to be
!C the median of the first, last, and middle values of the subsequence;
!C and the QuickSort is cut off when a subsequence has 9 or fewer
!C elements, and a straight insertion sort of the entire array is done
!C at the end.  The result is comparable to a pure insertion sort for
!C very short arrays, and very fast for very large arrays (of order 12
!C micro-sec/element on the 3081K for arrays of 10K elements).  It is
!C also not subject to the poor performance of the pure QuickSort on
!C partially ordered data.
!C
!C Created:  15 Jul 1986  Len Moss
!C
!C===================================================================
    
    INTEGER, INTENT(IN)::  N
    REAL(SP)    :: DATA(N)
    INTEGER :: INDX(N)
    
    INTEGER  :: LSTK(31),RSTK(31),ISTK
    INTEGER  :: L,R,I,J,P,INDXP,INDXT
    REAL     :: DATAP
    
!C     QuickSort Cutoff
!C
!C     Quit QuickSort-ing when a subsequence contains M or fewer
!C     elements and finish off at end with straight insertion sort.
!C     According to Knuth, V.3, the optimum value of M is around 9.
 
    INTEGER, PARAMETER ::   M = 9
    
!C===================================================================
!C
!C     Make initial guess for INDX
 
    DO 50 i=1,n
       indx(i)=i
50     CONTINUE
 
!C     If array is short, skip QuickSort and go directly to
!C     the straight insertion sort.
 
       IF (n.LE.m) GOTO 900
       
!C===================================================================
!C
!C     QuickSort
!C
!C     The "Qn:"s correspond roughly to steps in Algorithm Q,
!C     Knuth, V.3, PP.116-117, modified to select the median
!C     of the first, last, and middle elements as the "pivot
!C     key" (in Knuth's notation, "K").  Also modified to leave
!C     data in place and produce an INDX array.  To simplify
!C     comments, let DATA[I]=DATA(INDX(I)).
 
!C Q1: Initialize
       istk=0
       l=1
       r=n
       
200    CONTINUE
 
!C Q2: Sort the subsequence DATA[L]..DATA[R].
!C
!C     At this point, DATA[l] &lt;= DATA[m] &lt;= DATA[r] for all l &lt; L,
!C     r &gt; R, and L &lt;= m &lt;= R.  (First time through, there is no
!C     DATA for l &lt; L or r &gt; R.)
 
       i=l
       j=r
 
!C Q2.5: Select pivot key
!C
!C     Let the pivot, P, be the midpoint of this subsequence,
!C     P=(L+R)/2; then rearrange INDX(L), INDX(P), and INDX(R)
!C     so the corresponding DATA values are in increasing order.
!C     The pivot key, DATAP, is then DATA[P].
 
       p=(l+r)/2
       indxp=indx(p)
       datap=DATA(indxp)
       
       IF (DATA(indx(l)) .GT. datap) THEN
          indx(p)=indx(l)
          indx(l)=indxp
          indxp=indx(p)
          datap=DATA(indxp)
       ENDIF
       
       IF (datap .GT. DATA(indx(r))) THEN
          IF (DATA(indx(l)) .GT. DATA(indx(r))) THEN
             indx(p)=indx(l)
             indx(l)=indx(r)
          ELSE
             indx(p)=indx(r)
          ENDIF
          indx(r)=indxp
          indxp=indx(p)
          datap=DATA(indxp)
       ENDIF
 
!C     Now we swap values between the right and left sides and/or
!C     move DATAP until all smaller values are on the left and all
!C     larger values are on the right.  Neither the left or right
!C     side will be internally ordered yet; however, DATAP will be
!C     in its final position.
 
300    CONTINUE
 
!C Q3: Search for datum on left &gt;= DATAP
!C
!C     At this point, DATA[L] &lt;= DATAP.  We can therefore start scanning
!C     up from L, looking for a value &gt;= DATAP (this scan is guaranteed
!C     to terminate since we initially placed DATAP near the middle of
!C     the subsequence).
 
       i=i+1
       IF (DATA(indx(i)).LT.datap) GOTO 300
       
400    CONTINUE
       
!C Q4: Search for datum on right &lt;= DATAP
!C
!C     At this point, DATA[R] &gt;= DATAP.  We can therefore start scanning
!C     down from R, looking for a value &lt;= DATAP (this scan is guaranteed
!C     to terminate since we initially placed DATAP near the middle of
!C     the subsequence).
 
       j=j-1
       IF (DATA(indx(j)).GT.datap) GOTO 400
 
!C Q5: Have the two scans collided?
 
       IF (i.LT.j) THEN
 
!C Q6: No, interchange DATA[I] &lt;--&gt; DATA[J] and continue
          
          indxt=indx(i)
          indx(i)=indx(j)
          indx(j)=indxt
          GOTO 300
       ELSE
          
!C Q7: Yes, select next subsequence to sort
!C
!C    At this point, I &gt;= J and DATA[l] &lt;= DATA[I] == DATAP &lt;= DATA[r],
!C     for all L &lt;= l &lt; I and J &lt; r &lt;= R.  If both subsequences are
!C     more than M elements long, push the longer one on the stack and
!C     go back to QuickSort the shorter; if only one is more than M
!C     elements long, go back and QuickSort it; otherwise, pop a
!C     subsequence off the stack and QuickSort it.
 
          IF (r-j .GE. i-l .AND. i-l .GT. m) THEN
             istk=istk+1
             lstk(istk)=j+1
             rstk(istk)=r
             r=i-1
          ELSE IF (i-l .GT. r-j .AND. r-j .GT. m) THEN
             istk=istk+1
             lstk(istk)=l
             rstk(istk)=i-1
             l=j+1
          ELSE IF (r-j .GT. m) THEN
             l=j+1
          ELSE IF (i-l .GT. m) THEN
             r=i-1
          ELSE
!C Q8: Pop the stack, or terminate QuickSort if empty
             IF (istk.LT.1) GOTO 900
             l=lstk(istk)
             r=rstk(istk)
             istk=istk-1
          ENDIF
          GOTO 200
       ENDIF
 
900    CONTINUE
 
!C===================================================================
!C
!C Q9: Straight Insertion sort
 
       DO 950 i=2,n
          !         if(mod(I,10000).eq.0) print*,'i=',i
          IF (DATA(indx(i-1)) .GT. DATA(indx(i))) THEN
             indxp=indx(i)
             datap=DATA(indxp)
             p=i-1
920          CONTINUE
             indx(p+1) = indx(p)
             p=p-1
             IF (p.GT.0) THEN
                IF (DATA(indx(p)).GT.datap) GOTO 920
             ENDIF
             indx(p+1) = indxp
          ENDIF
950       CONTINUE
          
!C===================================================================
!C
!C     All done
          RETURN
!       END DO
!    END DO
  END SUBROUTINE sortrx
!******************************


  SUBROUTINE interp_quad_a(zin,Weights,zout)
    IMPLICIT NONE
    REAL(SP), ALLOCATABLE, TARGET, INTENT(IN) :: Zin(:)
    TYPE(interp_weights), INTENT(in) :: WEIGHTS
    REAL(SP), ALLOCATABLE, TARGET :: Zout(:)

    REAL(SP), POINTER :: ZinP(:),ZoutP(:)

    IF (allocated(zin)) THEN
       zinp => zin
    ELSE
       CALL fatal_error &
          &("INTERP: ZIN IS NOT ALLOCATED")
    END IF

    IF (allocated(zout)) THEN
       zoutp => zout   !(lbound(zout):ubound(zout))
    ELSE
       CALL fatal_error &
          &("INTERP: ZOUT IS NOT ALLOCATED")
    END IF
 
    CALL interp_quad_p(zinp,weights,zoutp)

  END SUBROUTINE interp_quad_a
  
  SUBROUTINE interp_quad_p(zin,Weights,zout)
     REAL(SP), POINTER :: Zin(:)
     TYPE(interp_weights), INTENT(in) :: WEIGHTS
     REAL(SP), POINTER :: Zout(:)

     REAL(SP), ALLOCATABLE :: Zvec(:)
     
     integer :: I,status, msze,nsze, lb, ub

     IF (.NOT. ASSOCIATED(zin))CALL fatal_error &
          &("INTERP: ZIN IS NOT ALLOCATED")

     ! CAN'T CHECK ALLOCATION STATUS OF INTENT OUT VARIABLE!

     
     IF (weights%Nin /= ubound(zin,1)) CALL fatal_error &
          &("INTERP: THE DATA INPUT SIZE DOES NOT MATCH THE WEIGHT!")

     ! MAKE VECTOR IN RIGHT ORDER...
     ALLOCATE(zvec(0:weights%Nin),stat=status)
     if(status /= 0) CALL fatal_error("INPTERP: COULD NOT ALLOCATE SPACE")
     zvec = 0.0_sp

     DO i = 1,weights%Nin
        zvec(i) = zin(weights%INDEX(i))
     END DO

     DO i = 1,weights%Nout

        CALL interp_weigh(zvec,weights%PTW(i),zout(i))

     END DO

     !DEALLOCATE IS AUTOMATIC FOR ALLOCATABLES

   end subroutine interp_quad_p

! VECTORIZE THIS - THERE IS NO REASON TO USE AN ARRAY TO STORE THE WEIGHTS!
   subroutine interp_weigh(Zvec,PTW,zval)
     implicit none
     real(SP), INTENT(OUT) :: ZVAL
     TYPE(r2pts), INTENT(IN) :: PTW
     real(SP), ALLOCATABLE , INTENT(IN):: zvec(:)
     integer  :: i,j,n

     ! WEIGHTS ARE NORMALIZED TO ONE ALREADY: DO NOT DIVIDE BY THE SUM!

     zval=0.0_sp
     do i=1,4
        do j=1,2
           n=ptw%quad_num(i,j)
           zval=zvec(n)*ptw%weights(i,j)+zval
        end do
     end do

   end subroutine interp_weigh
   

END MODULE mod_interp