nesting.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE nesting_mod
 
#ifdef NESTING
!
!git $Id$
!=======================================================================
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license           Hernan G. Arango   !
!    See License_ROMS.md                              John C. Warner   !
!=======================================================================
!                                                                      !
!  This module contains several routines  to process the connectivity  !
!  between nested grids. It process the contact region points between  !
!  data donor and data receiver grids.                                 !
!                                                                      !
!  The locations of the linear interpolation weights in the donor      !
!  grid with respect the receiver grid contact region at contact       !
!  point x(Irg,Jrg,Krg) are:                                           !
!                                                                      !
!                       8___________7   (Idg+1,Jdg+1,Kdg)              !
!                      /.          /|                                  !
!                     / .         / |                                  !
!  (Idg,Jdg+1,Kdg)  5/___________/6 |                                  !
!                    |  .        |  |                                  !
!                    |  .   x    |  |                                  !
!                    | 4.........|..|3  (Idg+1,Jdg+1,Kdg-1)            !
!                    | .         |  /                                  !
!                    |.          | /                                   !
!                    |___________|/                                    !
!  (Idg,Jdg,Kdg-1)   1           2                                     !
!                                                                      !
!                                        Suffix:   dg = donor grid     !
!                                                  rg = receiver grid  !
!                                                                      !
!  Routines:                                                           !
!  ========                                                            !
!                                                                      !
!  nesting          Public interface to time-stepping kernel           !
!                                                                      !
!  get_composite    Composite  grid, extract contact points donor data !
!  get_metrics      Extract grid spacing metrics, on_u and om_v.       !
!  get_refine       Refinement grid, extract contact points donor data !
!  put_composite    Composite  grid, fill contact points (interpolate) !
!  put_refine       Refinement grid, fill contact points (interpolate) !
!                                                                      !
!  bry_fluxes       Extracts horizontat advective fluxes the contact   !
!                     boundary of donor and receiver grids.            !
# ifdef NESTING_DEBUG
!  check_massflux   If refinement, check mass fluxes between coarse    !
!                     and fine grids for volume conservation. It is    !
!                     use only for debugging and diagnostics.          !
# endif
!  correct_tracer   Correct coarse grid tracer at the refinement grid  !
!                     boundary with the refined accumulated fluxes     !
!  do_twoway        Logical function to determine at which time-step   !
!                     the two-way exchange takes place.                !
!  fill_contact     Used to Fill grid metrics at contact points.       !
!  fine2coarse      Replace coarse grid state variables with the       !
!                     averaged fine grid values (two-way nesting)      !
!                                                                      !
!  get_contact2d    Get 2D field donor grid cell holding contact point !
!  get_contact3d    Get 3D field donor grid cell holding contact point !
!  get_persisted2d  Get 2D field persisted values on contact points    !
!  mask_hweights    Scale horizontal interpolation weights with masking!
!  put_contact2d    Set 2D field contact points, spatial interpolation !
!  put_contact3d    Set 3D field contact points, spatial interpolation !
!                                                                      !
!  put_refine2d     Interpolate (space-time) 2D state variables        !
!  put_refine3d     Interpolate (space-time) 3D state variables        !
!                                                                      !
!  z_weights        Set donor grid vertical indices (cell holding      !
!                     contact point) and vertical interpolation        !
!                     weights                                          !
!                                                                      !
!  WARNINGS:                                                           !
!  ========                                                            !
!                                                                      !
!  All the routines contained in this module are inside of a parallel  !
!  region, except the main driver routine "nesting",  which is called  !
!  serially several times from main2d or main3d to perform  different  !
!  tasks.  Notice that the calls to private "get_***"  and  "put_***"  !
!  routines need to be in separated  parallel loops because of serial  !
!  with partitions and  shared-memory rules.  Furthermore,  the donor  !
!  and receiver grids may have different tile partitions. There is no  !
!  I/O management inside the nesting routines.                         !
!                                                                      !
!  The connectivity between donor and receiver grids can be  complex.  !
!  The horizontal mapping between grids is static and done outside of  !
!  ROMS.  Only the time-dependent  vertical interpolation weights are  !
!  computed here.  The contact region points  I- and  J-cell  indices  !
!  between donor and receiver grids, and the horizontal interpolation  !
!  weights are read from  the input nesting connectivity NetCDF file.  !
!  It makes the nesting efficient and greatly simplifies parallelism.  !
!                                                                      !
!=======================================================================
!
      implicit none
!
      PUBLIC  :: nesting
      PUBLIC  :: bry_fluxes
# ifndef ONE_WAY
      PUBLIC  :: do_twoway
# endif
      PUBLIC  :: fill_contact
# ifdef NESTING_DEBUG
      PUBLIC  :: check_massflux
# endif
# ifdef SOLVE3D
      PRIVATE :: correct_tracer
      PRIVATE :: correct_tracer_tile
# endif
      PRIVATE :: fine2coarse
      PUBLIC  :: fine2coarse2d
# ifdef SOLVE3D
      PUBLIC  :: fine2coarse3d
# endif
      PUBLIC  :: get_contact2d
# ifdef SOLVE3D
      PUBLIC  :: get_contact3d
# endif
      PRIVATE :: get_composite
      PUBLIC  :: get_metrics
      PUBLIC  :: get_persisted2d
      PRIVATE :: get_refine
# if defined MASKING || defined WET_DRY
      PUBLIC  :: mask_hweights
# endif
      PRIVATE :: put_composite
      PRIVATE :: put_refine
      PRIVATE :: put_refine2d
# ifdef SOLVE3D
      PRIVATE :: put_refine3d
      PUBLIC  :: z_weights
# endif
!
      CONTAINS
!
      SUBROUTINE nesting (ng, model, isection)
!
!=======================================================================
!                                                                      !
!  This routine process the contact region points between composite    !
!  grids.  In composite grids, it is possible to have more than one    !
!  contact region.                                                     !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Data receiver grid number (integer)                   !
!     model      Calling model identifier (integer)                    !
!     isection   Governing equations time-stepping section in          !
!                  main2d or main3d indicating which state             !
!                  variables to process (integer)                      !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_ncparam
      USE mod_nesting
      USE mod_scalars
!
# ifdef SOLVE3D
      USE set_depth_mod, ONLY : set_depth
# endif
      USE strings_mod,   ONLY : founderror
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model, isection
!
!  Local variable declarations.
!
      logical :: lputfsur
!
      integer :: subs, tile, thread
      integer :: ngc
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__
 
# ifdef PROFILE
!
!-----------------------------------------------------------------------
!  Turn on time clocks.
!-----------------------------------------------------------------------
!
      CALL wclock_on (ng, model, 36, __line__, myfile)
# endif
# ifdef SOLVE3D
!
!-----------------------------------------------------------------------
!  Process vertical indices and interpolation weigths associated with
!  depth. Currently, vertical weights are used in composite grids
!  because their grids are not coincident. They are not needed in
!  refinement grids because the donor and receiver grids have the same
!  number of vertical levels and matching bathymetry.  However, in
!  the future, it is possible to have configurations that require
!  vertical weights in refinement. The switch "get_Vweights" controls
!  if such weights are computed or not.  If false, it will accelerate
!  computations because of less distributed-memory communications.
!-----------------------------------------------------------------------
!
      IF ((isection.eq.nzwgt).and.get_vweights) THEN
        DO tile=last_tile(ng),first_tile(ng),-1
          CALL z_weights (ng, model, tile)
        END DO
!$OMP BARRIER
        RETURN
      END IF
# endif
 
# if defined MASKING || defined WET_DRY
!
!-----------------------------------------------------------------------
!  If Land/Sea masking, scale horizontal interpolation weights to
!  account for land contact points. If wetting and drying, the scaling
!  is done at every time-step because masking is time dependent.
!-----------------------------------------------------------------------
!
      IF (isection.eq.nmask) THEN
        DO tile=last_tile(ng),first_tile(ng),-1
          CALL mask_hweights (ng, model, tile)
        END DO
!$OMP BARRIER
        RETURN
      END IF
# endif
!
!-----------------------------------------------------------------------
!  If refinement grid, process contact points.
!-----------------------------------------------------------------------
!
      IF (refinedgrid(ng)) THEN
!
!  Extract grid spacing metrics (on_u and om_v) and load then to
!  REFINE(:) structure.  These metrics are needed to impose mass
!  flux at the finer grid physical boundaries.  It need to be done
!  separately because parallelism partions between all nested grid.
!
        IF (isection.eq.ndxdy) THEN
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL get_metrics (ng, model, tile)
          END DO
!$OMP BARRIER
!
!  Extract and store donor grid data at contact points.
!
        ELSE IF (isection.eq.ngetd) THEN
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL get_refine (ng, model, tile)
          END DO
!$OMP BARRIER
!
!  Fill refinement grid contact points variables by interpolating
!  (space, time) from extracted donor grid data.  The free-surface
!  needs to be processed first and in a separate parallel region
!  because of shared-memory applications.
!
        ELSE IF (isection.eq.nputd) THEN
          lputfsur=.true.
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL put_refine (ng, model, tile, lputfsur)
          END DO
!$OMP BARRIER
!
          lputfsur=.false.
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL put_refine (ng, model, tile, lputfsur)
          END DO
!$OMP BARRIER
 
# ifdef NESTING_DEBUG
!
!  If refinement, check mass flux conservation between coarser and
!  finer grids.
!
        ELSE IF (isection.eq.nmflx) THEN
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL check_massflux (ng, model, tile)
          END DO
# endif
 
# ifndef ONE_WAY
!
!  Fine to coarse coupling: two-way nesting.
!
        ELSE IF (isection.eq.n2way) THEN
 
          ngc=coarserdonor(ng)               ! coarse grid number
 
#  if defined SOLVE3D && !defined NO_CORRECT_TRACER
!
!  Correct coarse grid tracer values at the refinement grid, ng,
!  boundary with the refined accumulated fluxes (Hz*u*T/n, Hz*v*T/m).
!
          DO tile=first_tile(ngc),last_tile(ngc),+1
            CALL correct_tracer (ngc, ng, model, tile)
          END DO
!$OMP BARRIER
#  endif
!
!  Replace coarse grid 2D state variables with the averaged fine grid
!  values (two-way coupling).
!
          DO tile=last_tile(ngc),first_tile(ngc),-1
            CALL fine2coarse (ng, model, r2dvar, tile)
          END DO
!$OMP BARRIER
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#  ifdef SOLVE3D
!
!  Update coarse grid depth variables. We have a new coarse grid
!  adjusted free-surface, Zt_avg1.
!
          DO tile=first_tile(ngc),last_tile(ngc),+1
            CALL set_depth (ngc, tile, model)
          END DO
!$OMP BARRIER
!
!  Replace coarse grid 3D state variables with the averaged fine grid
!  values (two-way coupling).
!
          DO tile=last_tile(ngc),first_tile(ngc),-1
            CALL fine2coarse (ng, model, r3dvar, tile)
          END DO
!$OMP BARRIER
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
# else
!
!  Fine to coarse coupling (two-way nesting) is not activated!
!
        ELSE IF (isection.eq.n2way) THEN
# endif
        END IF
!
!-----------------------------------------------------------------------
!  Otherwise, process contact points in composite grid.
!-----------------------------------------------------------------------
!
      ELSE
!
!  Get composite grid contact points data from donor grid. It extracts
!  the donor grid cell data necessary to interpolate state variables
!  at each contact point.
!
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL get_composite (ng, model, isection, tile)
        END DO
!$OMP BARRIER
!
!  Fill composite grid contact points variables by interpolating from
!  extracted donor grid data.
!
        DO tile=last_tile(ng),first_tile(ng),-1
          CALL put_composite (ng, model, isection, tile)
        END DO
!$OMP BARRIER
 
      END IF
 
# ifdef PROFILE
!
!-----------------------------------------------------------------------
!  Turn off time clocks.
!-----------------------------------------------------------------------
!
      CALL wclock_off (ng, model, 36, __line__, myfile)
# endif
!
      RETURN
      END SUBROUTINE nesting
!
      SUBROUTINE bry_fluxes (dg, rg, cr, model, tile,                   &
     &                       IminS, ImaxS, JminS, JmaxS,                &
     &                       ILB, IUB, JLB, JUB,                        &
     &                       scale, FX, FE,                             &
     &                       F_west, F_east, F_south, F_north)
!
!=======================================================================
!                                                                      !
!  This routine extracts tracer horizontal advective fluxes (Hz*u*T/n, !
!  Hz*v*T/m) at the grid contact boundary (physical domain perimeter). !
!  The data source is either the coarse or finer grid.  These fluxes   !
!  are used for in two-way nesting.         b                          !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     dg         Donor grid number (integer)                           !
!     rg         Receiver grid number (integer)                        !
!     cr         Contact region number to process (integer)            !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!     scale      Advective flux scale (floating-point)                 !
!     IminS      Advective flux, I-dimension Lower bound (integer)     !
!     ImaxS      Advective flux, I-dimension Upper bound (integer)     !
!     JminS      Advective flux, J-dimension Lower bound (integer)     !
!     JmaxS      Advective flux, J-dimension Upper bound (integer)     !
!     ILB        Western/Eastern   boundary flux Lower bound (integer) !
!     IUB        Western/Eastern   boundary flux Upper bound (integer) !
!     JLB        Southern/Northern boundary flux Lower bound (integer) !
!     JUB        Southern/Northern boundary flux Lower bound (integer) !
!     FX         Horizontal advetive flux in the XI-direction (array)  !
!     FE         Horizontal advetive flux in the ETA-direction (array) !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     F_west     Western  boundary advective flux (1D array)           !
!     F_east     Eastern  boundary advective flux (1D array)           !
!     F_south    Southern boundary advective flux (1D array)           !
!     F_north    Northerb boundary advective flux (1D array)           !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_nesting
      USE mod_scalars
!
# ifdef DISTRIBUTE
 
      USE distribute_mod, ONLY : mp_assemble
# endif
      USE strings_mod,    ONLY : founderror
!
!  Imported variable declarations.
!
      integer, intent(in) :: dg, rg, cr, model, tile
      integer, intent(in) :: imins, imaxs, jmins, jmaxs
      integer, intent(in) :: ilb, iub, jlb, jub
 
      real(dp), intent(in) :: scale
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: fx(imins:,jmins:)
      real(r8), intent(in) :: fe(imins:,jmins:)
 
      real(r8), intent(inout) :: f_west (jlb:)
      real(r8), intent(inout) :: f_east (jlb:)
      real(r8), intent(inout) :: f_south(ilb:)
      real(r8), intent(inout) :: f_north(ilb:)
# else
      real(r8), intent(in) :: fx(imins:imaxs,jmins:jmaxs)
      real(r8), intent(in) :: fe(imins:imaxs,jmins:jmaxs)
 
      real(r8), intent(inout) :: f_west (jlb:jub)
      real(r8), intent(inout) :: f_east (jlb:jub)
      real(r8), intent(inout) :: f_south(ilb:iub)
      real(r8), intent(inout) :: f_north(ilb:iub)
# endif
!
!  Local variable declarations.
!
      integer :: istr, iend, jstr, jend
      integer :: ib_east, ib_west, jb_north, jb_south
      integer :: i, j, m
 
# ifdef DISTRIBUTE
      integer :: nptswe, nptssn
# endif
!
      real(r8), parameter :: fspv = 0.0_r8
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", bry_fluxes"
!
!-----------------------------------------------------------------------
!  Initialize local variables.
!-----------------------------------------------------------------------
!
!  Set tile starting and ending indices.
!
      istr=bounds(rg)%Istr(tile)
      iend=bounds(rg)%Iend(tile)
      jstr=bounds(rg)%Jstr(tile)
      jend=bounds(rg)%Jend(tile)
 
# ifdef DISTRIBUTE
!
!  Initialize arrays to facilitate collective communications.
!
      nptswe=jub-jlb+1
      nptssn=iub-ilb+1
# endif
!
!-----------------------------------------------------------------------
!  If "rg" is the finer grid, extract advective tracer flux at its
!  physical domain boundaries (grid perimeter).
!-----------------------------------------------------------------------
!
!  Receiver finer grid number is greater than donor coaser grid number
!  because of refinement nesting layers.
!
      IF (rg.gt.dg) THEN
!
! Reset fluxes to zero for first entry of receiver finer grid.
!
        IF (mod(iic(rg)-1,refinesteps(rg)).eq.0) THEN
          f_west =fspv
          f_east =fspv
          f_south=fspv
          f_north=fspv
        END IF
!
!  Western boundary.
!
        IF (domain(dg)%Western_Edge(tile)) THEN
          DO j=jstr,jend
            f_west(j)=f_west(j)+fx(istr,j)*scale
          END DO
          DO j=jlb,jstr-1
            f_west(j)=fspv
          END DO
          DO j=jend+1,jub
            f_west(j)=fspv
          END DO
        ELSE
          f_west=fspv
        END IF
!
!  Eastern boundary.
!
        IF (domain(dg)%Eastern_Edge(tile)) THEN
          DO j=jstr,jend
            f_east(j)=f_east(j)+fx(iend+1,j)*scale
          END DO
          DO j=jlb,jstr-1
            f_east(j)=fspv
          END DO
          DO j=jend+1,jub
            f_east(j)=fspv
          END DO
        ELSE
          f_east=fspv
        END IF
!
!  Southern boundary.
!
        IF (domain(dg)%Southern_Edge(tile)) THEN
          DO i=istr,iend
            f_south(i)=f_south(i)+fe(i,jstr)*scale
          END DO
          DO i=ilb,istr-1
            f_south(i)=fspv
          END DO
          DO i=iend+1,iub
            f_south(i)=fspv
          END DO
        ELSE
          f_south=fspv
        END IF
!
!  Northern boundary.
!
        IF (domain(dg)%Northern_Edge(tile)) THEN
          DO i=istr,iend
            f_north(i)=f_north(i)+fe(i,jend+1)*scale
          END DO
          DO i=ilb,istr-1
            f_north(i)=fspv
          END DO
          DO i=iend+1,iub
            f_north(i)=fspv
          END DO
        ELSE
          f_north=fspv
        END IF
!
!-----------------------------------------------------------------------
!  If "rg" is the coarser grid, extract coarser grid advective tracer
!  flux at the location of the finer grid physical domain boundaries
!  (grid perimeter).
!-----------------------------------------------------------------------
!
!  Receiver coarser grid number is smaller than donor finer grid number
!  because of refinement nesting layers.
!
      ELSE IF (rg.lt.dg) THEN
!
!  Western/Eastern boundaries.
!
        f_west =fspv
        f_east =fspv
        f_south=fspv
        f_north=fspv
!
        ib_west=i_left(dg)
        ib_east=i_right(dg)
        DO j=jstr,jend
          IF ((istr.le.ib_west).and.(ib_west.le.iend)) THEN
            f_west(j)=fx(ib_west,j)*scale
          END IF
!
          IF ((istr.le.ib_east).and.(ib_east.le.iend)) THEN
            f_east(j)=fx(ib_east,j)*scale
          END IF
        END DO
!
!  Southern/Northern boundaries.
!
        jb_south=j_bottom(dg)
        jb_north=j_top(dg)
        DO i=istr,iend
          IF ((jstr.le.jb_south).and.(jb_south.le.jend)) THEN
            f_south(i)=fe(i,jb_south)*scale
          END IF
!
          IF ((jstr.le.jb_north).and.(jb_north.le.jend)) THEN
            f_north(i)=fe(i,jb_north)*scale
          END IF
        END DO
      END IF
 
# ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Gather and broadcast data from all nodes.
!-----------------------------------------------------------------------
!
      CALL mp_assemble (dg, model, nptswe, fspv, f_west(jlb:))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      CALL mp_assemble (dg, model, nptswe, fspv, f_east(jlb:))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      CALL mp_assemble (dg, model, nptssn, fspv, f_south(ilb:))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      CALL mp_assemble (dg, model, nptssn, fspv, f_north(ilb:))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
!
      RETURN
      END SUBROUTINE bry_fluxes
 
# ifdef NESTING_DEBUG
!
      SUBROUTINE check_massflux (ngf, model, tile)
!
!=======================================================================
!                                                                      !
!  If refinement, this routine check mass fluxes between coarse and    !
!  fine grids for mass and volume conservation. It is only used for    !
!  diagnostic purposes.                                                !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ngf          Finer grid number (integer)                         !
!     model        Calling model identifier (integer)                  !
!     tile         Domain tile partition (integer)                     !
!                                                                      !
!  On Output:    (mod_nesting)                                         !
!                                                                      !
!     BRY_CONTACT  Updated Mflux in structure.                         !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_nesting
      USE mod_scalars
 
#  ifdef DISTRIBUTE
!
      USE distribute_mod, ONLY : mp_assemble
#  endif
!
!  Imported variable declarations.
!
      integer, intent(in) :: ngf, model, tile
!
!  Local variable declarations.
!
#  ifdef DISTRIBUTE
      integer :: ilb, iub, jlb, jub, nptssn, nptswe, my_tile
#  endif
      integer :: iedge, ibc, ibc_min, ibc_max, ibf, io
      integer :: jedge, jbc, jbc_min, jbc_max, jbf, jo
      integer :: istr, iend, jstr, jend
      integer :: cjcr, cr, dg, half, icr, isum, jsum, m, rg
      integer :: tnew, told
 
#  ifdef DISTRIBUTE
      real(r8), parameter :: spv = 0.0_r8
#  endif
      real(r8) :: eastsum, northsum, southsum, westsum
      real(r8) :: mfratio
!
!-----------------------------------------------------------------------
!  Check mass and volume conservation during refinement between coarse
!  and fine grids.
!-----------------------------------------------------------------------
!
      DO cr=1,ncontact
!
!  Get data donor and data receiver grid numbers.
!
        dg=rcontact(cr)%donor_grid
        rg=rcontact(cr)%receiver_grid
!
!  Process only contact region data for requested nested finer grid
!  "ngf". Notice that the donor grid is coarser than receiver grid.
!
        IF ((rg.eq.ngf).and.(dxmax(dg).gt.dxmax(rg))) THEN
!
!  Set tile starting and ending indices for donor coarser grid.
!
          istr=bounds(dg)%Istr(tile)
          iend=bounds(dg)%Iend(tile)
          jstr=bounds(dg)%Jstr(tile)
          jend=bounds(dg)%Jend(tile)
!
!  Set time rolling indices and conjugate region where the coarser
!  donor grid becomes the receiver grid.
!
          told=3-rollingindex(cr)
          tnew=rollingindex(cr)
          DO icr=1,ncontact
            IF ((rg.eq.rcontact(icr)%donor_grid).and.                   &
     &          (dg.eq.rcontact(icr)%receiver_grid)) THEN
              cjcr=icr
              EXIT
            END IF
          END DO
 
#  ifdef DISTRIBUTE
!
!  Set global size of boundary edges for coarse grid (donor index).
!
          my_tile=-1
          ilb=bounds(dg)%LBi(my_tile)
          iub=bounds(dg)%UBi(my_tile)
          jlb=bounds(dg)%LBj(my_tile)
          jub=bounds(dg)%UBj(my_tile)
          nptswe=jub-jlb+1
          nptssn=iub-ilb+1
!
!  If distributed-memory, initialize arrays used to check mass flux
!  conservation with special value (zero) to facilitate the global
!  reduction when collecting data between all nodes.
!
          bry_contact(iwest ,cjcr)%Mflux=spv
          bry_contact(ieast ,cjcr)%Mflux=spv
          bry_contact(isouth,cjcr)%Mflux=spv
          bry_contact(inorth,cjcr)%Mflux=spv
#  endif
!
!  Set finer grid center (half) and offset indices (Io and Jo) for
!  coarser grid (I,J) coordinates.
!
          half=(refinescale(ngf)-1)/2
          io=half+1
          jo=half+1
!
!-----------------------------------------------------------------------
!  Average finer grid western boundary mass fluxes and load them to the
!  BRY_CONTACT structure.
!-----------------------------------------------------------------------
!
          ibc=i_left(ngf)
          jbc_min=j_bottom(ngf)
          jbc_max=j_top(ngf)-1            ! interior points, no top
!                                           left corner
          IF (domain(ngf)%SouthWest_Test(tile)) THEN
            IF (master) THEN
              WRITE (300,10) 'Western Boundary Mass Fluxes: ',          &
     &                       cr, dg, rg, iif(rg), iic(rg), int(time(rg))
              FLUSH (300)
            END IF
          END IF
!
          DO jbc=jstr,jend
            IF (((istr.le.ibc).and.(ibc.le.iend)).and.                  &
     &          ((jbc_min.le.jbc).and.(jbc.le.jbc_max))) THEN
!
!  Sum finer grid western boundary mass fluxes within coarser grid cell.
!
              westsum=0.0_r8
              jedge=jo+(jbc-jbc_min)*refinescale(ngf)
              DO jsum=-half,half
                jbf=jedge+jsum
                westsum=westsum+bry_contact(iwest,cr)%Mflux(jbf)
              END DO
              m=bry_contact(iwest,cr)%C2Bindex(jbf)     ! pick last one
!
!  Load coarser grid western boundary mass flux that have been averaged
!  from finer grid. These values can be compared with the coarser grid
!  values REFINED(cr)%U2d_flux to check if the mass flux between coarser
!  and finer grid is conserved.
!
              bry_contact(iwest,cjcr)%Mflux(jbc)=westsum
              IF (westsum.ne.0) THEN
                mfratio=refined(cr)%U2d_flux(1,m,tnew)/westsum
              ELSE
                mfratio=1.0_r8
              END IF
              WRITE (300,30) jbc, refined(cr)%U2d_flux(1,m,tnew),        &
     &                       westsum, mfratio
              FLUSH (300)
            END IF
          END DO
!
!-----------------------------------------------------------------------
!  Average finer grid eastern boundary mass fluxes and load them to the
!  BRY_CONTACT structure.
!-----------------------------------------------------------------------
!
          ibc=i_right(ngf)
          jbc_min=j_bottom(ngf)
          jbc_max=j_top(ngf)-1            ! interior points, no top
!                                           right corner
          IF (domain(ngf)%SouthWest_Test(tile)) THEN
            IF (master) THEN
              WRITE (300,10) 'Eastern Boundary Mass Fluxes: ',          &
     &                       cr, dg, rg, iif(rg), iic(rg), int(time(rg))
              FLUSH (300)
            END IF
          END IF
!
          DO jbc=jstr,jend
            IF (((istr.le.ibc).and.(ibc.le.iend)).and.                  &
     &          ((jbc_min.le.jbc).and.(jbc.le.jbc_max))) THEN
!
!  Sum finer grid eastern boundary mass fluxes within coarser grid cell.
!
              eastsum=0.0_r8
              jedge=jo+(jbc-jbc_min)*refinescale(ngf)
              DO jsum=-half,half
                jbf=jedge+jsum
                eastsum=eastsum+bry_contact(ieast,cr)%Mflux(jbf)
              END DO
              m=bry_contact(ieast,cr)%C2Bindex(jbf)     ! pick last one
!
!  Load coarser grid eastern boundary mass flux that have been averaged
!  from finer grid. These values can be compared with the coarser grid
!  values REFINED(cr)%U2d_flux to check if the mass flux between coarser
!  and finer grid is conserved.
!
              bry_contact(ieast,cjcr)%Mflux(jbc)=eastsum
              IF (eastsum.ne.0) THEN
                mfratio=refined(cr)%U2d_flux(1,m,tnew)/eastsum
              ELSE
                mfratio=1.0_r8
              END IF
              WRITE (300,30) jbc, refined(cr)%U2d_flux(1,m,tnew),        &
     &                       eastsum, mfratio
              FLUSH (300)
            END IF
          END DO
!
!-----------------------------------------------------------------------
!  Average finer grid southern boundary mass fluxes and load them to the
!  BRY_CONTACT structure.
!-----------------------------------------------------------------------
!
          jbc=j_bottom(ngf)
          ibc_min=i_left(ngf)
          ibc_max=i_right(ngf)-1          ! interior points, no bottom
!                                           right corner
          IF (domain(ngf)%SouthWest_Test(tile)) THEN
            IF (master) THEN
              WRITE (300,20) 'Southern Boundary Mass Fluxes: ',         &
     &                       cr, dg, rg, iif(rg), iic(rg), int(time(rg))
              FLUSH (300)
            END IF
          END IF
!
          DO ibc=istr,iend
            IF (((ibc_min.le.ibc).and.(ibc.le.ibc_max)).and.            &
     &          ((jstr.le.jbc).and.(jbc.le.jend))) THEN
!
!  Sum finer grid southern boundary mass fluxes within coarser grid
!  cell.
!
              southsum=0.0_r8
              iedge=io+(ibc-ibc_min)*refinescale(ngf)
              DO isum=-half,half
                ibf=iedge+isum
                southsum=southsum+bry_contact(isouth,cr)%Mflux(ibf)
              END DO
              m=bry_contact(isouth,cr)%C2Bindex(ibf)    ! pick last one
!
!  Load coarser grid southern boundary mass flux that have been averaged
!  from finer grid. These values can be compared with the coarser grid
!  values REFINED(cr)%V2d_flux to check if the mass flux between coarser
!  and finer grid is conserved.
!
              bry_contact(isouth,cjcr)%Mflux(ibc)=southsum
              IF (southsum.ne.0) THEN
                mfratio=refined(cr)%V2d_flux(1,m,tnew)/southsum
              ELSE
                mfratio=1.0_r8
              END IF
              WRITE (300,30) ibc, refined(cr)%V2d_flux(1,m,tnew),        &
     &                       southsum, mfratio
              FLUSH (300)
            END IF
          END DO
!
!-----------------------------------------------------------------------
!  Average finer grid northern boundary mass fluxes and load them to the
!  BRY_CONTACT structure.
!-----------------------------------------------------------------------
!
          jbc=j_top(ngf)
          ibc_min=i_left(ngf)
          ibc_max=i_right(ngf)-1          ! interior points, no top
!                                           right corner
          IF (domain(ngf)%SouthWest_Test(tile)) THEN
            IF (master) THEN
              WRITE (300,20) 'Northern Boundary Mass Fluxes: ',         &
     &                       cr, dg, rg, iif(rg), iic(rg), int(time(rg))
              FLUSH (300)
            END IF
          END IF
!
          DO ibc=istr,iend
            IF (((ibc_min.le.ibc).and.(ibc.le.ibc_max)).and.            &
     &          ((jstr.le.jbc).and.(jbc.le.jend))) THEN
!
!  Sum finer grid northern boundary mass fluxes within coarser grid
!  cell.
!
              northsum=0.0_r8
              iedge=io+(ibc-ibc_min)*refinescale(ngf)
              DO isum=-half,half
                ibf=iedge+isum
                northsum=northsum+bry_contact(inorth,cr)%Mflux(ibf)
              END DO
              m=bry_contact(inorth,cr)%C2Bindex(ibf)    ! pick last one
!
!  Load coarser grid northern boundary mass flux that have been averaged
!  from finer grid. These values can be compared with the coarser grid
!  values REFINED(cr)%V2d_flux to check if the mass flux between coarser
!  and finer grid is conserved.
!
              bry_contact(inorth,cjcr)%Mflux(ibc)=northsum
              IF (northsum.ne.0) THEN
                mfratio=refined(cr)%V2d_flux(1,m,tnew)/northsum
              ELSE
                mfratio=1.0_r8
              END IF
              WRITE (300,30) ibc, refined(cr)%V2d_flux(1,m,tnew),        &
     &                       northsum, mfratio
            END IF
          END DO
 
#  ifdef DISTRIBUTE
!
!  Collect data from all nodes.
!
          CALL mp_assemble (dg, model, nptswe, spv,                     &
     &                      bry_contact(iwest ,cjcr)%Mflux(jlb:))
          CALL mp_assemble (dg, model, nptswe, spv,                     &
     &                      bry_contact(ieast ,cjcr)%Mflux(jlb:))
          CALL mp_assemble (dg, model, nptssn, spv,                     &
     &                      bry_contact(isouth,cjcr)%Mflux(ilb:))
          CALL mp_assemble (dg, model, nptssn, spv,                     &
     &                      bry_contact(inorth,cjcr)%Mflux(ilb:))
#  endif
        END IF
      END DO
!
      FLUSH (300)
!
  10  FORMAT (/,1x,a,/,/,4x,'cr = ',i2.2,4x,'dg = ',i2.2,4x,'rg = ',    &
     &        i2.2,4x,'iif(rg) = ',i3.3,4x,'iic(rg) = ',i7.7,4x,        &
     &        'time(rg) = ',i7.7,/,/,2x,'Coarse',6x,'Coarse Grid',8x,   &
     &        'Fine Grid',11x,'Ratio',/,4x,'Jb',9x,'U2d_flux',9x,       &
     &        'SUM(U2d_flux)',/)
  20  FORMAT (/,1x,a,/,/,4x,'cr = ',i2.2,4x,'dg = ',i2.2,4x,'rg = ',    &
     &        i2.2,4x,'iif(rg) = ',i3.3,4x,'iic(rg) = ',i7.7,4x,        &
     &        'time(rg) = ',i7.7,/,/,2x,'Coarse',6x,'Coarse Grid',8x,   &
     &        'Fine Grid',11x,'Ratio',/,4x,'Ib',9x,'V2d_flux',9x,       &
     &        'SUM(V2d_flux)',/)
  30  FORMAT (4x,i4.4,3(3x,1p,e15.8))
!
      RETURN
      END SUBROUTINE check_massflux
# endif
 
# ifndef ONE_WAY
!
      FUNCTION do_twoway (model, nl, il, ng, istep) RESULT (doit)
!
!=======================================================================
!                                                                      !
!  This function determines if the two-way exchange between finer to   !
!  coarser grid is appropriate. In complex nesting applications with   !
!  telescoping grids, grid with RefineScale > 0 including finer grids  !
!  inside, the two-way feedback between finer and coarse grids needs   !
!  to be in the correct sequence.                                      !
!                                                                      !
!  This function is called from either main2d or main3d.               !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     nl         Latest time-stepped nested layer index (integer)      !
!     il         Current nested layer index in the two-way DO loop     !
!                  call from main2d or main3d (integer)                !
!     ng         Finer grid number to process in the fine to coarse    !
!                  feedback (integer)                                  !
!     istep      Current rime-step counter from STEP_LOOP in main2d    !
!                  or main3d (integer)                                 !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     doit       The value of the result is TRUE/FALSE if the exchage  !
!                  is required or not (logical)                        !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_nesting
      USE mod_scalars
!
!  Imported variable declarations.
!
      integer, intent(in) :: model, nl, il, ng, istep
!
!  Local variable declarations.
!
      logical :: doit
 
      integer :: dgc, dgf
!
!-----------------------------------------------------------------------
!  Determine if two-way exchage is required at the corrent time-step.
!-----------------------------------------------------------------------
!
      doit=.false.
!
!  Chech if ng is a refined grid, RefineScale(ng) > 0.
!
      IF (refinedgrid(ng).and.(refinescale(ng).gt.0)) THEN
!
!  If telescoping grid, the sequence of exchanges when all the nested
!  grid reach the coarser grid (ng=1) time needs to be done in the
!  correct sequence: finest grid exchage first.
!
        IF (telescoping(ng)) THEN
          dgc=coarserdonor(ng)
          dgf=finerdonor(ng)
          IF (twowayinterval(dgf).eq.dt(dgc)) THEN
            doit=.true.
          END IF
!
!  If not telescoping, the exchange is in terms of the "istep" counter
!  but also at the right sequence.
!
        ELSE
          IF ((istep.eq.refinesteps(ng)).and.(il.eq.nl)) THEN
            doit=.true.
          END IF
        END IF
      END IF
!
!  Update two-way exchange counter between fine and coarse grids.
!
      IF (doit) THEN
        IF (model.eq.iadm) THEN
          twowayinterval(ng)=twowayinterval(ng)-                        &
     &                       real(refinesteps(ng),r8)*dt(ng)
        ELSE
          twowayinterval(ng)=twowayinterval(ng)+                        &
     &                       real(refinesteps(ng),r8)*dt(ng)
        END IF
      END IF
 
      RETURN
      END FUNCTION do_twoway
# endif
!
      SUBROUTINE fill_contact (rg, model, tile,                         &
     &                         cr, Npoints, contact,                    &
     &                         gtype, mvname, SpValCheck,               &
     &                         LBi, UBi, LBj, UBj,                      &
     &                         Ac, Ar)
!
!=======================================================================
!                                                                      !
!  This routine is used during initialization to fill the contact      !
!  points of a specified grid metric array. We need to have metric     !
!  data in all the extended computational points of the grid. No       !
!  attempt is done here to interpolate such values since the are       !
!  read in "set_contact" from input contact points NetCDF file.        !
!  This routine just unpack data into global arrays and check if       !
!  all needed values are filled.                                       !
!                                                                      !
!  During allocation these special metric grid arrays are initialized  !
!  to "spval" to avoid resetting those values processed already from   !
!  the regular Grid NetCDF file.  That is, only those contact points   !
!  outside the physical grid are processed here.  This is a good way   !
!  to check if all the extra numerical points have been processed.     !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     rg         Receiver grid number (integer)                        !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!     cr         Contact region number to process (integer)            !
!     Npoints    Number of points in the contact region (integer)      !
!     contact    Contact region information variables (T_NGC structure)!
!     gtype      C-grid variable type (integer)                        !
!     mvname     Metrics variable name (string)                        !
!     SpValCheck Special value used to check if the contact point      !
!                  needs to be processed.                              !
!     LBi        Receiver grid, I-dimension Lower bound (integer)      !
!     UBi        Receiver grid, I-dimension Upper bound (integer)      !
!     LBj        Receiver grid, J-dimension Lower bound (integer)      !
!     UBj        Receiver grid, J-dimension Upper bound (integer)      !
!     Ac         Metric data at Contact points.                        !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     Ar         Updated metric grid array                             !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_ncparam
      USE mod_nesting
!
!  Imported variable declarations.
!
      integer, intent(in) :: rg, model, tile
      integer, intent(in) :: cr, gtype, npoints
      integer, intent(in) :: lbi, ubi, lbj, ubj
!
      real(dp), intent(in) :: spvalcheck
!
      character(len=*), intent(in) :: mvname
!
      TYPE (t_ngc), intent(in) :: contact(:)
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: ac(:)
      real(r8), intent(inout) :: ar(lbi:,lbj:)
# else
      real(r8), intent(in) :: ac(npoints)
      real(r8), intent(inout) :: ar(lbi:ubi,lbj:ubj)
# endif
!
!  Local variable declarations.
!
      integer :: i, j, m
      integer :: istr, iend, jstr, jend
!
!-----------------------------------------------------------------------
!  Interpolate 2D data from donor grid to receiver grid contact points.
!-----------------------------------------------------------------------
!
!
!  Set starting and ending tile indices for the receiver grid.
!
      SELECT CASE (gtype)
        CASE (p2dvar)
          istr=bounds(rg) % IstrP(tile)
          iend=bounds(rg) % IendP(tile)
          jstr=bounds(rg) % JstrP(tile)
          jend=bounds(rg) % JendP(tile)
        CASE (r2dvar)
          istr=bounds(rg) % IstrT(tile)
          iend=bounds(rg) % IendT(tile)
          jstr=bounds(rg) % JstrT(tile)
          jend=bounds(rg) % JendT(tile)
        CASE (u2dvar)
          istr=bounds(rg) % IstrP(tile)
          iend=bounds(rg) % IendT(tile)
          jstr=bounds(rg) % JstrT(tile)
          jend=bounds(rg) % JendT(tile)
        CASE (v2dvar)
          istr=bounds(rg) % IstrT(tile)
          iend=bounds(rg) % IendT(tile)
          jstr=bounds(rg) % JstrP(tile)
          jend=bounds(rg) % JendT(tile)
      END SELECT
!
!  Interpolate.
!
      DO m=1,npoints
        i=contact(cr)%Irg(m)
        j=contact(cr)%Jrg(m)
        IF (((istr.le.i).and.(i.le.iend)).and.                          &
     &      ((jstr.le.j).and.(j.le.jend))) THEN
!!        IF (Ar(i,j).gt.SpValCheck) THEN        ! Only process contact
            ar(i,j)=ac(m)                        ! points outside in the
!!        END IF                                 ! regular grid
        END IF
      END DO
!
      RETURN
      END SUBROUTINE fill_contact
 
# if defined MASKING || defined WET_DRY
!
      SUBROUTINE mask_hweights (ng, model, tile)
!
!=======================================================================
!                                                                      !
!  This routine scales the horizontal interpolation weights to account !
!  for Land/Sea masking in  the nested contact region.  If wetting and !
!  drying,  the scaling is done at every time step since the  Land/Sea !
!  masking is time dependent.                                          !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     rg         Receiver grid number (integer)                        !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     Lweight    Updated linear interpolation weights                  !
#  ifdef QUADRATIC_WEIGHTS
!     Qweight    Updated quadratic interpolation weights               !
#  endif
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_grid
      USE mod_nesting
      USE mod_scalars
!
#  ifdef DISTRIBUTE
      USE distribute_mod, ONLY : mp_assemble
#  endif
      USE strings_mod,    ONLY : founderror
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model, tile
!
!  Local variable declarations.
!
      integer :: cr, dg,  m, rg
      integer :: istr, iend, jstr, jend
      integer :: idg, idgp1, jdg, jdgp1
      integer :: npointsr, npointsu, npointsv
#  ifdef DISTRIBUTE
      integer :: lpoints, qpoints
#  endif
!
      real(r8) :: cff
      real(r8) :: lwsum, lwsum_check, masklwsum
#  ifdef QUADRATIC_WEIGHTS
      real(r8) :: qwsum, qwsum_check, maskqwsum
#  endif
#  ifdef DISTRIBUTE
      real(r8), parameter :: spv = 0.0_r8
#  endif
      real(r8), dimension(4) :: lweight
#  ifdef QUADRATIC_WEIGHTS
      real(r8), dimension(9) :: qweight
#  endif
#  ifdef DISTRIBUTE
      real(r8), allocatable :: lw(:,:)
#   ifdef QUADRATIC_WEIGHTS
      real(r8), allocatable :: qw(:,:)
#   endif
#  endif
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", mask_hweights"
!
!=======================================================================
!  If appropriate, scale horizontal interpolation weights to account
!  for Land/Sea masking.
!=======================================================================
!
      DO cr=1,ncontact
!
!  Get data donor and data receiver grid numbers.
!
        dg=rcontact(cr)%donor_grid
        rg=rcontact(cr)%receiver_grid
!
        scale_weights : IF (refinedgrid(rg).and.                        &
     &                      (dg.eq.ng).and.(dg.lt.rg)) THEN
!
!  Scale interpolation weigths for RHO-contact points.
!
          istr=bounds(dg) % Istr(tile)
          iend=bounds(dg) % Iend(tile)
          jstr=bounds(dg) % Jstr(tile)
          jend=bounds(dg) % Jend(tile)
!
!-----------------------------------------------------------------------
!  Scale horizontal interpolation weigths for RHO-contact points.
!-----------------------------------------------------------------------
!
          npointsr=rcontact(cr)%Npoints
 
#  ifdef DISTRIBUTE
!
!  If distributed-memory, allocate and initialize working arrays
!  with special value (zero) to facilitate the global reduction
!  when collecting data between all nodes.
!
          lpoints=4*npointsr
          IF (.not.allocated(lw)) THEN
            allocate ( lw(4,npointsr) )
          END IF
          lw=spv
 
#   ifdef QUADRATIC_WEIGHTS
          qpoints=9*npointsr
          IF (.not.allocated(qw)) THEN
            allocate ( qw(9,npointsr) )
          END IF
          qw=spv
#   endif
#  endif
!
!  Scale interpolation weights.
!
          DO m=1,npointsr
            idg  =rcontact(cr)%Idg(m)
            idgp1=min(idg+1, bounds(dg)%UBi(-1))
            jdg  =rcontact(cr)%Jdg(m)
            jdgp1=min(jdg+1, bounds(dg)%UBj(-1))
            IF (((istr.le.idg).and.(idg.le.iend)).and.                  &
     &          ((jstr.le.jdg).and.(jdg.le.jend))) THEN
!
!  Linear interpolation weights.
!
              masklwsum=grid(dg)%rmask(idg  ,jdg  )+                    &
     &                  grid(dg)%rmask(idgp1,jdg  )+                    &
     &                  grid(dg)%rmask(idgp1,jdgp1)+                    &
     &                  grid(dg)%rmask(idg  ,jdgp1)
              IF (masklwsum.lt.4) THEN
#  ifdef WET_DRY
                lweight(1)=rcontact(cr)%LweightUnmasked(1,m)*           &
     &                     grid(dg)%rmask_full(idg  ,jdg  )
                lweight(2)=rcontact(cr)%LweightUnmasked(2,m)*           &
     &                     grid(dg)%rmask_full(idgp1,jdg  )
                lweight(3)=rcontact(cr)%LweightUnmasked(3,m)*           &
     &                     grid(dg)%rmask_full(idgp1,jdgp1)
                lweight(4)=rcontact(cr)%LweightUnmasked(4,m)*           &
     &                     grid(dg)%rmask_full(idg  ,jdgp1)
#  else
                lweight(1)=rcontact(cr)%Lweight(1,m)*                   &
     &                     grid(dg)%rmask_full(idg  ,jdg  )
                lweight(2)=rcontact(cr)%Lweight(2,m)*                   &
     &                     grid(dg)%rmask_full(idgp1,jdg  )
                lweight(3)=rcontact(cr)%Lweight(3,m)*                   &
     &                     grid(dg)%rmask_full(idgp1,jdgp1)
                lweight(4)=rcontact(cr)%Lweight(4,m)*                   &
     &                     grid(dg)%rmask_full(idg  ,jdgp1)
#  endif
                lwsum=sum(lweight)
                IF (lwsum.gt.0) THEN
                  cff=1.0_r8/lwsum
                  lweight(1)=cff*lweight(1)    ! using only water points
                  lweight(2)=cff*lweight(2)    ! in the interpolation of
                  lweight(3)=cff*lweight(3)    ! of contact points
                  lweight(4)=cff*lweight(4)
                ELSE
                  lweight=0.0_r8               ! all the donor grid
                END IF                         ! corners are on land
                lwsum_check=sum(lweight)
#  ifdef DISTRIBUTE
                lw(1,m)=lweight(1)
                lw(2,m)=lweight(2)
                lw(3,m)=lweight(3)
                lw(4,m)=lweight(4)
              ELSE
                lw(1,m)=rcontact(cr)%Lweight(1,m)   ! we need to load
                lw(2,m)=rcontact(cr)%Lweight(2,m)   ! unscaled values
                lw(3,m)=rcontact(cr)%Lweight(3,m)   ! to facilitate
                lw(4,m)=rcontact(cr)%Lweight(4,m)   ! parallel reduction
#  else
                rcontact(cr)%Lweight(1,m)=lweight(1)
                rcontact(cr)%Lweight(2,m)=lweight(2)
                rcontact(cr)%Lweight(3,m)=lweight(3)
                rcontact(cr)%Lweight(4,m)=lweight(4)
#  endif
              END IF
 
#  ifdef QUADRATIC_WEIGHTS
!
!  Quadratic interpolation weights.
!
              maskqwsum=grid(dg)%rmask(idg-1, jdg-1)+                   &
     &                  grid(dg)%rmask(idg  , jdg-1)+                   &
     &                  grid(dg)%rmask(idgp1, jdg-1)+                   &
     &                  grid(dg)%rmask(idg-1, jdg  )+                   &
     &                  grid(dg)%rmask(idg  , jdg  )+                   &
     &                  grid(dg)%rmask(idgp1, jdg  )+                   &
     &                  grid(dg)%rmask(idg-1, jdgp1)+                   &
     &                  grid(dg)%rmask(idg  , jdgp1)+                   &
     &                  grid(dg)%rmask(idgp1, jdgp1)
              IF (maskqwsum.lt.9) THEN
#   ifdef WET_DRY
                qweight(1)=rcontact(cr)%QweightUnmasked(1,m)*           &
     &                     grid(dg)%rmask_full(idg-1,jdg-1)
                qweight(2)=rcontact(cr)%QweightUnmasked(2,m)*           &
     &                     grid(dg)%rmask_full(idg  ,jdg-1)
                qweight(3)=rcontact(cr)%QweightUnmasked(3,m)*           &
     &                     grid(dg)%rmask_full(idgp1,jdg-1)
                qweight(4)=rcontact(cr)%QweightUnmasked(4,m)*           &
     &                     grid(dg)%rmask_full(idg-1,jdg  )
                qweight(5)=rcontact(cr)%QweightUnmasked(5,m)*           &
     &                     grid(dg)%rmask_full(idg  ,jdg  )
                qweight(6)=rcontact(cr)%QweightUnmasked(6,m)*           &
     &                     grid(dg)%rmask_full(idgp1,jdg  )
                qweight(7)=rcontact(cr)%QweightUnmasked(7,m)*           &
     &                     grid(dg)%rmask_full(idg-1,jdgp1)
                qweight(8)=rcontact(cr)%QweightUnmasked(8,m)*           &
     &                     grid(dg)%rmask_full(idg  ,jdgp1)
                qweight(9)=rcontact(cr)%QweightUnmasked(9,m)*           &
     &                     grid(dg)%rmask_full(idgp1,jdgp1)
#   else
                qweight(1)=rcontact(cr)%Qweight(1,m)*                   &
     &                     grid(dg)%rmask_full(idg-1,jdg-1)
                qweight(2)=rcontact(cr)%Qweight(2,m)*                   &
     &                     grid(dg)%rmask_full(idg  ,jdg-1)
                qweight(3)=rcontact(cr)%Qweight(3,m)*                   &
     &                     grid(dg)%rmask_full(idgp1,jdg-1)
                qweight(4)=rcontact(cr)%Qweight(4,m)*                   &
     &                     grid(dg)%rmask_full(idg-1,jdg  )
                qweight(5)=rcontact(cr)%Qweight(5,m)*                   &
     &                     grid(dg)%rmask_full(idg  ,jdg  )
                qweight(6)=rcontact(cr)%Qweight(6,m)*                   &
     &                     grid(dg)%rmask_full(idgp1,jdg  )
                qweight(7)=rcontact(cr)%Qweight(7,m)*                   &
     &                     grid(dg)%rmask_full(idg-1,jdgp1)
                qweight(8)=rcontact(cr)%Qweight(8,m)*                   &
     &                     grid(dg)%rmask_full(idg  ,jdgp1)
                qweight(9)=rcontact(cr)%Qweight(9,m)*                   &
     &                     grid(dg)%rmask_full(idgp1,jdgp1)
#   endif
                qwsum=sum(qweight)
                IF (qwsum.gt.0) THEN
                  cff=1.0_r8/qwsum
                  qweight(1)=cff*qweight(1)    ! using only water points
                  qweight(2)=cff*qweight(2)    ! in the interpolation of
                  qweight(3)=cff*qweight(3)    ! of contact points
                  qweight(4)=cff*qweight(4)
                  qweight(5)=cff*qweight(5)
                  qweight(6)=cff*qweight(6)
                  qweight(7)=cff*qweight(7)
                  qweight(8)=cff*qweight(8)
                  qweight(9)=cff*qweight(9)
                ELSE
                  qweight=0.0_r8               ! all the donor grid
                END IF                         ! corners are on land
                qwsum_check=sum(qweight)
#   ifdef DISTRIBUTE
                qw(1,m)=qweight(1)
                qw(2,m)=qweight(2)
                qw(3,m)=qweight(3)
                qw(4,m)=qweight(4)
                qw(5,m)=qweight(5)
                qw(6,m)=qweight(6)
                qw(7,m)=qweight(7)
                qw(8,m)=qweight(8)
                qw(9,m)=qweight(9)
              ELSE
                qw(1,m)=rcontact(cr)%Qweight(1,m)   ! we need to load
                qw(2,m)=rcontact(cr)%Qweight(2,m)   ! unscaled values
                qw(3,m)=rcontact(cr)%Qweight(3,m)   ! to facilitate
                qw(4,m)=rcontact(cr)%Qweight(4,m)   ! parallel reduction
                qw(5,m)=rcontact(cr)%Qweight(5,m)
                qw(6,m)=rcontact(cr)%Qweight(6,m)
                qw(7,m)=rcontact(cr)%Qweight(7,m)
                qw(8,m)=rcontact(cr)%Qweight(8,m)
                qw(9,m)=rcontact(cr)%Qweight(9,m)
#   else
                rcontact(cr)%Qweight(1,m)=qweight(1)
                rcontact(cr)%Qweight(2,m)=qweight(2)
                rcontact(cr)%Qweight(3,m)=qweight(3)
                rcontact(cr)%Qweight(4,m)=qweight(4)
                rcontact(cr)%Qweight(5,m)=qweight(5)
                rcontact(cr)%Qweight(6,m)=qweight(6)
                rcontact(cr)%Qweight(7,m)=qweight(7)
                rcontact(cr)%Qweight(8,m)=qweight(8)
                rcontact(cr)%Qweight(9,m)=qweight(9)
#   endif
              END IF
#  endif
            END IF
          END DO
 
#  ifdef DISTRIBUTE
!
!  Exchange data between all parallel nodes.
!
          CALL mp_assemble (dg, model, lpoints, spv, lw)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#   ifdef QUADRATIC_WEIGHTS
          CALL mp_assemble (dg, model, qpoints, spv, qw)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#   endif
!
!  Load exchanged weights.
!
          rcontact(cr)%Lweight(1:4,1:npointsr)=lw(1:4,1:npointsr)
#   ifdef QUADRATIC_WEIGHTS
          rcontact(cr)%Qweight(1:9,1:npointsr)=qw(1:9,1:npointsr)
#   endif
!
!  Deallocate local work arrays.
!
          IF (allocated(lw)) THEN
            deallocate (lw)
          END IF
 
#   ifdef QUADRATIC_WEIGHTS
          IF (allocated(qw)) THEN
            deallocate (qw)
          END IF
#   endif
#  endif
!
!-----------------------------------------------------------------------
!  Scale interpolation weigths for U-contact points.
!-----------------------------------------------------------------------
!
          npointsu=ucontact(cr)%Npoints
 
#  ifdef DISTRIBUTE
!
!  If distributed-memory, allocate and initialize working arrays
!  with special value (zero) to facilitate the global reduction
!  when collecting data between all nodes.
!
          lpoints=4*npointsu
          IF (.not.allocated(lw)) THEN
            allocate ( lw(4,npointsu) )
          END IF
          lw=spv
 
#   ifdef QUADRATIC_WEIGHTS
          qpoints=9*npointsu
          IF (.not.allocated(qw)) THEN
            allocate ( qw(9,npointsu) )
          END IF
          qw=spv
#   endif
#  endif
!
!  Scale interpolation weights.
!
          DO m=1,npointsu
            idg  =ucontact(cr)%Idg(m)
            idgp1=min(idg+1, bounds(dg)%UBi(-1))
            jdg  =ucontact(cr)%Jdg(m)
            jdgp1=min(jdg+1, bounds(dg)%UBj(-1))
            IF (((istr.le.idg).and.(idg.le.iend)).and.                  &
     &          ((jstr.le.jdg).and.(jdg.le.jend))) THEN
!
!  Linear interpolation weights.
!
              masklwsum=grid(dg)%umask_full(idg  ,jdg  )+               &
     &                  grid(dg)%umask_full(idgp1,jdg  )+               &
     &                  grid(dg)%umask_full(idgp1,jdgp1)+               &
     &                  grid(dg)%umask_full(idg  ,jdgp1)
              IF (masklwsum.lt.4) THEN
#  ifdef WET_DRY
                lweight(1)=ucontact(cr)%LweightUnmasked(1,m)*           &
     &                     grid(dg)%umask_full(idg  ,jdg  )
                lweight(2)=ucontact(cr)%LweightUnmasked(2,m)*           &
     &                     grid(dg)%umask_full(idgp1,jdg  )
                lweight(3)=ucontact(cr)%LweightUnmasked(3,m)*           &
     &                     grid(dg)%umask_full(idgp1,jdgp1)
                lweight(4)=ucontact(cr)%LweightUnmasked(4,m)*           &
     &                     grid(dg)%umask_full(idg  ,jdgp1)
#  else
                lweight(1)=ucontact(cr)%Lweight(1,m)*                   &
     &                     grid(dg)%umask_full(idg  ,jdg  )
                lweight(2)=ucontact(cr)%Lweight(2,m)*                   &
     &                     grid(dg)%umask_full(idgp1,jdg  )
                lweight(3)=ucontact(cr)%Lweight(3,m)*                   &
     &                     grid(dg)%umask_full(idgp1,jdgp1)
                lweight(4)=ucontact(cr)%Lweight(4,m)*                   &
     &                     grid(dg)%umask_full(idg  ,jdgp1)
#  endif
                lwsum=sum(lweight)
                IF (lwsum.gt.0) THEN
                  cff=1.0_r8/lwsum
                  lweight(1)=cff*lweight(1)    ! using only water points
                  lweight(2)=cff*lweight(2)    ! in the interpolation of
                  lweight(3)=cff*lweight(3)    ! of contact points
                  lweight(4)=cff*lweight(4)
                ELSE
                  lweight=0.0_r8               ! All the donor grid
                END IF                         ! corners are on land
                lwsum_check=sum(lweight)
#  ifdef DISTRIBUTE
                lw(1,m)=lweight(1)
                lw(2,m)=lweight(2)
                lw(3,m)=lweight(3)
                lw(4,m)=lweight(4)
              ELSE
                lw(1,m)=ucontact(cr)%Lweight(1,m)   ! we need to load
                lw(2,m)=ucontact(cr)%Lweight(2,m)   ! unscaled values
                lw(3,m)=ucontact(cr)%Lweight(3,m)   ! to facilitate
                lw(4,m)=ucontact(cr)%Lweight(4,m)   ! parallel reduction
#  else
                ucontact(cr)%Lweight(1,m)=lweight(1)
                ucontact(cr)%Lweight(2,m)=lweight(2)
                ucontact(cr)%Lweight(3,m)=lweight(3)
                ucontact(cr)%Lweight(4,m)=lweight(4)
#  endif
              END IF
 
#  ifdef QUADRATIC_WEIGHTS
!
!  Quadratic interpolation weights.
!
              maskqwsum=grid(dg)%umask_full(idg-1, jdg-1)+              &
     &                  grid(dg)%umask_full(idg  , jdg-1)+              &
     &                  grid(dg)%umask_full(idgp1, jdg-1)+              &
     &                  grid(dg)%umask_full(idg-1, jdg  )+              &
     &                  grid(dg)%umask_full(idg  , jdg  )+              &
     &                  grid(dg)%umask_full(idgp1, jdg  )+              &
     &                  grid(dg)%umask_full(idg-1, jdgp1)+              &
     &                  grid(dg)%umask_full(idg  , jdgp1)+              &
     &                  grid(dg)%umask_full(idgp1, jdgp1)
              IF (maskqwsum.lt.9) THEN
#   ifdef WET_DRY
                qweight(1)=ucontact(cr)%QweightUnmasked(1,m)*           &
     &                     grid(dg)%umask_full(idg-1,jdg-1)
                qweight(2)=ucontact(cr)%QweightUnmasked(2,m)*           &
     &                     grid(dg)%umask_full(idg  ,jdg-1)
                qweight(3)=ucontact(cr)%QweightUnmasked(3,m)*           &
     &                     grid(dg)%umask_full(idgp1,jdg-1)
                qweight(4)=ucontact(cr)%QweightUnmasked(4,m)*           &
     &                     grid(dg)%umask_full(idg-1,jdg  )
                qweight(5)=ucontact(cr)%QweightUnmasked(5,m)*           &
     &                     grid(dg)%umask_full(idg  ,jdg  )
                qweight(6)=ucontact(cr)%QweightUnmasked(6,m)*           &
     &                     grid(dg)%umask_full(idgp1,jdg  )
                qweight(7)=ucontact(cr)%QweightUnmasked(7,m)*           &
     &                     grid(dg)%umask_full(idg-1,jdgp1)
                qweight(8)=ucontact(cr)%QweightUnmasked(8,m)*           &
     &                     grid(dg)%umask_full(idg  ,jdgp1)
                qweight(9)=ucontact(cr)%QweightUnmasked(9,m)*           &
     &                     grid(dg)%umask_full(idgp1,jdgp1)
#   else
                qweight(1)=ucontact(cr)%Qweight(1,m)*                   &
     &                     grid(dg)%umask_full(idg-1,jdg-1)
                qweight(2)=ucontact(cr)%Qweight(2,m)*                   &
     &                     grid(dg)%umask_full(idg  ,jdg-1)
                qweight(3)=ucontact(cr)%Qweight(3,m)*                   &
     &                     grid(dg)%umask_full(idgp1,jdg-1)
                qweight(4)=ucontact(cr)%Qweight(4,m)*                   &
     &                     grid(dg)%umask_full(idg-1,jdg  )
                qweight(5)=ucontact(cr)%Qweight(5,m)*                   &
     &                     grid(dg)%umask_full(idg  ,jdg  )
                qweight(6)=ucontact(cr)%Qweight(6,m)*                   &
     &                     grid(dg)%umask_full(idgp1,jdg  )
                qweight(7)=ucontact(cr)%Qweight(7,m)*                   &
     &                     grid(dg)%umask_full(idg-1,jdgp1)
                qweight(8)=ucontact(cr)%Qweight(8,m)*                   &
     &                     grid(dg)%umask_full(idg  ,jdgp1)
                qweight(9)=ucontact(cr)%Qweight(9,m)*                   &
     &                     grid(dg)%umask_full(idgp1,jdgp1)
#   endif
                qwsum=sum(qweight)
                IF (qwsum.gt.0) THEN
                  cff=1.0_r8/qwsum
                  qweight(1)=cff*qweight(1)    ! using only water points
                  qweight(2)=cff*qweight(2)    ! in the interpolation of
                  qweight(3)=cff*qweight(3)    ! of contact points
                  qweight(4)=cff*qweight(4)
                  qweight(5)=cff*qweight(5)
                  qweight(6)=cff*qweight(6)
                  qweight(7)=cff*qweight(7)
                  qweight(8)=cff*qweight(8)
                  qweight(9)=cff*qweight(9)
                ELSE
                  qweight=0.0_r8               ! All the donor grid
                END IF                         ! corners are on land
                qwsum_check=sum(qweight)
#   ifdef DISTRIBUTE
                qw(1,m)=qweight(1)
                qw(2,m)=qweight(2)
                qw(3,m)=qweight(3)
                qw(4,m)=qweight(4)
                qw(5,m)=qweight(5)
                qw(6,m)=qweight(6)
                qw(7,m)=qweight(7)
                qw(8,m)=qweight(8)
                qw(9,m)=qweight(9)
              ELSE
                qw(1,m)=ucontact(cr)%Qweight(1,m)   ! we need to load
                qw(2,m)=ucontact(cr)%Qweight(2,m)   ! unscaled values
                qw(3,m)=ucontact(cr)%Qweight(3,m)   ! to facilitate
                qw(4,m)=ucontact(cr)%Qweight(4,m)   ! parallel reduction
                qw(5,m)=ucontact(cr)%Qweight(5,m)
                qw(6,m)=ucontact(cr)%Qweight(6,m)
                qw(7,m)=ucontact(cr)%Qweight(7,m)
                qw(8,m)=ucontact(cr)%Qweight(8,m)
                qw(9,m)=ucontact(cr)%Qweight(9,m)
#   else
                ucontact(cr)%Qweight(1,m)=qweight(1)
                ucontact(cr)%Qweight(2,m)=qweight(2)
                ucontact(cr)%Qweight(3,m)=qweight(3)
                ucontact(cr)%Qweight(4,m)=qweight(4)
                ucontact(cr)%Qweight(5,m)=qweight(5)
                ucontact(cr)%Qweight(6,m)=qweight(6)
                ucontact(cr)%Qweight(7,m)=qweight(7)
                ucontact(cr)%Qweight(8,m)=qweight(8)
                ucontact(cr)%Qweight(9,m)=qweight(9)
#   endif
              END IF
#  endif
            END IF
          END DO
 
#  ifdef DISTRIBUTE
!
!  Exchange data between all parallel nodes.
!
          CALL mp_assemble (dg, model, lpoints, spv, lw)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#   ifdef QUADRATIC_WEIGHTS
          CALL mp_assemble (dg, model, qpoints, spv, qw)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#   endif
!
!  Load exchanged weights.
!
          ucontact(cr)%Lweight(1:4,1:npointsu)=lw(1:4,1:npointsu)
#   ifdef QUADRATIC_WEIGHTS
          ucontact(cr)%Qweight(1:9,1:npointsu)=qw(1:9,1:npointsu)
#   endif
!
!  Deallocate local work arrays.
!
          IF (allocated(lw)) THEN
            deallocate (lw)
          END IF
 
#   ifdef QUADRATIC_WEIGHTS
          IF (allocated(qw)) THEN
            deallocate (qw)
          END IF
#   endif
#  endif
!
!-----------------------------------------------------------------------
!  Scale interpolation weigths for V-contact points.
!-----------------------------------------------------------------------
!
          npointsv=vcontact(cr)%Npoints
 
#  ifdef DISTRIBUTE
!
!  If distributed-memory, allocate and initialize working arrays
!  with special value (zero) to facilitate the global reduction
!  when collecting data between all nodes.
!
          lpoints=4*npointsv
          IF (.not.allocated(lw)) THEN
            allocate ( lw(4,npointsv) )
          END IF
          lw=spv
 
#   ifdef QUADRATIC_WEIGHTS
          lpoints=9*npointsv
          IF (.not.allocated(qw)) THEN
            allocate ( qw(9,npointsv) )
          END IF
          qw=spv
#   endif
#  endif
!
!  Scale interpolation weights.
!
          DO m=1,npointsv
            idg  =vcontact(cr)%Idg(m)
            idgp1=min(idg+1, bounds(dg)%UBi(-1))
            jdg  =vcontact(cr)%Jdg(m)
            jdgp1=min(jdg+1, bounds(dg)%UBj(-1))
            IF (((istr.le.idg).and.(idg.le.iend)).and.                  &
     &          ((jstr.le.jdg).and.(jdg.le.jend))) THEN
!
!  Linear interpolation weights.
!
              masklwsum=grid(dg)%vmask_full(idg  ,jdg  )+               &
     &                  grid(dg)%vmask_full(idgp1,jdg  )+               &
     &                  grid(dg)%vmask_full(idgp1,jdgp1)+               &
     &                  grid(dg)%vmask_full(idg  ,jdgp1)
              IF (masklwsum.lt.4) THEN
#  ifdef WET_DRY
                lweight(1)=vcontact(cr)%LweightUnmasked(1,m)*           &
     &                     grid(dg)%vmask_full(idg  ,jdg  )
                lweight(2)=vcontact(cr)%LweightUnmasked(2,m)*           &
     &                     grid(dg)%vmask_full(idgp1,jdg  )
                lweight(3)=vcontact(cr)%LweightUnmasked(3,m)*           &
     &                     grid(dg)%vmask_full(idgp1,jdgp1)
                lweight(4)=vcontact(cr)%LweightUnmasked(4,m)*           &
     &                     grid(dg)%vmask_full(idg  ,jdgp1)
#  else
                lweight(1)=vcontact(cr)%Lweight(1,m)*                   &
     &                     grid(dg)%vmask_full(idg  ,jdg  )
                lweight(2)=vcontact(cr)%Lweight(2,m)*                   &
     &                     grid(dg)%vmask_full(idgp1,jdg  )
                lweight(3)=vcontact(cr)%Lweight(3,m)*                   &
     &                     grid(dg)%vmask_full(idgp1,jdgp1)
                lweight(4)=vcontact(cr)%Lweight(4,m)*                   &
     &                     grid(dg)%vmask_full(idg  ,jdgp1)
#  endif
                lwsum=sum(lweight)
                IF (lwsum.gt.0) THEN
                  cff=1.0_r8/lwsum
                  lweight(1)=cff*lweight(1)    ! using only water points
                  lweight(2)=cff*lweight(2)    ! in the interpolation of
                  lweight(3)=cff*lweight(3)    ! of contact points
                  lweight(4)=cff*lweight(4)
                ELSE
                  lweight=0.0_r8               ! All the donor grid
                END IF                         ! corners are on land
                lwsum_check=sum(lweight)
#  ifdef DISTRIBUTE
                lw(1,m)=lweight(1)
                lw(2,m)=lweight(2)
                lw(3,m)=lweight(3)
                lw(4,m)=lweight(4)
              ELSE
                lw(1,m)=vcontact(cr)%Lweight(1,m)   ! we need to load
                lw(2,m)=vcontact(cr)%Lweight(2,m)   ! unscaled values
                lw(3,m)=vcontact(cr)%Lweight(3,m)   ! to facilitate
                lw(4,m)=vcontact(cr)%Lweight(4,m)   ! parallel reduction
#  else
                vcontact(cr)%Lweight(1,m)=lweight(1)
                vcontact(cr)%Lweight(2,m)=lweight(2)
                vcontact(cr)%Lweight(3,m)=lweight(3)
                vcontact(cr)%Lweight(4,m)=lweight(4)
#  endif
              END IF
 
#  ifdef QUADRATIC_WEIGHTS
!
!  Quadratic interpolation weights.
!
              maskqwsum=grid(dg)%vmask_full(idg-1, jdg-1)+              &
     &                  grid(dg)%vmask_full(idg  , jdg-1)+              &
     &                  grid(dg)%vmask_full(idgp1, jdg-1)+              &
     &                  grid(dg)%vmask_full(idg-1, jdg  )+              &
     &                  grid(dg)%vmask_full(idg  , jdg  )+              &
     &                  grid(dg)%vmask_full(idgp1, jdg  )+              &
     &                  grid(dg)%vmask_full(idg-1, jdgp1)+              &
     &                  grid(dg)%vmask_full(idg  , jdgp1)+              &
     &                  grid(dg)%vmask_full(idgp1, jdgp1)
              IF (maskqwsum.lt.9) THEN
#   ifdef WET_DRY
                qweight(1)=vcontact(cr)%QweightUnmasked(1,m)*           &
     &                     grid(dg)%vmask_full(idg-1,jdg-1)
                qweight(2)=vcontact(cr)%QweightUnmasked(2,m)*           &
     &                     grid(dg)%vmask_full(idg  ,jdg-1)
                qweight(3)=vcontact(cr)%QweightUnmasked(3,m)*           &
     &                     grid(dg)%vmask_full(idgp1,jdg-1)
                qweight(4)=vcontact(cr)%QweightUnmasked(4,m)*           &
     &                     grid(dg)%vmask_full(idg-1,jdg  )
                qweight(5)=vcontact(cr)%QweightUnmasked(5,m)*           &
     &                     grid(dg)%vmask_full(idg  ,jdg  )
                qweight(6)=vcontact(cr)%QweightUnmasked(6,m)*           &
     &                     grid(dg)%vmask_full(idgp1,jdg  )
                qweight(7)=vcontact(cr)%QweightUnmasked(7,m)*           &
     &                     grid(dg)%vmask_full(idg-1,jdgp1)
                qweight(8)=vcontact(cr)%QweightUnmasked(8,m)*           &
     &                     grid(dg)%vmask_full(idg  ,jdgp1)
                qweight(9)=vcontact(cr)%QweightUnmasked(9,m)*           &
     &                     grid(dg)%vmask_full(idgp1,jdgp1)
#   else
                qweight(1)=vcontact(cr)%Qweight(1,m)*                   &
     &                     grid(dg)%vmask_full(idg-1,jdg-1)
                qweight(2)=vcontact(cr)%Qweight(2,m)*                   &
     &                     grid(dg)%vmask_full(idg  ,jdg-1)
                qweight(3)=vcontact(cr)%Qweight(3,m)*                   &
     &                     grid(dg)%vmask_full(idgp1,jdg-1)
                qweight(4)=vcontact(cr)%Qweight(4,m)*                   &
     &                     grid(dg)%vmask_full(idg-1,jdg  )
                qweight(5)=vcontact(cr)%Qweight(5,m)*                   &
     &                     grid(dg)%vmask_full(idg  ,jdg  )
                qweight(6)=vcontact(cr)%Qweight(6,m)*                   &
     &                     grid(dg)%vmask_full(idgp1,jdg  )
                qweight(7)=vcontact(cr)%Qweight(7,m)*                   &
     &                     grid(dg)%vmask_full(idg-1,jdgp1)
                qweight(8)=vcontact(cr)%Qweight(8,m)*                   &
     &                     grid(dg)%vmask_full(idg  ,jdgp1)
                qweight(9)=vcontact(cr)%Qweight(9,m)*                   &
     &                     grid(dg)%vmask_full(idgp1,jdgp1)
#   endif
                qwsum=sum(qweight)
                IF (qwsum.gt.0) THEN
                  cff=1.0_r8/qwsum
                  qweight(1)=cff*qweight(1)    ! using only water points
                  qweight(2)=cff*qweight(2)    ! in the interpolation of
                  qweight(3)=cff*qweight(3)    ! of contact points
                  qweight(4)=cff*qweight(4)
                  qweight(5)=cff*qweight(5)
                  qweight(6)=cff*qweight(6)
                  qweight(7)=cff*qweight(7)
                  qweight(8)=cff*qweight(8)
                  qweight(9)=cff*qweight(9)
                ELSE
                  qweight=0.0_r8               ! All the donor grid
                END IF                         ! corners are on land
                END IF
                qwsum_check=sum(qweight)
#   ifdef DISTRIBUTE
                qw(1,m)=qweight(1)
                qw(2,m)=qweight(2)
                qw(3,m)=qweight(3)
                qw(4,m)=qweight(4)
                qw(5,m)=qweight(5)
                qw(6,m)=qweight(6)
                qw(7,m)=qweight(7)
                qw(8,m)=qweight(8)
                qw(9,m)=qweight(9)
              ELSE
                qw(1,m)=vcontact(cr)%Qweight(1,m)   ! we need to load
                qw(2,m)=vcontact(cr)%Qweight(2,m)   ! unscaled values
                qw(3,m)=vcontact(cr)%Qweight(3,m)   ! to facilitate
                qw(4,m)=vcontact(cr)%Qweight(4,m)   ! parallel reduction
                qw(5,m)=vcontact(cr)%Qweight(5,m)
                qw(6,m)=vcontact(cr)%Qweight(6,m)
                qw(7,m)=vcontact(cr)%Qweight(7,m)
                qw(8,m)=vcontact(cr)%Qweight(8,m)
                qw(9,m)=vcontact(cr)%Qweight(9,m)
#   else
                vcontact(cr)%Qweight(1,m)=qweight(1)
                vcontact(cr)%Qweight(2,m)=qweight(2)
                vcontact(cr)%Qweight(3,m)=qweight(3)
                vcontact(cr)%Qweight(4,m)=qweight(4)
                vcontact(cr)%Qweight(5,m)=qweight(5)
                vcontact(cr)%Qweight(6,m)=qweight(6)
                vcontact(cr)%Qweight(7,m)=qweight(7)
                vcontact(cr)%Qweight(8,m)=qweight(8)
                vcontact(cr)%Qweight(9,m)=qweight(9)
#   endif
              END IF
#  endif
            END IF
          END DO
 
#  ifdef DISTRIBUTE
!
!  Exchange data between all parallel nodes.
!
          CALL mp_assemble (dg, model, lpoints, spv, lw)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#   ifdef QUADRATIC_WEIGHTS
          CALL mp_assemble (dg, model, qpoints, spv, qw)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#   endif
!
!  Load exchanged weights.
!
          vcontact(cr)%Lweight(1:4,1:npointsv)=lw(1:4,1:npointsv)
#   ifdef QUADRATIC_WEIGHTS
          vcontact(cr)%Qweight(1:9,1:npointsv)=qw(1:9,1:npointsv)
#   endif
!
!  Deallocate local work arrays.
!
          IF (allocated(lw)) THEN
            deallocate (lw)
          END IF
 
#   ifdef QUADRATIC_WEIGHTS
          IF (allocated(qw)) THEN
            deallocate (qw)
          END IF
#   endif
#  endif
 
        END IF scale_weights
      END DO
!
      RETURN
      END SUBROUTINE mask_hweights
# endif
!
      SUBROUTINE get_composite (ng, model, isection, tile)
!
!=======================================================================
!                                                                      !
!  This routine gets the donor grid data required to process the       !
!  contact points of the current composite grid. It extracts the       !
!  donor cell points containing each contact point. In composite       !
!  grids, it is possible to have more than one contact region.         !
!                                                                      !
!  The interpolation of composite grid contact points from donor       !
!  grid data is carried out in a different parallel region using       !
!  'put_composite'.                                                    !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Composite grid number (integer)                       !
!     model      Calling model identifier (integer)                    !
!     isection   Governing equations time-stepping section in          !
!                  main2d or main3d indicating which state             !
!                  variables to process (integer)                      !
!     tile       Domain tile partition (integer)                       !
!                                                                      !
!  On Output:    (mod_nesting)                                         !
!                                                                      !
!     COMPOSITE  Updated contact points structure.                     !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_coupling
      USE mod_forces
      USE mod_grid
      USE mod_ncparam
      USE mod_nesting
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model, isection, tile
!
!  Local variable declarations.
!
      integer :: cr, dg, rg, nrec, rec
# ifdef SOLVE3D
      integer :: itrc
# endif
      integer :: lbi, ubi, lbj, ubj
      integer :: tindex
!
!-----------------------------------------------------------------------
!  Get donor grid data needed to process composite grid contact points.
!  Only process those variables associated with the governing equation
!  time-stepping section.
!-----------------------------------------------------------------------
!
      DO cr=1,ncontact
!
!  Get data donor and data receiver grid numbers.
!
        dg=rcontact(cr)%donor_grid
        rg=rcontact(cr)%receiver_grid
!
!  Process only contact region data for requested nested grid "ng".
!
        IF (rg.eq.ng) THEN
!
!  Set donor grid lower and upper array indices.
!
          lbi=bounds(dg)%LBi(tile)
          ubi=bounds(dg)%UBi(tile)
          lbj=bounds(dg)%LBj(tile)
          ubj=bounds(dg)%UBj(tile)
!
!  Process bottom stress (bustr, bvstr).
!
          IF (isection.eq.nbstr) THEN
            CALL get_contact2d (dg, model, tile,                        &
     &                          u2dvar, vname(1,idubms),                &
     &                          cr, ucontact(cr)%Npoints, ucontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          forces(dg) % bustr,                     &
     &                          composite(cr) % bustr)
            CALL get_contact2d (dg, model, tile,                        &
     &                          v2dvar, vname(1,idvbms),                &
     &                          cr, vcontact(cr)%Npoints, vcontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          forces(dg) % bvstr,                     &
     &                          composite(cr) % bvstr)
          END IF
!
!  Process free-surface (zeta) at the appropriate time index.
!
          IF ((isection.eq.nfsic).or.                                   &
     &        (isection.eq.nzeta).or.                                   &
     &        (isection.eq.n2dps).or.                                   &
     &        (isection.eq.n2dcs)) THEN
            IF (isection.eq.nzeta) THEN
              nrec=2                   ! process time records 1 and 2
            ELSE
              nrec=1                   ! process knew record
            END IF
            DO rec=1,nrec
              IF (isection.eq.nzeta) THEN
                tindex=rec
              ELSE
                tindex=knew(dg)
              END IF
              CALL get_contact2d (dg, model, tile,                      &
     &                            r2dvar, vname(1,idfsur),              &
     &                            cr, rcontact(cr)%Npoints, rcontact,   &
     &                            lbi, ubi, lbj, ubj,                   &
     &                            ocean(dg) % zeta(:,:,tindex),         &
     &                            composite(cr) % zeta(:,:,rec))
            END DO
          END IF
 
# if !(defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3)
!
!  Process free-surface equation rigth-hand-side (rzeta) term.
!
          IF (isection.eq.n2dps) THEN
            tindex=1
            CALL get_contact2d (dg, model, tile,                        &
     &                          r2dvar, vname(1,idrzet),                &
     &                          cr, rcontact(cr)%Npoints, rcontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          ocean(dg) % rzeta(:,:,tindex),          &
     &                          composite(cr) % rzeta)
          END IF
# endif
!
!  Process 2D momentum components (ubar,vbar) at the appropriate time
!  index.
!
          IF ((isection.eq.n2dic).or.                                   &
     &        (isection.eq.n2dps).or.                                   &
     &        (isection.eq.n2dcs).or.                                   &
     &        (isection.eq.n3duv)) THEN
            IF (isection.eq.n3duv) THEN
              nrec=2                   ! process time records 1 and 2
            ELSE
              nrec=1                   ! process knew record
            END IF
            DO rec=1,nrec
              IF (isection.eq.n3duv) THEN
                tindex=rec
              ELSE
                tindex=knew(dg)
              END IF
              CALL get_contact2d (dg, model, tile,                      &
     &                            u2dvar, vname(1,idubar),              &
     &                            cr, ucontact(cr)%Npoints, ucontact,   &
     &                            lbi, ubi, lbj, ubj,                   &
     &                            ocean(dg) % ubar(:,:,tindex),         &
     &                            composite(cr) % ubar(:,:,rec))
              CALL get_contact2d (dg, model, tile,                      &
     &                            v2dvar, vname(1,idvbar),              &
     &                            cr, vcontact(cr)%Npoints, vcontact,   &
     &                            lbi, ubi, lbj, ubj,                   &
     &                            ocean(dg) % vbar(:,:,tindex),         &
     &                            composite(cr) % vbar(:,:,rec))
            END DO
          END IF
 
# ifdef SOLVE3D
!
!  Process time averaged free-surface (Zt_avg1) and 2D momentum fluxes
!  (DU_avg1, DV_avg1).
!
          IF (isection.eq.n2dfx) THEN
            CALL get_contact2d (dg, model, tile,                        &
     &                          r2dvar, 'Zt_avg1',                      &
     &                          cr, rcontact(cr)%Npoints, rcontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          coupling(dg) % Zt_avg1,                 &
     &                          composite(cr) % Zt_avg1)
            CALL get_contact2d (dg, model, tile,                        &
     &                          u2dvar, 'DU_avg1',                      &
     &                          cr, ucontact(cr)%Npoints, ucontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          coupling(dg) % DU_avg1,                 &
     &                          composite(cr) % DU_avg1)
            CALL get_contact2d (dg, model, tile,                        &
     &                          v2dvar, 'DV_avg1',                      &
     &                          cr, vcontact(cr)%Npoints, vcontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          coupling(dg) % DV_avg1,                 &
     &                          composite(cr) % DV_avg1)
          END IF
 
#  if !defined TS_FIXED
!
!  Process tracer variables (t) at the appropriate time index.
!
          IF ((isection.eq.ntvic).or.                                   &
     &        (isection.eq.nrhst).or.                                   &
     &        (isection.eq.n3dtv)) THEN
            DO itrc=1,nt(ng)
              IF (isection.eq.nrhst) THEN
                tindex=3
              ELSE
                tindex=nnew(dg)
              END IF
              CALL get_contact3d (dg, model, tile,                      &
     &                            r3dvar, vname(1,idtvar(itrc)),        &
     &                            cr, rcontact(cr)%Npoints, rcontact,   &
     &                            lbi, ubi, lbj, ubj, 1, n(dg),         &
     &                            ocean(dg) % t(:,:,:,tindex,itrc),     &
     &                            composite(cr) % t(:,:,:,itrc))
            END DO
          END IF
#  endif
!
!  Process 3D momentum (u, v) at the appropriate time-index.
!
          IF ((isection.eq.n3dic).or.                                   &
     &        (isection.eq.n3duv)) THEN
            tindex=nnew(dg)
            CALL get_contact3d (dg, model, tile,                        &
     &                          u3dvar, vname(1,iduvel),                &
     &                          cr, ucontact(cr)%Npoints, ucontact,     &
     &                          lbi, ubi, lbj, ubj, 1, n(dg),           &
     &                          ocean(dg) % u(:,:,:,tindex),            &
     &                          composite(cr) % u)
            CALL get_contact3d (dg, model, tile,                        &
     &                          v3dvar, vname(1,idvvel),                &
     &                          cr, vcontact(cr)%Npoints, vcontact,     &
     &                          lbi, ubi, lbj, ubj, 1, n(dg),           &
     &                          ocean(dg) % v(:,:,:,tindex),            &
     &                          composite(cr) % v)
          END IF
!
!  Process 3D momentum fluxes (Huon, Hvom).
!
          IF (isection.eq.n3duv) THEN
            CALL get_contact3d (dg, model, tile,                        &
     &                          u3dvar, 'Huon',                         &
     &                          cr, ucontact(cr)%Npoints, ucontact,     &
     &                          lbi, ubi, lbj, ubj, 1, n(dg),           &
     &                          grid(dg) % Huon,                        &
     &                          composite(cr) % Huon)
            CALL get_contact3d (dg, model, tile,                        &
     &                          v3dvar, 'Hvom',                         &
     &                          cr, vcontact(cr)%Npoints, vcontact,     &
     &                          lbi, ubi, lbj, ubj, 1, n(dg),           &
     &                          grid(dg) % Hvom,                        &
     &                          composite(cr) % Hvom)
          END IF
# endif
 
        END IF
      END DO
!
      RETURN
      END SUBROUTINE get_composite
!
      SUBROUTINE get_metrics (ng, model, tile)
!
!=======================================================================
!                                                                      !
!  This routine extracts grid spacing metrics "on_u" and "om_v",       !
!  which are used to impose mass flux at the finer grid physical       !
!  boundaries in refinement applications.                              !
!                                                                      !
!  The extracted metrics is stored in REFINED structure and are        !
!  needed in 'put_refine2d'.                                           !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Refinement grid number (integer)                      !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!                                                                      !
!  On Output:    (mod_nesting)                                         !
!                                                                      !
!     REFINED    Updated contact points structure.                     !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_grid
      USE mod_nesting
      USE mod_scalars
!
#  ifdef DISTRIBUTE
      USE distribute_mod, ONLY : mp_assemble
#  endif
      USE strings_mod,    ONLY : founderror
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model, tile
!
!  Local variable declarations.
!
      integer :: cr, dg, i, j, m
      integer :: istr, iend, jstr, jend
      integer :: npointsu, npointsv
!
      real(r8), parameter :: spv = 0.0_r8
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", get_metrics"
!
!-----------------------------------------------------------------------
!  Extract grid spacing metrics.
!-----------------------------------------------------------------------
!
      DO cr=1,ncontact
!
!  Get data donor grid number.
!
        dg=rcontact(cr)%donor_grid
!
!  Extract grid spacing at U-points.
!
        istr=bounds(dg) % IstrP(tile)
        iend=bounds(dg) % IendT(tile)
        jstr=bounds(dg) % JstrT(tile)
        jend=bounds(dg) % JendT(tile)
!
        npointsu=ucontact(cr) % Npoints
        refined(cr) % on_u(1:npointsu) = spv
!
        DO m=1,npointsu
          i=ucontact(cr) % Idg(m)
          j=ucontact(cr) % Jdg(m)
          IF (((istr.le.i).and.(i.le.iend)).and.                        &
     &        ((jstr.le.j).and.(j.le.jend))) THEN
            refined(cr) % on_u(m) = grid(dg) % on_u(i,j)
          END IF
        END DO
!
!  Extract grid spacing at V-points.
!
        istr=bounds(dg) % IstrT(tile)
        iend=bounds(dg) % IendT(tile)
        jstr=bounds(dg) % JstrP(tile)
        jend=bounds(dg) % JendT(tile)
!
        npointsv=vcontact(cr) % Npoints
        refined(cr) % om_v(1:npointsv) = spv
!
        DO m=1,npointsv
          i=vcontact(cr) % Idg(m)
          j=vcontact(cr) % Jdg(m)
          IF (((istr.le.i).and.(i.le.iend)).and.                        &
     &        ((jstr.le.j).and.(j.le.jend))) THEN
            refined(cr) % om_v(m) = grid(dg) % om_v(i,j)
          END IF
        END DO
 
# ifdef DISTRIBUTE
!
!  Exchange data between all nodes.
!
        CALL mp_assemble (dg, model, npointsu, spv, refined(cr) % on_u)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL mp_assemble (dg, model, npointsv, spv, refined(cr) % om_v)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
 
      END DO
!
      RETURN
      END SUBROUTINE get_metrics
!
      SUBROUTINE get_refine (ng, model, tile)
!
!=======================================================================
!                                                                      !
!  This routine gets the donor grid data required to process the       !
!  contact points of the  current  refinement  grid. It extracts       !
!  the donor cell points containing each contact point.                !
!                                                                      !
!  The extracted data is stored in two-time rolling records which      !
!  are needed for the space and time interpolation in 'put_refine'.    !
!                                                                      !
!  Except for initialization, this routine is called at the bottom     !
!  of the donor grid time step so all the values are updated for the   !
!  time(dg)+dt(dg). That is, in 2D applications it is called after     !
!  "step2d" corrector step and in 3D applications it is called after   !
!  "step3d_t". This is done to have the coarser grid snapshots at      !
!  time(dg) and time(dg)+dt(dg) to bound the interpolation of the      !
!  finer grid contact points.                                          !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Refinement grid number (integer)                      !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!                                                                      !
!  On Output:    (mod_nesting)                                         !
!                                                                      !
!     REFINED    Updated contact points structure.                     !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_coupling
      USE mod_ncparam
      USE mod_nesting
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model, tile
!
!  Local variable declarations.
!
# ifdef NESTING_DEBUG
      logical, save :: first = .true.
# endif
      integer :: tindex2d, cr, dg, ir, rg, tnew
# ifdef SOLVE3D
      integer :: tindex3d, itrc
# endif
      integer :: lbi, ubi, lbj, ubj
!
!-----------------------------------------------------------------------
!  Get donor grid data needed to process refinement grid contact points.
!  The extracted contact point data is stored in two time records to
!  facilitate the space-time interpolation elsewhere.
!-----------------------------------------------------------------------
!
      DO cr=1,ncontact
!
!  Get data donor and data receiver grid numbers.
!
        dg=rcontact(cr)%donor_grid
        rg=rcontact(cr)%receiver_grid
!
!  Process only contact region data for requested nested grid "ng".
!
        IF ((dg.eq.coarserdonor(rg)).and.(dg.eq.ng)) THEN
!
!  Set donor grid lower and upper array indices.
!
          lbi=bounds(dg)%LBi(tile)
          ubi=bounds(dg)%UBi(tile)
          lbj=bounds(dg)%LBj(tile)
          ubj=bounds(dg)%UBj(tile)
!
!  Update rolling time indices. The contact data is stored in two time
!  levels.  We need a special case for ROMS initialization in "main2d"
!  or "main3d" after the processing "ini_fields".  Notice that a dt(dg)
!  is added because this routine is called after the end of the time
!  step.
!
          IF (rollingindex(cr).eq.0) THEN
            tnew=1                                  ! ROMS
            rollingindex(cr)=tnew                   ! initialization
            rollingtime(tnew,cr)=time(dg)           ! step
          ELSE
            tnew=3-rollingindex(cr)
            rollingindex(cr)=tnew
            rollingtime(tnew,cr)=time(dg)+dt(dg)
          END IF
!
!  Set donor grid time index to process. In 3D applications, the 2D
!  record index to use can be either 1 or 2 since both ubar(:,:,1:2)
!  and vbar(:,:,1:2) are set to its time-averaged values in "step3d_uv".
!  That is, we can use Tindex2d=kstp(dg) or Tindex2d=knew(dg). However,
!  in 2D applications we need to use Tindex2d=knew(dg).
!
          tindex2d=knew(dg)
# ifdef SOLVE3D
          tindex3d=nnew(dg)
# endif
 
# ifdef NESTING_DEBUG
!
!  If debugging, write information into Fortran unit 100 to check the
!  logic of processing donor grid data.
!
          IF (domain(ng)%SouthWest_Test(tile)) THEN
            IF (master) THEN
              IF (first) THEN
                first=.false.
                WRITE (100,10)
              END IF
              WRITE (100,20) ng, cr, dg, rg, iic(dg), iic(rg),          &
#  ifdef SOLVE3D
     &                       3-tnew, tnew, tindex2d, tindex3d,          &
#  else
     &                       3-tnew, tnew, tindex2d, tindex2d,          &
#  endif
     &                       int(rollingtime(3-tnew,cr)),               &
     &                       int(rollingtime(tnew,cr)),                 &
     &                       int(time(ng))
 10           FORMAT (2x,'ng',2x,'cr',2x,'dg',2x,'rg',5x,'iic',4x,'iic',&
     &                3x,'told',3x,'tnew',2x,'Tindex',1x,'Tindex',      &
     &                2x,'time',3x,'time',3x,'time',/,20x,'(dg)',       &
     &                3x,'(rg)',18x,'2D',5x,'3D',4x,'told',3x,'tnew',   &
     &                3x,'(ng)',/)
 20           FORMAT (4i4,9i7)
              FLUSH (100)
            END IF
          END IF
# endif
!
!  Extract free-surface.
!
# ifdef SOLVE3D
          CALL get_contact2d (dg, model, tile,                          &
     &                        r2dvar, 'Zt_avg1',                        &
     &                        cr, rcontact(cr)%Npoints, rcontact,       &
     &                        lbi, ubi, lbj, ubj,                       &
     &                        coupling(dg) % Zt_avg1,                   &
     &                        refined(cr) % zeta(:,:,tnew))
# else
          CALL get_contact2d (dg, model, tile,                          &
     &                        r2dvar, 'zeta',                           &
     &                        cr, rcontact(cr)%Npoints, rcontact,       &
     &                        lbi, ubi, lbj, ubj,                       &
     &                        ocean(dg) % zeta(:,:,tindex2d),           &
     &                        refined(cr) % zeta(:,:,tnew))
# endif
!
!  Extract 2D momentum components (ubar, vbar).
!
          CALL get_contact2d (dg, model, tile,                          &
     &                        u2dvar, vname(1,idubar),                  &
     &                        cr, ucontact(cr)%Npoints, ucontact,       &
     &                        lbi, ubi, lbj, ubj,                       &
     &                        ocean(dg) % ubar(:,:,tindex2d),           &
     &                        refined(cr) % ubar(:,:,tnew))
 
          CALL get_contact2d (dg, model, tile,                          &
     &                        v2dvar, vname(1,idvbar),                  &
     &                        cr, vcontact(cr)%Npoints, vcontact,       &
     &                        lbi, ubi, lbj, ubj,                       &
     &                        ocean(dg) % vbar(:,:,tindex2d),           &
     &                        refined(cr) % vbar(:,:,tnew))
!
# ifdef SOLVE3D
!  Interpolate time-averaged fluxes (U2d_flux, V2d_flux) at contact
# else
!  Interpolate 2D momentum fluxes (U2d_flux, V2d_flux) at contact
# endif
!  points. They will be used later to impose mass flux conservation
!  at the finer grid boundary (see routine "put_refine2d").
!
          CALL get_persisted2d (dg, rg, model, tile,                    &
     &                          u2dvar, 'U2d_flux',                     &
     &                          cr, ucontact(cr)%Npoints, ucontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
# ifdef SOLVE3D
     &                          coupling(dg) % DU_avg2,                 &
# else
     &                          ocean(dg) % DU_flux,                    &
# endif
     &                          refined(cr) % U2d_flux(:,:,tnew))
 
          CALL get_persisted2d (dg, rg, model, tile,                    &
     &                          v2dvar, 'V2d_flux',                     &
     &                          cr, vcontact(cr)%Npoints, vcontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
# ifdef SOLVE3D
     &                          coupling(dg) % DV_avg2,                 &
# else
     &                          ocean(dg) % DV_flux,                    &
# endif
     &                          refined(cr) % V2d_flux(:,:,tnew))
 
# ifdef SOLVE3D
!
!  Tracer-type variables.
!
          DO itrc=1,nt(dg)
            CALL get_contact3d (dg, model, tile,                        &
     &                          r3dvar, vname(1,idtvar(itrc)),          &
     &                          cr, rcontact(cr)%Npoints, rcontact,     &
     &                          lbi, ubi, lbj, ubj, 1, n(dg),           &
     &                          ocean(dg) % t(:,:,:,tindex3d,itrc),     &
     &                          refined(cr) % t(:,:,:,tnew,itrc))
          END DO
!
!  Extract 3D momentum components (u, v).
!
          CALL get_contact3d (dg, model, tile,                          &
     &                        u3dvar, vname(1,iduvel),                  &
     &                        cr, ucontact(cr)%Npoints, ucontact,       &
     &                        lbi, ubi, lbj, ubj, 1, n(dg),             &
     &                        ocean(dg) % u(:,:,:,tindex3d),            &
     &                        refined(cr) % u(:,:,:,tnew))
 
          CALL get_contact3d (dg, model, tile,                          &
     &                        v3dvar, vname(1,idvvel),                  &
     &                        cr, vcontact(cr)%Npoints, vcontact,       &
     &                        lbi, ubi, lbj, ubj, 1, n(dg),             &
     &                        ocean(dg) % v(:,:,:,tindex3d),            &
     &                        refined(cr) % v(:,:,:,tnew))
# endif
        END IF
      END DO
!
      RETURN
      END SUBROUTINE get_refine
!
      SUBROUTINE put_composite (ng, model, isection, tile)
!
!=======================================================================
!                                                                      !
!  This routine interpolates composite grid contact points from donor  !
!  grid data extracted in routine 'get_composite'.                     !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Composite grid number (integer)                       !
!     model      Calling model identifier (integer)                    !
!     isection   Governing equations time-stepping section in          !
!                  main2d or main3d indicating which state             !
!                  variables to process (integer)                      !
!     tile       Domain tile partition (integer)                       !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_coupling
      USE mod_forces
      USE mod_grid
      USE mod_ncparam
      USE mod_nesting
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
 
# ifdef DISTRIBUTE
!
      USE mp_exchange_mod, ONLY : mp_exchange2d
#  ifdef SOLVE3D
      USE mp_exchange_mod, ONLY : mp_exchange3d, mp_exchange4d
#  endif
# endif
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model, isection, tile
!
!  Local variable declarations.
!
      integer :: dg, rg, cr, nrec, rec
# ifdef SOLVE3D
      integer :: itrc
# endif
      integer :: lbi, ubi, lbj, ubj
      integer :: tindex
!
!-----------------------------------------------------------------------
!  Interpolate composite grid contact points from donor grid data.
!  Only process those variables associated with the governing equation
!  time-stepping section.
!-----------------------------------------------------------------------
!
      cr_loop : DO cr=1,ncontact
!
!  Get data donor and data receiver grid numbers.
!
        dg=rcontact(cr)%donor_grid
        rg=rcontact(cr)%receiver_grid
!
!  Process only contact region data for requested nested grid "ng".
!
        IF (rg.eq.ng) THEN
!
!  Set receiver grid lower and upper array indices.
!
          lbi=bounds(rg)%LBi(tile)
          ubi=bounds(rg)%UBi(tile)
          lbj=bounds(rg)%LBj(tile)
          ubj=bounds(rg)%UBj(tile)
!
!  Process bottom stress (bustr, bvstr).
!
          IF (isection.eq.nbstr) THEN
            CALL put_contact2d (rg, model, tile,                        &
     &                          u2dvar, vname(1,idubms),                &
     &                          cr, ucontact(cr)%Npoints, ucontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
# ifdef MASKING
     &                          grid(rg) % umask,                       &
# endif
     &                          composite(cr) % bustr,                  &
     &                          forces(rg) % bustr)
            CALL put_contact2d (rg, model, tile,                        &
     &                          v2dvar, vname(1,idvbms),                &
     &                          cr, vcontact(cr)%Npoints, vcontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
# ifdef MASKING
     &                          grid(rg) % vmask,                       &
# endif
     &                          composite(cr) % bvstr,                  &
     &                          forces(rg) % bvstr)
# ifdef DISTRIBUTE
            CALL mp_exchange2d (rg, tile, model, 2,                     &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          nghostpoints,                           &
     &                          ewperiodic(rg), nsperiodic(rg),         &
     &                          forces(rg) % bustr,                     &
     &                          forces(rg) % bvstr)
# endif
          END IF
!
!  Process free-surface (zeta) at the appropriate time index.
!
          IF ((isection.eq.nfsic).or.                                   &
     &        (isection.eq.nzeta).or.                                   &
     &        (isection.eq.n2dps).or.                                   &
     &        (isection.eq.n2dcs)) THEN
            IF (isection.eq.nzeta) THEN
              nrec=2                   ! process time records 1 and 2
            ELSE
              nrec=1                   ! process knew record
            END IF
            DO rec=1,nrec
              IF (isection.eq.nzeta) THEN
                tindex=rec
              ELSE
                tindex=knew(rg)
              END IF
              CALL put_contact2d (rg, model, tile,                      &
     &                            r2dvar, vname(1,idfsur),              &
     &                            cr, rcontact(cr)%Npoints, rcontact,   &
     &                            lbi, ubi, lbj, ubj,                   &
# ifdef MASKING
     &                            grid(rg) % rmask,                     &
# endif
     &                            composite(cr) % zeta(:,:,rec),        &
     &                            ocean(rg) % zeta(:,:,tindex))
# ifdef DISTRIBUTE
              CALL mp_exchange2d (rg, tile, model, 1,                   &
     &                            lbi, ubi, lbj, ubj,                   &
     &                            nghostpoints,                         &
     &                            ewperiodic(rg), nsperiodic(rg),       &
     &                            ocean(rg) % zeta(:,:,tindex))
# endif
            END DO
          END IF
 
# if !(defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3)
!
!  Process free-surface equation rigth-hand-side (rzeta) term.
!
          IF (isection.eq.n2dps) THEN
            tindex=1
            CALL put_contact2d (rg, model, tile,                        &
     &                          r2dvar, vname(1,idrzet),                &
     &                          cr, rcontact(cr)%Npoints, rcontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
#  ifdef MASKING
     &                          grid(rg) % rmask,                       &
#  endif
     &                          composite(cr) % rzeta,                  &
     &                          ocean(rg) % rzeta(:,:,tindex))
#  ifdef DISTRIBUTE
            CALL mp_exchange2d (rg, tile, model, 1,                     &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          nghostpoints,                           &
     &                          ewperiodic(rg), nsperiodic(rg),         &
     &                          ocean(rg) % rzeta(:,:,tindex))
#  endif
          END IF
# endif
!
!  Process 2D momentum components (ubar,vbar) at the appropriate time
!  index.
!
          IF ((isection.eq.n2dic).or.                                   &
     &        (isection.eq.n2dps).or.                                   &
     &        (isection.eq.n2dcs).or.                                   &
     &        (isection.eq.n3duv)) THEN
            IF (isection.eq.n3duv) THEN
              nrec=2                   ! process time records 1 and 2
            ELSE
              nrec=1                   ! process KNEW record
            END IF
            DO rec=1,nrec
              IF (isection.eq.n3duv) THEN
                tindex=rec
              ELSE
                tindex=knew(rg)
              END IF
              CALL put_contact2d (rg, model, tile,                      &
     &                            u2dvar, vname(1,idubar),              &
     &                            cr, ucontact(cr)%Npoints, ucontact,   &
     &                            lbi, ubi, lbj, ubj,                   &
# ifdef MASKING
     &                            grid(rg) % umask,                     &
# endif
     &                            composite(cr) % ubar(:,:,rec),        &
     &                            ocean(rg) % ubar(:,:,tindex))
              CALL put_contact2d (rg, model, tile,                      &
     &                            v2dvar, vname(1,idvbar),              &
     &                            cr, vcontact(cr)%Npoints, vcontact,   &
     &                            lbi, ubi, lbj, ubj,                   &
# ifdef MASKING
     &                            grid(rg) % vmask,                     &
# endif
     &                            composite(cr) % vbar(:,:,rec),        &
     &                            ocean(rg) % vbar(:,:,tindex))
# ifdef DISTRIBUTE
              CALL mp_exchange2d (rg, tile, model, 2,                   &
     &                            lbi, ubi, lbj, ubj,                   &
     &                            nghostpoints,                         &
     &                            ewperiodic(rg), nsperiodic(rg),       &
     &                            ocean(rg) % ubar(:,:,tindex),         &
     &                            ocean(rg) % vbar(:,:,tindex))
# endif
            END DO
          END IF
 
# ifdef SOLVE3D
!
!  Process time averaged free-surface (Zt_avg1) and 2D momentum fluxes
!  (DU_avg1, DV_avg1).
!
          IF (isection.eq.n2dfx) THEN
            CALL put_contact2d (rg, model, tile,                        &
     &                          r2dvar, 'Zt_avg1',                      &
     &                          cr, rcontact(cr)%Npoints, rcontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
#  ifdef MASKING
     &                          grid(rg) % rmask,                       &
#  endif
     &                          composite(cr) % Zt_avg1,                &
     &                          coupling(rg) % Zt_avg1)
            CALL put_contact2d (rg, model, tile,                        &
     &                          u2dvar, vname(1,idufx1),                &
     &                          cr, ucontact(cr)%Npoints, ucontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
#  ifdef MASKING
     &                          grid(rg) % umask,                       &
#  endif
     &                          composite(cr) % DU_avg1,                &
     &                          coupling(rg) % DU_avg1)
            CALL put_contact2d (rg, model, tile,                        &
     &                          v2dvar, vname(1,idvfx1),                &
     &                          cr, vcontact(cr)%Npoints, vcontact,     &
     &                          lbi, ubi, lbj, ubj,                     &
#  ifdef MASKING
     &                          grid(rg) % vmask,                       &
#  endif
     &                          composite(cr) % DV_avg1,                &
     &                          coupling(rg) % DV_avg1)
# ifdef DISTRIBUTE
            CALL mp_exchange2d (rg, tile, model, 3,                     &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          nghostpoints,                           &
     &                          ewperiodic(rg), nsperiodic(rg),         &
     &                          coupling(rg) % Zt_avg1,                 &
     &                          coupling(rg) % DU_avg1,                 &
     &                          coupling(rg) % DV_avg1)
# endif
          END IF
 
#  if !defined TS_FIXED
!
!  Process tracer variables (t) at the appropriate time index.
!
          IF ((isection.eq.ntvic).or.                                   &
     &        (isection.eq.nrhst).or.                                   &
     &        (isection.eq.n3dtv)) THEN
            DO itrc=1,nt(ng)
              IF (isection.eq.nrhst) THEN
                tindex=3
              ELSE
                tindex=nnew(rg)
              END IF
              CALL put_contact3d (rg, model, tile,                      &
     &                            r3dvar, vname(1,idtvar(itrc)),        &
     &                            cr, rcontact(cr)%Npoints, rcontact,   &
     &                            lbi, ubi, lbj, ubj, 1, n(rg),         &
#   ifdef MASKING
     &                            grid(rg) % rmask,                     &
#   endif
     &                            composite(cr) % t(:,:,:,itrc),        &
     &                            ocean(rg) % t(:,:,:,tindex,itrc))
            END DO
#   ifdef DISTRIBUTE
            CALL mp_exchange4d (rg, tile, model, 1,                     &
     &                          lbi, ubi, lbj, ubj, 1, n(rg), 1, nt(rg),&
     &                          nghostpoints,                           &
     &                          ewperiodic(rg), nsperiodic(rg),         &
     &                          ocean(rg) % t(:,:,:,tindex,:))
#   endif
          END IF
#  endif
!
!  Process 3D momentum (u, v) at the appropriate time-index.
!
          IF ((isection.eq.n3dic).or.                                   &
     &        (isection.eq.n3duv)) THEN
            tindex=nnew(rg)
            CALL put_contact3d (rg, model, tile,                        &
     &                          u3dvar, vname(1,iduvel),                &
     &                          cr, ucontact(cr)%Npoints, ucontact,     &
     &                          lbi, ubi, lbj, ubj, 1, n(rg),           &
#  ifdef MASKING
     &                          grid(rg) % umask,                       &
#  endif
     &                          composite(cr) % u,                      &
     &                          ocean(rg) % u(:,:,:,tindex))
            CALL put_contact3d (rg, model, tile,                        &
     &                          v3dvar, vname(1,idvvel),                &
     &                          cr, vcontact(cr)%Npoints, vcontact,     &
     &                          lbi, ubi, lbj, ubj, 1, n(rg),           &
#  ifdef MASKING
     &                          grid(rg) % vmask,                       &
#  endif
     &                          composite(cr) % v,                      &
     &                          ocean(rg) % v(:,:,:,tindex))
#  ifdef DISTRIBUTE
            CALL mp_exchange3d (rg, tile, model, 2,                     &
     &                          lbi, ubi, lbj, ubj, 1, n(rg),           &
     &                          nghostpoints,                           &
     &                          ewperiodic(rg), nsperiodic(rg),         &
     &                          ocean(rg) % u(:,:,:,tindex),            &
     &                          ocean(rg) % v(:,:,:,tindex))
#  endif
          END IF
!
!  Process 3D momentum fluxes (Huon, Hvom).
!
          IF (isection.eq.n3duv) THEN
            CALL put_contact3d (rg, model, tile,                        &
     &                          u3dvar, 'Huon',                         &
     &                          cr, ucontact(cr)%Npoints, ucontact,     &
     &                          lbi, ubi, lbj, ubj, 1, n(rg),           &
#  ifdef MASKING
     &                          grid(rg) % umask,                       &
#  endif
     &                          composite(cr) % Huon,                   &
     &                          grid(rg) % Huon)
            CALL put_contact3d (rg, model, tile,                        &
     &                          v3dvar, 'Hvom',                         &
     &                          cr, vcontact(cr)%Npoints, vcontact,     &
     &                          lbi, ubi, lbj, ubj, 1, n(rg),           &
#  ifdef MASKING
     &                          grid(rg) % vmask,                       &
#  endif
     &                          composite(cr) % Hvom,                   &
     &                          grid(rg) % Hvom)
#  ifdef DISTRIBUTE
            CALL mp_exchange3d (rg, tile, model, 2,                     &
     &                          lbi, ubi, lbj, ubj, 1, n(rg),           &
     &                          nghostpoints,                           &
     &                          ewperiodic(rg), nsperiodic(rg),         &
     &                          grid(rg) % Huon,                        &
     &                          grid(rg) % Hvom)
#  endif
          END IF
# endif
 
        END IF
      END DO cr_loop
!
      RETURN
      END SUBROUTINE put_composite
!
      SUBROUTINE put_refine (ng, model, tile, LputFsur)
!
!=======================================================================
!                                                                      !
!  This routine interpolates refinement grid contact points from donor !
!  grid data extracted in routine 'get_refine'. Notice that because of !
!  shared-memory parallelism,  the free-surface is processed first and !
!  in a different parallel region.
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Refinement grid number (integer)                      !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!     LputFsur   Switch to process or not free-surface (logical)       !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_coupling
      USE mod_forces
      USE mod_grid
      USE mod_ncparam
      USE mod_nesting
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
!
!  Imported variable declarations.
!
      logical, intent(in) :: lputfsur
      integer, intent(in) :: ng, model, tile
!
!  Local variable declarations.
!
      integer :: dg, rg, cr, nrec, rec
# ifdef SOLVE3D
      integer :: itrc
# endif
      integer :: lbi, ubi, lbj, ubj
      integer :: tindex
!
!-----------------------------------------------------------------------
!  Interpolate refinement grid contact points from donor grid data
!  (space-time interpolation)
!-----------------------------------------------------------------------
!
      DO cr=1,ncontact
!
!  Get data donor and data receiver grid numbers.
!
        dg=rcontact(cr)%donor_grid
        rg=rcontact(cr)%receiver_grid
!
!  Process only contact region data for requested nested grid "ng", if
!  donor grid is coarser than receiver grid.  That is, we are only
!  processing external contact points areas.
!
        IF ((rg.eq.ng).and.(dxmax(dg).gt.dxmax(rg))) THEN
!
!  Set receiver grid lower and upper array indices.
!
          lbi=bounds(rg)%LBi(tile)
          ubi=bounds(rg)%UBi(tile)
          lbj=bounds(rg)%LBj(tile)
          ubj=bounds(rg)%UBj(tile)
!
!  Fill free-surface separatelly.
!
          IF (lputfsur) THEN
            CALL put_refine2d (ng, dg, cr, model, tile, lputfsur,       &
     &                         lbi, ubi, lbj, ubj)
          ELSE
!
!  Fill other 2D state variables (like momentum) contact points.
!
            CALL put_refine2d (ng, dg, cr, model, tile, lputfsur,       &
     &                         lbi, ubi, lbj, ubj)
 
# ifdef SOLVE3D
!
!  Fill 3D state variables contact points.
!
            CALL put_refine3d (ng, dg, cr, model, tile,                 &
     &                         lbi, ubi, lbj, ubj)
# endif
          END IF
        END IF
      END DO
!
      RETURN
      END SUBROUTINE put_refine
 
# ifdef SOLVE3D
!
      SUBROUTINE correct_tracer (ng, ngf, model, tile)
!
!=======================================================================
!                                                                      !
!  This routine corrects the tracer values in the coarser grid at the  !
!  location of the finer grid physical domain perimeter by comparing   !
!  vertically accumulated horizontal tracer flux (Hz*u*T/n, Hz*v*T/m)  !
!  in two-way nesting refinement:                                      !
!                                                                      !
!  coarse grid,  t(:,jb,:,nstp,:) = t(:,jb,:,nstp,:) - FacJ    (west,  !
!                                                               east)  !
!                t(ib,:,:,nstp,:) = t(ib,:,:,nstp,:) - FacI    (south, !
!                                                               north) !
!  where                                                               !
!                                                                      !
!                FacJ = (TFF(jb,itrc) - TFC(jb,itrc)) *                !
!                       pm(:,jb) * pn(:,jb) / D(:,jb)                  !
!                                                                      !
!                TFF(ib,itrc) = SUM[SUM[Tflux(ib,k,itrc)]]     finer   !
!                                                              grid    !
!                               for  k=1:N, 1:RefineScale      flux    !
!                                                                      !
!                TFC(ib,itrc) = SUM[Tflux(ib,k,itrc)]          coarser !
!                                                              grid    !
!                               for  k=1:N                     flux    !
!                                                                      !
!  Similarly, for the southern and northern tracer fluxes.             !
!                                                                      !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ngc        Coarser grid number (integer)                         !
!     ngf        Finer grid number (integer)                           !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!                                                                      !
!  On Output:    (mod_ocean)                                           !
!                                                                      !
!     t          Updated coarse grid tracer values at finer grid       !
!                perimeter                                             !
!                                                                      !
!=======================================================================
!
      USE mod_param
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, ngf, model, tile
!
!  Local variable declarations.
!
#  include "tile.h"
!
      CALL correct_tracer_tile (ng, ngf, model, tile,                   &
                                lbi, ubi, lbj, ubj,                     &
     &                          imins, imaxs, jmins, jmaxs)
!
      RETURN
      END SUBROUTINE correct_tracer
!
!***********************************************************************
      SUBROUTINE correct_tracer_tile (ngc, ngf, model, tile,            &
     &                                LBi, UBi, LBj, UBj,               &
     &                                IminS, ImaxS, JminS, JmaxS)
!***********************************************************************
!
      USE mod_param
      USE mod_clima
      USE mod_grid
      USE mod_ocean
      USE mod_nesting
      USE mod_scalars
      USE mod_stepping
 
#  ifdef DISTRIBUTE
!
      USE mp_exchange_mod, ONLY : mp_exchange4d
#  endif
!
!  Imported variable declarations.
!
      integer, intent(in) :: ngc, ngf, model, tile
      integer, intent(in) :: lbi, ubi, lbj, ubj
      integer, intent(in) :: imins, imaxs, jmins, jmaxs
!
!  Local variable declarations.
!
      integer :: iedge, ibc, ibc_min, ibc_max, ibf, io
      integer :: jedge, jbc, jbc_min, jbc_max, jbf, jo
      integer :: istr, iend, jstr, jend
      integer :: istrm2, iendp2, jstrm2, jendp2
      integer :: tindex, i, ic, isum, itrc, j, jsum, k, half
      integer :: cr, dg, dgcr, rg, rgcr
 
      real(r8) :: tfc, tff, tvalue, cff
 
      real(r8) :: dinv(imins:imaxs,jmins:jmaxs)
!
!-----------------------------------------------------------------------
!  Correct coarser grid tracer values at finer grid perimeter.
!-----------------------------------------------------------------------
!
!  Determine contact regions where coarse grid is the donor and coarse
!  grid is the receiver..
!
      DO cr=1,ncontact
        dg=donor_grid(cr)
        rg=receiver_grid(cr)
        IF ((ngc.eq.dg).and.(ngf.eq.rg)) THEN
          dgcr=cr                                   ! coarse is donor
        ELSE IF ((ngc.eq.rg).and.(ngf.eq.dg)) THEN
          rgcr=cr                                   ! coarse is receiver
        END IF
      END DO
!
!  Set tile starting and ending indices for coarser grid.
!
      istr  =bounds(ngc)%Istr  (tile)
      iend  =bounds(ngc)%Iend  (tile)
      jstr  =bounds(ngc)%Jstr  (tile)
      jend  =bounds(ngc)%Jend  (tile)
!
      istrm2=bounds(ngc)%Istrm2(tile)
      iendp2=bounds(ngc)%Iendp2(tile)
      jstrm2=bounds(ngc)%Jstrm2(tile)
      jendp2=bounds(ngc)%Jendp2(tile)
!
!  Compute coarser grid inverse water colunm thickness.
!
      DO j=jstrm2,jendp2
        DO i=istrm2,iendp2
          cff=grid(ngc)%Hz(i,j,1)
          DO k=2,n(rg)
            cff=cff+grid(ngc)%Hz(i,j,k)
          END DO
          dinv(i,j)=1.0_r8/cff
        END DO
      END DO
!
!  Set finer grid center (half) and offset indices (Io and Jo) for
!  coarser grid (I,J) coordinates.
!
      half=(refinescale(ngf)-1)/2
      io=half+1
      jo=half+1
!
!  Set coarse grid tracer index to correct. Since the exchange of data
!  is done at the bottom of main3d, we need to use the newest time
!  index, I think.
!
!!    Tindex=nstp(ngc)  ! HGA: Why this index is stable? newest is nosy
      tindex=nnew(ngc)
!
!=======================================================================
!  Compute vertically integrated horizontal advective tracer flux for
!  coarser at the finer grid physical boundary.  Then, correct coarser
!  grid tracer values at that boundary.
!=======================================================================
!
!  Initialize tracer counter index. The "tclm" array is only allocated
!  to the NTCLM fields that need to be processed. This is done to
!  reduce memory.
!
      ic=0
!
      t_loop : DO itrc=1,nt(ngc)
        ic=ic+1
!
!-----------------------------------------------------------------------
!  Finer grid western boundary.
!-----------------------------------------------------------------------
!
        ibc=i_left(ngf)
        jbc_min=j_bottom(ngf)
        jbc_max=j_top(ngf)-1              ! interior points, no top
!                                           left corner
        DO jbc=jstr,jend
          IF (((istr.le.ibc-1).and.(ibc-1.le.iend)).and.                &
     &        ((jbc_min.le.jbc).and.(jbc.le.jbc_max))) THEN
!
!  Sum vertically coarse grid horizontal advective tracer flux,
!  Hz*u*T/n, from last time-step.
!
            tfc=0.0_r8
            DO k=1,n(ngc)
              tfc=tfc+bry_contact(iwest,rgcr)%Tflux(jbc,k,itrc)
            END DO
!
!  Sum vertically and horizontally finer grid advective tracer flux.
!  This is a vertical and horizontal J-integral because "RefineScale"
!  sub-divisions are done in the finer grid in each single coarse grid
!  at the J-edge.
!
            tff=0.0_r8
            jedge=jo+(jbc-jbc_min)*refinescale(ngf)
            DO jsum=-half,half
              jbf=jedge+jsum
              DO k=1,n(ngf)
                tff=tff+bry_contact(iwest,dgcr)%Tflux(jbf,k,itrc)
              END DO
            END DO
!
!  Zeroth order correction to fine grid time integral (RIL, 2016).
!  Correct coarse grid tracer at the finer grid western boundary.
!
            cff=grid(ngc)%pm(ibc-1,jbc)*                                &
     &          grid(ngc)%pn(ibc-1,jbc)*                                &
     &          dinv(ibc-1,jbc)
            DO k=1,n(ngc)
              tvalue=max(0.0_r8,                                        &
     &                   ocean(ngc)%t(ibc-1,jbc,k,tindex,itrc)-         &
     &                   cff*(tff-tfc))
              IF (ltracerclm(itrc,ngc).and.lnudgetclm(itrc,ngc)) THEN
                tvalue=tvalue+                                          &
     &                 dt(ngc)*clima(ngc)%Tnudgcof(ibc-1,jbc,k,itrc)*   &
     &                 (clima(ngc)%tclm(ibc-1,jbc,k,ic)-tvalue)
              END IF
#  ifdef MASKING
              tvalue=tvalue*grid(ngc)%rmask(ibc-1,jbc)
#  endif
              ocean(ngc)%t(ibc-1,jbc,k,tindex,itrc)=tvalue
            END DO
          END IF
        END DO
!
!-----------------------------------------------------------------------
!  Finer grid eastern boundary.
!-----------------------------------------------------------------------
!
        ibc=i_right(ngf)
        jbc_min=j_bottom(ngf)
        jbc_max=j_top(ngf)-1              ! interior points, no top
!                                           right corner
        DO jbc=jstr,jend
          IF (((istr.le.ibc).and.(ibc.le.iend)).and.                    &
     &        ((jbc_min.le.jbc).and.(jbc.le.jbc_max))) THEN
!
!  Sum vertically coarse grid horizontal advective tracer flux,
!  Hz*u*T/n, from last time-step.
!
            tfc=0.0_r8
            DO k=1,n(ngc)
              tfc=tfc+bry_contact(ieast,rgcr)%Tflux(jbc,k,itrc)
            END DO
!
!  Sum vertically and horizontally finer grid advective tracer flux.
!  This is a vertical and horizontal J-integral because "RefineScale"
!  sub-divisions are done in the finer grid in each single coarse grid
!  at the J-edge.
!
            tff=0.0_r8
            jedge=jo+(jbc-jbc_min)*refinescale(ngf)
            DO jsum=-half,half
              jbf=jedge+jsum
              DO k=1,n(ngf)
                tff=tff+bry_contact(ieast,dgcr)%Tflux(jbf,k,itrc)
              END DO
            END DO
!
!  Zeroth order correction to fine grid time integral (RIL, 2016).
!  Correct coarse grid tracer at the finer grid eastern boundary.
!
            cff=grid(ngc)%pm(ibc,jbc)*                                  &
     &          grid(ngc)%pn(ibc,jbc)*                                  &
     &          dinv(ibc,jbc)
            DO k=1,n(ngc)
              tvalue=max(0.0_r8,                                        &
     &                   ocean(ngc)%t(ibc,jbc,k,tindex,itrc)-           &
     &                   cff*(tff-tfc))
              IF (ltracerclm(itrc,ngc).and.lnudgetclm(itrc,ngc)) THEN
                tvalue=tvalue+                                          &
     &                 dt(ngc)*clima(ngc)%Tnudgcof(ibc,jbc,k,itrc)*     &
     &                 (clima(ngc)%tclm(ibc,jbc,k,ic)-tvalue)
              END IF
#  ifdef MASKING
              tvalue=tvalue*grid(ngc)%rmask(ibc,jbc)
#  endif
              ocean(ngc)%t(ibc,jbc,k,tindex,itrc)=tvalue
            END DO
          END IF
        END DO
!
!-----------------------------------------------------------------------
!  Finer grid southern boundary.
!-----------------------------------------------------------------------
!
        jbc=j_bottom(ngf)
        ibc_min=i_left(ngf)
        ibc_max=i_right(ngf)-1            ! interior points, no bottom
!                                           right corner
        DO ibc=istr,iend
          IF (((ibc_min.le.ibc).and.(ibc.le.ibc_max)).and.              &
     &        ((jstr.le.jbc-1).and.(jbc-1.le.jend))) THEN
!
!  Sum vertically coarse grid horizontal advective tracer flux,
!  Hz*v*T/m, from last time-step.
!
            tfc=0.0_r8
            DO k=1,n(ngc)
              tfc=tfc+bry_contact(isouth,rgcr)%Tflux(ibc,k,itrc)
            END DO
!
!  Sum vertically and horizontally finer grid advective tracer flux.
!  This is a vertical and horizontal I-integral because "RefineScale"
!  sub-divisions are done in the finer grid in each single coarse grid
!  at the I-edge.
!
            tff=0.0_r8
            iedge=io+(ibc-ibc_min)*refinescale(ngf)
            DO isum=-half,half
              ibf=iedge+isum
              DO k=1,n(ngf)
                tff=tff+bry_contact(isouth,dgcr)%Tflux(ibf,k,itrc)
              END DO
            END DO
!
!  Zeroth order correction to fine grid time integral (RIL, 2016).
!  Correct coarse grid tracer at the finer grid southern boundary.
!
            cff=grid(ngc)%pm(ibc,jbc-1)*                                &
     &          grid(ngc)%pn(ibc,jbc-1)*                                &
     &          dinv(ibc,jbc-1)
            DO k=1,n(ngc)
              tvalue=max(0.0_r8,                                        &
     &                   ocean(ngc)%t(ibc,jbc-1,k,tindex,itrc)-         &
     &                   cff*(tff-tfc))
              IF (ltracerclm(itrc,ngc).and.lnudgetclm(itrc,ngc)) THEN
                tvalue=tvalue+                                          &
     &                 dt(ngc)*clima(ngc)%Tnudgcof(ibc,jbc-1,k,itrc)*   &
     &                 (clima(ngc)%tclm(ibc,jbc-1,k,ic)-tvalue)
              END IF
#  ifdef MASKING
              tvalue=tvalue*grid(ngc)%rmask(ibc,jbc-1)
#  endif
              ocean(ngc)%t(ibc,jbc-1,k,tindex,itrc)=tvalue
            END DO
          END IF
        END DO
!
!-----------------------------------------------------------------------
!  Finer grid northern boundary.
!-----------------------------------------------------------------------
!
        jbc=j_top(ngf)
        ibc_min=i_left(ngf)
        ibc_max=i_right(ngf)-1            ! interior points, no top
!                                           right corner
        DO ibc=istr,iend
          IF (((ibc_min.le.ibc).and.(ibc.le.ibc_max)).and.              &
     &        ((jstr.le.jbc).and.(jbc.le.jend))) THEN
!
!  Sum vertically coarse grid horizontal advective tracer flux,
!  Hz*v*T/m, from last time-step.
!
            tfc=0.0_r8
            DO k=1,n(ngc)
              tfc=tfc+bry_contact(inorth,rgcr)%Tflux(ibc,k,itrc)
            END DO
!
!  Sum vertically and horizontally finer grid advective tracer flux.
!  This is a vertical and horizontal I-integral because "RefineScale"
!  sub-divisions are done in the finer grid in each single coarse grid
!  at the I-edge.
!
            tff=0.0_r8
            iedge=io+(ibc-ibc_min)*refinescale(ngf)
            DO isum=-half,half
              ibf=iedge+isum
              DO k=1,n(ngf)
                tff=tff+bry_contact(inorth,dgcr)%Tflux(ibf,k,itrc)
              END DO
            END DO
!
!  Zeroth order correction to fine grid time integral.
!  Correct coarse grid tracer at the finer grid northern boundary.
!
            cff=grid(ngc)%pm(ibc,jbc)*                                  &
     &          grid(ngc)%pn(ibc,jbc)*                                  &
     &          dinv(ibc,jbc)
            DO k=1,n(ngc)
              tvalue=max(0.0_r8,                                        &
     &                   ocean(ngc)%t(ibc,jbc,k,tindex,itrc)-           &
     &                   cff*(tff-tfc))
              IF (ltracerclm(itrc,ngc).and.lnudgetclm(itrc,ngc)) THEN
                tvalue=tvalue+                                          &
     &                 dt(ngc)*clima(ngc)%Tnudgcof(ibc,jbc,k,itrc)*     &
     &                 (clima(ngc)%tclm(ibc,jbc,k,ic)-tvalue)
              END IF
#  ifdef MASKING
              tvalue=tvalue*grid(ngc)%rmask(ibc,jbc)
#  endif
              ocean(ngc)%t(ibc,jbc,k,tindex,itrc)=tvalue
            END DO
          END IF
        END DO
      END DO t_loop
 
#  ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Exchange boundary data.
!-----------------------------------------------------------------------
!
      CALL mp_exchange4d (ngc, tile, model, 1,                          &
     &                    lbi, ubi, lbj, ubj, 1, n(ngc),                &
     &                    1, nt(ngc),                                   &
     &                    nghostpoints,                                 &
     &                    ewperiodic(ngc), nsperiodic(ngc),             &
     &                    ocean(ngc)%t(:,:,:,tindex,:))
#  endif
!
      RETURN
      END SUBROUTINE correct_tracer_tile
# endif
!
      SUBROUTINE fine2coarse (ng, model, vtype, tile)
!
!=======================================================================
!                                                                      !
!  This routine replaces interior coarse grid data with the refined    !
!  averaged values: two-way nesting.                                   !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Refinement grid number (integer)                      !
!     model      Calling model identifier (integer)                    !
!     vtype      State variables to process (integer):                 !
!                  vtype = r2dvar      2D state variables              !
!                  vtype = r3dvar      3D state variables              !
!     tile       Domain tile partition (integer)                       !
!                                                                      !
!  On Output:    (mod_coupling, mod_ocean)                             !
!                                                                      !
!                Updated state variable with average refined grid      !
!                  solution                                            !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_coupling
      USE mod_forces
      USE mod_grid
      USE mod_iounits
      USE mod_ncparam
      USE mod_nesting
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
!
      USE exchange_2d_mod
# ifdef SOLVE3D
      USE exchange_3d_mod
# endif
# ifdef DISTRIBUTE
      USE mp_exchange_mod, ONLY : mp_exchange2d
#  ifdef SOLVE3D
      USE mp_exchange_mod, ONLY : mp_exchange3d, mp_exchange4d
#  endif
# endif
      USE strings_mod,     ONLY : founderror
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model, vtype, tile
!
!  Local variable declarations.
!
      logical :: areaavg
!
      integer :: lbid, ubid, lbjd, ubjd
      integer :: lbir, ubir, lbjr, ubjr
      integer :: dindex2d, rindex2d
# ifdef SOLVE3D
      integer :: dindex3d, rindex3d
# endif
      integer :: cr, dg, k, rg, nrec, rec
# ifdef SOLVE3D
      integer :: itrc
# endif
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", fine2coarse"
!
!-----------------------------------------------------------------------
!  Average interior fine grid state variable data to the coarse grid
!  location. Then, replace coarse grid values with averaged data.
!-----------------------------------------------------------------------
!
      DO cr=1,ncontact
!
!  Get data donor and data receiver grid numbers.
!
        dg=rcontact(cr)%donor_grid
        rg=rcontact(cr)%receiver_grid
!
!  Process contact region if the current refinement grid "ng" is the
!  donor grid.  The coarse grid "rg" is the receiver grid and the
!  contact structure has all the information necessary for fine to
!  coarse coupling. The donor grid size is always smaller than the
!  receiver coarser grid.
!
        IF ((ng.eq.dg).and.(dxmax(dg).lt.dxmax(rg))) THEN
!
!  Set donor and receiver grids lower and upper array indices.
!
          lbid=bounds(dg)%LBi(tile)
          ubid=bounds(dg)%UBi(tile)
          lbjd=bounds(dg)%LBj(tile)
          ubjd=bounds(dg)%UBj(tile)
!
          lbir=bounds(rg)%LBi(tile)
          ubir=bounds(rg)%UBi(tile)
          lbjr=bounds(rg)%LBj(tile)
          ubjr=bounds(rg)%UBj(tile)
!
!  Report.
!
          IF (domain(ng)%SouthWest_Test(tile)) THEN
            IF (master.and.(vtype.eq.r2dvar)) THEN
              WRITE (stdout,10) dg, rg, cr
  10          FORMAT (6x,'FINE2COARSE - exchanging data between grids:',&
     &                ' dg = ',i2.2,' and rg = ',i2.2,'  at cr = ',i2.2)
            END IF
          END IF
!
!  Set state variable indices to process for donor and receiver grids.
!  Since the exchange of data is done at the bottom of main2d/main3d,
!  we need to use the newest time indices.
!
          dindex2d=knew(dg)         ! Donor    2D variables index
          rindex2d=knew(rg)         ! Receiver 3D variables index
# ifdef SOLVE3D
          dindex3d=nnew(dg)         ! Donor    3D variables index
          rindex3d=nnew(rg)         ! Receiver 3D variables index
# endif
!
!-----------------------------------------------------------------------
!  Process 2D state variables.
!-----------------------------------------------------------------------
!
          IF (vtype.eq.r2dvar) THEN
!
!  Free-surface.
!
            areaavg=.false.
# ifdef SOLVE3D
            CALL fine2coarse2d (rg, dg, model, tile,                    &
     &                          r2dvar, 'Zt_avg1',                      &
     &                          areaavg, refinescale(dg),               &
     &                          cr, rcontact(cr)%Npoints, rcontact,     &
     &                          lbid, ubid, lbjd, ubjd,                 &
     &                          lbir, ubir, lbjr, ubjr,                 &
     &                          grid(dg)%om_r,                          &
     &                          grid(dg)%on_r,                          &
     &                          grid(rg)%pm,                            &
     &                          grid(rg)%pn,                            &
     &                          grid(rg)%h,                             &
#  ifdef MASKING
     &                          grid(dg)%rmask_full,                    &
     &                          grid(rg)%rmask,                         &
#  endif
     &                          coupling(dg)%Zt_avg1,                   &
     &                          coupling(rg)%Zt_avg1)
# else
            CALL fine2coarse2d (rg, dg, model, tile,                    &
     &                          r2dvar, vname(1,idfsur),                &
     &                          areaavg, refinescale(dg),               &
     &                          cr, rcontact(cr)%Npoints, rcontact,     &
     &                          lbid, ubid, lbjd, ubjd,                 &
     &                          lbir, ubir, lbjr, ubjr,                 &
     &                          grid(dg)%om_r,                          &
     &                          grid(dg)%on_r,                          &
     &                          grid(rg)%pm,                            &
     &                          grid(rg)%pn,                            &
     &                          grid(rg)%h,                             &
#  ifdef MASKING
     &                          grid(dg)%rmask,                         &
     &                          grid(rg)%rmask,                         &
#  endif
     &                          ocean(dg)%zeta(:,:,dindex2d),           &
     &                          ocean(rg)%zeta(:,:,rindex2d))
# endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Process 2D momentum components (ubar,vbar).
!
            areaavg=.false.
            CALL fine2coarse2d (rg, dg, model, tile,                    &
     &                          u2dvar, vname(1,idubar),                &
     &                          areaavg, refinescale(dg),               &
     &                          cr, ucontact(cr)%Npoints, ucontact,     &
     &                          lbid, ubid, lbjd, ubjd,                 &
     &                          lbir, ubir, lbjr, ubjr,                 &
     &                          grid(dg)%om_u,                          &
     &                          grid(dg)%on_u,                          &
     &                          grid(rg)%pm,                            &
     &                          grid(rg)%pn,                            &
     &                          grid(rg)%h,                             &
# ifdef MASKING
     &                          grid(dg)%umask_full,                    &
     &                          grid(rg)%umask_full,                    &
# endif
     &                          ocean(dg)%ubar(:,:,dindex2d),           &
# ifdef SOLVE3D
     &                          ocean(rg)%ubar(:,:,1),                  &
     &                          ocean(rg)%ubar(:,:,2))
# else
     &                          ocean(rg)%ubar(:,:,rindex2d))
# endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
            CALL fine2coarse2d (rg, dg, model, tile,                    &
     &                          v2dvar, vname(1,idvbar),                &
     &                          areaavg, refinescale(dg),               &
     &                          cr, vcontact(cr)%Npoints, vcontact,     &
     &                          lbid, ubid, lbjd, ubjd,                 &
     &                          lbir, ubir, lbjr, ubjr,                 &
     &                          grid(dg)%om_v,                          &
     &                          grid(dg)%on_v,                          &
     &                          grid(rg)%pm,                            &
     &                          grid(rg)%pn,                            &
     &                          grid(rg)%h,                             &
# ifdef MASKING
     &                          grid(dg)%vmask_full,                    &
     &                          grid(rg)%vmask_full,                    &
# endif
     &                          ocean(dg)%vbar(:,:,dindex2d),           &
# ifdef SOLVE3D
     &                          ocean(rg)%vbar(:,:,1),                  &
     &                          ocean(rg)%vbar(:,:,2))
# else
     &                          ocean(rg)%vbar(:,:,rindex2d))
# endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# ifdef SOLVE3D
!
!-----------------------------------------------------------------------
!  Process 3D state variables.
!-----------------------------------------------------------------------
!
          ELSE IF (vtype.eq.r3dvar) THEN
!
!  Tracer type-variables.
!
            areaavg=.false.
            DO itrc=1,nt(rg)
              CALL fine2coarse3d (rg, dg, model, tile,                  &
     &                            r3dvar, vname(1,idtvar(itrc)),        &
     &                            areaavg, refinescale(dg),             &
     &                            cr, rcontact(cr)%Npoints, rcontact,   &
     &                            lbid, ubid, lbjd, ubjd, 1, n(dg),     &
     &                            lbir, ubir, lbjr, ubjr, 1, n(rg),     &
     &                            grid(dg)%om_r,                        &
     &                            grid(dg)%on_r,                        &
     &                            grid(rg)%pm,                          &
     &                            grid(rg)%pn,                          &
#  ifdef MASKING
     &                            grid(dg)%rmask_full,                  &
     &                            grid(rg)%rmask,                       &
#  endif
     &                            ocean(dg)%t(:,:,:,dindex3d,itrc),     &
     &                            ocean(rg)%t(:,:,:,rindex3d,itrc))
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
!
!  Process 3D momentum components (u, v).
!
            areaavg=.false.
            CALL fine2coarse3d (rg, dg, model, tile,                    &
     &                          u3dvar, vname(1,iduvel),                &
     &                          areaavg, refinescale(dg),               &
     &                          cr, ucontact(cr)%Npoints, ucontact,     &
     &                          lbid, ubid, lbjd, ubjd, 1, n(dg),       &
     &                          lbir, ubir, lbjr, ubjr, 1, n(rg),       &
     &                          grid(dg)%om_u,                          &
     &                          grid(dg)%on_u,                          &
     &                          grid(rg)%pm,                            &
     &                          grid(rg)%pn,                            &
#  ifdef MASKING
     &                          grid(dg)%umask_full,                    &
     &                          grid(rg)%umask_full,                    &
#  endif
     &                          ocean(dg)%u(:,:,:,dindex3d),            &
     &                          ocean(rg)%u(:,:,:,rindex3d))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
            CALL fine2coarse3d (rg, dg, model, tile,                    &
     &                          v3dvar, vname(1,idvvel),                &
     &                          areaavg, refinescale(dg),               &
     &                          cr, vcontact(cr)%Npoints, vcontact,     &
     &                          lbid, ubid, lbjd, ubjd, 1, n(dg),       &
     &                          lbir, ubir, lbjr, ubjr, 1, n(rg),       &
     &                          grid(dg)%om_v,                          &
     &                          grid(dg)%on_v,                          &
     &                          grid(rg)%pm,                            &
     &                          grid(rg)%pn,                            &
#  ifdef MASKING
     &                          grid(dg)%vmask_full,                    &
     &                          grid(rg)%vmask_full,                    &
#  endif
     &                          ocean(dg)%v(:,:,:,dindex3d),            &
     &                          ocean(rg)%v(:,:,:,rindex3d))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
          END IF
!
!-----------------------------------------------------------------------
!  Exchange boundary data.
!-----------------------------------------------------------------------
!
          IF (ewperiodic(rg).or.nsperiodic(rg)) THEN
            IF (vtype.eq.r2dvar) THEN
# ifdef SOLVE3D
              CALL exchange_r2d_tile (rg, tile,                         &
     &                                lbir, ubir, lbjr, ubjr,           &
     &                                coupling(rg)%Zt_avg1)
              DO k=1,2
                CALL exchange_u2d_tile (rg, tile,                       &
     &                                  lbir, ubir, lbjr, ubjr,         &
     &                                  ocean(rg)%ubar(:,:,k))
                CALL exchange_v2d_tile (rg, tile,                       &
     &                                  lbir, ubir, lbjr, ubjr,         &
     &                                  ocean(rg)%vbar(:,:,k))
              END DO
# else
              CALL exchange_r2d_tile (rg, tile,                         &
     &                                lbir, ubir, lbjr, ubjr,           &
     &                                ocean(rg)%zeta(:,:,rindex2d))
              CALL exchange_u2d_tile (rg, tile,                         &
     &                                lbir, ubir, lbjr, ubjr,           &
     &                                ocean(rg)%ubar(:,:,rindex2d))
              CALL exchange_v2d_tile (rg, tile,                         &
     &                                lbir, ubir, lbjr, ubjr,           &
     &                                ocean(rg)%vbar(:,:,rindex2d))
 
# endif
# ifdef SOLVE3D
            ELSE IF (vtype.eq.r3dvar) THEN
              CALL exchange_u3d_tile (rg, tile,                         &
     &                                lbir, ubir, lbjr, ubjr, 1, n(rg), &
     &                                ocean(rg)%u(:,:,:,rindex3d))
              CALL exchange_v3d_tile (rg, tile,                         &
     &                                lbir, ubir, lbjr, ubjr, 1, n(rg), &
     &                                ocean(rg)%v(:,:,:,rindex3d))
              DO itrc=1,nt(rg)
                CALL exchange_r3d_tile (rg, tile,                       &
     &                                  lbir, ubir, lbjr, ubjr,         &
     &                                  1, n(rg),                       &
     &                                ocean(rg)%t(:,:,:,rindex3d,itrc))
              END DO
# endif
            END IF
          END IF
 
# ifdef DISTRIBUTE
!
          IF (vtype.eq.r2dvar) THEN
#  ifdef SOLVE3D
            CALL mp_exchange2d (rg, tile, model, 1,                     &
     &                          lbir, ubir, lbjr, ubjr,                 &
     &                          nghostpoints,                           &
     &                          ewperiodic(rg), nsperiodic(rg),         &
     &                          coupling(rg)%Zt_avg1)
            CALL mp_exchange2d (rg, tile, model, 4,                     &
     &                          lbir, ubir, lbjr, ubjr,                 &
     &                          nghostpoints,                           &
     &                          ewperiodic(rg), nsperiodic(rg),         &
     &                          ocean(rg)%ubar(:,:,1),                  &
     &                          ocean(rg)%vbar(:,:,1),                  &
     &                          ocean(rg)%ubar(:,:,2),                  &
     &                          ocean(rg)%vbar(:,:,2))
#  else
            CALL mp_exchange2d (rg, tile, model, 3,                     &
     &                          lbir, ubir, lbjr, ubjr,                 &
     &                          nghostpoints,                           &
     &                          ewperiodic(rg), nsperiodic(rg),         &
     &                          ocean(rg)%zeta(:,:,rindex2d),           &
     &                          ocean(rg)%ubar(:,:,rindex2d),           &
     &                          ocean(rg)%vbar(:,:,rindex2d))
#  endif
#  ifdef SOLVE3D
          ELSE IF (vtype.eq.r3dvar) THEN
            CALL mp_exchange3d (rg, tile, model, 2,                     &
     &                          lbir, ubir, lbjr, ubjr, 1, n(rg),       &
     &                          nghostpoints,                           &
     &                          ewperiodic(rg), nsperiodic(rg),         &
     &                          ocean(rg)%u(:,:,:,rindex3d),            &
     &                          ocean(rg)%v(:,:,:,rindex3d))
            CALL mp_exchange4d (rg, tile, model, 1,                     &
     &                          lbir, ubir, lbjr, ubjr, 1, n(rg),       &
     &                          1, nt(rg),                              &
     &                          nghostpoints,                           &
     &                          ewperiodic(rg), nsperiodic(rg),         &
     &                          ocean(rg)%t(:,:,:,rindex3d,:))
#  endif
          END IF
# endif
        END IF
      END DO
!
      RETURN
      END SUBROUTINE fine2coarse
!
      SUBROUTINE fine2coarse2d (ng, dg, model, tile,                    &
     &                          gtype, svname,                          &
     &                          AreaAvg, Rscale,                        &
     &                          cr, Npoints, contact,                   &
     &                          LBiF, UBiF, LBjF, UBjF,                 &
     &                          LBiC, UBiC, LBjC, UBjC,                 &
# ifdef DISTRIBUTE
     &                          Adx, Ady,                               &
# else
     &                          dxF, dyF,                               &
# endif
     &                          pmC, pnC, hhC,                          &
# ifdef MASKING
#  ifdef DISTRIBUTE
     &                          Amsk,                                   &
#  else
     &                          Fmsk,                                   &
#  endif
     &                          Cmsk,                                   &
# endif
# ifdef DISTRIBUTE
     &                          A,                                      &
# else
     &                          F,                                      &
# endif
     &                          C1, C2)
!
!=======================================================================
!                                                                      !
!  This routine replaces the coarse grid data inside the refinement    !
!  grid interior for a 2D state variable with its refined averaged     !
!  values: two-way nesting.                                            !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Coarser grid number (integer)                         !
!     dg         Finer grid number (integer)                           !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!     gtype      C-grid variable type (integer)                        !
!     svname     State variable name (string)                          !
!     AreaAvg    Switch for area averaging (logical)                   !
!     Rscale     Refinement grid scale (integer)                       !
!     cr         Contact region number to process (integer)            !
!     Npoints    Number of points in the contact zone (integer)        !
!     contact    Contact zone information variables (T_NGC structure)  !
!     LBiF       Finer grid, I-dimension Lower bound (integer)         !
!     UBiF       Finer grid, I-dimension Upper bound (integer)         !
!     LBjF       Finer grid, J-dimension Lower bound (integer)         !
!     UBjF       Finer grid, J-dimension Upper bound (integer)         !
!     LBiC       Coarser grid, I-dimension Lower bound (integer)       !
!     UBiC       Coarser grid, I-dimension Upper bound (integer)       !
!     LBjC       Coarser grid, J-dimension Lower bound (integer)       !
!     UBjC       Coarser grid, J-dimension Upper bound (integer)       !
# ifdef DISTRIBUTE
!     Adx        Finer grid, X-grid spacing (1/pm: om_r, om_u, om_v)   !
!     Ady        Finer grid, Y-grid spacing (1/pn: on_r, on_u, on_v)   !
# else
!     dxF        Finer grid, X-grid spacing (1/pm: om_r, om_u, om_v)   !
!     dyF        Finer grid, Y-grid spacing (1/pn: on_r, on_u, on_v)   !
# endif
!     pmC        Coarser grid, inverse X-grid spacing (1/dx) at RHO    !
!     pnC        Coarser grid, inverse Y-grid spacing (1/dy) at RHO    !
!     hhC        Coarser grid, bathymetry at RHO                       !
# ifdef MASKING
#  ifdef DISTRIBUTE
!     Amsk       Finer grid land/sea masking (2D array)                !
#  else
!     Fmsk       Finer grid land/sea masking (2D array)                !
#  endif
!     Cmsk       Coarser grid land/sea masking (2D array)              !
# endif
# ifdef DISTRIBUTE
!     A          Finer grid 2D data                                    !
# else
!     F          Finer grid 2D data                                    !
# endif
!     C1         Coarser grid 2D data, record 1                        !
!     C2         Coarser grid 2D data, record 2 (OPTIONAL)             !
!                                                                      !
!  On Output:    (mod_nesting)                                         !
!                                                                      !
!     C1         Updated Coarser grid 2D data, record 1                !
!     C2         Uodated Coarser grid 2D data, record 2 (OPTIONAL)     !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_ncparam
      USE mod_nesting
      USE mod_scalars
!
# ifdef DISTRIBUTE
      USE distribute_mod, ONLY : mp_aggregate2d
# endif
      USE strings_mod,    ONLY : founderror
!
!  Imported variable declarations.
!
      logical, intent(in) :: areaavg
      integer, intent(in) :: ng, dg, model, tile
      integer, intent(in) :: gtype, cr, npoints, rscale
      integer, intent(in) :: lbif, ubif, lbjf, ubjf
      integer, intent(in) :: lbic, ubic, lbjc, ubjc
!
      character(len=*), intent(in) :: svname
!
      TYPE (t_ngc), intent(in) :: contact(:)
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: pmc(lbic:,lbjc:)
      real(r8), intent(in) :: pnc(lbic:,lbjc:)
      real(r8), intent(in) :: hhc(lbic:,lbjc:)
#  ifdef MASKING
      real(r8), intent(in) :: cmsk(lbic:,lbjc:)
#   ifdef DISTRIBUTE
      real(r8), intent(in) :: amsk(lbif:,lbjf:)
#   else
      real(r8), intent(in) :: fmsk(lbif:,lbjf:)
#   endif
#  endif
#  ifdef DISTRIBUTE
      real(r8), intent(in) :: a(lbif:,lbjf:)
      real(r8), intent(in) :: adx(lbif:,lbjf:)
      real(r8), intent(in) :: ady(lbif:,lbjf:)
#  else
      real(r8), intent(in) :: f(lbif:,lbjf:)
      real(r8), intent(in) :: dxf(lbif:,lbjf:)
      real(r8), intent(in) :: dyf(lbif:,lbjf:)
#  endif
      real(r8), intent(inout) :: c1(lbic:,lbjc:)
      real(r8), intent(inout), optional :: c2(lbic:,lbjc:)
# else
      real(r8), intent(in) :: pmc(lbic:ubic,lbjc:ubjc)
      real(r8), intent(in) :: pnc(lbic:ubic,lbjc:ubjc)
      real(r8), intent(in) :: hhc(lbic:ubic,lbjc:ubjc)
#  ifdef MASKING
      real(r8), intent(in) :: cmsk(lbic:ubic,lbjc:ubjc)
#   ifdef DISTRIBUTE
      real(r8), intent(in) :: amsk(lbif:ubif,lbjf:ubjf)
#   else
      real(r8), intent(in) :: fmsk(lbif:ubif,lbjf:ubjf)
#   endif
#  endif
#  ifdef DISTRIBUTE
      real(r8), intent(in) :: a(lbif:ubif,lbjf:ubjf)
      real(r8), intent(in) :: adx(lbif:ubif,lbjf:ubjf)
      real(r8), intent(in) :: adx(lbif:ubif,lbjf:ubjf)
#  else
      real(r8), intent(in) :: f(lbif:ubif,lbjf:ubjf)
      real(r8), intent(in) :: dxf(lbif:ubif,lbjf:ubjf)
      real(r8), intent(in) :: dyf(lbif:ubif,lbjf:ubjf)
#  endif
      real(r8), intent(inout) :: c1(lbic:ubic,lbjc:ubjc)
      real(r8), intent(inout), optional :: c2(lbic:ubic,lbjc:ubjc)
# endif
!
!  Local variable declarations.
!
      integer :: iadd, ic, jadd, jc, half, i, j, m
      integer :: ib_east, ib_west, jb_north, jb_south
# ifdef DISTRIBUTE
      integer :: lbi, ubi, lbj, ubj
# endif
!
      real(r8) :: areac_inv, my_area, my_areasum, ratio
      real(r8) :: my_avg, my_count, my_sum
 
# ifdef DISTRIBUTE
      real(r8), allocatable :: f(:,:)
      real(r8), allocatable :: dxf(:,:)
      real(r8), allocatable :: dyf(:,:)
#  ifdef MASKING
      real(r8), allocatable :: fmsk(:,:)
#  endif
# endif
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", fine2coarse2d"
 
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
!  Average interior fine grid state variable data to the coarse grid
!  location. Then, replace coarse grid values with averaged data.
!-----------------------------------------------------------------------
 
# ifdef DISTRIBUTE
!
!  Allocate global work array(s).
!
      lbi=bounds(dg)%LBi(-1)
      ubi=bounds(dg)%UBi(-1)
      lbj=bounds(dg)%LBj(-1)
      ubj=bounds(dg)%UBj(-1)
      IF (.not.allocated(f)) THEN
        allocate ( f(lbi:ubi,lbj:ubj) )
      END IF
      IF (areaavg) THEN
        IF (.not.allocated(dxf)) THEN
          allocate ( dxf(lbi:ubi,lbj:ubj) )
        END IF
        IF (.not.allocated(dyf)) THEN
          allocate ( dyf(lbi:ubi,lbj:ubj) )
        END IF
      END IF
#  ifdef MASKING
      IF (.not.allocated(fmsk)) THEN
        allocate ( fmsk(lbi:ubi,lbj:ubj) )
      END IF
#  endif
!
!  Gather finer grid data from all nodes in the group to build a global
!  array.
!
      CALL mp_aggregate2d (dg, model, gtype,                            &
     &                     lbif, ubif, lbjf, ubjf,                      &
     &                     lbi,  ubi,  lbj,  ubj,                       &
     &                     a, f)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      IF (areaavg) THEN
        CALL mp_aggregate2d (dg, model, gtype,                          &
     &                       lbif, ubif, lbjf, ubjf,                    &
     &                       lbi,  ubi,  lbj,  ubj,                     &
     &                       adx, dxf)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL mp_aggregate2d (dg, model, gtype,                          &
     &                       lbif, ubif, lbjf, ubjf,                    &
     &                       lbi,  ubi,  lbj,  ubj,                     &
     &                       ady, dyf)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#  ifdef MASKING
!
      CALL mp_aggregate2d (dg, model, gtype,                            &
     &                     lbif, ubif, lbjf, ubjf,                      &
     &                     lbi,  ubi,  lbj,  ubj,                       &
     &                     amsk, fmsk)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
# endif
!
!  Average finer grid data to coarse grid according to the refinement
!  ratio.
!
      half=(rscale-1)/2
      IF (areaavg) THEN               ! area averaging
        DO m=1,npoints
          i=contact(cr)%Idg(m)
          j=contact(cr)%Jdg(m)
          ic=contact(cr)%Irg(m)
          jc=contact(cr)%Jrg(m)
          IF (((istr.le.ic).and.(ic.le.iend)).and.                      &
     &        ((jstr.le.jc).and.(jc.le.jend))) THEN
            my_count=0.0_r8
            my_sum=0.0_r8
            my_areasum=0.0_r8
            DO jadd=-half,half
              DO iadd=-half,half
                my_area=dxf(i+iadd,j+jadd)*dyf(i+iadd,j+jadd)
                my_areasum=my_areasum+my_area
# ifdef MASKING
                my_sum=my_sum+                                          &
     &                 f(i+iadd,j+jadd)*my_area*                        &
     &                 min(1.0_r8,fmsk(i+iadd,j+jadd))
                my_count=my_count+min(1.0_r8,fmsk(i+iadd,j+jadd))
# else
                my_sum=my_sum+                                          &
     &                 f(i+iadd,j+jadd)*my_area
# endif
              END DO
            END DO
            SELECT CASE (gtype)              ! coarse grid inverse area
              CASE (r2dvar)
                areac_inv=pmc(ic,jc)*pnc(ic,jc)
              CASE (u2dvar)
                areac_inv=0.25_r8*(pmc(ic-1,jc)+pmc(ic,jc))*            &
     &                            (pnc(ic-1,jc)+pnc(ic,jc))
              CASE (v2dvar)
                areac_inv=0.25_r8*(pmc(ic,jc-1)+pmc(ic,jc))*            &
     &                            (pnc(ic,jc-1)+pnc(ic,jc))
              CASE DEFAULT
                areac_inv=pmc(ic,jc)*pnc(ic,jc)
            END SELECT
            ratio=my_areasum*areac_inv       ! for debugging purposes
            my_avg=my_sum*areac_inv
# ifdef MASKING
            IF (my_count.gt.0.0_r8) THEN
              my_avg=my_avg*rscale*rscale/my_count
            END IF
            my_avg=my_avg*cmsk(ic,jc)
# endif
            c1(ic,jc)=my_avg
            IF (PRESENT(c2)) THEN
              c2(ic,jc)=my_avg
            END IF
          END IF
        END DO
      ELSE                            ! simple averaging
        DO m=1,npoints
          i=contact(cr)%Idg(m)
          j=contact(cr)%Jdg(m)
          ic=contact(cr)%Irg(m)
          jc=contact(cr)%Jrg(m)
          IF (((istr.le.ic).and.(ic.le.iend)).and.                      &
     &        ((jstr.le.jc).and.(jc.le.jend))) THEN
            my_count=0.0_r8
            my_avg=0.0_r8
            my_sum=0.0_r8
            DO jadd=-half,half
              DO iadd=-half,half
# ifdef MASKING
                my_sum=my_sum+                                          &
     &                 f(i+iadd,j+jadd)*fmsk(i+iadd,j+jadd)
                my_count=my_count+min(1.0_r8,fmsk(i+iadd,j+jadd))
# else
                my_sum=my_sum+                                          &
     &                 f(i+iadd,j+jadd)
                my_count=my_count+1.0_r8
# endif
              END DO
            END DO
            IF (my_count.gt.0.0_r8) my_avg=my_sum/my_count
# ifdef MASKING
            my_avg=my_avg*cmsk(ic,jc)
#  ifdef WET_DRY
            IF (gtype.eq.r2dvar) THEN
              IF (my_avg.le.(dcrit(ng)-hhc(ic,jc))) THEN
                my_avg=dcrit(ng)-hhc(ic,jc)
              END IF
            END IF
#  endif
# endif
!
            c1(ic,jc)=my_avg
            IF (PRESENT(c2)) THEN
              c2(ic,jc)=my_avg
            END IF
          END IF
        END DO
      END IF
# ifdef REFINE_BOUNDARY
!
!-----------------------------------------------------------------------
!  Average finer grid BOUNDARY data to coarse grid U-type and V-type
!  variables according to the refinement ratio.
!-----------------------------------------------------------------------
!
!  U-type variables finer grid eastern and western boundaries.
!
      IF (gtype.eq.u2dvar) THEN
!
!  Get indices of coarser grid (ng) corresponding to the corners of
!  the finer grid (dg).
!
        ib_west=i_left(dg)
        ib_east=i_right(dg)
        jb_south=j_bottom(dg)
        jb_north=j_top(dg)
!
!  Eastern boundary.
!
        i=lm(dg)+1                            ! donor finer grid
        ic=ib_east                            ! receiver coarser grid
        DO jc=jb_south,jb_north-1
          j=(jc-jb_south)*rscale+half+1
          IF (((istr.le.ic).and.(ic.le.iend)).and.                      &
     &        ((jstr.le.jc).and.(jc.le.jend))) THEN
            my_count=0.0_r8
            my_avg=0.0_r8
            my_sum=0.0_r8
            DO jadd=-half,half
              my_sum=my_sum+                                            &
     &               f(i,j+jadd)
              my_count=my_count+1.0_r8
            END DO
            my_avg=my_sum/my_count
            c1(ic,jc)=my_avg
          END IF
        END DO
!
!  Western boundary.
!
        i=1                                   ! donor finer grid
        ic=ib_west                            ! receiver coarser grid
        DO jc=jb_south,jb_north-1
          j=(jc-jb_south)*rscale+half+1
          IF (((istr.le.ic).and.(ic.le.iend)).and.                      &
     &        ((jstr.le.jc).and.(jc.le.jend))) THEN
            my_count=0.0_r8
            my_avg=0.0_r8
            my_sum=0.0_r8
            DO jadd=-half,half
              my_sum=my_sum+                                            &
     &               f(i,j+jadd)
              my_count=my_count+1.0_r8
            END DO
            my_avg=my_sum/my_count
            c1(ic,jc)=my_avg
          END IF
        END DO
      END IF
!
!  V-type variables finer grid northern and southern boundaries.
!
      IF (gtype.eq.v2dvar) THEN
!
!  Get indices of coarser grid (ng) corresponding to the corners of
!  the finer grid (dg).
!
        ib_west=i_left(dg)
        ib_east=i_right(dg)
        jb_south=j_bottom(dg)
        jb_north=j_top(dg)
!
!  Northern boundary.
!
        j=mm(dg)+1                            ! donor finer grid
        jc=jb_north                           ! receiver coarser grid
        DO ic=ib_west,ib_east-1
          i=(ic-ib_west)*rscale+half+1
          IF (((istr.le.ic).and.(ic.le.iend)).and.                      &
     &        ((jstr.le.jc).and.(jc.le.jend))) THEN
            my_count=0.0_r8
            my_avg=0.0_r8
            my_sum=0.0_r8
            DO iadd=-half,half
              my_sum=my_sum+                                            &
     &               f(i+iadd,j)
              my_count=my_count+1.0_r8
            END DO
            my_avg=my_sum/my_count
            c1(ic,jc)=my_avg
          END IF
        END DO
!
!  Southern boundary.
!
        j=1                                   ! donor finer grid
        jc=jb_south                           ! receiver coarser grid
        DO ic=ib_west,ib_east-1
          i=(ic-ib_west)*rscale+half+1
          IF (((istr.le.ic).and.(ic.le.iend)).and.                      &
     &        ((jstr.le.jc).and.(jc.le.jend))) THEN
            my_count=0.0_r8
            my_avg=0.0_r8
            my_sum=0.0_r8
            DO iadd=-half,half
              my_sum=my_sum+                                            &
     &               f(i+iadd,j)
              my_count=my_count+1.0_r8
            END DO
            my_avg=my_sum/my_count
            c1(ic,jc)=my_avg
          END IF
        END DO
      END IF
# endif
# ifdef DISTRIBUTE
!
!  Deallocate work array.
!
      IF (allocated(f)) THEN
        deallocate (f)
      END IF
      IF (areaavg) THEN
        IF (allocated(dxf)) THEN
          deallocate (dxf)
        END IF
        IF (allocated(dyf)) THEN
          deallocate (dyf)
        END IF
      END IF
#  ifdef MASKING
      IF (allocated(fmsk)) THEN
        deallocate (fmsk)
      END IF
#  endif
# endif
!
      RETURN
      END SUBROUTINE fine2coarse2d
!
# ifdef SOLVE3D
      SUBROUTINE fine2coarse3d (ng, dg, model, tile,                    &
     &                          gtype, svname,                          &
     &                          AreaAvg, Rscale,                        &
     &                          cr, Npoints, contact,                   &
     &                          LBiF, UBiF, LBjF, UBjF, LBkF, UBkF,     &
     &                          LBiC, UBiC, LBjC, UBjC, LBkC, UBkC,     &
#  ifdef DISTRIBUTE
     &                          Adx, Ady,                               &
#  else
     &                          dxF, dyF,                               &
#  endif
     &                          pmC, pnC,                               &
#  ifdef MASKING
#   ifdef DISTRIBUTE
     &                          Amsk,                                   &
#   else
     &                          Fmsk,                                   &
#   endif
     &                          Cmsk,                                   &
#  endif
#  ifdef DISTRIBUTE
     &                          A,                                      &
#  else
     &                          F,                                      &
#  endif
     &                          C)
!
!=======================================================================
!                                                                      !
!  This routine replaces the coarse grid data inside the refinement    !
!  grid interior for a 3D state variable with its refined averaged     !
!  values: two-way nesting.                                            !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Coarser grid number (integer)                         !
!     dg         Finer grid number (integer)                           !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!     gtype      C-grid variable type (integer)                        !
!     svname     State variable name (string)                          !
!     AreaAvg    Switch for area averaging (logical)                   !
!     Rscale     Refinement grid scale (integer)                       !
!     cr         Contact region number to process (integer)            !
!     Npoints    Number of points in the contact zone (integer)        !
!     contact    Contact zone information variables (T_NGC structure)  !
!     LBiF       Finer grid, I-dimension Lower bound (integer)         !
!     UBiF       Finer grid, I-dimension Upper bound (integer)         !
!     LBjF       Finer grid, J-dimension Lower bound (integer)         !
!     UBjF       Finer grid, J-dimension Upper bound (integer)         !
!     LBkF       Finer grid, K-dimension Lower bound (integer)         !
!     UBkF       Finer grid, K-dimension Upper bound (integer)         !
!     LBiC       Coarser grid, I-dimension Lower bound (integer)       !
!     UBiC       Coarser grid, I-dimension Upper bound (integer)       !
!     LBjC       Coarser grid, J-dimension Lower bound (integer)       !
!     UBjC       Coarser grid, J-dimension Upper bound (integer)       !
!     LBkC       Coarser grid, K-dimension Lower bound (integer)       !
!     UBkC       Coarser grid, K-dimension Upper bound (integer)       !
#  ifdef DISTRIBUTE
!     Adx        Finer grid, X-grid spacing (1/pm: om_r, om_u, om_v)   !
!     Ady        Finer grid, Y-grid spacing (1/pn: on_r, on_u, on_v)   !
#  else
!     dxF        Finer grid, X-grid spacing (1/pm: om_r, om_u, om_v)   !
!     dyF        Finer grid, Y-grid spacing (1/pn: on_r, on_u, on_v)   !
#  endif
!     pmC        Coarser grid, inverse X-grid spacing (1/dx) at RHO    !
!     pnC        Coarser grid, inverse Y-grid spacing (1/dy) at RHO    !
#  ifdef MASKING
#   ifdef DISTRIBUTE
!     Amsk       Finer grid land/sea masking (2D array)                !
#   else
!     Fmsk       Finer grid land/sea masking (2D array)                !
#   endif
!     Cmsk       Coarser grid land/sea masking (2D array)              !
#  endif
#  ifdef DISTRIBUTE
!     A          Finer grid 2D data                                    !
#  else
!     F          Finer grid 2D data                                    !
#  endif
!     C          Coarser grid 3D data                                  !
!                                                                      !
!  On Output:    (mod_nesting)                                         !
!                                                                      !
!     C          Updated Coarser grid 3D data                          !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_ncparam
      USE mod_nesting
      USE mod_scalars
!
#  ifdef DISTRIBUTE
      USE distribute_mod, ONLY : mp_aggregate2d
      USE distribute_mod, ONLY : mp_aggregate3d
#  endif
      USE strings_mod,    ONLY : founderror
!
!  Imported variable declarations.
!
      logical, intent(in) :: areaavg
      integer, intent(in) :: ng, dg, model, tile
      integer, intent(in) :: gtype, cr, npoints, rscale
      integer, intent(in) :: lbif, ubif, lbjf, ubjf, lbkf, ubkf
      integer, intent(in) :: lbic, ubic, lbjc, ubjc, lbkc, ubkc
!
      character(len=*), intent(in) :: svname
!
      TYPE (t_ngc), intent(in) :: contact(:)
!
#  ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: pmc(lbic:,lbjc:)
      real(r8), intent(in) :: pnc(lbic:,lbjc:)
#   ifdef MASKING
      real(r8), intent(in) :: cmsk(lbic:,lbjc:)
#    ifdef DISTRIBUTE
      real(r8), intent(in) :: amsk(lbif:,lbjf:)
#    else
      real(r8), intent(in) :: fmsk(lbif:,lbjf:)
#    endif
#   endif
#   ifdef DISTRIBUTE
      real(r8), intent(in) :: a(lbif:,lbjf:,lbkf:)
      real(r8), intent(in) :: adx(lbif:,lbjf:)
      real(r8), intent(in) :: ady(lbif:,lbjf:)
#   else
      real(r8), intent(in) :: f(lbif:,lbjf:,lbkf:)
      real(r8), intent(in) :: dxf(lbif:,lbjf:)
      real(r8), intent(in) :: dyf(lbif:,lbjf:)
#   endif
      real(r8), intent(inout) :: c(lbic:,lbjc:,lbkc:)
#  else
      real(r8), intent(in) :: pmc(lbic:ubic,lbjc:ubjc)
      real(r8), intent(in) :: pnc(lbic:ubic,lbjc:ubjc)
#   ifdef MASKING
      real(r8), intent(in) :: cmsk(lbic:ubic,lbjc:ubjc)
#    ifdef DISTRIBUTE
      real(r8), intent(in) :: amsk(lbif:ubif,lbjf:ubjf)
#    else
      real(r8), intent(in) :: fmsk(lbif:ubif,lbjf:ubjf)
#    endif
#   endif
#   ifdef DISTRIBUTE
      real(r8), intent(in) :: a(lbif:ubif,lbjf:ubjf,lbkf:ubkf)
      real(r8), intent(in) :: adx(lbif:ubif,lbjf:ubjf)
      real(r8), intent(in) :: adx(lbif:ubif,lbjf:ubjf)
#   else
      real(r8), intent(in) :: f(lbif:ubif,lbjf:ubjf,lbkf:ubkf)
      real(r8), intent(in) :: dxf(lbif:ubif,lbjf:ubjf)
      real(r8), intent(in) :: dyf(lbif:ubif,lbjf:ubjf)
#   endif
      real(r8), intent(inout) :: c(lbic:ubic,lbjc:ubjc,lbkc:ubkc)
#  endif
!
!  Local variable declarations.
!
      integer :: iadd, ic, jadd, jc, half, i, j, k, m
      integer :: ib_east, ib_west, jb_north, jb_south
#  ifdef DISTRIBUTE
      integer :: lbi, ubi, lbj, ubj
#  endif
!
      real(r8) :: areac_inv, my_area, my_areasum, ratio
      real(r8) :: my_avg, my_count, my_sum
 
#  ifdef DISTRIBUTE
      real(r8), allocatable :: f(:,:,:)
      real(r8), allocatable :: dxf(:,:)
      real(r8), allocatable :: dyf(:,:)
#   ifdef MASKING
      real(r8), allocatable :: fmsk(:,:)
#   endif
#  endif
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", fine2coarse3d"
 
#  include "set_bounds.h"
!
!-----------------------------------------------------------------------
!  Average interior fine grid state variable data to the coarse grid
!  location. Then, replace coarse grid values with averaged data.
!-----------------------------------------------------------------------
 
#  ifdef DISTRIBUTE
!
!  Allocate global work array(s).
!
      lbi=bounds(dg)%LBi(-1)
      ubi=bounds(dg)%UBi(-1)
      lbj=bounds(dg)%LBj(-1)
      ubj=bounds(dg)%UBj(-1)
      IF (.not.allocated(f)) THEN
        allocate ( f(lbi:ubi,lbj:ubj,lbkf:ubkf) )
      END IF
      IF (areaavg) THEN
        IF (.not.allocated(dxf)) THEN
          allocate ( dxf(lbi:ubi,lbj:ubj) )
        END IF
        IF (.not.allocated(dyf)) THEN
          allocate ( dyf(lbi:ubi,lbj:ubj) )
        END IF
      END IF
#   ifdef MASKING
      IF (.not.allocated(fmsk)) THEN
        allocate ( fmsk(lbi:ubi,lbj:ubj) )
      END IF
#   endif
!
!  Gather finer grid data from all nodes in the group to build a global
!  array.
!
      CALL mp_aggregate3d (dg, model, gtype,                            &
     &                     lbif, ubif, lbjf, ubjf,                      &
     &                     lbi,  ubi,  lbj,  ubj,                       &
     &                     lbkf, ubkf,                                  &
     &                     a, f)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      IF (areaavg) THEN
        CALL mp_aggregate2d (dg, model, gtype,                          &
     &                       lbif, ubif, lbjf, ubjf,                    &
     &                       lbi,  ubi,  lbj,  ubj,                     &
     &                       adx, dxf)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL mp_aggregate2d (dg, model, gtype,                          &
     &                       lbif, ubif, lbjf, ubjf,                    &
     &                       lbi,  ubi,  lbj,  ubj,                     &
     &                       ady, dyf)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#   ifdef MASKING
!
      CALL mp_aggregate2d (dg, model, gtype,                            &
     &                     lbif, ubif, lbjf, ubjf,                      &
     &                     lbi,  ubi,  lbj,  ubj,                       &
     &                     amsk, fmsk)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#   endif
#  endif
!
!-----------------------------------------------------------------------
!  Average finer grid INTERIOR data to coarse grid according to the
!  refinement ratio.
!-----------------------------------------------------------------------
!
      half=(rscale-1)/2
      IF (areaavg) THEN               ! area averaging
        DO k=lbkc,ubkc
          DO m=1,npoints
            i=contact(cr)%Idg(m)
            j=contact(cr)%Jdg(m)
            ic=contact(cr)%Irg(m)
            jc=contact(cr)%Jrg(m)
            IF (((istr.le.ic).and.(ic.le.iend)).and.                    &
     &          ((jstr.le.jc).and.(jc.le.jend))) THEN
              my_count=0.0_r8
              my_sum=0.0_r8
              my_areasum=0.0_r8
              DO jadd=-half,half
                DO iadd=-half,half
                  my_area=dxf(i+iadd,j+jadd)*dyf(i+iadd,j+jadd)
                  my_areasum=my_areasum+my_area
#  ifdef MASKING
                  my_sum=my_sum+                                        &
     &                   f(i+iadd,j+jadd,k)*my_area*                    &
     &                   min(1.0_r8,fmsk(i+iadd,j+jadd))
                  my_count=my_count+min(1.0_r8,fmsk(i+iadd,j+jadd))
#  else
                  my_sum=my_sum+                                        &
     &                   f(i+iadd,j+jadd,k)*my_area
#  endif
                END DO
              END DO
              SELECT CASE (gtype)            ! coarse grid inverse area
                CASE (r3dvar)
                  areac_inv=pmc(ic,jc)*pnc(ic,jc)
                CASE (u3dvar)
                  areac_inv=0.25_r8*(pmc(ic-1,jc)+pmc(ic,jc))*          &
     &                              (pnc(ic-1,jc)+pnc(ic,jc))
                CASE (v3dvar)
                  areac_inv=0.25_r8*(pmc(ic,jc-1)+pmc(ic,jc))*          &
     &                              (pnc(ic,jc-1)+pnc(ic,jc))
                CASE DEFAULT
                  areac_inv=pmc(ic,jc)*pnc(ic,jc)
              END SELECT
              ratio=my_areasum*areac_inv     ! for debugging purposes
              my_avg=my_sum*areac_inv
#  ifdef MASKING
              IF (my_count.gt.0.0_r8) THEN
                my_avg=my_avg*rscale*rscale/my_count
              END IF
              my_avg=my_avg*cmsk(ic,jc)
#  endif
              c(ic,jc,k)=my_avg
            END IF
          END DO
        END DO
      ELSE                            ! simple averaging
        DO k=lbkc,ubkc
          DO m=1,npoints
            i=contact(cr)%Idg(m)
            j=contact(cr)%Jdg(m)
            ic=contact(cr)%Irg(m)
            jc=contact(cr)%Jrg(m)
            IF (((istr.le.ic).and.(ic.le.iend)).and.                    &
     &          ((jstr.le.jc).and.(jc.le.jend))) THEN
              my_count=0.0_r8
              my_avg=0.0_r8
              my_sum=0.0_r8
              DO jadd=-half,half
                DO iadd=-half,half
#  ifdef MASKING
                  my_sum=my_sum+                                        &
     &                   f(i+iadd,j+jadd,k)*fmsk(i+iadd,j+jadd)
                  my_count=my_count+min(1.0_r8,fmsk(i+iadd,j+jadd))
#  else
                  my_sum=my_sum+                                        &
     &                   f(i+iadd,j+jadd,k)
                  my_count=my_count+1.0_r8
#  endif
                END DO
              END DO
              IF (my_count.gt.0.0_r8) my_avg=my_sum/my_count
#  ifdef MASKING
              my_avg=my_avg*cmsk(ic,jc)
#  endif
              c(ic,jc,k)=my_avg
            END IF
          END DO
        END DO
      END IF
#  ifdef REFINE_BOUNDARY
!
!-----------------------------------------------------------------------
!  Average finer grid BOUNDARY data to coarse grid U-type and V-type
!  variables according to the refinement ratio.
!-----------------------------------------------------------------------
!
!  U-type variables finer grid eastern and western boundaries.
!
      IF (gtype.eq.u3dvar) THEN
!
!  Get indices of coarser grid (ng) corresponding to the corners of
!  the finer grid (dg).
!
        ib_west=i_left(dg)
        ib_east=i_right(dg)
        jb_south=j_bottom(dg)
        jb_north=j_top(dg)
!
!  Eastern boundary.
!
        DO k=lbkc,ubkc
          i=lm(dg)+1                            ! donor finer grid
          ic=ib_east                            ! receiver coarser grid
          DO jc=jb_south,jb_north-1
            j=(jc-jb_south)*rscale+half+1
            IF (((istr.le.ic).and.(ic.le.iend)).and.                    &
     &          ((jstr.le.jc).and.(jc.le.jend))) THEN
              my_count=0.0_r8
              my_avg=0.0_r8
              my_sum=0.0_r8
              DO jadd=-half,half
                my_sum=my_sum+                                          &
     &                 f(i,j+jadd,k)
                my_count=my_count+1.0_r8
              END DO
              my_avg=my_sum/my_count
              c(ic,jc,k)=my_avg
            END IF
          END DO
        END DO
!
!  Western boundary.
!
        DO k=lbkc,ubkc
          i=1                                   ! donor finer grid
          ic=ib_west                            ! receiver coarser grid
          DO jc=jb_south,jb_north-1
            j=(jc-jb_south)*rscale+half+1
            IF (((istr.le.ic).and.(ic.le.iend)).and.                    &
     &          ((jstr.le.jc).and.(jc.le.jend))) THEN
              my_count=0.0_r8
              my_avg=0.0_r8
              my_sum=0.0_r8
              DO jadd=-half,half
                my_sum=my_sum+                                          &
     &                 f(i,j+jadd,k)
                my_count=my_count+1.0_r8
              END DO
              my_avg=my_sum/my_count
              c(ic,jc,k)=my_avg
            END IF
          END DO
        END DO
      END IF
!
!  V-type variables finer grid northern and southern boundaries.
!
      IF (gtype.eq.v3dvar) THEN
!
!  Get indices of coarser grid (ng) corresponding to the corners of
!  the finer grid (dg).
!
        ib_west=i_left(dg)
        ib_east=i_right(dg)
        jb_south=j_bottom(dg)
        jb_north=j_top(dg)
!
!  Northern boundary.
!
        DO k=lbkc,ubkc
          j=mm(dg)+1                            ! donor finer grid
          jc=jb_north                           ! receiver coarser grid
          DO ic=ib_west,ib_east-1
            i=(ic-ib_west)*rscale+half+1
            IF (((istr.le.ic).and.(ic.le.iend)).and.                    &
     &          ((jstr.le.jc).and.(jc.le.jend))) THEN
              my_count=0.0_r8
              my_avg=0.0_r8
              my_sum=0.0_r8
              DO iadd=-half,half
                my_sum=my_sum+                                          &
     &                 f(i+iadd,j,k)
                my_count=my_count+1.0_r8
              END DO
              my_avg=my_sum/my_count
              c(ic,jc,k)=my_avg
            END IF
          END DO
        END DO
!
!  Southern boundary.
!
        DO k=lbkc,ubkc
          j=1                                   ! donor finer grid
          jc=jb_south                           ! receiver coarser grid
          DO ic=ib_west,ib_east-1
            i=(ic-ib_west)*rscale+half+1
            IF (((istr.le.ic).and.(ic.le.iend)).and.                    &
     &          ((jstr.le.jc).and.(jc.le.jend))) THEN
              my_count=0.0_r8
              my_avg=0.0_r8
              my_sum=0.0_r8
              DO iadd=-half,half
                my_sum=my_sum+                                          &
     &                 f(i+iadd,j,k)
                my_count=my_count+1.0_r8
              END DO
              my_avg=my_sum/my_count
              c(ic,jc,k)=my_avg
            END IF
          END DO
        END DO
      END IF
#  endif
#  ifdef DISTRIBUTE
!
!  Deallocate work array.
!
      IF (allocated(f)) THEN
        deallocate (f)
      END IF
      IF (areaavg) THEN
        IF (allocated(dxf)) THEN
          deallocate (dxf)
        END IF
        IF (allocated(dyf)) THEN
          deallocate (dyf)
        END IF
      END IF
#   ifdef MASKING
      IF (allocated(fmsk)) THEN
        deallocate (fmsk)
      END IF
#   endif
#  endif
!
      RETURN
      END SUBROUTINE fine2coarse3d
# endif
!
      SUBROUTINE get_contact2d (dg, model, tile,                        &
     &                          gtype, svname,                          &
     &                          cr, Npoints, contact,                   &
     &                          LBi, UBi, LBj, UBj,                     &
     &                          Ad, Ac)
!
!=======================================================================
!                                                                      !
!  This routine gets the donor grid data (Ac) necessary  to process    !
!  the contact points for a 2D state variable (Ad). It extracts the    !
!  donor cell points containing each contact point, Ac(1:4,:).         !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     dg         Donor grid number (integer)                           !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!     gtype      C-grid variable type (integer)                        !
!     svname     State variable name (string)                          !
!     cr         Contact region number to process (integer)            !
!     Npoints    Number of points in the contact region (integer)      !
!     contact    Contact region information variables (T_NGC structure)!
!     LBi        Donor grid, I-dimension Lower bound (integer)         !
!     UBi        Donor grid, I-dimension Upper bound (integer)         !
!     LBj        Donor grid, J-dimension Lower bound (integer)         !
!     UBj        Donor grid, J-dimension Upper bound (integer)         !
!     Ad         Donor grid data (2D array)                            !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     Ac         2D state variable contact point data                  !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_ncparam
      USE mod_nesting
 
# ifdef DISTRIBUTE
!
      USE distribute_mod, ONLY : mp_assemble
# endif
!
!  Imported variable declarations.
!
      integer, intent(in) :: dg, model, tile
      integer, intent(in) :: gtype, cr, npoints
      integer, intent(in) :: lbi, ubi, lbj, ubj
!
      character(len=*), intent(in) :: svname
!
      TYPE (t_ngc), intent(in) :: contact(:)
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: ad(lbi:,lbj:)
      real(r8), intent(inout) :: ac(:,:)
# else
      real(r8), intent(in) :: ad(lbi:ubi,lbj:ubj)
      real(r8), intent(inout) :: ac(npoints,4)
# endif
!
!  Local variable declarations.
!
      integer :: i, ip1, j, jp1, m
      integer :: imin, imax, jmin, jmax
      integer :: istr, iend, jstr, jend
# ifdef DISTRIBUTE
      integer :: npts
# endif
 
      real(r8), parameter :: aspv = 0.0_r8
!
!-----------------------------------------------------------------------
!  Initialize.
!-----------------------------------------------------------------------
 
# ifdef DISTRIBUTE
!
!  Initialize contact points array to special value to facilite
!  distribute-memory data collection from all nodes.
!
      npts=4*npoints
      DO m=1,npoints
        ac(1,m)=aspv
        ac(2,m)=aspv
        ac(3,m)=aspv
        ac(4,m)=aspv
      END DO
# endif
!
!  Set starting and ending tile indices for the donor grids.
!
      SELECT CASE (gtype)
        CASE (r2dvar)
          imin=bounds(dg) % IstrT(-1)    ! full RHO-grid range
          imax=bounds(dg) % IendT(-1)
          jmin=bounds(dg) % JstrT(-1)
          jmax=bounds(dg) % JendT(-1)
!
          istr=bounds(dg) % IstrT(tile)  ! domain partition range
          iend=bounds(dg) % IendT(tile)
          jstr=bounds(dg) % JstrT(tile)
          jend=bounds(dg) % JendT(tile)
        CASE (u2dvar)
          imin=bounds(dg) % IstrP(-1)    ! full U-grid range
          imax=bounds(dg) % IendT(-1)
          jmin=bounds(dg) % JstrT(-1)
          jmax=bounds(dg) % JendT(-1)
!
          istr=bounds(dg) % IstrP(tile)  ! domain partition range
          iend=bounds(dg) % IendT(tile)
          jstr=bounds(dg) % JstrT(tile)
          jend=bounds(dg) % JendT(tile)
        CASE (v2dvar)
          imin=bounds(dg) % IstrT(-1)    ! full V-grid range
          imax=bounds(dg) % IendT(-1)
          jmin=bounds(dg) % JstrP(-1)
          jmax=bounds(dg) % JendT(-1)
!
          istr=bounds(dg) % IstrT(tile)  ! domain partition range
          iend=bounds(dg) % IendT(tile)
          jstr=bounds(dg) % JstrP(tile)
          jend=bounds(dg) % JendT(tile)
      END SELECT
!
!-----------------------------------------------------------------------
!  Extract donor grid data at contact points.
!-----------------------------------------------------------------------
!
!  Notice that the indices i+1 and j+1 are bounded the maximum values
!  of the grid. This implies that contact point lies on the grid
!  boundary.
!
      DO m=1,npoints
        i=contact(cr)%Idg(m)
        j=contact(cr)%Jdg(m)
        ip1=min(i+1,imax)
        jp1=min(j+1,jmax)
        IF (((istr.le.i).and.(i.le.iend)).and.                          &
     &      ((jstr.le.j).and.(j.le.jend))) THEN
          ac(1,m)=ad(i  ,j  )
          ac(2,m)=ad(ip1,j  )
          ac(3,m)=ad(ip1,jp1)
          ac(4,m)=ad(i  ,jp1)
        END IF
      END DO
 
# ifdef DISTRIBUTE
!
!  Gather and broadcast data from all nodes.
!
      CALL mp_assemble (dg, model, npts, aspv, ac)
# endif
!
      RETURN
      END SUBROUTINE get_contact2d
 
# ifdef SOLVE3D
!
      SUBROUTINE get_contact3d (dg, model, tile,                        &
     &                          gtype, svname,                          &
     &                          cr, Npoints, contact,                   &
     &                          LBi, UBi, LBj, UBj, LBk, UBk,           &
     &                          Ad, Ac)
!
!=======================================================================
!                                                                      !
!  This routine gets the donor grid data (Ac) necessary  to process    !
!  the contact points for a 3D state variable (Ad). It extracts the    !
!  donor cell points containing each contact point, Ac(1:4,k,:).       !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     dg         Donor grid number (integer)                           !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!     gtype      C-grid variable type (integer)                        !
!     svname     State variable name (string)                          !
!     cr         Contact region number to process (integer)            !
!     Npoints    Number of points in the contact region (integer)      !
!     contact    Contact region information variables (T_NGC structure)!
!     LBi        Donor grid, I-dimension Lower bound (integer)         !
!     UBi        Donor grid, I-dimension Upper bound (integer)         !
!     LBj        Donor grid, J-dimension Lower bound (integer)         !
!     UBj        Donor grid, J-dimension Upper bound (integer)         !
!     Ad         Donor grid data (3D array)                            !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     Ac         3D state variable contact point data                  !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_ncparam
      USE mod_nesting
 
#  ifdef DISTRIBUTE
!
      USE distribute_mod, ONLY : mp_assemble
#  endif
!
!  Imported variable declarations.
!
      integer, intent(in) :: dg, model, tile
      integer, intent(in) :: gtype, cr, npoints
      integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk
!
      character(len=*), intent(in) :: svname
!
      TYPE (t_ngc), intent(in) :: contact(:)
!
#  ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: ad(lbi:,lbj:,lbk:)
      real(r8), intent(inout) :: ac(:,lbk:,:)
#  else
      real(r8), intent(in) :: ad(lbi:ubi,lbj:ubj,lbk:ubk)
      real(r8), intent(inout) :: ac(4,lbk:ubk,npoints)
#  endif
!
!  Local variable declarations.
!
      integer :: i, ip1, j, jp1, k, m
      integer :: imin, imax, jmin, jmax
      integer :: istr, iend, jstr, jend
#  ifdef DISTRIBUTE
      integer :: npts
#  endif
 
      real(r8), parameter :: aspv = 0.0_r8
!
!-----------------------------------------------------------------------
!  Initialize.
!-----------------------------------------------------------------------
 
#  ifdef DISTRIBUTE
!
!  Initialize contact points array to special value to facilite
!  distribute-memory data collection from all nodes.
!
      npts=4*(ubk-lbk+1)*npoints
      DO k=lbk,ubk
        DO m=1,npoints
          ac(1,k,m)=aspv
          ac(2,k,m)=aspv
          ac(3,k,m)=aspv
          ac(4,k,m)=aspv
        END DO
      END DO
#  endif
!
!  Set starting and ending tile indices for the donor grid.
!
      SELECT CASE (gtype)
        CASE (r3dvar)
          imin=bounds(dg) % IstrT(-1)    ! full RHO-grid range
          imax=bounds(dg) % IendT(-1)
          jmin=bounds(dg) % JstrT(-1)
          jmax=bounds(dg) % JendT(-1)
!
          istr=bounds(dg) % IstrT(tile)  ! domain partition range
          iend=bounds(dg) % IendT(tile)
          jstr=bounds(dg) % JstrT(tile)
          jend=bounds(dg) % JendT(tile)
        CASE (u3dvar)
          imin=bounds(dg) % IstrP(-1)    ! full U-grid range
          imax=bounds(dg) % IendT(-1)
          jmin=bounds(dg) % JstrT(-1)
          jmax=bounds(dg) % JendT(-1)
!
          istr=bounds(dg) % IstrP(tile)  ! domain partition range
          iend=bounds(dg) % IendT(tile)
          jstr=bounds(dg) % JstrT(tile)
          jend=bounds(dg) % JendT(tile)
        CASE (v3dvar)
          imin=bounds(dg) % IstrT(-1)    ! full V-grid range
          imax=bounds(dg) % IendT(-1)
          jmin=bounds(dg) % JstrP(-1)
          jmax=bounds(dg) % JendT(-1)
!
          istr=bounds(dg) % IstrT(tile)  ! domain partition range
          iend=bounds(dg) % IendT(tile)
          jstr=bounds(dg) % JstrP(tile)
          jend=bounds(dg) % JendT(tile)
      END SELECT
!
!-----------------------------------------------------------------------
!  Extract donor grid data at contact points.
!-----------------------------------------------------------------------
!
!  Notice that the indices i+1 and j+1 are bounded the maximum values
!  of the grid. This implies that contact point lies on the grid
!  boundary.
!
      DO k=lbk,ubk
        DO m=1,npoints
          i=contact(cr)%Idg(m)
          j=contact(cr)%Jdg(m)
          ip1=min(i+1,imax)
          jp1=min(j+1,jmax)
          IF (((istr.le.i).and.(i.le.iend)).and.                        &
     &        ((jstr.le.j).and.(j.le.jend))) THEN
            ac(1,k,m)=ad(i  ,j  ,k)
            ac(2,k,m)=ad(ip1,j  ,k)
            ac(3,k,m)=ad(ip1,jp1,k)
            ac(4,k,m)=ad(i  ,jp1,k)
          END IF
        END DO
      END DO
 
#  ifdef DISTRIBUTE
!
!  Gather and broadcast data from all nodes.
!
      CALL mp_assemble (dg, model, npts, aspv, ac(:,lbk:,:))
#  endif
!
      RETURN
      END SUBROUTINE get_contact3d
# endif
!
      SUBROUTINE get_persisted2d (dg, rg, model, tile,                  &
     &                            gtype, svname,                        &
     &                            cr, Npoints, contact,                 &
     &                            LBi, UBi, LBj, UBj,                   &
     &                            Ad, Ac)
!
!=======================================================================
!                                                                      !
!  This routine gets the donor grid data (Ac) necessary  to process    !
!  the contact points for a 2D flux variable (Ad). It extracts the     !
!  donor cell points containing each contact point, Ac(1:4,:).         !
!                                                                      !
!  This routine is different that 'get_contact2d'.  It is used in      !
!  refinement to impose the appropriate coarser grid flux to insure    !
!  volume and mass conservation.  The value of the coarse grid cell    !
!  is presisted over the refined grid points along its physical        !
!  boundary.  This will facilitate that the sum of all the refined     !
!  grid point is the same as that of the coarse grid containing such   !
!  points.  The spatial interpolation as set in 'get_contact2d' will   !
!  not conserve volume and mass.                                       !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     dg         Donor grid number (integer)                           !
!     rg         Receiver grid number (integer)                        !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!     gtype      C-grid variable type (integer)                        !
!     svname     State variable name (string)                          !
!     cr         Contact region number to process (integer)            !
!     Npoints    Number of points in the contact region (integer)      !
!     contact    Contact region information variables (T_NGC structure)!
!     LBi        Donor grid, I-dimension Lower bound (integer)         !
!     UBi        Donor grid, I-dimension Upper bound (integer)         !
!     LBj        Donor grid, J-dimension Lower bound (integer)         !
!     UBj        Donor grid, J-dimension Upper bound (integer)         !
!     Ad         Donor grid data (2D array)                            !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     Ac         2D flux variable contact point data                   !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_ncparam
      USE mod_nesting
      USE mod_scalars
!
# ifdef DISTRIBUTE
      USE distribute_mod, ONLY : mp_assemble
# endif
      USE strings_mod,    ONLY : founderror
!
!  Imported variable declarations.
!
      integer, intent(in) :: dg, rg, model, tile
      integer, intent(in) :: gtype, cr, npoints
      integer, intent(in) :: lbi, ubi, lbj, ubj
!
      character(len=*), intent(in) :: svname
!
      TYPE (t_ngc), intent(in) :: contact(:)
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: ad(lbi:,lbj:)
      real(r8), intent(inout) :: ac(:,:)
# else
      real(r8), intent(in) :: ad(lbi:ubi,lbj:ubj)
      real(r8), intent(inout) :: ac(npoints,4)
# endif
!
!  Local variable declarations.
!
      integer :: idg, ip1, irg, jdg, jp1, jrg
      integer :: imin, imax, jmin, jmax
      integer :: istr, iend, jstr, jend
      integer :: i, i_add, j, j_add, m, m_add
# ifdef DISTRIBUTE
      integer :: npts
# endif
!
      real(r8), parameter :: aspv = 0.0_r8
      real(r8):: rscale
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", get_persisted2d"
!
!-----------------------------------------------------------------------
!  Initialize.
!-----------------------------------------------------------------------
 
# ifdef DISTRIBUTE
!
!  Initialize contact points array to special value to facilite
!  distribute-memory data collection from all nodes.
!
      npts=4*npoints
      DO m=1,npoints
        ac(1,m)=aspv
        ac(2,m)=aspv
        ac(3,m)=aspv
        ac(4,m)=aspv
      END DO
# endif
!
!  Set starting and ending tile indices for the donor grids.
!
      SELECT CASE (gtype)
        CASE (r2dvar)
          imin=bounds(dg) % IstrT(-1)    ! full RHO-grid range
          imax=bounds(dg) % IendT(-1)
          jmin=bounds(dg) % JstrT(-1)
          jmax=bounds(dg) % JendT(-1)
!
          istr=bounds(dg) % IstrT(tile)  ! domain partition range
          iend=bounds(dg) % IendT(tile)
          jstr=bounds(dg) % JstrT(tile)
          jend=bounds(dg) % JendT(tile)
!
          m_add=nstrr(cr)-1
        CASE (u2dvar)
          imin=bounds(dg) % IstrP(-1)    ! full U-grid range
          imax=bounds(dg) % IendT(-1)
          jmin=bounds(dg) % JstrT(-1)
          jmax=bounds(dg) % JendT(-1)
!
          istr=bounds(dg) % IstrP(tile)  ! domain partition range
          iend=bounds(dg) % IendT(tile)
          jstr=bounds(dg) % JstrT(tile)
          jend=bounds(dg) % JendT(tile)
!
          m_add=nstru(cr)-1
        CASE (v2dvar)
          imin=bounds(dg) % IstrT(-1)    ! full V-grid range
          imax=bounds(dg) % IendT(-1)
          jmin=bounds(dg) % JstrP(-1)
          jmax=bounds(dg) % JendT(-1)
!
          istr=bounds(dg) % IstrT(tile)  ! domain partition range
          iend=bounds(dg) % IendT(tile)
          jstr=bounds(dg) % JstrP(tile)
          jend=bounds(dg) % JendT(tile)
!
          m_add=nstrv(cr)-1
      END SELECT
!
!-----------------------------------------------------------------------
!  Extract donor grid data at contact points.
!-----------------------------------------------------------------------
!
!  Notice that the indices i+1 and j+1 are bounded the maximum values
!  of the grid. This implies that contact point lies on the grid
!  boundary.
!
      rscale=1.0_r8/real(refinescale(rg))
      DO m=1,npoints
        idg=contact(cr)%Idg(m)
        jdg=contact(cr)%Jdg(m)
        irg=contact(cr)%Irg(m)
        jrg=contact(cr)%Jrg(m)
        ip1=min(idg+1,imax)
        jp1=min(jdg+1,jmax)
        IF (((istr.le.idg).and.(idg.le.iend)).and.                      &
     &      ((jstr.le.jdg).and.(jdg.le.jend))) THEN
          IF (on_boundary(m+m_add).gt.0) THEN
            IF ((on_boundary(m+m_add).eq.1).or.                         &
     &          (on_boundary(m+m_add).eq.3)) THEN       ! western and
              j_add=int(real(jrg-1,r8)*rscale)          ! eastern edges
              j=j_bottom(rg)+j_add
              ac(1,m)=ad(idg,j)
              ac(2,m)=ad(idg,j)
              ac(3,m)=ad(idg,j)
              ac(4,m)=ad(idg,j)
            ELSE IF ((on_boundary(m+m_add).eq.2).or.                    &
     &               (on_boundary(m+m_add).eq.4)) THEN  ! southern and
              i_add=int(real(irg-1,r8)*rscale)          ! northern edges
              i=i_left(rg)+i_add
              ac(1,m)=ad(i,jdg)
              ac(2,m)=ad(i,jdg)
              ac(3,m)=ad(i,jdg)
              ac(4,m)=ad(i,jdg)
            END IF
          ELSE
            ac(1,m)=ad(idg,jdg)               ! contact point is not at
            ac(2,m)=ad(ip1,jdg)               ! physical boundary, just
            ac(3,m)=ad(ip1,jp1)               ! set values for spatial
            ac(4,m)=ad(idg,jp1)               ! interpolation (not used)
          END IF
        END IF
      END DO
 
# ifdef DISTRIBUTE
!
!  Gather and broadcast data from all nodes.
!
      CALL mp_assemble (dg, model, npts, aspv, ac)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
!
      RETURN
      END SUBROUTINE get_persisted2d
!
      SUBROUTINE put_contact2d (rg, model, tile,                        &
     &                          gtype, svname,                          &
     &                          cr, Npoints, contact,                   &
     &                          LBi, UBi, LBj, UBj,                     &
# ifdef MASKING
     &                          Amask,                                  &
# endif
     &                          Ac, Ar)
!
!=======================================================================
!                                                                      !
!  This routine uses extracted donor grid data (Ac) to spatially       !
!  interpolate a 2D state variable  at the receiver grid contact       !
!  points.  If the donor and receiver grids are coincident,  the       !
!  Lweight(1,:) is unity and Lweight(2:4,:) are zero.                  !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     rg         Receiver grid number (integer)                        !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!     gtype      C-grid variable type (integer)                        !
!     svname     State variable name (string)                          !
!     cr         Contact region number to process (integer)            !
!     Npoints    Number of points in the contact region (integer)      !
!     contact    Contact region information variables (T_NGC structure)!
!     LBi        Receiver grid, I-dimension Lower bound (integer)      !
!     UBi        Receiver grid, I-dimension Upper bound (integer)      !
!     LBj        Receiver grid, J-dimension Lower bound (integer)      !
!     UBj        Receiver grid, J-dimension Upper bound (integer)      !
# ifdef MASKING
!     Amask      Receiver grid land/sea masking                        !
# endif
!     Ac         Contact point data extracted from donor grid          !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     Ar         Updated receiver grid 2D state array                  !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_ncparam
      USE mod_nesting
!
!  Imported variable declarations.
!
      integer, intent(in) :: rg, model, tile
      integer, intent(in) :: gtype, cr, Npoints
      integer, intent(in) :: LBi, UBi, LBj, UBj
!
      character(len=*), intent(in) :: svname
!
      TYPE (T_NGC), intent(in) :: contact(:)
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: Ac(:,:)
#  ifdef MASKING
      real(r8), intent(in) :: Amask(LBi:,LBj:)
#  endif
      real(r8), intent(inout) :: Ar(LBi:,LBj:)
# else
      real(r8), intent(in) :: Ac(4,Npoints)
#  ifdef MASKING
      real(r8), intent(in) :: Amask(LBi:UBi,LBj:UBj)
#  endif
      real(r8), intent(inout) :: Ar(LBi:UBi,LBj:UBj)
# endif
!
!  Local variable declarations.
!
      integer :: i, j, m
      integer :: Istr, Iend, Jstr, Jend
!
!-----------------------------------------------------------------------
!  Interpolate 2D data from donor grid to receiver grid contact points.
!-----------------------------------------------------------------------
!
!  Set starting and ending tile indices for the receiver grid.
!
      SELECT CASE (gtype)
        CASE (r2dvar)
          istr=bounds(rg) % IstrT(tile)
          iend=bounds(rg) % IendT(tile)
          jstr=bounds(rg) % JstrT(tile)
          jend=bounds(rg) % JendT(tile)
        CASE (u2dvar)
          istr=bounds(rg) % IstrP(tile)
          iend=bounds(rg) % IendT(tile)
          jstr=bounds(rg) % JstrT(tile)
          jend=bounds(rg) % JendT(tile)
        CASE (v2dvar)
          istr=bounds(rg) % IstrT(tile)
          iend=bounds(rg) % IendT(tile)
          jstr=bounds(rg) % JstrP(tile)
          jend=bounds(rg) % JendT(tile)
      END SELECT
!
!  Interpolate.
!
      DO m=1,npoints
        i=contact(cr)%Irg(m)
        j=contact(cr)%Jrg(m)
        IF (((istr.le.i).and.(i.le.iend)).and.                          &
     &      ((jstr.le.j).and.(j.le.jend))) THEN
          ar(i,j)=contact(cr)%Lweight(1,m)*ac(1,m)+                     &
     &            contact(cr)%Lweight(2,m)*ac(2,m)+                     &
     &            contact(cr)%Lweight(3,m)*ac(3,m)+                     &
     &            contact(cr)%Lweight(4,m)*ac(4,m)
# ifdef MASKING
          ar(i,j)=ar(i,j)*amask(i,j)
# endif
        END IF
      END DO
!
      RETURN
      END SUBROUTINE put_contact2d
 
# ifdef SOLVE3D
!
      SUBROUTINE put_contact3d (rg, model, tile,                        &
     &                          gtype, svname,                          &
     &                          cr, Npoints, contact,                   &
     &                          LBi, UBi, LBj, UBj, LBk, UBk,           &
#  ifdef MASKING
     &                          Amask,                                  &
#  endif
     &                          Ac, Ar)
!
!=======================================================================
!                                                                      !
!  This routine uses extracted donor grid data (Ac) to spatially       !
!  interpolate a 3D state variable  at the receiver grid contact       !
!  points.  If the donor and receiver grids  are concident,  the       !
!  Lweight(1,:) is unity and Lweight(2:4,:) are zero.                  !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     rg         Receiver grid number (integer)                        !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!     gtype      C-grid variable type (integer)                        !
!     svname     State variable name (string)                          !
!     cr         Contact region number to process (integer)            !
!     Npoints    Number of points in the contact region (integer)      !
!     contact    Contact region information variables (T_NGC structure)!
!     LBi        Receiver grid, I-dimension Lower bound (integer)      !
!     UBi        Receiver grid, I-dimension Upper bound (integer)      !
!     LBj        Receiver grid, J-dimension Lower bound (integer)      !
!     UBj        Receiver grid, J-dimension Upper bound (integer)      !
!     LBk        Receiver grid, K-dimension Lower bound (integer)      !
!     UBk        Receiver grid, K-dimension Upper bound (integer)      !
#  ifdef MASKING
!     Amask      Receiver grid land/sea masking                        !
#  endif
!     Ac         Contact point data extracted from donor grid          !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     Ar         Updated receiver grid 3D state array                  !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_ncparam
      USE mod_nesting
!
!  Imported variable declarations.
!
      integer, intent(in) :: rg, model, tile
      integer, intent(in) :: gtype, cr, Npoints
      integer, intent(in) :: LBi, UBi, LBj, UBj, LBk, UBk
!
      character(len=*), intent(in) :: svname
!
      TYPE (T_NGC), intent(in) :: contact(:)
!
#  ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: Ac(:,:,:)
#   ifdef MASKING
      real(r8), intent(in) :: Amask(LBi:,LBj:)
#   endif
      real(r8), intent(inout) :: Ar(LBi:,LBj:,LBk:)
#  else
      real(r8), intent(in) :: Ac(Npoints,LBk:UBk,4)
#   ifdef MASKING
      real(r8), intent(in) :: Amask(LBi:UBi,LBj:UBj)
#   endif
      real(r8), intent(inout) :: Ar(LBi:UBi,LBj:UBj,LBk:UBk)
#  endif
!
!  Local variable declarations.
!
      integer :: i, j, k, kdg, kdgm1, m
      integer :: Istr, Iend, Jstr, Jend, Kmin
 
      real(r8), dimension(8) :: cff
!
!-----------------------------------------------------------------------
!  Interpolate 3D data from donor grid to receiver grid contact points.
!-----------------------------------------------------------------------
!
!  Set starting and ending tile indices for the receiver grid.
!
      SELECT CASE (gtype)
        CASE (r3dvar)
          istr=bounds(rg) % IstrT(tile)
          iend=bounds(rg) % IendT(tile)
          jstr=bounds(rg) % JstrT(tile)
          jend=bounds(rg) % JendT(tile)
          kmin=1
        CASE (u3dvar)
          istr=bounds(rg) % IstrP(tile)
          iend=bounds(rg) % IendT(tile)
          jstr=bounds(rg) % JstrT(tile)
          jend=bounds(rg) % JendT(tile)
          kmin=1
        CASE (v3dvar)
          istr=bounds(rg) % IstrT(tile)
          iend=bounds(rg) % IendT(tile)
          jstr=bounds(rg) % JstrP(tile)
          jend=bounds(rg) % JendT(tile)
          kmin=1
        CASE (w3dvar)
          istr=bounds(rg) % IstrT(tile)
          iend=bounds(rg) % IendT(tile)
          jstr=bounds(rg) % JstrT(tile)
          jend=bounds(rg) % JendT(tile)
          kmin=0
      END SELECT
!
!  Interpolate.
!
      DO k=lbk,ubk
        DO m=1,npoints
          i=contact(cr)%Irg(m)
          j=contact(cr)%Jrg(m)
          kdg=contact(cr)%Kdg(k,m)
          kdgm1=max(kdg-1,kmin)
          IF (((istr.le.i).and.(i.le.iend)).and.                        &
     &        ((jstr.le.j).and.(j.le.jend))) THEN
            cff(1)=contact(cr)%Lweight(1,m)*contact(cr)%Vweight(1,k,m)
            cff(2)=contact(cr)%Lweight(2,m)*contact(cr)%Vweight(1,k,m)
            cff(3)=contact(cr)%Lweight(3,m)*contact(cr)%Vweight(1,k,m)
            cff(4)=contact(cr)%Lweight(4,m)*contact(cr)%Vweight(1,k,m)
            cff(5)=contact(cr)%Lweight(1,m)*contact(cr)%Vweight(2,k,m)
            cff(6)=contact(cr)%Lweight(2,m)*contact(cr)%Vweight(2,k,m)
            cff(7)=contact(cr)%Lweight(3,m)*contact(cr)%Vweight(2,k,m)
            cff(8)=contact(cr)%Lweight(4,m)*contact(cr)%Vweight(2,k,m)
            ar(i,j,k)=cff(1)*ac(1,kdgm1,m)+                             &
     &                cff(2)*ac(2,kdgm1,m)+                             &
     &                cff(3)*ac(3,kdgm1,m)+                             &
     &                cff(4)*ac(4,kdgm1,m)+                             &
     &                cff(5)*ac(1,kdg  ,m)+                             &
     &                cff(6)*ac(2,kdg  ,m)+                             &
     &                cff(7)*ac(3,kdg  ,m)+                             &
     &                cff(8)*ac(4,kdg  ,m)
#  ifdef MASKING
            ar(i,j,k)=ar(i,j,k)*amask(i,j)
#  endif
          END IF
        END DO
      END DO
!
      RETURN
      END SUBROUTINE put_contact3d
# endif
!
      SUBROUTINE put_refine2d (ng, dg, cr, model, tile, LputFsur,       &
     &                         LBi, UBi, LBj, UBj)
!
!=======================================================================
!                                                                      !
!  This routine interpolates (space, time) refinement grid 2D state    !
!  variables contact points using data from the donor grid.  Notice    !
!  that because of shared-memory parallelism,  the  free-surface is    !
!  processed first and in a different parallel region.                 !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Refinement (receiver) grid number (integer)           !
!     dg         Donor grid number (integer)                           !
!     cr         Contact region number to process (integer)            !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!     LputFsur   Switch to process or not free-surface (logical)       !
!     LBi        Receiver grid, I-dimension Lower bound (integer)      !
!     UBi        Receiver grid, I-dimension Upper bound (integer)      !
!     LBj        Receiver grid, J-dimension Lower bound (integer)      !
!     UBj        Receiver grid, J-dimension Upper bound (integer)      !
!                                                                      !
!  On Output:    OCEAN(ng) structure                                   !
!                                                                      !
!     zeta       Updated free-surface                                  !
!     ubar       Updated 2D momentum in the XI-direction               !
!     vbar       Updated 2D momentum in the ETA-direction              !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_coupling
      USE mod_grid
      USE mod_nesting
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
      USE mod_iounits
 
# ifdef DISTRIBUTE
!
      USE distribute_mod,  ONLY : mp_assemble
      USE mp_exchange_mod, ONLY : mp_exchange2d
# endif
      USE strings_mod,     ONLY : founderror
!
!  Imported variable declarations.
!
      logical, intent(in) :: lputfsur
      integer, intent(in) :: ng, dg, cr, model, tile
      integer, intent(in) :: lbi, ubi, lbj, ubj
!
!  Local variable declarations.
!
      logical :: uboundary, vboundary
!
# ifdef DISTRIBUTE
      integer :: ilb, iub, jlb, jub, nptssn, nptswe, my_tile
# endif
      integer :: nsub, i, irec, j, kindex, m, tnew, told
      integer :: idg, jdg
!
# ifdef DISTRIBUTE
      real(r8), parameter :: spv = 0.0_r8
# endif
      real(dp) :: wnew, wold, secscale, fac
      real(r8) :: cff, cff1
      real(r8) :: my_value
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", put_refined2d"
 
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
!  Interpolate (space, time) refinement grid contact points for 2D state
!  variables from donor grid.
!-----------------------------------------------------------------------
 
# ifdef DISTRIBUTE
!
!  Set global size of boundary edges.
!
      IF (.not.lputfsur) THEN
        my_tile=-1
        ilb=bounds(ng)%LBi(my_tile)
        iub=bounds(ng)%UBi(my_tile)
        jlb=bounds(ng)%LBj(my_tile)
        jub=bounds(ng)%UBj(my_tile)
        nptswe=jub-jlb+1
        nptssn=iub-ilb+1
 
#  ifdef NESTING_DEBUG
!
!  If distributed-memory, initialize arrays used to check mass flux
!  conservation with special value (zero) to facilitate the global
!  reduction when collecting data between all nodes.
!
        bry_contact(iwest ,cr)%Mflux=spv
        bry_contact(ieast ,cr)%Mflux=spv
        bry_contact(isouth,cr)%Mflux=spv
        bry_contact(inorth,cr)%Mflux=spv
#  endif
      END IF
# endif
!
!  Set time snapshot indices for the donor grid data.
!
      told=3-rollingindex(cr)
      tnew=rollingindex(cr)
!
!  Set linear time interpolation weights. Fractional seconds are
!  rounded to the nearest milliseconds integer towards zero in the
!  time interpolation weights.
!
      secscale=1000.0_dp              ! seconds to milliseconds
!
      wold=anint((rollingtime(tnew,cr)-time(ng))*secscale,dp)
      wnew=anint((time(ng)-rollingtime(told,cr))*secscale,dp)
      fac=1.0_dp/(wold+wnew)
      wold=fac*wold
      wnew=fac*wnew
!
      IF (((wold*wnew).lt.0.0_dp).or.((wold+wnew).le.0.0_dp)) THEN
        IF (domain(ng)%SouthWest_Test(tile)) THEN
          IF (master) THEN
            WRITE (stdout,10) cr, dg, ng,                               &
     &                        iic(dg), told, tnew,                      &
     &                        iic(ng), wold, wnew,                      &
     &                        int(time(ng)),                            &
     &                        int(rollingtime(told,cr)),                &
     &                        int(rollingtime(tnew,cr))
          END IF
          exit_flag=8
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
      END IF
!
!-----------------------------------------------------------------------
!  Process free-surface.
!-----------------------------------------------------------------------
!
      free_surface : IF (lputfsur) THEN
        DO m=1,rcontact(cr)%Npoints
          i=rcontact(cr)%Irg(m)
          j=rcontact(cr)%Jrg(m)
          IF (((istrt.le.i).and.(i.le.iendt)).and.                      &
     &        ((jstrt.le.j).and.(j.le.jendt))) THEN
            my_value=wold*(rcontact(cr)%Lweight(1,m)*                   &
     &                                  refined(cr)%zeta(1,m,told)+     &
     &                     rcontact(cr)%Lweight(2,m)*                   &
     &                                  refined(cr)%zeta(2,m,told)+     &
     &                     rcontact(cr)%Lweight(3,m)*                   &
     &                                  refined(cr)%zeta(3,m,told)+     &
     &                     rcontact(cr)%Lweight(4,m)*                   &
     &                                  refined(cr)%zeta(4,m,told))+    &
     &               wnew*(rcontact(cr)%Lweight(1,m)*                   &
     &                                  refined(cr)%zeta(1,m,tnew)+     &
     &                     rcontact(cr)%Lweight(2,m)*                   &
     &                                  refined(cr)%zeta(2,m,tnew)+     &
     &                     rcontact(cr)%Lweight(3,m)*                   &
     &                                  refined(cr)%zeta(3,m,tnew)+     &
     &                     rcontact(cr)%Lweight(4,m)*                   &
     &                                  refined(cr)%zeta(4,m,tnew))
# ifdef MASKING
            my_value=my_value*grid(ng)%rmask(i,j)
# endif
# ifdef WET_DRY
            IF (my_value.le.(dcrit(ng)-grid(ng)%h(i,j))) THEN
              my_value=dcrit(ng)-grid(ng)%h(i,j)
            END IF
# endif
# ifdef SOLVE3D
            ocean(ng)%zeta(i,j,1)=my_value
            ocean(ng)%zeta(i,j,2)=my_value
            ocean(ng)%zeta(i,j,3)=my_value
            coupling(ng)%Zt_avg1(i,j)=my_value
# else
            ocean(ng)%zeta(i,j,knew(ng))=my_value
# endif
          END IF
        END DO
 
      ELSE
!
!-----------------------------------------------------------------------
!  Process 2D momentum.
!-----------------------------------------------------------------------
!
!  Notice that contact points at the domain western, eastern, southern
!  and northern physical boundaries are avoided for the "kindex" time
!  record. They are  assigned in the mass flux computations below.
!  This exception is done for adjoint correctness.
!
# ifdef SOLVE3D
        kindex=indx1(ng)
# else
        kindex=knew(ng)
# endif
!
!  2D momentum in the XI-direction.
!
        DO m=1,ucontact(cr)%Npoints
          i=ucontact(cr)%Irg(m)
          j=ucontact(cr)%Jrg(m)
          IF (((istrp.le.i).and.(i.le.iendt)).and.                      &
     &        ((jstrt.le.j).and.(j.le.jendt))) THEN
            my_value=wold*(ucontact(cr)%Lweight(1,m)*                   &
     &                                  refined(cr)%ubar(1,m,told)+     &
     &                     ucontact(cr)%Lweight(2,m)*                   &
     &                                  refined(cr)%ubar(2,m,told)+     &
     &                     ucontact(cr)%Lweight(3,m)*                   &
     &                                  refined(cr)%ubar(3,m,told)+     &
     &                     ucontact(cr)%Lweight(4,m)*                   &
     &                                  refined(cr)%ubar(4,m,told))+    &
     &               wnew*(ucontact(cr)%Lweight(1,m)*                   &
     &                                  refined(cr)%ubar(1,m,tnew)+     &
     &                     ucontact(cr)%Lweight(2,m)*                   &
     &                                  refined(cr)%ubar(2,m,tnew)+     &
     &                     ucontact(cr)%Lweight(3,m)*                   &
     &                                  refined(cr)%ubar(3,m,tnew)+     &
     &                     ucontact(cr)%Lweight(4,m)*                   &
     &                                  refined(cr)%ubar(4,m,tnew))
# ifdef MASKING
            my_value=my_value*grid(ng)%umask(i,j)
# endif
# ifdef WET_DRY
            my_value=my_value*grid(ng)%umask_wet(i,j)
# endif
            uboundary=(m.eq.bry_contact(iwest,cr)%C2Bindex(j)).or.    &
     &                (m.eq.bry_contact(ieast,cr)%C2Bindex(j))
# ifdef SOLVE3D
            DO irec=1,3
              IF(.not.(uboundary.and.(irec.eq.kindex))) THEN
                ocean(ng)%ubar(i,j,irec)=my_value
!!            ELSE                                   ! for debugging
!!              OCEAN(ng)%ubar(i,j,irec)=0.0_r8      ! purposes
              END IF
            END DO
# else
            IF (.not.uboundary) THEN
              ocean(ng)%ubar(i,j,knew(ng))=my_value
            END IF
# endif
          END IF
        END DO
!
!  2D momentum in the ETA-direction.
!
        DO m=1,vcontact(cr)%Npoints
          i=vcontact(cr)%Irg(m)
          j=vcontact(cr)%Jrg(m)
          IF (((istrt.le.i).and.(i.le.iendt)).and.                      &
     &        ((jstrp.le.j).and.(j.le.jendt))) THEN
            my_value=wold*(vcontact(cr)%Lweight(1,m)*                   &
     &                                  refined(cr)%vbar(1,m,told)+     &
     &                     vcontact(cr)%Lweight(2,m)*                   &
     &                                  refined(cr)%vbar(2,m,told)+     &
     &                     vcontact(cr)%Lweight(3,m)*                   &
     &                                  refined(cr)%vbar(3,m,told)+     &
     &                     vcontact(cr)%Lweight(4,m)*                   &
     &                                  refined(cr)%vbar(4,m,told))+    &
     &               wnew*(vcontact(cr)%Lweight(1,m)*                   &
     &                                  refined(cr)%vbar(1,m,tnew)+     &
     &                     vcontact(cr)%Lweight(2,m)*                   &
     &                                  refined(cr)%vbar(2,m,tnew)+     &
     &                     vcontact(cr)%Lweight(3,m)*                   &
     &                                  refined(cr)%vbar(3,m,tnew)+     &
     &                     vcontact(cr)%Lweight(4,m)*                   &
     &                                  refined(cr)%vbar(4,m,tnew))
# ifdef MASKING
            my_value=my_value*grid(ng)%vmask(i,j)
# endif
# ifdef WET_DRY
            my_value=my_value*grid(ng)%vmask_wet(i,j)
# endif
            vboundary=(m.eq.bry_contact(isouth,cr)%C2Bindex(i)).or.     &
     &                (m.eq.bry_contact(inorth,cr)%C2Bindex(i))
# ifdef SOLVE3D
            DO irec=1,3
              IF (.not.(vboundary.and.(irec.eq.kindex))) THEN
                ocean(ng)%vbar(i,j,irec)=my_value
!!            ELSE                                   ! for debugging
!!              OCEAN(ng)%vbar(i,j,irec)=0.0_r8      ! purposes
              END IF
            END DO
# else
            IF (.not.vboundary) THEN
              ocean(ng)%vbar(i,j,knew(ng))=my_value
            END IF
# endif
          END IF
        END DO
!
!-----------------------------------------------------------------------
!  Impose mass flux at the finer grid physical boundaries. This is only
!  done for "kindex" time record.
!
!  Western/Eastern boundary:
!
!    ubar(Ibry,:,kindex) = U2d_flux(Ibry,:) * pn(Ibry,:) / D(Ibry,:)
!
!  Southern/Northern boundary:
!
!    vbar(:,Jbry,kindex) = V2d_flux(:,Jbry) * pm(:,Jbry) / D(:,Jbry)
!
# ifdef SOLVE3D
!  Notice that in 3D applications, "REFINED(cr)%U2d_flux" and
!  "REFINED(cr)%V2d_flux" are computed from "DU_avg2" and "DV_avg2d"
!  state variables, repectively.
# else
!  Notice that in 2D applications, "REFINED(cr)%U2d_flux" and
!  "REFINED(cr)%V2d_flux" are computed from "DU_flux" and "DV_flux"
!  state variables, repectively.
# endif
!
!  Use the latest coarse grid mass flux REFINED(cr)%U2D_flux(1,:,tnew)
!  with a linear variation (cff1) to ensure that the sum of the refined
!  grid fluxes equals the coarse grid flux.
!-----------------------------------------------------------------------
!
!  Western edge.
!
        IF (domain(ng)%Western_Edge(tile)) THEN
          DO j=jstr,jend
            m=bry_contact(iwest,cr)%C2Bindex(j)
            idg=ucontact(cr)%Idg(m)                 ! for debugging
            jdg=ucontact(cr)%Jdg(m)                 ! purposes
            cff=0.5_r8*grid(ng)%on_u(istr,j)*                           &
                (grid(ng)%h(istr-1,j)+                                  &
     &           ocean(ng)%zeta(istr-1,j,kindex)+                       &
     &           grid(ng)%h(istr  ,j)+                                  &
     &           ocean(ng)%zeta(istr  ,j,kindex))
            cff1=grid(ng)%on_u(istr,j)/refined(cr)%on_u(m)
# ifdef TIME_INTERP_FLUX
            my_value=cff1*(wold*refined(cr)%U2d_flux(1,m,told)+         &
     &                     wnew*refined(cr)%U2d_flux(1,m,tnew))/cff
# else
            my_value=cff1*refined(cr)%U2d_flux(1,m,tnew)/cff
# endif
# ifdef WEC
            my_value=my_value-ocean(ng)%ubar_stokes(istr,j)
# endif
# ifdef MASKING
            my_value=my_value*grid(ng)%umask(istr,j)
# endif
# ifdef WET_DRY
            my_value=my_value*grid(ng)%umask_wet(istr,j)
# endif
# ifdef NESTING_DEBUG
            bry_contact(iwest,cr)%Mflux(j)=cff*my_value
# endif
            ocean(ng)%ubar(istr,j,kindex)=my_value
          END DO
        END IF
!
!  Eastern edge.
!
        IF (domain(ng)%Eastern_Edge(tile)) THEN
          DO j=jstr,jend
            m=bry_contact(ieast,cr)%C2Bindex(j)
            idg=ucontact(cr)%Idg(m)                 ! for debugging
            jdg=ucontact(cr)%Jdg(m)                 ! purposes
            cff=0.5_r8*grid(ng)%on_u(iend+1,j)*                         &
     &          (grid(ng)%h(iend+1,j)+                                  &
     &           ocean(ng)%zeta(iend+1,j,kindex)+                       &
     &           grid(ng)%h(iend  ,j)+                                  &
     &           ocean(ng)%zeta(iend  ,j,kindex))
            cff1=grid(ng)%on_u(iend+1,j)/refined(cr)%on_u(m)
# ifdef TIME_INTERP_FLUX
            my_value=cff1*(wold*refined(cr)%U2d_flux(1,m,told)+         &
     &                     wnew*refined(cr)%U2d_flux(1,m,tnew))/cff
# else
            my_value=cff1*refined(cr)%U2d_flux(1,m,tnew)/cff
# endif
# ifdef WEC
            my_value=my_value-ocean(ng)%ubar_stokes(iend+1,j)
# endif
# ifdef MASKING
            my_value=my_value*grid(ng)%umask(iend+1,j)
# endif
# ifdef WET_DRY
            my_value=my_value*grid(ng)%umask_wet(iend+1,j)
# endif
# ifdef NESTING_DEBUG
            bry_contact(ieast,cr)%Mflux(j)=cff*my_value
# endif
            ocean(ng)%ubar(iend+1,j,kindex)=my_value
          END DO
        END IF
!
!  Southern edge.
!
        IF (domain(ng)%Southern_Edge(tile)) THEN
          DO i=istr,iend
            m=bry_contact(isouth,cr)%C2Bindex(i)
            idg=vcontact(cr)%Idg(m)                 ! for debugging
            jdg=vcontact(cr)%Jdg(m)                 ! purposes
            cff=0.5_r8*grid(ng)%om_v(i,jstr)*                           &
     &          (grid(ng)%h(i,jstr-1)+                                  &
     &           ocean(ng)%zeta(i,jstr-1,kindex)+                       &
     &           grid(ng)%h(i,jstr  )+                                  &
     &           ocean(ng)%zeta(i,jstr  ,kindex))
            cff1=grid(ng)%om_v(i,jstr)/refined(cr)%om_v(m)
# ifdef TIME_INTERP_FLUX
            my_value=cff1*(wold*refined(cr)%V2d_flux(1,m,told)+         &
     &                     wnew*refined(cr)%V2d_flux(1,m,tnew))/cff
# else
            my_value=cff1*refined(cr)%V2d_flux(1,m,tnew)/cff
# endif
# ifdef WEC
            my_value=my_value-ocean(ng)%vbar_stokes(i,jstr)
# endif
# ifdef MASKING
            my_value=my_value*grid(ng)%vmask(i,jstr)
# endif
# ifdef WET_DRY
            my_value=my_value*grid(ng)%vmask_wet(i,jstr)
# endif
# ifdef NESTING_DEBUG
            bry_contact(isouth,cr)%Mflux(i)=cff*my_value
# endif
            ocean(ng)%vbar(i,jstr,kindex)=my_value
          END DO
        END IF
!
!  Northern edge.
!
        IF (domain(ng)%Northern_Edge(tile)) THEN
          DO i=istr,iend
            m=bry_contact(inorth,cr)%C2Bindex(i)
            idg=vcontact(cr)%Idg(m)                 ! for debugging
            jdg=vcontact(cr)%Jdg(m)                 ! purposes
            cff=0.5_r8*grid(ng)%om_v(i,jend+1)*                         &
     &          (grid(ng)%h(i,jend+1)+                                  &
     &           ocean(ng)%zeta(i,jend+1,kindex)+                       &
     &           grid(ng)%h(i,jend  )+                                  &
     &           ocean(ng)%zeta(i,jend  ,kindex))
            cff1=grid(ng)%om_v(i,jend+1)/refined(cr)%om_v(m)
# ifdef TIME_INTERP_FLUX
            my_value=cff1*(wold*refined(cr)%V2d_flux(1,m,told)+         &
     &                     wnew*refined(cr)%V2d_flux(1,m,tnew))/cff
# else
            my_value=cff1*refined(cr)%V2d_flux(1,m,tnew)/cff
# endif
# ifdef WEC
            my_value=my_value-ocean(ng)%vbar_stokes(i,jend+1)
# endif
# ifdef MASKING
            my_value=my_value*grid(ng)%vmask(i,jend+1)
# endif
# ifdef WET_DRY
            my_value=my_value*grid(ng)%vmask_wet(i,jend+1)
# endif
# ifdef NESTING_DEBUG
            bry_contact(inorth,cr)%Mflux(i)=cff*my_value
# endif
            ocean(ng)%vbar(i,jend+1,kindex)=my_value
          END DO
        END IF
 
      END IF free_surface
 
# ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Exchange tile information.
!-----------------------------------------------------------------------
!
!  Free-surface.
!
      IF (lputfsur) THEN
#  ifdef SOLVE3D
        CALL mp_exchange2d (ng, tile, model, 4,                         &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      nghostpoints,                               &
     &                      ewperiodic(ng), nsperiodic(ng),             &
     &                      coupling(ng)%Zt_avg1,                       &
     &                      ocean(ng)%zeta(:,:,1),                      &
     &                      ocean(ng)%zeta(:,:,2),                      &
     &                      ocean(ng)%zeta(:,:,3))
#  else
        CALL mp_exchange2d (ng, tile, model, 1,                         &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      nghostpoints,                               &
     &                      ewperiodic(ng), nsperiodic(ng),             &
     &                      ocean(ng)%zeta(:,:,knew(ng)))
#  endif
!
!  2D momentum.
!
      ELSE
#  ifdef SOLVE3D
        CALL mp_exchange2d (ng, tile, model, 3,                         &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      nghostpoints,                               &
     &                      ewperiodic(ng), nsperiodic(ng),             &
     &                      ocean(ng)%ubar(:,:,1),                      &
     &                      ocean(ng)%ubar(:,:,2),                      &
     &                      ocean(ng)%ubar(:,:,3))
 
        CALL mp_exchange2d (ng, tile, model, 3,                         &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      nghostpoints,                               &
     &                      ewperiodic(ng), nsperiodic(ng),             &
     &                      ocean(ng)%vbar(:,:,1),                      &
     &                      ocean(ng)%vbar(:,:,2),                      &
     &                      ocean(ng)%vbar(:,:,3))
#  else
        CALL mp_exchange2d (ng, tile, model, 2,                         &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      nghostpoints,                               &
     &                      ewperiodic(ng), nsperiodic(ng),             &
     &                      ocean(ng)%ubar(:,:,knew(ng)),               &
     &                      ocean(ng)%vbar(:,:,knew(ng)))
#  endif
 
#  ifdef NESTING_DEBUG
!
        CALL mp_assemble (ng, model, nptswe, spv,                       &
     &                    bry_contact(iwest ,cr)%Mflux(jlb:))
        CALL mp_assemble (ng, model, nptswe, spv,                       &
     &                    bry_contact(ieast ,cr)%Mflux(jlb:))
        CALL mp_assemble (ng, model, nptssn, spv,                       &
     &                    bry_contact(isouth,cr)%Mflux(ilb:))
        CALL mp_assemble (ng, model, nptssn, spv,                       &
     &                    bry_contact(inorth,cr)%Mflux(ilb:))
#  endif
      END IF
# endif
!
  10  FORMAT (/,' PUT_REFINE2D - unbounded contact points temporal: ',  &
     &          ' interpolation:',                                      &
     &        /,2x,  'cr       = ',i2.2,                                &
     &          8x,'dg         = ',i2.2,                                &
     &          8x,'ng         = ',i2.2,                                &
     &        /,2x,  'iic(dg)  = ',i7.7,                                &
     &          3x,'told       = ',i1,                                  &
     &          9x,'tnew       = ',i1,                                  &
     &        /,2x,  'iic(ng)  = ',i7.7,                                &
     &          3x,'Wold       = ',f8.5,                                &
     &          2x,'Wnew       = ',f8.5,                                &
     &        /,2x,  'time(ng) = ',i7.7,                                &
     &          3x,'time(told) = ',i7.7,                                &
     &          3x,'time(tnew) = ',i7.7)
!
      RETURN
      END SUBROUTINE put_refine2d
 
# ifdef SOLVE3D
!
      SUBROUTINE put_refine3d (ng, dg, cr, model, tile,                 &
     &                         LBi, UBi, LBj, UBj)
!
!=======================================================================
!                                                                      !
!  This routine interpolates (space, time) refinement grid 3D state    !
!  variables contact points using data from the donor grid.            !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Refinement (receiver) grid number (integer)           !
!     dg         Donor grid number (integer)                           !
!     cr         Contact region number to process (integer)            !
!     model      Calling model identifier (integer)                    !
!     tile       Domain tile partition (integer)                       !
!     LBi        Receiver grid, I-dimension Lower bound (integer)      !
!     UBi        Receiver grid, I-dimension Upper bound (integer)      !
!     LBj        Receiver grid, J-dimension Lower bound (integer)      !
!     UBj        Receiver grid, J-dimension Upper bound (integer)      !
!                                                                      !
!  On Output:    OCEAN(ng) structure                                   !
!                                                                      !
!     t          Updated tracer-type variables                         !
!     u          Updated 3D momentum in the XI-direction               !
!     v          Updated 3D momentum in the ETA-direction              !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_grid
      USE mod_nesting
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
      USE mod_iounits
 
#  ifdef DISTRIBUTE
!
      USE mp_exchange_mod, ONLY : mp_exchange3d, mp_exchange4d
#  endif
      USE strings_mod,     ONLY : founderror
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, dg, cr, model, tile
      integer, intent(in) :: lbi, ubi, lbj, ubj
!
!  Local variable declarations.
!
# ifdef NESTING_DEBUG
      logical, save :: first = .true.
!
# endif
      integer :: i, itrc, j, k, m, tnew, told
!
      real(dp) :: wnew, wold, secscale, fac
      real(r8) :: my_value
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", put_refined3d"
 
#  include "set_bounds.h"
!
!-----------------------------------------------------------------------
!  Interpolate (space, time) refinement grid contact points for 2D state
!  variables from donor grid.
!-----------------------------------------------------------------------
!
!  Set time snapshot indices for the donor grid data.
!
      told=3-rollingindex(cr)
      tnew=rollingindex(cr)
!
!  Set linear time interpolation weights. Fractional seconds are
!  rounded to the nearest milliseconds integer towards zero in the
!  time interpolation weights.
!
      secscale=1000.0_dp              ! seconds to milliseconds
!
      wold=anint((rollingtime(tnew,cr)-time(ng))*secscale,dp)
      wnew=anint((time(ng)-rollingtime(told,cr))*secscale,dp)
      fac=1.0_dp/(wold+wnew)
      wold=fac*wold
      wnew=fac*wnew
!
      IF (((wold*wnew).lt.0.0_dp).or.((wold+wnew).le.0.0_dp)) THEN
        IF (domain(ng)%SouthWest_Test(tile)) THEN
          IF (master) THEN
            WRITE (stdout,10) cr, dg, ng,                               &
     &                        iic(dg), told, tnew,                      &
     &                        iic(ng), wold, wnew,                      &
     &                        int(time(ng)),                            &
     &                        int(rollingtime(told,cr)),                &
     &                        int(rollingtime(tnew,cr))
          END IF
          exit_flag=8
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
      END IF
 
# ifdef NESTING_DEBUG
!
!  If debugging, write information into Fortran unit 200 to check the
!  logic of interpolating from donor grid data.
!
      IF (domain(ng)%SouthWest_Test(tile)) THEN
        IF (master) THEN
          IF (first) THEN
            first=.false.
            WRITE (200,20)
          END IF
          WRITE (200,30) cr, dg, ng, iic(dg), iic(ng), told,  tnew,     &
     &                   int(time(ng)),                                 &
     &                   int(rollingtime(told,cr)),                     &
     &                   int(rollingtime(tnew,cr)),                     &
     &                   wold, wnew
 20       FORMAT (3x,'cr',3x,'dg',3x,'ng',3x,'iic',2x,'iic',1x,'told',  &
     &            1x,'tnew',3x,'time',3x,'time',3x,'time',5x,'Wold',    &
     &            7x,'Wnew',/,17x,'(dg)',1x,'(ng)',13x,'(ng)',3x,       &
     &            'told',3x,'tnew',/)
 30       FORMAT (7i5,3i7,2f11.4)
          FLUSH (200)
        END IF
      END IF
# endif
!
!  Tracer-type variables.
!
      DO m=1,rcontact(cr)%Npoints
        i=rcontact(cr)%Irg(m)
        j=rcontact(cr)%Jrg(m)
        IF (((istrt.le.i).and.(i.le.iendt)).and.                        &
     &      ((jstrt.le.j).and.(j.le.jendt))) THEN
          DO itrc=1,nt(ng)
            DO k=1,n(ng)
              my_value=wold*                                            &
     &                 (rcontact(cr)%Lweight(1,m)*                      &
     &                               refined(cr)%t(1,k,m,told,itrc)+    &
     &                  rcontact(cr)%Lweight(2,m)*                      &
     &                               refined(cr)%t(2,k,m,told,itrc)+    &
     &                  rcontact(cr)%Lweight(3,m)*                      &
     &                               refined(cr)%t(3,k,m,told,itrc)+    &
     &                  rcontact(cr)%Lweight(4,m)*                      &
     &                               refined(cr)%t(4,k,m,told,itrc))+   &
     &                 wnew*                                            &
     &                 (rcontact(cr)%Lweight(1,m)*                      &
     &                               refined(cr)%t(1,k,m,tnew,itrc)+    &
     &                  rcontact(cr)%Lweight(2,m)*                      &
     &                               refined(cr)%t(2,k,m,tnew,itrc)+    &
     &                  rcontact(cr)%Lweight(3,m)*                      &
     &                               refined(cr)%t(3,k,m,tnew,itrc)+    &
     &                  rcontact(cr)%Lweight(4,m)*                      &
     &                               refined(cr)%t(4,k,m,tnew,itrc))
#  ifdef MASKING
              my_value=my_value*grid(ng)%rmask(i,j)
#  endif
              ocean(ng)%t(i,j,k,1,itrc)=my_value
              ocean(ng)%t(i,j,k,2,itrc)=my_value
              ocean(ng)%t(i,j,k,3,itrc)=my_value
            END DO
          END DO
        END IF
      END DO
!
!  3D momentum in the XI-direction.
!
      DO m=1,ucontact(cr)%Npoints
        i=ucontact(cr)%Irg(m)
        j=ucontact(cr)%Jrg(m)
        IF (((istrp.le.i).and.(i.le.iendt)).and.                        &
     &      ((jstrt.le.j).and.(j.le.jendt))) THEN
          DO k=1,n(ng)
            my_value=wold*                                              &
     &               (ucontact(cr)%Lweight(1,m)*                        &
     &                             refined(cr)%u(1,k,m,told)+           &
     &                ucontact(cr)%Lweight(2,m)*                        &
     &                             refined(cr)%u(2,k,m,told)+           &
     &                ucontact(cr)%Lweight(3,m)*                        &
     &                             refined(cr)%u(3,k,m,told)+           &
     &                ucontact(cr)%Lweight(4,m)*                        &
     &                             refined(cr)%u(4,k,m,told))+          &
     &               wnew*                                              &
     &               (ucontact(cr)%Lweight(1,m)*                        &
     &                             refined(cr)%u(1,k,m,tnew)+           &
     &                ucontact(cr)%Lweight(2,m)*                        &
     &                             refined(cr)%u(2,k,m,tnew)+           &
     &                ucontact(cr)%Lweight(3,m)*                        &
     &                             refined(cr)%u(3,k,m,tnew)+           &
     &                ucontact(cr)%Lweight(4,m)*                        &
     &                             refined(cr)%u(4,k,m,tnew))
#  ifdef MASKING
            my_value=my_value*grid(ng)%umask(i,j)
#  endif
            ocean(ng)%u(i,j,k,1)=my_value
            ocean(ng)%u(i,j,k,2)=my_value
          END DO
        END IF
      END DO
!
!  3D momentum in the ETA-direction.
!
      DO m=1,vcontact(cr)%Npoints
        i=vcontact(cr)%Irg(m)
        j=vcontact(cr)%Jrg(m)
        IF (((istrt.le.i).and.(i.le.iendt)).and.                        &
     &      ((jstrp.le.j).and.(j.le.jendt))) THEN
          DO k=1,n(ng)
            my_value=wold*                                              &
     &               (vcontact(cr)%Lweight(1,m)*                        &
     &                             refined(cr)%v(1,k,m,told)+           &
     &                vcontact(cr)%Lweight(2,m)*                        &
     &                             refined(cr)%v(2,k,m,told)+           &
     &                vcontact(cr)%Lweight(3,m)*                        &
     &                             refined(cr)%v(3,k,m,told)+           &
     &                vcontact(cr)%Lweight(4,m)*                        &
     &                             refined(cr)%v(4,k,m,told))+          &
     &               wnew*                                              &
     &               (vcontact(cr)%Lweight(1,m)*                        &
     &                             refined(cr)%v(1,k,m,tnew)+           &
     &                vcontact(cr)%Lweight(2,m)*                        &
     &                             refined(cr)%v(2,k,m,tnew)+           &
     &                vcontact(cr)%Lweight(3,m)*                        &
     &                             refined(cr)%v(3,k,m,tnew)+           &
     &                vcontact(cr)%Lweight(4,m)*                        &
     &                             refined(cr)%v(4,k,m,tnew))
#  ifdef MASKING
            my_value=my_value*grid(ng)%vmask(i,j)
#  endif
            ocean(ng)%v(i,j,k,1)=my_value
            ocean(ng)%v(i,j,k,2)=my_value
          END DO
        END IF
      END DO
 
#  ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Exchange tile information.
!-----------------------------------------------------------------------
!
      CALL mp_exchange4d (ng, tile, model, 3,                           &
     &                    lbi, ubi, lbj, ubj, 1, n(ng), 1, nt(ng),      &
     &                    nghostpoints,                                 &
     &                    ewperiodic(ng), nsperiodic(ng),               &
     &                    ocean(ng)%t(:,:,:,1,:),                       &
     &                    ocean(ng)%t(:,:,:,2,:),                       &
     &                    ocean(ng)%t(:,:,:,3,:))
 
      CALL mp_exchange3d (ng, tile, model, 4,                           &
     &                    lbi, ubi, lbj, ubj, 1, n(ng),                 &
     &                    nghostpoints,                                 &
     &                    ewperiodic(ng), nsperiodic(ng),               &
     &                    ocean(ng)%u(:,:,:,1),                         &
     &                    ocean(ng)%u(:,:,:,2),                         &
     &                    ocean(ng)%v(:,:,:,1),                         &
     &                    ocean(ng)%v(:,:,:,2))
#  endif
!
  10  FORMAT (/,' PUT_REFINE3D - unbounded contact points temporal: ',  &
     &          ' interpolation:',                                      &
     &        /,2x,  'cr       = ',i2.2,                                &
     &          8x,'dg         = ',i2.2,                                &
     &          8x,'ng         = ',i2.2,                                &
     &        /,2x,  'iic(dg)  = ',i7.7,                                &
     &          3x,'told       = ',i1,                                  &
     &          9x,'tnew       = ',i1,                                  &
     &        /,2x,  'iic(ng)  = ',i7.7,                                &
     &          3x,'Wold       = ',f8.5,                                &
     &          2x,'Wnew       = ',f8.5,                                &
     &        /,2x,  'time(ng) = ',i7.7,                                &
     &          3x,'time(told) = ',i7.7,                                &
     &          3x,'time(tnew) = ',i7.7)
!
      RETURN
      END SUBROUTINE put_refine3d
# endif
 
# ifdef SOLVE3D
!
      SUBROUTINE z_weights (ng, model, tile)
!
!=======================================================================
!                                                                      !
!  This routine determines the vertical indices and interpolation      !
!  weights associated with depth,  which are needed to process 3D      !
!  fields in the contact region.                                       !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     model      Calling model identifier (integer)                    !
!     tile       Domain partition for composite grid ng (integer)      !
!                                                                      !
!  On Output:    Updated T_NGC type structures in mod_param:           !
!                                                                      !
!     Rcontact   Updated values for Kdg(:,:) and Vweigths (:,:,:)      !
!     Ucontact   Updated values for Kdg(:,:) and Vweigths (:,:,:)      !
!     Vcontact   Updated values for Kdg(:,:) and Vweigths (:,:,:)      !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_grid
      USE mod_nesting
      USE mod_scalars
!
#  ifdef DISTRIBUTE
      USE distribute_mod, ONLY : mp_assemble
#  endif
      USE strings_mod,    ONLY : founderror
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model, tile
!
!  Local variable declarations.
!
      integer :: cr, dg, rg, i, j, k, m
      integer :: idg, jdg, kdg, imind, imaxd, jmind, jmaxd
      integer :: irg, jrg, krg, iminr, imaxr, jminr, jmaxr
      integer :: idgm1, idgp1, jdgm1, jdgp1
      integer :: npoints
#  ifdef DISTRIBUTE
      integer :: nkpts, nwpts, nzpts
 
      integer, parameter :: ispv = 0
#  endif
!
      real(r8), parameter :: spv = 0.0_r8
 
      real(r8) :: zbot, zr, ztop, dz, r1, r2
 
      real(r8), allocatable :: zd(:,:,:)
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", z_weights"
!
!=======================================================================
!  Compute vertical indices and weights for each contact region.
!=======================================================================
!
      DO cr=1,ncontact
!
!  Get donor and receiver grid numbers.
!
        dg=rcontact(cr)%donor_grid
        rg=rcontact(cr)%receiver_grid
!
!  Process only contact region data for requested nested grid "ng".
!
        IF (rg.eq.ng) THEN
!
!-----------------------------------------------------------------------
!  Process variables in structure Rcontact(cr).
!-----------------------------------------------------------------------
!
!  Get number of contact points to process.
!
          npoints=rcontact(cr)%Npoints
!
!  Set starting and ending tile indices for the donor and receiver
!  grids.
!
          imind=bounds(dg) % IstrT(tile)
          imaxd=bounds(dg) % IendT(tile)
          jmind=bounds(dg) % JstrT(tile)
          jmaxd=bounds(dg) % JendT(tile)
!
          iminr=bounds(rg) % IstrT(tile)
          imaxr=bounds(rg) % IendT(tile)
          jminr=bounds(rg) % JstrT(tile)
          jmaxr=bounds(rg) % JendT(tile)
 
#  ifdef DISTRIBUTE
!
!  If distributed-memory, initialize with special value (zero) to
!  facilitate the global reduction when collecting data between all
!  nodes.
!
          nkpts=n(rg)*npoints
          nwpts=2*nkpts
          nzpts=4*nkpts
 
          rcontact(cr)%Kdg(1:n(rg),1:npoints)=ispv
          rcontact(cr)%Vweight(1:2,1:n(rg),1:npoints)=spv
#  endif
!
!  If coincident grids and requested, avoid vertical interpolation.
!
          r_contact : IF (.not.rcontact(cr)%interpolate.and.            &
     &                    rcontact(cr)%coincident) THEN
            DO krg=1,n(rg)
              DO m=1,npoints
                irg=rcontact(cr)%Irg(m)
                jrg=rcontact(cr)%Jrg(m)
                IF (((iminr.le.irg).and.(irg.le.imaxr)).and.            &
     &              ((jminr.le.jrg).and.(jrg.le.jmaxr))) THEN
                  rcontact(cr)%Kdg(krg,m)=krg
                  rcontact(cr)%Vweight(1,krg,m)=1.0_r8
                  rcontact(cr)%Vweight(2,krg,m)=0.0_r8
                END IF
              END DO
            END DO
!
!  Otherwise, vertically interpolate because donor and receiver grids
!  are not coincident.
!
          ELSE
!
!  Allocate and initialize local working arrays.
!
            IF (.not.allocated(zd)) THEN
              allocate ( zd(4,n(dg),npoints) )
            END IF
            zd=spv
!
!  Extract donor grid depths for each cell containing the receiver grid
!  contact point.  Notice that indices i+1 and j+1 are bounded to the
!  maximum possible values in contact points at the edge of the contact
!  region.  In such cases, Lweight(1,m)=1 and Lweight(2:3,m)=0.  This is
!  done to avoid out of range errors. We need to take care of this in
!  the adjoint code.
!
            DO kdg=1,n(dg)
              DO m=1,npoints
                idg  =rcontact(cr)%Idg(m)
                idgp1=min(idg+1, bounds(dg)%UBi(-1))
                jdg  =rcontact(cr)%Jdg(m)
                jdgp1=min(jdg+1, bounds(dg)%UBj(-1))
                IF (((imind.le.idg).and.(idg.le.imaxd)).and.            &
     &              ((jmind.le.jdg).and.(jdg.le.jmaxd))) THEN
                  zd(1,kdg,m)=grid(dg)%z_r(idg  ,jdg  ,kdg)
                  zd(2,kdg,m)=grid(dg)%z_r(idgp1,jdg  ,kdg)
                  zd(3,kdg,m)=grid(dg)%z_r(idgp1,jdgp1,kdg)
                  zd(4,kdg,m)=grid(dg)%z_r(idg  ,jdgp1,kdg)
                END IF
              END DO
            END DO
 
#  ifdef DISTRIBUTE
!
!  Exchange data between all parallel nodes.
!
            CALL mp_assemble (dg, model, nzpts, spv, zd)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
!
!  Determine donor grid vertical indices (Kdg) and weights (Vweight)
!  needed for the interpolation of data at the receiver grid contact
!  points.
!
            DO krg=1,n(rg)
              DO m=1,npoints
                irg=rcontact(cr)%Irg(m)
                jrg=rcontact(cr)%Jrg(m)
                IF (((iminr.le.irg).and.(irg.le.imaxr)).and.            &
     &              ((jminr.le.jrg).and.(jrg.le.jmaxr))) THEN
                  ztop=rcontact(cr)%Lweight(1,m)*zd(1,n(dg),m)+         &
     &                 rcontact(cr)%Lweight(2,m)*zd(2,n(dg),m)+         &
     &                 rcontact(cr)%Lweight(3,m)*zd(3,n(dg),m)+         &
     &                 rcontact(cr)%Lweight(4,m)*zd(4,n(dg),m)
                  zbot=rcontact(cr)%Lweight(1,m)*zd(1,1    ,m)+         &
     &                 rcontact(cr)%Lweight(2,m)*zd(2,1    ,m)+         &
     &                 rcontact(cr)%Lweight(3,m)*zd(3,1    ,m)+         &
     &                 rcontact(cr)%Lweight(4,m)*zd(4,1    ,m)
                  zr=grid(rg)%z_r(irg,jrg,krg)
                  IF (zr.ge.ztop) THEN           ! If shallower, use top
                    rcontact(cr)%Kdg(krg,m)=n(dg)! donor grid cell value
                    rcontact(cr)%Vweight(1,krg,m)=0.0_r8
                    rcontact(cr)%Vweight(2,krg,m)=1.0_r8
                  ELSE IF (zbot.ge.zr) THEN      ! If deeper, use bottom
                    rcontact(cr)%Kdg(krg,m)=1    ! donor grid cell value
                    rcontact(cr)%Vweight(1,krg,m)=0.0_r8
                    rcontact(cr)%Vweight(2,krg,m)=1.0_r8
                  ELSE                           ! bounded, interpolate
                    DO kdg=n(dg),2,-1
                      ztop=rcontact(cr)%Lweight(1,m)*zd(1,kdg  ,m)+     &
     &                     rcontact(cr)%Lweight(2,m)*zd(2,kdg  ,m)+     &
     &                     rcontact(cr)%Lweight(3,m)*zd(3,kdg  ,m)+     &
     &                     rcontact(cr)%Lweight(4,m)*zd(4,kdg  ,m)
                      zbot=rcontact(cr)%Lweight(1,m)*zd(1,kdg-1,m)+     &
     &                     rcontact(cr)%Lweight(2,m)*zd(2,kdg-1,m)+     &
     &                     rcontact(cr)%Lweight(3,m)*zd(3,kdg-1,m)+     &
     &                     rcontact(cr)%Lweight(4,m)*zd(4,kdg-1,m)
                      IF ((ztop.gt.zr).and.(zr.ge.zbot)) THEN
                        dz=ztop-zbot
                        r2=(zr-zbot)/dz
                        r1=1.0_r8-r2
                        rcontact(cr)%Kdg(krg,m)=kdg
                        rcontact(cr)%Vweight(1,krg,m)=r1
                        rcontact(cr)%Vweight(2,krg,m)=r2
                      END IF
                    END DO
                  END IF
                END IF
              END DO
            END DO
          END IF r_contact
 
#  ifdef DISTRIBUTE
!
!  Exchange data between all parallel nodes.
!
          CALL mp_assemble (rg, model, nkpts, ispv, rcontact(cr)%Kdg)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          CALL mp_assemble (rg, model, nwpts, spv, rcontact(cr)%Vweight)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
!
!  Deallocate local work arrays.
!
          IF (allocated(zd)) THEN
            deallocate (zd)
          END IF
!
!-----------------------------------------------------------------------
!  Process variables in structure Ucontact(cr).
!-----------------------------------------------------------------------
!
!  Get number of contact points to process.
!
          npoints=ucontact(cr)%Npoints
!
!  Set starting and ending tile indices for the donor and receiver
!  grids.
!
          imind=bounds(dg) % IstrP(tile)
          imaxd=bounds(dg) % IendT(tile)
          jmind=bounds(dg) % JstrT(tile)
          jmaxd=bounds(dg) % JendT(tile)
!
          iminr=bounds(rg) % IstrP(tile)
          imaxr=bounds(rg) % IendT(tile)
          jminr=bounds(rg) % JstrT(tile)
          jmaxr=bounds(rg) % JendT(tile)
 
#  ifdef DISTRIBUTE
!
!  If distributed-memory, initialize with special value (zero) to
!  facilitate the global reduction when collecting data between all
!  nodes.
!
          nkpts=n(rg)*npoints
          nwpts=2*nkpts
          nzpts=4*nkpts
 
          ucontact(cr)%Kdg(1:n(rg),1:npoints)=ispv
          ucontact(cr)%Vweight(1:2,1:n(rg),1:npoints)=spv
#  endif
!
!  If coincident grids and requested, avoid vertical interpolation.
!
          u_contact : IF (.not.ucontact(cr)%interpolate.and.            &
     &                    ucontact(cr)%coincident) THEN
            DO krg=1,n(rg)
              DO m=1,npoints
                irg=ucontact(cr)%Irg(m)
                jrg=ucontact(cr)%Jrg(m)
                IF (((iminr.le.irg).and.(irg.le.imaxr)).and.            &
     &              ((jminr.le.jrg).and.(jrg.le.jmaxr))) THEN
                  ucontact(cr)%Kdg(krg,m)=krg
                  ucontact(cr)%Vweight(1,krg,m)=1.0_r8
                  ucontact(cr)%Vweight(2,krg,m)=0.0_r8
                END IF
              END DO
            END DO
!
!  Otherwise, vertically interpolate because donor and receiver grids
!  are not coincident.
!
          ELSE
!
!  Allocate and initialize local working arrays.
!
            IF (.not.allocated(zd)) THEN
              allocate (zd(4,n(dg),npoints))
            END IF
            zd=spv
!
!  Extract donor grid depths for each cell containing the receiver grid
!  contact point.  Notice that indices i-1, i+1 and j-1, j+1 are bounded
!  the minimum/maximum possible values in contact points at the edge of
!  the contact region.  In such cases, the interpolation weights
!  Lweight(1,m)=1 and Lweight(2:3,m)=0.  This is done to avoid out of
!  range errors. We need to take care of this in the adjoint code.
!
            DO kdg=1,n(dg)
              DO m=1,npoints
                idg  =ucontact(cr)%Idg(m)
                idgm1=max(idg-1, bounds(dg)%LBi(-1))
                idgp1=min(idg+1, bounds(dg)%UBi(-1))
                jdg  =ucontact(cr)%Jdg(m)
                jdgp1=min(jdg+1, bounds(dg)%UBj(-1))
                IF (((imind.le.idg).and.(idg.le.imaxd)).and.            &
     &              ((jmind.le.jdg).and.(jdg.le.jmaxd))) THEN
                  zd(1,kdg,m)=0.5_r8*(grid(dg)%z_r(idgm1,jdg  ,kdg)+    &
     &                                grid(dg)%z_r(idg  ,jdg  ,kdg))
                  zd(2,kdg,m)=0.5_r8*(grid(dg)%z_r(idg  ,jdg  ,kdg)+    &
     &                                grid(dg)%z_r(idgp1,jdg  ,kdg))
                  zd(3,kdg,m)=0.5_r8*(grid(dg)%z_r(idg  ,jdgp1,kdg)+    &
     &                                grid(dg)%z_r(idgp1,jdgp1,kdg))
                  zd(4,kdg,m)=0.5_r8*(grid(dg)%z_r(idgm1,jdgp1,kdg)+    &
     &                                grid(dg)%z_r(idg  ,jdgp1,kdg))
                END IF
              END DO
            END DO
 
#  ifdef DISTRIBUTE
!
!  Exchange data between all parallel nodes.
!
            CALL mp_assemble (dg, model, nzpts, spv, zd)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
!
!  Determine donor grid vertical indices (Kdg) and weights (Vweight)
!  needed for the interpolation of data at the receiver grid contact
!  points.
!
            DO krg=1,n(rg)
              DO m=1,npoints
                irg=ucontact(cr)%Irg(m)
                jrg=ucontact(cr)%Jrg(m)
                IF (((iminr.le.irg).and.(irg.le.imaxr)).and.            &
     &              ((jminr.le.jrg).and.(jrg.le.jmaxr))) THEN
                  ztop=ucontact(cr)%Lweight(1,m)*zd(1,n(dg),m)+         &
     &                 ucontact(cr)%Lweight(2,m)*zd(2,n(dg),m)+         &
     &                 ucontact(cr)%Lweight(3,m)*zd(3,n(dg),m)+         &
     &                 ucontact(cr)%Lweight(4,m)*zd(4,n(dg),m)
                  zbot=ucontact(cr)%Lweight(1,m)*zd(1,1    ,m)+         &
     &                 ucontact(cr)%Lweight(2,m)*zd(2,1    ,m)+         &
     &                 ucontact(cr)%Lweight(3,m)*zd(3,1    ,m)+         &
     &                 ucontact(cr)%Lweight(4,m)*zd(4,1    ,m)
                  zr=0.5_r8*(grid(rg)%z_r(irg  ,jrg,krg)+               &
     &                       grid(rg)%z_r(irg-1,jrg,krg))
                  IF (zr.ge.ztop) THEN           ! If shallower, use top
                    ucontact(cr)%Kdg(krg,m)=n(dg)! donor grid cell value
                    ucontact(cr)%Vweight(1,krg,m)=0.0_r8
                    ucontact(cr)%Vweight(2,krg,m)=1.0_r8
                  ELSE IF (zbot.ge.zr) THEN      ! If deeper, use bottom
                    ucontact(cr)%Kdg(krg,m)=1    ! donor grid cell value
                    ucontact(cr)%Vweight(1,krg,m)=0.0_r8
                    ucontact(cr)%Vweight(2,krg,m)=1.0_r8
                  ELSE                           ! bounded, interpolate
                    DO kdg=n(dg),2,-1
                      ztop=ucontact(cr)%Lweight(1,m)*zd(1,kdg  ,m)+     &
     &                     ucontact(cr)%Lweight(2,m)*zd(2,kdg  ,m)+     &
     &                     ucontact(cr)%Lweight(3,m)*zd(3,kdg  ,m)+     &
     &                     ucontact(cr)%Lweight(4,m)*zd(4,kdg  ,m)
                      zbot=ucontact(cr)%Lweight(1,m)*zd(1,kdg-1,m)+     &
     &                     ucontact(cr)%Lweight(2,m)*zd(2,kdg-1,m)+     &
     &                     ucontact(cr)%Lweight(3,m)*zd(3,kdg-1,m)+     &
     &                     ucontact(cr)%Lweight(4,m)*zd(4,kdg-1,m)
                      IF ((ztop.gt.zr).and.(zr.ge.zbot)) THEN
                        dz=ztop-zbot
                        r2=(zr-zbot)/dz
                        r1=1.0_r8-r2
                        ucontact(cr)%Kdg(krg,m)=kdg
                        ucontact(cr)%Vweight(1,krg,m)=r1
                        ucontact(cr)%Vweight(2,krg,m)=r2
                      END IF
                    END DO
                  END IF
                END IF
              END DO
            END DO
          END IF u_contact
 
#  ifdef DISTRIBUTE
!
!  Exchange data between all parallel nodes.
!
          CALL mp_assemble (rg, model, nkpts, ispv, ucontact(cr)%Kdg)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          CALL mp_assemble (rg, model, nwpts, spv, ucontact(cr)%Vweight)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
!
!  Deallocate local work arrays.
!
          IF (allocated(zd)) THEN
            deallocate (zd)
          END IF
!
!-----------------------------------------------------------------------
!  Process variables in structure Vcontact(cr).
!-----------------------------------------------------------------------
!
!  Get number of contact points to process.
!
          npoints=vcontact(cr)%Npoints
!
!  Set starting and ending tile indices for the donor and receiver
!  grids.
!
          imind=bounds(dg) % IstrT(tile)
          imaxd=bounds(dg) % IendT(tile)
          jmind=bounds(dg) % JstrP(tile)
          jmaxd=bounds(dg) % JendT(tile)
!
          iminr=bounds(rg) % IstrT(tile)
          imaxr=bounds(rg) % IendT(tile)
          jminr=bounds(rg) % JstrP(tile)
          jmaxr=bounds(rg) % JendT(tile)
 
#  ifdef DISTRIBUTE
!
!  If distributed-memory, initialize with special value (zero) to
!  facilitate the global reduction when collecting data between all
!  nodes.
!
          nkpts=n(rg)*npoints
          nwpts=2*nkpts
          nzpts=4*nkpts
 
          vcontact(cr)%Kdg(1:n(rg),1:npoints)=ispv
          vcontact(cr)%Vweight(1:2,1:n(rg),1:npoints)=spv
#  endif
!
!  If coincident grids and requested, avoid vertical interpolation.
!
          v_contact : IF (.not.vcontact(cr)%interpolate.and.            &
     &                    vcontact(cr)%coincident) THEN
            DO krg=1,n(rg)
              DO m=1,npoints
                irg=vcontact(cr)%Irg(m)
                jrg=vcontact(cr)%Jrg(m)
                IF (((iminr.le.irg).and.(irg.le.imaxr)).and.            &
     &              ((jminr.le.jrg).and.(jrg.le.jmaxr))) THEN
                  vcontact(cr)%Kdg(krg,m)=krg
                  vcontact(cr)%Vweight(1,krg,m)=1.0_r8
                  vcontact(cr)%Vweight(2,krg,m)=0.0_r8
                END IF
              END DO
            END DO
!
!  Otherwise, vertically interpolate because donor and receiver grids
!  are not coincident.
!
          ELSE
!
!  Allocate and initialize local working arrays.
!
            IF (.not.allocated(zd)) THEN
              allocate (zd(4,n(dg),npoints))
            END IF
            zd=spv
!
!  Extract donor grid depths for each cell containing the receiver grid
!  contact point.
!
            DO kdg=1,n(dg)
              DO m=1,npoints
                idg=vcontact(cr)%Idg(m)
                idgp1=min(idg+1, bounds(dg)%UBi(-1))
                jdg=vcontact(cr)%Jdg(m)
                jdgm1=max(jdg-1, bounds(dg)%LBj(-1))
                jdgp1=min(jdg+1, bounds(dg)%UBj(-1))
                IF (((imind.le.idg).and.(idg.le.imaxd)).and.            &
     &              ((jmind.le.jdg).and.(jdg.le.jmaxd))) THEN
                  zd(1,kdg,m)=0.5_r8*(grid(dg)%z_r(idg  ,jdgm1,kdg)+    &
     &                                grid(dg)%z_r(idg  ,jdg  ,kdg))
                  zd(2,kdg,m)=0.5_r8*(grid(dg)%z_r(idgp1,jdgm1,kdg)+    &
     &                                grid(dg)%z_r(idgp1,jdg  ,kdg))
                  zd(3,kdg,m)=0.5_r8*(grid(dg)%z_r(idgp1,jdg  ,kdg)+    &
     &                                grid(dg)%z_r(idgp1,jdgp1,kdg))
                  zd(4,kdg,m)=0.5_r8*(grid(dg)%z_r(idg  ,jdg  ,kdg)+    &
     &                                grid(dg)%z_r(idg  ,jdgp1,kdg))
                END IF
              END DO
            END DO
 
#  ifdef DISTRIBUTE
!
!  Exchange data between all parallel nodes.
!
            CALL mp_assemble (dg, model, nzpts, spv, zd)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
!
!  Determine donor grid vertical indices (Kdg) and weights (Vweight)
!  needed for the interpolation of data at the receiver grid contact
!  points.
!
            DO krg=1,n(rg)
              DO m=1,npoints
                irg=vcontact(cr)%Irg(m)
                jrg=vcontact(cr)%Jrg(m)
                IF (((iminr.le.irg).and.(irg.le.imaxr)).and.            &
     &              ((jminr.le.jrg).and.(jrg.le.jmaxr))) THEN
                  ztop=vcontact(cr)%Lweight(1,m)*zd(1,n(dg),m)+         &
     &                 vcontact(cr)%Lweight(2,m)*zd(2,n(dg),m)+         &
     &                 vcontact(cr)%Lweight(3,m)*zd(3,n(dg),m)+         &
     &                 vcontact(cr)%Lweight(4,m)*zd(4,n(dg),m)
                  zbot=vcontact(cr)%Lweight(1,m)*zd(1,1    ,m)+         &
     &                 vcontact(cr)%Lweight(2,m)*zd(2,1    ,m)+         &
     &                 vcontact(cr)%Lweight(3,m)*zd(3,1    ,m)+         &
     &                 vcontact(cr)%Lweight(4,m)*zd(4,1    ,m)
                  zr=0.5_r8*(grid(rg)%z_r(irg,jrg  ,krg)+               &
     &                       grid(rg)%z_r(irg,jrg-1,krg))
                  IF (zr.ge.ztop) THEN           ! If shallower, use top
                    vcontact(cr)%Kdg(krg,m)=n(dg)! donor grid cell value
                    vcontact(cr)%Vweight(1,krg,m)=0.0_r8
                    vcontact(cr)%Vweight(2,krg,m)=1.0_r8
                  ELSE IF (zbot.ge.zr) THEN      ! If deeper, use bottom
                    vcontact(cr)%Kdg(krg,m)=1    ! donor grid cell value
                    vcontact(cr)%Vweight(1,krg,m)=0.0_r8
                    vcontact(cr)%Vweight(2,krg,m)=1.0_r8
                  ELSE                           ! bounded, interpolate
                    DO kdg=n(dg),2,-1
                      ztop=vcontact(cr)%Lweight(1,m)*zd(1,kdg  ,m)+     &
     &                     vcontact(cr)%Lweight(2,m)*zd(2,kdg  ,m)+     &
     &                     vcontact(cr)%Lweight(3,m)*zd(3,kdg  ,m)+     &
     &                     vcontact(cr)%Lweight(4,m)*zd(4,kdg  ,m)
                      zbot=vcontact(cr)%Lweight(1,m)*zd(1,kdg-1,m)+     &
     &                     vcontact(cr)%Lweight(2,m)*zd(2,kdg-1,m)+     &
     &                     vcontact(cr)%Lweight(3,m)*zd(3,kdg-1,m)+     &
     &                     vcontact(cr)%Lweight(4,m)*zd(4,kdg-1,m)
                      IF ((ztop.gt.zr).and.(zr.ge.zbot)) THEN
                        dz=ztop-zbot
                        r2=(zr-zbot)/dz
                        r1=1.0_r8-r2
                        vcontact(cr)%Kdg(krg,m)=kdg
                        vcontact(cr)%Vweight(1,krg,m)=r1
                        vcontact(cr)%Vweight(2,krg,m)=r2
                      END IF
                    END DO
                  END IF
                END IF
              END DO
            END DO
          END IF v_contact
 
#  ifdef DISTRIBUTE
!
!  Exchange data between all parallel nodes.
!
          CALL mp_assemble (rg, model, nkpts, ispv, vcontact(cr)%Kdg)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          CALL mp_assemble (rg, model, nwpts, spv, vcontact(cr)%Vweight)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
!
!  Deallocate local work arrays.
!
          IF (allocated(zd)) THEN
            deallocate (zd)
          END IF
 
        END IF
      END DO
!
      RETURN
      END SUBROUTINE z_weights
# endif
#endif
      END MODULE nesting_mod