interpolate.F

Go to the documentation of this file.

#include "cppdefs.h"
 
#undef CUBIC_MASKED    /* Bicubic interpolation with mask, needs work */
 
      MODULE roms_interpolate_mod
!
!git $Id$
!=======================================================================
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                            Hernan G. Arango   !
!========================================== Alexander F. Shchepetkin ===
!                                                                      !
!  This module contains several all-purpose generic interpolations     !
!  routines:                                                           !
!                                                                      !
!  Routines:                                                           !
!                                                                      !
!    cinterp2d     Bicubic  interpolation for any 2D field.            !
!    linterp2d     Bilinear interpolation for any 2D field.            !
!    hindices      Finds model grid cell for any datum.                !
!    try_range     Binary search of model grid cell for any datum.     !
!    inside        Closed polygon datum search.                        !
!                                                                      !
!  Overloading routines using "roms_interp_type" structure, S, that    !
!  can be used in the ROMS-JEDI interface.                             !
!                                                                      !
!    - Interpolate a ROMS 2D state field to observation locations      !
!      (ObsVetting is an optional output argument)                     !
!                                                                      !
!      roms_datum_interp (S, fld2d, datum(:), method)                  !
!                                                                      !
!    - Interpolate a ROMS 3D state field to observation locations      !
!      (ObsVetting is an optional output argument)                     !
!                                                                      !
!      roms_datum_interp (S, fld3d, zfld3d, zlocs(:), datum(:), method)!
!                                                                      !
!    - Interpolate a ROMS 3D state field to observation locations      !
!      level-by-level and returns datum(1:nlevs,1:nlocs)               !
!                                                                      !
!      roms_datum_interp (S, fld3d, zfld3d, datum(:,:), method)        !
!                                                                      !
!    - Interpolate a ROMS 2D source field to a new destination         !
!      location (2D remapping)                                         !
!                                                                      !
!      roms_horiz_interp (S, fld2d_src, fld2d_dst, method)             !
!                                                                      !
!    - Interpolate a ROMS 3D source field to a new destination         !
!      locations level-by-level (3D remapping)                         !
!                                                                      !
!      roms_horiz_interp (S, fld2d_src, fld2d_dst, method)             !
!                                                                      !
!  Managing routines for "roms_interp_type" structure, S:              !
!                                                                      !
!    - Creates a ROMS interpolation object                             !
!                                                                      !
!      roms_interp_create (ng, S)                                      !
!                                                                      !
!    - Deallocates pointer arrays in the interpolation structure       !
!                                                                      !
!      roms_interp_delete (S)                                          !
!                                                                      !
!    - Computes the horizontal fractional coordinates S%x_dst and      !
!      S%y_dst of the source cells containing the destination values   !
!                                                                      !
!      roms_interp_fractional (S)                                      !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_grid
      USE mod_ncparam
      USE mod_scalars
!
#ifdef DISTRIBUTE
      USE distribute_mod, ONLY : mp_assemble
#endif
      USE mod_iounits,    ONLY : stdout
      USE strings_mod,    ONLY : founderror
!
      implicit none
!
!-----------------------------------------------------------------------
!  Set ROMS interpolation structure for Object Oriented Programming.
!-----------------------------------------------------------------------
!
      TYPE :: roms_interp_type
 
        logical :: rectangular = .false.            ! plaid grid switch
!
        integer :: ng                               ! nested grid number
        integer :: model                            ! model identifier
!
        real(r8):: spval                            ! unbounded value
!
!  Source grid array declaration bounds, tile partition range,
!  longitude/latitude, curvilinear angle, and land/sea mask
!
        integer :: lbi_src, ubi_src, lbj_src, ubj_src
        integer :: istr_src, iend_src, jstr_src, jend_src
!
        real(r8) :: min_src, max_src
!
        real(r8), allocatable :: lon_src(:,:)
        real(r8), allocatable :: lat_src(:,:)
 
        real(r8), allocatable :: angle_src(:,:)
        real(r8), allocatable :: mask_src(:,:)
!
!  Destination grid array declaration bounds, tile partition range,
!  longitude/latitude, fractional coordinates, and land/sea mask.
!
        integer :: lbi_dst, ubi_dst, lbj_dst, ubj_dst
        integer :: istr_dst, iend_dst, jstr_dst, jend_dst
!
        real(r8) :: min_dst, max_dst
!
        real(r8), allocatable :: lon_dst(:,:)
        real(r8), allocatable :: lat_dst(:,:)
 
        real(r8), allocatable :: mask_dst(:,:)
 
        real(r8), allocatable :: x_dst(:,:)
        real(r8), allocatable :: y_dst(:,:)
 
      END TYPE roms_interp_type
!
!  Interpolation methods available
!
      integer, parameter :: bilinearmethod = 0    ! bilinear
      integer, parameter :: bicubicmethod = 1     ! bicubic
!
!-----------------------------------------------------------------------
!  Interface for same name routine overloading.
!-----------------------------------------------------------------------
!
     INTERFACE roms_datum_interp
       MODULE PROCEDURE roms_datum_interp_2d
       MODULE PROCEDURE roms_datum_interp_3d
       MODULE PROCEDURE roms_datum_column_interp
     END INTERFACE roms_datum_interp
!
     INTERFACE roms_horiz_interp
       MODULE PROCEDURE roms_horiz_interp_2d
       MODULE PROCEDURE roms_horiz_interp_3d
     END INTERFACE roms_horiz_interp
!
!-----------------------------------------------------------------------
!  Interpolation routines.
!-----------------------------------------------------------------------
!
      PUBLIC  :: linterp2d
      PUBLIC  :: cinterp2d
      PUBLIC  :: hindices
!
      PUBLIC  :: roms_interp_delete
      PUBLIC  :: roms_interp_fractional
!
      PRIVATE :: inside
      PRIVATE :: try_range
!
!-----------------------------------------------------------------------
      CONTAINS
!-----------------------------------------------------------------------
!
      SUBROUTINE linterp2d (ng, LBx, UBx, LBy, UBy,                     &
     &                      Xinp, Yinp, Finp,                           &
     &                      LBi, UBi, LBj, UBj,                         &
     &                      Istr, Iend, Jstr, Jend,                     &
     &                      Iout, Jout,                                 &
     &                      Xout, Yout,                                 &
     &                      Fout, MinVal, MaxVal)
!
!=======================================================================
!                                                                      !
!  Given any gridded 2D field, Finp, this routine linearly interpolate !
!  to locations (Xout,Yout).  To facilitate the  interpolation  within !
!  any irregularly gridded 2D field,  the fractional grid cell indices !
!  (Iout,Jout) with respect Finp are needed at input.  Notice that the !
!  routine "hindices" can be used to compute these indices.            !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Nested grid number.                                   !
!     LBx        I-dimension lower bound of gridded field, Finp.       !
!     UBx        I-dimension upper bound of gridded field, Finp.       !
!     LBy        J-dimension lower bound of gridded field, Finp.       !
!     UBy        J-dimension upper bound of gridded field, Finp.       !
!     Xinp       X-locations of gridded field, Finp.                   !
!     Yinp       Y-locations of gridded field, Finp.                   !
!     Finp       2D field to interpolate from.                         !
!     LBi        I-dimension Lower bound of data to interpolate, Fout. !
!     UBi        I-dimension Upper bound of data to interpolate, Fout. !
!     LBj        J-dimension Lower bound of data to interpolate, Fout. !
!     UBj        J-dimension Upper bound of data to interpolate, Fout. !
!     Istr       Starting data I-index to interpolate, Fout.           !
!     Iend       Ending   data I-index to interpolate, Fout.           !
!     Jstr       Starting data J-index to interpolate, Fout.           !
!     Jend       Ending   data J-index to interpolate, Fout.           !
!     Iout       I-fractional Xinp grid cell containing Xout.          !
!     Jout       J-fractional Yinp grid cell containing Yout.          !
!     Xout       X-locations to interpolate, Fout.                     !
!     Yout       Y-locations to interpolate, Fout.                     !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     Fout       Interpolated 2D field.                                !
!     MinVal     Interpolated field minimum value (Optional).          !
!     MaxVal     Interpolated field maximum value (Optional).          !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, LBx, UBx, LBy, UBy
      integer, intent(in) :: LBi, UBi, LBj, UBj
      integer, intent(in) :: Istr, Iend, Jstr, Jend
!
      real(r8), intent(in) :: Xinp(LBx:UBx,LBy:UBy)
      real(r8), intent(in) :: Yinp(LBx:UBx,LBy:UBy)
      real(r8), intent(in) :: Finp(LBx:UBx,LBy:UBy)
 
      real(r8), intent(in) :: Iout(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Jout(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Xout(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Yout(LBi:UBi,LBj:UBj)
 
      real(r8), intent(out) :: Fout(LBi:UBi,LBj:UBj)
!
      real(r8), intent(out), optional :: MinVal
      real(r8), intent(out), optional :: MaxVal
!
!  Local variable declarations.
!
      integer :: i, i1, i2, j, j1, j2
!
      real(r8) :: p1, p2, q1, q2
      real(r8) :: Fmin, Fmax
!
!-----------------------------------------------------------------------
!  Linearly interpolate requested field
!-----------------------------------------------------------------------
!
      fmin=1.0e+35_r8
      fmax=-1.0e+35_r8
      DO j=jstr,jend
        DO i=istr,iend
          i1=int(iout(i,j))
          i2=i1+1
          j1=int(jout(i,j))
          j2=j1+1
          IF (((lbx.le.i1).and.(i1.le.ubx)).and.                        &
     &        ((lby.le.j1).and.(j1.le.uby))) THEN
            p2=real(i2-i1,r8)*(iout(i,j)-real(i1,r8))
            q2=real(j2-j1,r8)*(jout(i,j)-real(j1,r8))
            p1=1.0_r8-p2
            q1=1.0_r8-q2
            fout(i,j)=p1*q1*finp(i1,j1)+                                &
     &                p2*q1*finp(i2,j1)+                                &
     &                p2*q2*finp(i2,j2)+                                &
     &                p1*q2*finp(i1,j2)
            fmin=min(fmin,fout(i,j))
            fmax=max(fmax,fout(i,j))
          END IF
        END DO
      END DO
!
!  Return minimum and maximum values.
!
      IF (PRESENT(minval)) THEN
        minval=fmin
      END IF
      IF (PRESENT(maxval)) THEN
        maxval=fmax
      END IF
!
      RETURN
      END SUBROUTINE linterp2d
!
      SUBROUTINE cinterp2d (ng, LBx, UBx, LBy, UBy,                     &
     &                      Xinp, Yinp, Finp,                           &
     &                      LBi, UBi, LBj, UBj,                         &
     &                      Istr, Iend, Jstr, Jend,                     &
     &                      Iout, Jout,                                 &
     &                      Xout, Yout,                                 &
     &                      Fout, MinVal, MaxVal)
!
!=======================================================================
!                                                                      !
!  Given any gridded 2D field,  Finp, at locations (Xinp,Yinp) this    !
!  routine performs bicubic interpolation at locations (Xout,Yout).    !
!  To facilitate the interpolation within any  irregularly  gridded    !
!  field, the fractional grid cell indices (Iout,Jout) with respect    !
!  Finp are needed at input. Notice that the routine "hindices" can    !
!  be used to compute these indices.                                   !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Nested grid number.                                   !
!     LBx        I-dimension lower bound of gridded field, Finp.       !
!     UBx        I-dimension upper bound of gridded field, Finp.       !
!     LBy        J-dimension lower bound of gridded field, Finp.       !
!     UBy        J-dimension upper bound of gridded field, Finp.       !
!     Xinp       X-locations of gridded field, Finp.                   !
!     Yinp       Y-locations of gridded field, Finp.                   !
!     Finp       2D field to interpolate from.                         !
!     LBi        I-dimension Lower bound of data to interpolate, Fout. !
!     UBi        I-dimension Upper bound of data to interpolate, Fout. !
!     LBj        J-dimension Lower bound of data to interpolate, Fout. !
!     UBj        J-dimension Upper bound of data to interpolate, Fout. !
!     Istr       Starting data I-index to interpolate, Fout.           !
!     Iend       Ending   data I-index to interpolate, Fout.           !
!     Jstr       Starting data J-index to interpolate, Fout.           !
!     Jend       Ending   data J-index to interpolate, Fout.           !
!     Iout       I-fractional Xinp grid cell containing Xout.          !
!     Jout       J-fractional Yinp grid cell containing Yout.          !
!     Xout       X-locations to interpolate, Fout.                     !
!     Yout       Y-locations to interpolate, Fout.                     !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     Fout       Interpolated 2D field.                                !
!     MinVal     Interpolated field minimum value (OPTIONAL).          !
!     MaxVal     Interpolated field maximum value (OPTIONAL).          !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, LBx, UBx, LBy, UBy
      integer, intent(in) :: LBi, UBi, LBj, UBj
      integer, intent(in) :: Istr, Iend, Jstr, Jend
!
      real(r8), intent(in) :: Xinp(LBx:UBx,LBy:UBy)
      real(r8), intent(in) :: Yinp(LBx:UBx,LBy:UBy)
      real(r8), intent(in) :: Finp(LBx:UBx,LBy:UBy)
 
      real(r8), intent(in) :: Iout(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Jout(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Xout(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Yout(LBi:UBi,LBj:UBj)
 
      real(r8), intent(out) :: Fout(LBi:UBi,LBj:UBj)
!
      real(r8), intent(out), optional :: MinVal
      real(r8), intent(out), optional :: MaxVal
!
!  Local variable declarations.
!
      integer i, ic, iter, i1, i2, j, jc, j1, j2
!
      real(r8) :: a11, a12, a21, a22
      real(r8) :: e11, e12, e21, e22
      real(r8) :: cff, d1, d2, dfc, dx, dy, eta, xi, xy, yx
      real(r8) :: f0, fx, fxx, fxxx, fxxy, fxy, fxyy, fy, fyy, fyyy
!
      real(r8), parameter :: C01 = 1.0_r8/48.0_r8
      real(r8), parameter :: C02 = 1.0_r8/32.0_r8
      real(r8), parameter :: C03 = 0.0625_r8                  ! 1/16
      real(r8), parameter :: C04 = 1.0_r8/6.0_r8
      real(r8), parameter :: C05 = 0.25_r8
      real(r8), parameter :: C06 = 0.5_r8
      real(r8), parameter :: C07 = 0.3125_r8                  ! 5/16
      real(r8), parameter :: C08 = 0.625_r8                   ! 5/8
      real(r8), parameter :: C09 = 1.5_r8
      real(r8), parameter :: C10 = 13.0_r8/24.0_r8
 
      real(r8), parameter :: LIMTR = 3.0_r8
      real(r8), parameter :: spv = 0.0_r8            ! HGA need work
 
      real(r8) :: Fmin, Fmax
!
      real(r8), dimension(-1:2,-1:2) :: dfx, dfy, f
!
!-----------------------------------------------------------------------
!  Interpolates requested field locations (Xout,Yout).
!-----------------------------------------------------------------------
!
      fmin=1.0e+35_r8
      fmax=-1.0e+35_r8
      DO j=jstr,jend
        DO i=istr,iend
          i1=int(iout(i,j))
          i2=i1+1
          j1=int(jout(i,j))
          j2=j1+1
          IF (((lbx.le.i1).and.(i1.le.ubx)).and.                        &
     &        ((lby.le.j1).and.(j1.le.uby))) THEN
!
!  Determine local fractional coordinates (xi,eta) corresponding to
!  the target point (Xout,Yout) on the grid (Xinp,Yinp). Here, "xi"
!  and "eta" are defined, in such a way, that xi=eta=0 corresponds
!  to the middle of the cell (i1:i1+1,j1:j1+1), while xi=+/-1/2 and
!  eta=+/-1/2 (any combination +/- signs) corresponds to the four
!  corner points of the cell. Inside the cell it is assumed that
!  (Xout,Yout) are expressed via bi-linear functions of (xi,eta),
!  where term proportional to xi*eta does not vanish because
!  coordinate transformation may be at least weakly non-orthogonal
!  due to discretization errors. The associated non-linear system
!  is solved by iterative method of Newton.
!
            xy=xinp(i2,j2)-xinp(i1,j2)-xinp(i2,j1)+xinp(i1,j1)
            yx=yinp(i2,j2)-yinp(i1,j2)-yinp(i2,j1)+yinp(i1,j1)
            dx=xout(i,j)-0.25_r8*(xinp(i2,j2)+xinp(i1,j2)+              &
     &                            xinp(i2,j1)+xinp(i1,j1))
            dy=yout(i,j)-0.25_r8*(yinp(i2,j2)+yinp(i1,j2)+              &
     &                            yinp(i2,j1)+yinp(i1,j1))
!
!  The coordinate transformation matrix:
!
!           e11 e12
!           e21 e22
!
!  contains derivatives of (Xinp,Yinp) with respect to (xi,eta). Because
!  the coordinates may be non-orthogonal (at least due to discretization
!  errors), the nonlinear system
!
!           e11*xi+e12*eta+xy*xi*eta=dx
!           e21*xi+e22*eta+yx*xi*eta=dy
!
!  needs to be solved in order to retain symmetry.
!
            e11=0.5_r8*(xinp(i2,j2)-xinp(i1,j2)+xinp(i2,j1)-xinp(i1,j1))
            e12=0.5_r8*(xinp(i2,j2)+xinp(i1,j2)-xinp(i2,j1)-xinp(i1,j1))
            e21=0.5_r8*(yinp(i2,j2)-yinp(i1,j2)+yinp(i2,j1)-yinp(i1,j1))
            e22=0.5_r8*(yinp(i2,j2)+yinp(i1,j2)-yinp(i2,j1)-yinp(i1,j1))
!
            cff=1.0_r8/(e11*e22-e12*e21)
            xi=cff*(e22*dx-e12*dy)
            eta=cff*(e11*dy-e21*dx)
!
            DO iter=1,4
              d1=dx-e11*xi-e12*eta-xy*xi*eta
              d2=dy-e21*xi-e22*eta-yx*xi*eta
              a11=e11+xy*eta
              a12=e12+xy*xi
              a21=e21+yx*eta
              a22=e22+yx*xi
              cff=1.0_r8/(a11*a22-a12*a21)
              xi =xi +cff*(a22*d1-a12*d2)
              eta=eta+cff*(a11*d2-a21*d1)
            END DO
 
#ifndef CUBIC_MASKED
!
!  Genuinely two-dimensional, isotropic cubic interpolation scheme
!  using 12-point stencil.  In the code below the interpolated field,
!  Fout, is expanded into two-dimensional Taylor series of local
!  fractional coordinates "xi" and "eta", retaining all terms of
!  combined power up to third order (that is, xi, eta, xi^2, eta^2,
!  xi*eta, xi^3, eta^3, xi^2*eta, and xi*eta^2), with all
!  coefficients (i.e, derivatives) computed via           x  x
!  two-dimensional finite difference expressions          |  |
!  of "natural" order of accuracy: 4th-order for       x--x--x--x
!  the field itself and its first derivatives in          |  |
!  both directions; and 2nd-order for all higher-      x--x--x--x
!  order derivatives. The permissible range of            |  |
!  of coordinates is -1/2 < xi,eta < +1/2, which          x--x
!  covers the central cell on the stencil, while
!  xi=eta=0 corresponds to its center. This interpolation scheme has
!  the property that if xi,eta=+/-1/2 (any combination of +/- signs)
!  it reproduces exactly value of the function at the corresponding
!  corner of the central "working" cell. However, it does not pass
!  exactly through the  extreme points of the stencil, where either
!  xi=+/-3/2 or eta+/-3/2. And, unlike a split-directional scheme,
!  when interpolating along the line eta=+/-1/2 (similarly xi=+/-1/2),
!  it has non-zero contribution from points on the side from the line,
!  except if xi=-1/2; 0; +1/2 (similarly eta=-1/2; 0; +1/2).
!
            DO jc=-1,2
              DO ic=-1,2
                f(ic,jc)=finp(max(1,min(ubx,i1+ic)),                    &
     &                        max(1,min(uby,j1+jc)))
              END DO
            END DO
 
            f0=c07*(f(1,1)+f(1,0)+f(0,1)+f(0,0))-                       &
     &         c02*(f(2,0)+f(2,1)+f(1,2)+f(0,2)+                        &
     &              f(-1,1)+f(-1,0)+f(0,-1)+f(1,-1))
 
            fx=c08*(f(1,1)+f(1,0)-f(0,1)-f(0,0))-                       &
     &         c01*(f(2,1)+f(2,0)-f(-1,1)-f(-1,0))-                     &
     &         c03*(f(1,2)-f(0,2)+f(1,-1)-f(0,-1))
 
            fy=c08*(f(1,1)-f(1,0)+f(0,1)-f(0,0))-                       &
     &         c01*(f(1,2)+f(0,2)-f(1,-1)-f(0,-1))-                     &
     &         c03*(f(2,1)-f(2,0)+f(-1,1)-f(-1,0))
 
            fxy=f(1,1)-f(1,0)-f(0,1)+f(0,0)
 
            fxx=c05*(f(2,1)-f(1,1)-f(0,1)+f(-1,1)+                      &
     &               f(2,0)-f(1,0)-f(0,0)+f(-1,0))
 
            fyy=c05*(f(1,2)-f(1,1)-f(1,0)+f(1,-1)+                      &
     &               f(0,2)-f(0,1)-f(0,0)+f(0,-1))
 
            fxxx=c06*(f(2,1)+f(2,0)-f(-1,1)-f(-1,0))-                   &
     &           c09*(f(1,1)+f(1,0)-f(0,1)-f(0,0))
 
            fyyy=c06*(f(1,2)+f(0,2)-f(1,-1)-f(0,-1))-                   &
     &           c09*(f(1,1)-f(1,0)+f(0,1)-f(0,0))
 
            fxxy=c06*(f(2,1)-f(1,1)-f(0,1)+f(-1,1)-                     &
     &                f(2,0)+f(1,0)+f(0,0)-f(-1,0))
 
            fxyy=c06*(f(1,2)-f(1,1)-f(1,0)+f(1,-1)-                     &
     &                f(0,2)+f(0,1)+f(0,0)-f(0,-1))
#else
!
!  Algorithm below is equivalent to the one above, except that special
!  care is taken to avoid interpolation accross land. This is achieved
!  by shortening the stencil and reducing order of polynomial, if
!  extreme points of the stencil touch land. This is achieved by
!  expressing all f0,fx,fy,...,fxyy in terms of values of interpolated
!  field at the four corners of central cell (which already checked to
!  stay away from land), and eight one-sided differences dfx,dfy (see
!  below) in such a way that field values at the extreme points of the
!  12-point stencil do not participate directly into f0,fx,...,fxyy.
!  Should an extreme point of the stencil touch land, thus making it
!  impossible to compute the corresponding one-sided difference, this
!  difference is retracted toward the center of the stencil.
!
!  Optionally, a slope-limiting algorithm may be employed to prevent
!  spurious oscillations of the interpolant. This is a valuable property,
!  if dealing with rough data, however, as a side effect, it turns off
!  high-order interpolation in the vicinity of extrema.
!
!  The slope-limiting algorithm employed here checks that two consecutive
!  elementary differences, "dfx" and "dfc" have the same sign and differ
!  in magnitude by no more than factor of 3.
!
            f(0,0)=finp(i1,j1)
            f(1,0)=finp(i2,j1)
            f(0,1)=finp(i1,j2)
            f(1,1)=finp(i2,j2)
!
            dfc=f(1,1)-f(0,1)
            IF (i1+2.gt.ubx) THEN
              dfx(1,1)=dfc
            ELSE IF (finp(i1+2,j2).eq.spv) THEN
              dfx(1,1)=dfc
            ELSE
              dfx(1,1)=finp(i1+2,j2)-f(1,1)
# ifdef LIMTR
              IF ((dfx(1,1)*dfc).lt.0.0_r8) THEN
                dfx(1,1)=0.0_r8
              ELSE IF (abs(dfx(1,1)).gt.(limtr*abs(dfc))) THEN
                dfx(1,1)=limtr*dfc
              END IF
# endif
            END IF
!
            dfc=f(1,0)-f(0,0)
            IF ((i1+2).gt.ubx) THEN
              dfx(1,0)=dfc
            ELSE IF (finp(i1+2,j1).eq.spv) THEN
              dfx(1,0)=dfc
            ELSE
              dfx(1,0)=finp(i1+2,j1)-f(1,0)
# ifdef LIMTR
              IF ((dfx(1,0)*dfc).lt.0.0_r8) THEN
                dfx(1,0)=0.0_r8
              ELSE IF (abs(dfx(1,0)).gt.(limtr*abs(dfc))) THEN
                dfx(1,0)=limtr*dfc
              END IF
# endif
            END IF
!
            dfc=f(1,1)-f(0,1)
            IF (i1-1.lt.1) THEN
              dfx(0,1)=dfc
            ELSE IF (finp(i1-1,j2).eq.spv) THEN
              dfx(0,1)=dfc
            ELSE
              dfx(0,1)=f(0,1)-finp(i1-1,j2)
# ifdef LIMTR
              IF ((dfx(0,1)*dfc).lt.0.0_r8) THEN
                dfx(0,1)=0.0_r8
              ELSE IF (abs(dfx(0,1)).gt.(limtr*abs(dfc))) THEN
                dfx(0,1)=limtr*dfc
              END IF
# endif
            END IF
!
            dfc=f(1,0)-f(0,0)
            IF (i1-1.lt.1) THEN
              dfx(0,0)=dfc
            ELSE IF (finp(i1-1,j1).eq.spv) THEN
              dfx(0,0)=dfc
            ELSE
              dfx(0,0)=f(0,0)-finp(i1-1,j1)
# ifdef LIMTR
              IF ((dfx(0,0)*dfc).lt.0.0_r8) THEN
                dfx(0,0)=0.0_r8
              ELSE IF (abs(dfx(0,0)).gt.(limtr*abs(dfc))) THEN
                dfx(0,0)=limtr*dfc
              END IF
# endif
            END IF
!
            dfc=f(1,1)-f(1,0)
            IF (j1+2.gt.uby) THEN
              dfy(1,1)=dfc
            ELSE IF (finp(i2,j1+2).eq.spv) THEN
              dfy(1,1)=dfc
            ELSE
              dfy(1,1)=finp(i2,j1+2)-f(1,1)
# ifdef LIMTR
              IF ((dfy(1,1)*dfc).lt.0.0_r8) THEN
                dfy(1,1)=0.0_r8
              ELSEIF (abs(dfy(1,1)).gt.(limtr*abs(dfc))) THEN
                dfy(1,1)=limtr*dfc
              END IF
# endif
            END IF
!
            dfc=f(0,1)-f(0,0)
            IF (j1+2.gt.uby) THEN
              dfy(0,1)=dfc
            ELSE IF (finp(i1,j1+2).eq.spv) THEN
              dfy(0,1)=dfc
            ELSE
              dfy(0,1)=finp(i1,j1+2)-f(0,1)
# ifdef LIMTR
              IF ((dfy(0,1)*dfc).lt.0.0_r8) THEN
                dfy(0,1)=0.0_r8
              ELSE IF (abs(dfy(0,1)).gt.(limtr*abs(dfc))) THEN
                dfy(0,1)=limtr*dfc
              END IF
# endif
            END IF
!
            dfc=f(1,1)-f(1,0)
            IF (j1-1.lt.1) THEN
              dfy(1,0)=dfc
            ELSE IF (finp(i2,j1-1).eq.spv) THEN
              dfy(1,0)=dfc
            ELSE
              dfy(1,0)=f(1,0)-finp(i2,j1-1)
# ifdef LIMTR
              IF ((dfy(1,0)*dfc).lt.0.0_r8) THEN
                dfy(1,0)=0.0_r8
              ELSE IF (abs(dfy(1,0)).gt.(limtr*abs(dfc))) THEN
                dfy(1,0)=limtr*dfc
              END IF
# endif
            END IF
!
            dfc=f(0,1)-f(0,0)
            IF (j1-1.lt.1) THEN
              dfy(0,0)=dfc
            ELSE IF (finp(i1,j1-1).eq.spv) THEN
              dfy(0,0)=dfc
            ELSE
              dfy(0,0)=f(0,0)-finp(i1,j1-1)
# ifdef LIMTR
              IF ((dfy(0,0)*dfc).lt.0.0_r8) THEN
                dfy(0,0)=0.0_r8
              ELSEIF (abs(dfy(0,0)).gt.(limtr*abs(dfc))) THEN
                dfy(0,0)=limtr*dfc
              END IF
# endif
            END IF
!
            f0=c05*(f(1,1)+f(1,0)+f(0,1)+f(0,0))-                       &
     &         c02*(dfx(1,1)+dfx(1,0)-dfx(0,1)-dfx(0,0)+                &
     &              dfy(1,1)-dfy(1,0)+dfy(0,1)-dfy(0,0))
 
            fx=c10*(f(1,1)-f(0,1)+f(1,0)-f(0,0))-                       &
     &         c01*(dfx(1,1)+dfx(1,0)+dfx(0,1)+dfx(0,0))-               &
     &         c03*(dfy(1,1)-dfy(0,1)-dfy(1,0)+dfy(0,0))
 
            fy=c10*(f(1,1)-f(1,0)+f(0,1)-f(0,0))-                       &
     &         c01*(dfy(1,1)+dfy(0,1)+dfy(1,0)+dfy(0,0))-               &
     &         c03*(dfx(1,1)-dfx(1,0)-dfx(0,1)+dfx(0,0))
 
            fxy=f(1,1)-f(1,0)-f(0,1)+f(0,0)
 
            fxx=c05*(dfx(1,1)-dfx(0,1)+dfx(1,0)-dfx(0,0))
 
            fyy=c05*(dfy(1,1)-dfy(1,0)+dfy(0,1)-dfy(0,0))
 
            fxxx=c06*(dfx(1,1)+dfx(1,0)+dfx(0,1)+dfx(0,0))-             &
     &           f(1,1)+f(0,1)-f(1,0)+f(0,0)
 
            fyyy=c06*(dfy(1,1)+dfy(0,1)+dfy(1,0)+dfy(0,0))-             &
     &           f(1,1)+f(1,0)-f(0,1)+f(0,0)
 
            fxxy=c06*(dfx(1,1)-dfx(0,1)-dfx(1,0)+dfx(0,0))
 
            fxyy=c06*(dfy(1,1)-dfy(1,0)-dfy(0,1)+dfy(0,0))
#endif
            fout(i,j)=f0+                                               &
     &                fx*xi+                                            &
     &                fy*eta+                                           &
     &                c06*fxx*xi*xi+                                    &
     &                fxy*xi*eta+                                       &
     &                c06*fyy*eta*eta+                                  &
     &                c04*fxxx*xi*xi*xi+                                &
     &                c06*fxxy*xi*xi*eta+                               &
     &                c04*fyyy*eta*eta*eta+                             &
     &                c06*fxyy*xi*eta*eta
            fmin=min(fmin,fout(i,j))
            fmax=max(fmax,fout(i,j))
          END IF
        END DO
      END DO
!
!  Return minimum and maximum values.
!
      IF (PRESENT(minval)) THEN
        minval=fmin
      END IF
      IF (PRESENT(maxval)) THEN
        maxval=fmax
      END IF
!
      RETURN
      END SUBROUTINE cinterp2d
!
      SUBROUTINE hindices (ng, LBi, UBi, LBj, UBj,                      &
     &                     Is, Ie, Js, Je,                              &
     &                     angler, Xgrd, Ygrd,                          &
     &                     LBm, UBm, LBn, UBn,                          &
     &                     Ms, Me, Ns, Ne,                              &
     &                     Xpos, Ypos, Ipos, Jpos,                      &
     &                     IJspv, rectangular)
!
!=======================================================================
!                                                                      !
!  Given any geographical locations Xpos and Ypos, this routine finds  !
!  the corresponding array cell indices (Ipos, Jpos) of gridded  data  !
!  Xgrd and Ygrd containing each requested location. This indices are  !
!  usually used for interpolation.                                     !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng          Nested grid number.                                  !
!     LBi         I-dimension Lower bound of gridded data.             !
!     UBi         I-dimension Upper bound of gridded data.             !
!     LBj         J-dimension Lower bound of gridded data.             !
!     UBj         J-dimension Upper bound of gridded data.             !
!     Is          Starting gridded data I-index to search.             !
!     Ie          Ending   gridded data I-index to search.             !
!     Js          Starting gridded data J-index to search.             !
!     Je          Ending   gridded data J-index to search.             !
!     angler      Gridded data angle between X-axis and true EAST      !
!                   (radians).                                         !
!     Xgrd        Gridded data X-locations (usually, longitude).       !
!     Ygrd        Gridded data Y-locations (usually, latitude).        !
!     LBm         I-dimension Lower bound of requested locations.      !
!     UBm         I-dimension Upper bound of requested locations.      !
!     LBn         J-dimension Lower bound of requested locations.      !
!     UBn         J-dimension Upper bound of requested locations.      !
!     Ms          Starting requested locations I-index to search.      !
!     Me          Ending   requested locations I-index to search.      !
!     Ns          Starting requested locations J-index to search.      !
!     Ne          Ending   requested locations J-index to search.      !
!     Xpos        Requested X-locations to process (usually longitude).!
!     Ypos        Requested Y-locations to process (usually latitude). !
!     IJspv       Unbounded special value to assign.                   !
!     rectangular Logical switch indicating that gridded data has a    !
!                   plaid distribution.                                !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     Ipos       Fractional I-cell index containing locations in data. !
!     Jpos       Fractional J-cell index containing locations in data. !
!                                                                      !
!  Calls:    Try_Range                                                 !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      logical,  intent(in) :: rectangular
!
      integer,  intent(in) :: ng
      integer,  intent(in) :: LBi, UBi, LBj, UBj, Is, Ie, Js, Je
      integer,  intent(in) :: LBm, UBm, LBn, UBn, Ms, Me, Ns, Ne
 
      real(r8), intent(in) :: IJspv
!
      real(r8), intent(in) :: angler(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Xgrd(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Ygrd(LBi:UBi,LBj:UBj)
 
      real(r8), intent(in) :: Xpos(LBm:UBm,LBn:UBn)
      real(r8), intent(in) :: Ypos(LBm:UBm,LBn:UBn)
 
      real(r8), intent(out) :: Ipos(LBm:UBm,LBn:UBn)
      real(r8), intent(out) :: Jpos(LBm:UBm,LBn:UBn)
!
!  Local variable declarations.
!
      logical :: found, foundi, foundj
 
      integer :: Imax, Imin, Jmax, Jmin, i, i0, j, j0, mp, np
 
      real(r8) :: aa2, ang, bb2, diag2, dx, dy, phi
      real(r8) :: xfac, xpp, yfac, ypp
!
!-----------------------------------------------------------------------
!  Determine grid cell indices containing requested position points.
!  Then, interpolate to fractional cell position.
!-----------------------------------------------------------------------
!
      DO np=ns,ne
        DO mp=ms,me
          ipos(mp,np)=ijspv
          jpos(mp,np)=ijspv
!
!  The gridded data has a plaid distribution so the search is trivial.
!
          IF (rectangular) THEN
            foundi=.false.
            i_loop : DO i=lbi,ubi-1
              IF ((xgrd(i  ,1).le.xpos(mp,np)).and.                     &
     &            (xgrd(i+1,1).gt.xpos(mp,np))) THEN
                imin=i
                foundi=.true.
                EXIT i_loop
              END IF
            END DO i_loop
            foundj=.false.
            j_loop : DO j=lbj,ubj-1
              IF ((ygrd(1,j  ).le.ypos(mp,np)).and.                     &
     &            (ygrd(1,j+1).gt.ypos(mp,np))) THEN
                jmin=j
                foundj=.true.
                EXIT j_loop
              END IF
            END DO j_loop
            found=foundi.and.foundj
!
!  Check each position to find if it falls inside the whole domain.
!  Once it is stablished that it inside, find the exact cell to which
!  it belongs by successively dividing the domain by a half (binary
!  search).
!
          ELSE
            found=try_range(ng, lbi, ubi, lbj, ubj,                     &
     &                      xgrd, ygrd,                                 &
     &                      is, ie, js, je,                             &
     &                      xpos(mp,np), ypos(mp,np))
            IF (found) THEN
              imin=is
              imax=ie
              jmin=js
              jmax=je
              DO while (((imax-imin).gt.1).or.((jmax-jmin).gt.1))
                IF ((imax-imin).gt.1) THEN
                  i0=(imin+imax)/2
                  found=try_range(ng, lbi, ubi, lbj, ubj,               &
     &                            xgrd, ygrd,                           &
     &                            imin, i0, jmin, jmax,                 &
     &                            xpos(mp,np), ypos(mp,np))
                  IF (found) THEN
                    imax=i0
                  ELSE
                    imin=i0
                  END IF
                END IF
                IF ((jmax-jmin).gt.1) THEN
                  j0=(jmin+jmax)/2
                  found=try_range(ng, lbi, ubi, lbj, ubj,               &
     &                            xgrd, ygrd,                           &
     &                            imin, imax, jmin, j0,                 &
     &                            xpos(mp,np), ypos(mp,np))
                  IF (found) THEN
                    jmax=j0
                  ELSE
                    jmin=j0
                  END IF
                END IF
              END DO
              found=(is.le.imin).and.(imin.le.ie).and.                  &
     &              (is.le.imax).and.(imax.le.ie).and.                  &
     &              (js.le.jmin).and.(jmin.le.je).and.                  &
     &              (js.le.jmax).and.(jmax.le.je)
            END IF
          END IF
!
!  Knowing the correct cell, calculate the exact indices, accounting
!  for a possibly rotated grid.  If spherical, convert all positions
!  to meters first.
!
          IF (found) THEN
            IF (spherical) THEN
              yfac=eradius*deg2rad
              xfac=yfac*cos(ypos(mp,np)*deg2rad)
              xpp=(xpos(mp,np)-xgrd(imin,jmin))*xfac
              ypp=(ypos(mp,np)-ygrd(imin,jmin))*yfac
            ELSE
              xfac=1.0_r8
              yfac=1.0_r8
              xpp=xpos(mp,np)-xgrd(imin,jmin)
              ypp=ypos(mp,np)-ygrd(imin,jmin)
            END IF
!
!  Use Law of Cosines to get cell parallelogram "shear" angle.
!
            diag2=((xgrd(imin+1,jmin)-xgrd(imin,jmin+1))*xfac)**2+      &
     &            ((ygrd(imin+1,jmin)-ygrd(imin,jmin+1))*yfac)**2
            aa2=((xgrd(imin,jmin)-xgrd(imin+1,jmin))*xfac)**2+          &
     &          ((ygrd(imin,jmin)-ygrd(imin+1,jmin))*yfac)**2
            bb2=((xgrd(imin,jmin)-xgrd(imin,jmin+1))*xfac)**2+          &
     &          ((ygrd(imin,jmin)-ygrd(imin,jmin+1))*yfac)**2
            phi=asin((diag2-aa2-bb2)/(2.0_r8*sqrt(aa2*bb2)))
!
!  Transform float position into curvilinear coordinates. Assume the
!  cell is rectanglar, for now.
!
            ang=angler(imin,jmin)
            dx=xpp*cos(ang)+ypp*sin(ang)
            dy=ypp*cos(ang)-xpp*sin(ang)
!
!  Correct for parallelogram.
!
            dx=dx+dy*tan(phi)
            dy=dy/cos(phi)
!
!  Scale with cell side lengths to translate into cell indices.
!
            dx=min(max(0.0_r8,dx/sqrt(aa2)),1.0_r8)
            dy=min(max(0.0_r8,dy/sqrt(bb2)),1.0_r8)
            ipos(mp,np)=real(imin,r8)+dx
            jpos(mp,np)=real(jmin,r8)+dy
          END IF
        END DO
      END DO
!
      RETURN
      END SUBROUTINE hindices
!
      LOGICAL FUNCTION try_range (ng, LBi, UBi, LBj, UBj, Xgrd, Ygrd,   &
     &                            Imin, Imax, Jmin, Jmax, Xo, Yo)
!
!=======================================================================
!                                                                      !
!  Given a grided domain with matrix coordinates Xgrd and Ygrd, this   !
!  function finds if the point (Xo,Yo)  is inside the box defined by   !
!  the requested corners (Imin,Jmin) and (Imax,Jmax). It will return   !
!  logical switch  try_range=.TRUE.  if (Xo,Yo) is inside, otherwise   !
!  it will return false.                                               !
!                                                                      !
!  Calls:   inside                                                     !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer,  intent(in) :: ng, lbi, ubi, lbj, ubj
      integer,  intent(in) :: imin, imax, jmin, jmax
 
#ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: xgrd(lbi:,lbj:)
      real(r8), intent(in) :: ygrd(lbi:,lbj:)
#else
      real(r8), intent(in) :: xgrd(lbi:ubi,lbj:ubj)
      real(r8), intent(in) :: ygrd(lbi:ubi,lbj:ubj)
#endif
 
      real(r8), intent(in) :: xo, yo
!
!  Local variable declarations.
!
      integer ::  nb, i, j, shft, ic
 
      real(r8), dimension(2*(Jmax-Jmin+Imax-Imin)+1) :: xb, yb
!
!-----------------------------------------------------------------------
!  Define closed polygon.
!-----------------------------------------------------------------------
!
!  Note that the last point (Xb(Nb),Yb(Nb)) does not repeat first
!  point (Xb(1),Yb(1)).  Instead, in function inside, it is implied
!  that the closing segment is (Xb(Nb),Yb(Nb))-->(Xb(1),Yb(1)). In
!  fact, function inside sets Xb(Nb+1)=Xb(1) and Yb(Nb+1)=Yb(1).
!
      nb=2*(jmax-jmin+imax-imin)
      shft=1-imin
      DO i=imin,imax-1
        xb(i+shft)=xgrd(i,jmin)
        yb(i+shft)=ygrd(i,jmin)
      END DO
      shft=1-jmin+imax-imin
      DO j=jmin,jmax-1
        xb(j+shft)=xgrd(imax,j)
        yb(j+shft)=ygrd(imax,j)
      END DO
      shft=1+jmax-jmin+2*imax-imin
      DO i=imax,imin+1,-1
        xb(shft-i)=xgrd(i,jmax)
        yb(shft-i)=ygrd(i,jmax)
      END DO
      shft=1+2*jmax-jmin+2*(imax-imin)
      DO j=jmax,jmin+1,-1
        xb(shft-j)=xgrd(imin,j)
        yb(shft-j)=ygrd(imin,j)
      END DO
!
!-----------------------------------------------------------------------
!  Check if point (Xo,Yo) is inside of the defined polygon.
!-----------------------------------------------------------------------
!
      try_range=inside(nb, xb, yb, xo, yo)
!
      RETURN
      END FUNCTION try_range
!
      LOGICAL FUNCTION inside (Nb, Xb, Yb, Xo, Yo)
!
!=======================================================================
!                                                                      !
!  Given the vectors Xb and Yb of size Nb, defining the coordinates    !
!  of a closed polygon,  this function find if the point (Xo,Yo) is    !
!  inside the polygon.  If the point  (Xo,Yo)  falls exactly on the    !
!  boundary of the polygon, it still considered inside.                !
!                                                                      !
!  This algorithm does not rely on the setting of  Xb(Nb)=Xb(1) and    !
!  Yb(Nb)=Yb(1).  Instead, it assumes that the last closing segment    !
!  is (Xb(Nb),Yb(Nb)) --> (Xb(1),Yb(1)).                               !
!                                                                      !
!  Reference:                                                          !
!                                                                      !
!    Reid, C., 1969: A long way from Euclid. Oceanography EMR,         !
!      page 174.                                                       !
!                                                                      !
!  Algorithm:                                                          !
!                                                                      !
!  The decision whether the point is  inside or outside the polygon    !
!  is done by counting the number of crossings from the ray (Xo,Yo)    !
!  to (Xo,-infinity), hereafter called meridian, by the boundary of    !
!  the polygon.  In this counting procedure,  a crossing is counted    !
!  as +2 if the crossing happens from "left to right" or -2 if from    !
!  "right to left". If the counting adds up to zero, then the point    !
!  is outside.  Otherwise,  it is either inside or on the boundary.    !
!                                                                      !
!  This routine is a modified version of the Reid (1969) algorithm,    !
!  where all crossings were counted as positive and the decision is    !
!  made  based on  whether the  number of crossings is even or odd.    !
!  This new algorithm may produce different results  in cases where    !
!  Xo accidentally coinsides with one of the (Xb(k),k=1:Nb) points.    !
!  In this case, the crossing is counted here as +1 or -1 depending    !
!  of the sign of (Xb(k+1)-Xb(k)).  Crossings  are  not  counted if    !
!  Xo=Xb(k)=Xb(k+1).  Therefore, if Xo=Xb(k0) and Yo>Yb(k0), and if    !
!  Xb(k0-1) < Xb(k0) < Xb(k0+1),  the crossing is counted twice but    !
!  with weight +1 (for segments with k=k0-1 and k=k0). Similarly if    !
!  Xb(k0-1) > Xb(k0) > Xb(k0+1), the crossing is counted twice with    !
!  weight -1 each time.  If,  on the other hand,  the meridian only    !
!  touches the boundary, that is, for example, Xb(k0-1) < Xb(k0)=Xo    !
!  and Xb(k0+1) < Xb(k0)=Xo, then the crossing is counted as +1 for    !
!  segment k=k0-1 and -1 for segment k=k0, resulting in no crossing.   !
!                                                                      !
!  Note 1: (Explanation of the logical condition)                      !
!                                                                      !
!  Suppose  that there exist two points  (x1,y1)=(Xb(k),Yb(k))  and    !
!  (x2,y2)=(Xb(k+1),Yb(k+1)),  such that,  either (x1 < Xo < x2) or    !
!  (x1 > Xo > x2).  Therefore, meridian x=Xo intersects the segment    !
!  (x1,y1) -> (x2,x2) and the ordinate of the point of intersection    !
!  is:                                                                 !
!                                                                      !
!                 y1*(x2-Xo) + y2*(Xo-x1)                              !
!             y = -----------------------                              !
!                          x2-x1                                       !
!                                                                      !
!  The mathematical statement that point  (Xo,Yo)  either coinsides    !
!  with the point of intersection or lies to the north (Yo>=y) from    !
!  it is, therefore, equivalent to the statement:                      !
!                                                                      !
!         Yo*(x2-x1) >= y1*(x2-Xo) + y2*(Xo-x1),   if   x2-x1 > 0      !
!  or                                                                  !
!         Yo*(x2-x1) <= y1*(x2-Xo) + y2*(Xo-x1),   if   x2-x1 < 0      !
!                                                                      !
!  which, after noting that  Yo*(x2-x1) = Yo*(x2-Xo + Xo-x1) may be    !
!  rewritten as:                                                       !
!                                                                      !
!        (Yo-y1)*(x2-Xo) + (Yo-y2)*(Xo-x1) >= 0,   if   x2-x1 > 0      !
!  or                                                                  !
!        (Yo-y1)*(x2-Xo) + (Yo-y2)*(Xo-x1) <= 0,   if   x2-x1 < 0      !
!                                                                      !
!  and both versions can be merged into  essentially  the condition    !
!  that (Yo-y1)*(x2-Xo)+(Yo-y2)*(Xo-x1) has the same sign as x2-x1.    !
!  That is, the product of these two must be positive or zero.         !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer,  intent(in   ) :: nb
 
      real(r8), intent(in   ) :: xo, yo
 
#ifdef ASSUMED_SHAPE
      real(r8), intent(inout) :: xb(:), yb(:)
#else
      real(r8), intent(inout) :: xb(nb+1), yb(nb+1)
#endif
!
!  Local variable declarations.
!
      logical :: itisinside
!
      integer, parameter :: nstep =128
 
      integer :: crossings, i, inc, k, kk, nc
 
      integer, dimension(Nstep) :: sindex
!
      real(r8) :: dx1, dx2, dxy
!
!-----------------------------------------------------------------------
!  Find intersections.
!-----------------------------------------------------------------------
!
!  Set crossings counter and close the contour of the polygon.
!
      itisinside=.false.
      crossings=0
      xb(nb+1)=xb(1)
      yb(nb+1)=yb(1)
!
!  The search is optimized.  First select the indices of segments
!  where Xb(k) is different from Xb(k+1) and Xo falls between them.
!  Then, further investigate these segments in a separate loop.
!  Doing it in two stages takes less time because the first loop is
!  pipelined.
!
      DO kk=0,nb-1,nstep
        nc=0
        DO k=kk+1,min(kk+nstep,nb)
          IF (((xb(k+1)-xo)*(xo-xb(k)).ge.0.0_r8).and.                  &
     &        (xb(k).ne.xb(k+1))) THEN
            nc=nc+1
            sindex(nc)=k
          END IF
        END DO
        DO i=1,nc
          k=sindex(i)
          IF (xb(k).ne.xb(k+1)) THEN
            dx1=xo-xb(k)
            dx2=xb(k+1)-xo
            dxy=dx2*(yo-yb(k))-dx1*(yb(k+1)-yo)
            inc=0
            IF ((xb(k).eq.xo).and.(yb(k).eq.yo)) THEN
              crossings=1
              itisinside=.true.               ! terminate inner DO-loop
              EXIT
            ELSE IF (((dx1.eq.0.0_r8).and.(yo.ge.yb(k  ))).or.          &
     &              ((dx2.eq.0.0_r8).and.(yo.ge.yb(k+1)))) THEN
              inc=1
            ELSE IF ((dx1*dx2.gt.0.0_r8).and.                           &
     &              ((xb(k+1)-xb(k))*dxy.ge.0.0_r8)) THEN  ! see note 1
              inc=2
            END IF
            IF (xb(k+1).gt.xb(k)) THEN
              crossings=crossings+inc
            ELSE
              crossings=crossings-inc
            END IF
          END IF
        END DO
        IF (itisinside) EXIT                  ! terminate outer DO-loop
      END DO
!
!  Determine if point (Xo,Yo) is inside of closed polygon.
!
      IF (crossings.eq.0) THEN
        inside=.false.
      ELSE
        inside=.true.
      END IF
!
      RETURN
      END FUNCTION inside
!
      SUBROUTINE roms_datum_interp_2d (S, fld_src, datum, method,       &
     &                                 ObsVetting)
!
!=======================================================================
!                                                                      !
!  Interpolates a ROMS 2D field to observation locations.              !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     S          Interpolation structure, TYPE(roms_interp_type),      !
!                  initialized elsewhere                               !
!     fld_src    2D source field data (real)                           !
!     method     Interpolation method flag (integer)                   !
!                  BilinearMethod => bilinear interpolation            !
!                  BicubicMethod  => bicubic  interpolation            !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     datum      Interpolated values at datum locations (vector).      !
!     ObsVetting Observation screenning flag, 0: reject or 1: accept   !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer,  intent(in   ) :: method
!
      real(r8), intent(in   ) :: fld_src(:,:)
      real(r8), intent(out  ) :: datum(:)
 
      real(r8), intent(inout), optional :: ObsVetting(:)
!
      TYPE (roms_interp_type), intent(inout) :: S
!
!  Local variable declarations.
!
      integer :: ic, iobs, i1, i2, j1, j2, nlocs
!
      real(r8) :: p1, p2, q1, q2, wsum
 
      real(r8), dimension(4) :: Hmat
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", roms_datum_interp_2d"
!
!-----------------------------------------------------------------------
!  Interpolates a 2D ROMS state field to observation locations.
!-----------------------------------------------------------------------
!
      nlocs=SIZE(datum, dim=1)
!
      DO iobs=1,nlocs
        i1=int(s%x_dst(iobs,1))
        j1=int(s%y_dst(iobs,1))
        IF (((s%Istr_src.le.i1).and.(i1.lt.s%Iend_src)).and.            &
     &      ((s%Jstr_src.le.j1).and.(j1.lt.s%Jend_src))) THEN
          i2=i1+1
          j2=j1+1
          p2=real(i2-i1,r8)*(s%x_dst(iobs,1)-real(i1,r8))
          q2=real(j2-j1,r8)*(s%y_dst(iobs,1)-real(j1,r8))
          p1=1.0_r8-p2
          q1=1.0_r8-q2
          hmat(1)=p1*q1
          hmat(2)=p2*q1
          hmat(3)=p2*q2
          hmat(4)=p1*q2
#ifdef MASKING
          hmat(1)=hmat(1)*s%mask_src(i1,j1)
          hmat(2)=hmat(2)*s%mask_src(i2,j1)
          hmat(3)=hmat(3)*s%mask_src(i2,j2)
          hmat(4)=hmat(4)*s%mask_src(i1,j2)
          wsum=0.0_r8
          DO ic=1,4
            wsum=wsum+hmat(ic)
          END DO
          IF (wsum.gt.0.0_r8) THEN
            wsum=1.0_r8/wsum
            DO ic=1,4
              hmat(ic)=hmat(ic)*wsum
            END DO
          END IF
#endif
          datum(iobs)=hmat(1)*fld_src(i1,j1)+                          &
     &                hmat(2)*fld_src(i2,j1)+                          &
     &                hmat(3)*fld_src(i2,j2)+                          &
     &                hmat(4)*fld_src(i1,j2)
!
!  Set observation screening flag, 0: reject or 1: accept.
!
          IF (PRESENT(obsvetting)) THEN
#ifdef MASKING
            IF (wsum.gt.0.0_r8) obsvetting(iobs)=1.0_r8
#else
            obsvetting(iobs)=1.0_r8
#endif
          END IF
        END IF
      END DO
!
      RETURN
      END SUBROUTINE roms_datum_interp_2d
!
      SUBROUTINE roms_datum_interp_3d (S, fld_src, z_src, z_locs,      &
                                       datum, method, obsVetting)
!
!=======================================================================
!                                                                      !
!  Interpolates a 3D ROMS field to observation horizontal and vertical !
!  locations.                                                          !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     S          Interpolation structure, TYPE(roms_interp_type),      !
!                  initialized elsewhere                               !
!     fld_src    3D source field data (real)                           !
!     z_src      3D source field depth (real; meter; negative)         !
!     z_locs     Datum depth locations (real vector). Negative values  !
!                  are depths in meters while positive values are in   !
!                  vertical fractional coordinates of 1 to N.          !
!     method     Interpolation method flag (integer)                   !
!                  BilinearMethod => bilinear interpolation            !
!                  BicubicMethod  => bicubic  interpolation            !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     datum      Interpolated values at observation locations (vector) !
!     ObsVetting Observation screenning flag, 0: reject or 1: accept   !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer,  intent( in) :: method
!
      real(r8), intent( in) :: fld_src(:,:,:)
      real(r8), intent( in) :: z_src(:,:,:)
      real(r8), intent( in) :: z_locs(:)
      real(r8), intent(out) :: datum(:)
 
      real(r8), intent(inout), optional :: ObsVetting(:)
!
      TYPE (roms_interp_type), intent(inout) :: S
!
!  Local variable declarations.
!
      integer :: ic, iobs, i1, i2, j1, j2, k, k1, k2, nlevs, nlocs
!
      real(r8) :: Zbot, Ztop, dz, p1, p2, q1, q2, r1, r2
      real(r8) :: w11, w12, w21, w22, wsum
 
      real(r8), dimension(8) :: Hmat
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", roms_datum_interp_3d"
!
!-----------------------------------------------------------------------
!  Interpolate from 3D ROMS state field at observation locations.
!-----------------------------------------------------------------------
!
      nlevs=SIZE(fld_src, dim=3)
      nlocs=SIZE(datum, dim=1)
!
      DO iobs=1,nlocs
        i1=int(s%x_dst(iobs,1))
        j1=int(s%y_dst(iobs,1))
        IF (((s%Istr_src.le.i1).and.(i1.lt.s%Iend_src)).and.            &
     &      ((s%Jstr_src.le.j1).and.(j1.lt.s%Jend_src))) THEN
          i2=i1+1
          j2=j1+1
          p2=real(i2-i1,r8)*(s%x_dst(iobs,1)-real(i1,r8))
          q2=real(j2-j1,r8)*(s%y_dst(iobs,1)-real(j1,r8))
          p1=1.0_r8-p2
          q1=1.0_r8-q2
          w11=p1*q1
          w21=p2*q1
          w22=p2*q2
          w12=p1*q2
          IF (z_locs(iobs).gt.0.0_r8) THEN
            k1=max(1,int(z_locs(iobs)))        ! Positions in fractional
            k2=min(int(z_locs(iobs))+1,nlevs)  ! levels
            r2=real(k2-k1,r8)*(z_locs(iobs)-real(k1,r8))
            r1=1.0_r8-r2
          ELSE
            ztop=z_src(i1,j1,nlevs)
            zbot=z_src(i1,j1,1    )
            IF (z_locs(iobs).ge.ztop) THEN
              r1=0.0_r8                        ! If shallower, ignore.
              r2=0.0_r8
              IF (PRESENT(obsvetting)) obsvetting(iobs)=0.0_r8
            ELSE IF (zbot.ge.z_locs(iobs)) THEN
              r1=0.0_r8                        ! If deeper, ignore.
              r2=0.0_r8
              IF (PRESENT(obsvetting)) obsvetting(iobs)=0.0_r8
            ELSE
              DO k=nlevs,2,-1                  ! Otherwise, interpolate
                ztop=z_src(i1,j1,k  )          ! to fractional level
                zbot=z_src(i1,j1,k-1)
                IF ((ztop.gt.z_locs(iobs)).and.                         &
                    (z_locs(iobs).ge.zbot)) THEN
                  k1=k-1
                  k2=k
                END IF
              END DO
              dz=z_src(i1,j1,k2)-z_src(i1,j1,k1)
              r2=(z_locs(iobs)-z_src(i1,j1,k1))/dz
              r1=1.0_r8-r2
            END IF
          END IF
          IF ((r1+r2).gt.0.0_r8) THEN
            hmat(1)=w11*r1
            hmat(2)=w21*r1
            hmat(3)=w22*r1
            hmat(4)=w12*r1
            hmat(5)=w11*r2
            hmat(6)=w21*r2
            hmat(7)=w22*r2
            hmat(8)=w12*r2
#ifdef MASKING
            hmat(1)=hmat(1)*s%mask_src(i1,j1)
            hmat(2)=hmat(2)*s%mask_src(i2,j1)
            hmat(3)=hmat(3)*s%mask_src(i2,j2)
            hmat(4)=hmat(4)*s%mask_src(i1,j2)
            hmat(5)=hmat(5)*s%mask_src(i1,j1)
            hmat(6)=hmat(6)*s%mask_src(i2,j1)
            hmat(7)=hmat(7)*s%mask_src(i2,j2)
            hmat(8)=hmat(8)*s%mask_src(i1,j2)
            wsum=0.0_r8
            DO ic=1,8
              wsum=wsum+hmat(ic)
            END DO
            IF (wsum.gt.0.0_r8) THEN
              wsum=1.0_r8/wsum
              DO ic=1,8
                hmat(ic)=hmat(ic)*wsum
              END DO
            END IF
#endif
            datum(iobs)=hmat(1)*fld_src(i1,j1,k1)+                      &
     &                  hmat(2)*fld_src(i2,j1,k1)+                      &
     &                  hmat(3)*fld_src(i2,j2,k1)+                      &
     &                  hmat(4)*fld_src(i1,j2,k1)+                      &
     &                  hmat(5)*fld_src(i1,j1,k2)+                      &
     &                  hmat(6)*fld_src(i2,j1,k2)+                      &
     &                  hmat(7)*fld_src(i2,j2,k2)+                      &
     &                  hmat(8)*fld_src(i1,j2,k2)
!
!  Set observation screening flag, 0: reject or 1: accept.
!
            IF (PRESENT(obsvetting)) THEN
#ifdef MASKING
              IF (wsum.gt.0.0_r8) obsvetting(iobs)=1.0_r8
#else
              obsvetting(iobs)=1.0_r8
#endif
#ifndef ALLOW_BOTTOM_OBS
!
!  Reject observations that lie in the lower bottom grid cell (k=1) to
!  avoid clustering due shallowing of bathymetry during smoothing and
!  coarse level half-thickness (-h < z_locs < z_src(:,:,1)) in deep
!  water.
!
              IF ((z_locs(iobs).gt.0.0_r8).and.                         &
                  (z_locs(iobs).le.1.0_r8)) THEN
                obsvetting(iobs)=0.0_r8
              END IF
#endif
            END IF
          END IF
        END IF
      END DO
!
      RETURN
      END SUBROUTINE roms_datum_interp_3d
!
      SUBROUTINE roms_datum_column_interp (S, fld_src, datum, method)
!
!=======================================================================
!                                                                      !
!  Horizontally interpolates a 3D ROMS field to datum locations        !
!  level-by-level. The vertical interpolation will be carried out      !
!  elsewhere from the full column of values.  This is what is          !
!  required by the JEDI-UFO operator.                                  !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     S          Interpolation structure, TYPE(roms_interp_type),      !
!                  initialized elsewhere.                              !
!     fld_src    3D field source data (real).                          !
!     method     Interpolation method flag (integer)                   !
!                  BilinearMethod => bilinear interpolation            !
!                  BicubicMethod  => bicubic  interpolation            !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     datum      Interpolated values at datum locations (matrix):      !
!                  datum(nlocs,nlevs)                                  !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      TYPE (roms_interp_type), intent(inout) :: S
!
      real(r8), intent( in) :: fld_src(:,:,:)
      real(r8), intent(out) :: datum(:,:)
!
      integer,  intent( in) :: method
!
!  Local variable declarations.
!
      integer :: ic, iobs, i1, i2, j1, j2, k, nlevs, nlocs
!
      real(r8) :: p1, p2, q1, q2, wsum
 
      real(r8), dimension(4) :: Hmat
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", roms_datum_interp_2d"
!
!-----------------------------------------------------------------------
!  Interpolate from 2D state field to datum locations.
!-----------------------------------------------------------------------
!
      nlevs=SIZE(fld_src, dim=3)
      nlocs=SIZE(datum, dim=1)
!
      DO iobs=1,nlocs
        i1=int(s%x_dst(iobs,1))
        j1=int(s%y_dst(iobs,1))
        IF (((s%Istr_src.le.i1).and.(i1.lt.s%Iend_src)).and.            &
     &      ((s%Jstr_src.le.j1).and.(j1.lt.s%Jend_src))) THEN
          i2=i1+1
          j2=j1+1
          p2=real(i2-i1,r8)*(s%x_dst(iobs,1)-real(i1,r8))
          q2=real(j2-j1,r8)*(s%y_dst(iobs,1)-real(j1,r8))
          p1=1.0_r8-p2
          q1=1.0_r8-q2
          hmat(1)=p1*q1
          hmat(2)=p2*q1
          hmat(3)=p2*q2
          hmat(4)=p1*q2
#ifdef MASKING
          hmat(1)=hmat(1)*s%mask_src(i1,j1)
          hmat(2)=hmat(2)*s%mask_src(i2,j1)
          hmat(3)=hmat(3)*s%mask_src(i2,j2)
          hmat(4)=hmat(4)*s%mask_src(i1,j2)
          wsum=0.0_r8
          DO ic=1,4
            wsum=wsum+hmat(ic)
          END DO
          IF (wsum.gt.0.0_r8) THEN
            wsum=1.0_r8/wsum
            DO ic=1,4
              hmat(ic)=hmat(ic)*wsum
            END DO
          END IF
#endif
          DO k=1,nlevs
            datum(iobs,k)=hmat(1)*fld_src(i1,j1,k)+                    &
     &                    hmat(2)*fld_src(i2,j1,k)+                    &
     &                    hmat(3)*fld_src(i2,j2,k)+                    &
     &                    hmat(4)*fld_src(i1,j2,k)
          END DO
        END IF
      END DO
!
      RETURN
      END SUBROUTINE roms_datum_column_interp
!
      SUBROUTINE roms_horiz_interp_2d (S, fld_src, fld_dst, method)
!
!=======================================================================
!                                                                      !
!  Interpolate a ROMS 2D source field to a new destination locations.  !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     S          Interpolation structure, TYPE(roms_interp_type),      !
!                  initialized elsewhere.                              !
!     fld_src    2D field source data (real).                          !
!     method     Interpolation method flag (integer)                   !
!                  BilinearMethod => bilinear interpolation            !
!                  BicubicMethod  => bicubic  interpolation            !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     fld_dst    2D interpolated field data (real).                    !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      TYPE (roms_interp_type), intent(inout) :: S
!
      real(r8), intent(in ) :: fld_src(:,:)
      real(r8), intent(out) :: fld_dst(:,:)
!
      integer,  intent(in ) :: method
!
!  Local variable declarations.
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", roms_horiz_interp_2d"
!
!-----------------------------------------------------------------------
!  Interpolate 2D field.
!-----------------------------------------------------------------------
!
      SELECT CASE (method)
        CASE (bilinearmethod)
          CALL linterp2d (s%ng,                                         &
     &                    s%LBi_src, s%UBi_src, s%LBj_src, s%UBj_src,   &
     &                    s%lon_src, s%lat_src, fld_src,                &
     &                    s%LBi_dst, s%UBi_dst, s%LBj_dst, s%UBj_dst,   &
     &                    s%Istr_dst, s%Iend_dst,                       &
     &                    s%Jstr_dst, s%Jend_dst,                       &
     &                    s%x_dst, s%y_dst, s%lon_dst, s%lat_dst,       &
     &                    fld_dst,                                      &
     &                    minval = s%min_dst,                           &
     &                    maxval = s%max_dst)
        CASE (bicubicmethod)
          CALL cinterp2d (s%ng,                                         &
     &                    s%LBi_src, s%UBi_src, s%LBj_src, s%UBj_src,   &
     &                    s%lon_src, s%lat_src, fld_src,                &
     &                    s%LBi_dst, s%UBi_dst, s%LBj_dst, s%UBj_dst,   &
     &                    s%Istr_dst, s%Iend_dst,                       &
     &                    s%Jstr_dst, s%Jend_dst,                       &
     &                    s%x_dst, s%y_dst, s%lon_dst, s%lat_dst,       &
     &                    fld_dst,                                      &
     &                    minval = s%min_dst,                           &
     &                    maxval = s%max_dst)
        CASE DEFAULT
          IF (master) WRITE(stdout,10) method
          exit_flag=5
      END SELECT
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
  10  FORMAT (' ROMS_HORIZ_INTERP_2D - Illegal interpolation method =', &
     &        i0)
!
      RETURN
      END SUBROUTINE roms_horiz_interp_2d
!
      SUBROUTINE roms_horiz_interp_3d (S, fld_src, fld_dst, method)
!
!=======================================================================
!                                                                      !
!  Interpolate a ROMS 3D source field to a new destination locations   !
!  level-by-level. There is no vertical interpolation.                 !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     S          Interpolation structure, TYPE(roms_interp_type),      !
!                  initialized elsewhere.                              !
!     fld_src    3D field source data (real).                          !
!     method     Interpolation method flag (integer)                   !
!                  BilinearMethod => bilinear interpolation            !
!                  BicubicMethod  => bicubic  interpolation            !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     fld_dst    3D interpolated field data (real).                    !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      TYPE (roms_interp_type), intent(inout) :: S
!
      real(r8), intent( in) :: fld_src(:,:,:)
      real(r8), intent(out) :: fld_dst(:,:,:)
!
      integer,  intent( in) :: method
!
!  Local variable declarations.
!
      integer :: k, nk
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", roms_horiz_interp_3d"
!
!-----------------------------------------------------------------------
!  Interpolate 3D field level-by-level.
!-----------------------------------------------------------------------
!
      nk=SIZE(fld_src, dim=3)!
!
      SELECT CASE (method)
        CASE (bilinearmethod)
          DO k=1,nk
            CALL linterp2d (s%ng,                                       &
     &                      s%LBi_src, s%UBi_src, s%LBj_src, s%UBj_src, &
     &                      s%lon_src, s%lat_src, fld_src(:,:,k),       &
     &                      s%LBi_dst, s%UBi_dst, s%LBj_dst, s%UBj_dst, &
     &                      s%Istr_dst, s%Iend_dst,                     &
     &                      s%Jstr_dst, s%Jend_dst,                     &
     &                      s%x_dst, s%y_dst, s%lon_dst, s%lat_dst,     &
     &                      fld_dst(:,:,k),                             &
     &                      minval = s%min_dst,                         &
     &                      maxval = s%max_dst)
          END DO
        CASE (bicubicmethod)
          DO k=1,nk
            CALL cinterp2d (s%ng,                                       &
     &                      s%LBi_src, s%UBi_src, s%LBj_src, s%UBj_src, &
     &                      s%lon_src, s%lat_src, fld_src(:,:,k),       &
     &                      s%LBi_dst, s%UBi_dst, s%LBj_dst, s%UBj_dst, &
     &                      s%Istr_dst, s%Iend_dst,                     &
     &                      s%Jstr_dst, s%Jend_dst,                     &
     &                      s%x_dst, s%y_dst, s%lon_dst, s%lat_dst,     &
     &                      fld_dst(:,:,k),                             &
     &                      minval = s%min_dst,                         &
     &                      maxval = s%max_dst)
          END DO
        CASE DEFAULT
          IF (master) WRITE(stdout,10) method
          exit_flag=5
      END SELECT
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
  10  FORMAT (' ROMS_HORIZ_INTERP_3D - Illegal interpolation method =', &
     &        i0)
!
      RETURN
      END SUBROUTINE roms_horiz_interp_3d
!
      SUBROUTINE roms_interp_create (ng, tile, gtype, S)
!
!=======================================================================
!                                                                      !
!  Creates ROMS intepolation Object.  The interpolation structure has  !
!  two components: source and destination. Here, the source component  !
!  is initialized. The destination component is initialized elsewhere  !
!  for flexibility.                                                    !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, gtype
!
      TYPE (roms_interp_type), intent(inout) :: S
!
!  Local variable declarations.
!
      integer :: LBi, UBi, LBj, UBj
!
!-----------------------------------------------------------------------
!  Creates source component of interpolation object.
!-----------------------------------------------------------------------
!
!  If applicable, destroy interpolation object.
!
      IF (allocated(s%lon_src)) THEN
        CALL roms_interp_delete (s)
      END IF
!
!  Allocate source component arrays.
!
      lbi=bounds(ng)%LBi(tile)
      ubi=bounds(ng)%UBi(tile)
      lbj=bounds(ng)%LBj(tile)
      ubj=bounds(ng)%UBj(tile)
!
      allocate ( s%lon_src(lbi:ubi,lbj:ubj) )
      allocate ( s%lat_src(lbi:ubi,lbj:ubj) )
      allocate ( s%angle_src(lbi:ubi,lbj:ubj) )
#ifdef MASKING
      allocate ( s%mask_src(lbi:ubi,lbj:ubj) )
#endif
!
!  Associate and assign values to source component.
!
      s%ng = ng
      s%model = inlm
      s%spval = 0.0_r8
!
      SELECT CASE (gtype)
        CASE (r2dvar, r3dvar)
          s%lon_src   = grid(ng)%lonr
          s%lat_src   = grid(ng)%latr
          s%angle_src = grid(ng)%angler
#ifdef MASKING
          s%mask_src  = grid(ng)%rmask
#endif
          s%Istr_src  = bounds(ng)%IstrR(tile)
          s%Iend_src  = bounds(ng)%IendR(tile)
          s%Jstr_src  = bounds(ng)%JstrR(tile)
          s%Jend_src  = bounds(ng)%JendR(tile)
        CASE (u2dvar, u3dvar)
          s%lon_src   = grid(ng)%lonu
          s%lat_src   = grid(ng)%latu
          s%angle_src = grid(ng)%angler
#ifdef MASKING
          s%mask_src  = grid(ng)%umask
#endif
          s%Istr_src  = bounds(ng)%Istr (tile)
          s%Iend_src  = bounds(ng)%IendR(tile)
          s%Jstr_src  = bounds(ng)%JstrR(tile)
          s%Jend_src  = bounds(ng)%JendR(tile)
        CASE (v2dvar, v3dvar)
          s%lon_src   = grid(ng)%lonv
          s%lat_src   = grid(ng)%latv
          s%angle_src = grid(ng)%angler
#ifdef MASKING
          s%mask_src  = grid(ng)%vmask
#endif
          s%Istr_src  = bounds(ng)%IstrR(tile)
          s%Iend_src  = bounds(ng)%IendR(tile)
          s%Jstr_src  = bounds(ng)%Jstr (tile)
          s%Jend_src  = bounds(ng)%JendR(tile)
      END SELECT
!
!
      s%LBi_src = lbi
      s%UBi_src = ubi
      s%LBj_src = lbj
      s%UBj_src = ubj
!
      RETURN
      END SUBROUTINE roms_interp_create
!
      SUBROUTINE roms_interp_delete (S)
!
!=======================================================================
!                                                                      !
!  Deallocates pointer arrays in the intepolation structure, S.        !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      TYPE (roms_interp_type), intent(inout) :: S
!
!-----------------------------------------------------------------------
!  Deallocate interpolation structure.
!-----------------------------------------------------------------------
!
      IF (allocated(s%angle_src)) deallocate (s%angle_src)
      IF (allocated(s%lon_src))   deallocate (s%lon_src)
      IF (allocated(s%lat_src))   deallocate (s%lat_src)
      IF (allocated(s%mask_src))  deallocate (s%mask_src)
 
      IF (allocated(s%lon_dst))   deallocate (s%lon_dst)
      IF (allocated(s%lat_dst))   deallocate (s%lat_dst)
      IF (allocated(s%mask_dst))  deallocate (s%mask_dst)
      IF (allocated(s%x_dst))     deallocate (s%x_dst)
      IF (allocated(s%y_dst))     deallocate (s%y_dst)
!
      RETURN
      END SUBROUTINE roms_interp_delete
!
      SUBROUTINE roms_interp_fractional (S)
!
!=======================================================================
!                                                                      !
!  Computes the horizontal fractional coordinates (S%x_dst, S%y_dst)   !
!  of the source cells containing the destination values.              !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      TYPE (roms_interp_type), intent(inout) :: S
!
!  Local variable declarations.
!
      integer :: Npts
!
!-----------------------------------------------------------------------
!  Determine the horizontal fractional coordinates: S%x_dst and S%y_dst.
!-----------------------------------------------------------------------
!
      CALL hindices (s%ng,                                              &
     &               s%LBi_src,   s%UBi_src,  s%LBj_src,  s%UBj_src,    &
     &               s%Istr_src , s%Iend_src, s%Jstr_src, s%Jend_src,   &
     &               s%angle_src, s%lon_src,  s%lat_src,                &
     &               s%LBi_dst,   s%UBi_dst,  s%LBj_dst,  s%UBj_dst,    &
     &               s%Istr_dst,  s%Iend_dst, s%Jstr_dst, s%Jend_dst,   &
     &               s%lon_dst,   s%lat_dst,  s%x_dst,    s%y_dst,      &
     &               s%spval,     s%rectangular)
 
#ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Collect fractional coordinates.
!-----------------------------------------------------------------------
!
      npts=SIZE(s%x_dst)
      CALL mp_assemble (s%ng, s%model, npts, s%spval, s%x_dst)
      CALL mp_assemble (s%ng, s%model, npts, s%spval, s%y_dst)
#endif
!
      RETURN
      END SUBROUTINE roms_interp_fractional
!
      END MODULE roms_interpolate_mod