mod_nesting.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE mod_nesting
 
#ifdef NESTING
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This module defines structures for composite and refinement grids.  !
!                                                                      !
!  Composite Grids Structure: Donor grid data at contact points        !
!  =========================                                           !
!                                                                      !
!  bustr     Kinematic bottom momentum flux (bottom stress) in the     !
!              XI-direction (m2/s2)  at U-points.                      !
!  bvstr     Kinematic bottom momentum flux (bottom stress) in the     !
!              ETA-direction (m2/s2) at V-points.                      !
!  rzeta     Right-hand-side of free surface equation (m3/s).          !
!  ubar      Vertically integrated U-momentum component (m/s).         !
!  vbar      Vertically integrated V-momentum component (m/s).         !
!  zeta      Free surface (m).                                         !
# ifdef SOLVE3D
!                                                                      !
!  DU_avg1   Time averaged U-flux for 2D equations (m3/s).             !
!  DV_avg1   Time averaged V-flux for 2D equations (m3/s).             !
!  Huon      Total U-momentum flux term, Hz*u/pn.                      !
!  Hvom      Total V-momentum flux term, Hz*v/pm.                      !
!  Zt_avg1   Free-surface averaged over all short time-steps (m).      !
!  t         Tracer type variables (active and passive).               !
!  u         3D U-momentum component (m/s).                            !
!  v         3D U-momentum component (m/s).                            !
# endif
!                                                                      !
!  REFINED Grids Structure: Donor grid data at contact points          !
!  =======================  (two-time rolling snapshots)               !
!                                                                      !
!  ubar      Vertically integrated U-momentum component (m/s).         !
!  vbar      Vertically integrated V-momentum component (m/s).         !
!  zeta      Free surface (m).                                         !
# ifndef SOLVE3D
!  U2d_flux  U-flux for 2D equations transport (m3/s).                 !
!  V2d_flux  V-flux for 2D equations transport (m3/s).                 !
# else
!                                                                      !
!  U2d_flux  Time averaged U-flux for 3D equations coupling (m3/s).    !
!  V2d_flux  Time averaged V-flux for 3D equations coupling (m3/s).    !
!  t         Tracer type variables (active and passive).               !
!  u         3D U-momentum component (m/s).                            !
!  v         3D U-momentum component (m/s).                            !
# endif
!                                                                      !
!=======================================================================
!
      USE mod_kinds
!
      implicit none
!
      PUBLIC :: allocate_nesting
      PUBLIC :: deallocate_nesting
      PUBLIC :: initialize_nesting
!
!-----------------------------------------------------------------------
!  Nesting identification index of variables to process.
!-----------------------------------------------------------------------
!
!  The following identification indices are used in "initial" or
!  "main2d/main3d" to specify the variables that are processed in
!  each sub-timestep section.  Negative indices are used in grid
!  refinement whereas positive indices are used in composite grids.
!
      integer, parameter :: nmflx = -6    ! check mass flux conservation
      integer, parameter :: ndxdy = -5    ! extract on_u and om_v
      integer, parameter :: ngetd = -4    ! extract donor grid data
      integer, parameter :: nmask = -3    ! scale interpolation weights
      integer, parameter :: nputd = -2    ! fill contact points
      integer, parameter :: n2way = -1    ! fine to course coupling
!
      integer, parameter :: nfsic =  1    ! free surface initialization
      integer, parameter :: n2dic =  2    ! 2D momentum initialization
      integer, parameter :: n3dic =  3    ! 3D momentum initialization
      integer, parameter :: ntvic =  4    ! tracers initialization
      integer, parameter :: nbstr =  5    ! bottom stress (bustr,bvstr)
      integer, parameter :: nrhst =  6    ! RHS terms (tracers)
      integer, parameter :: nzeta =  7    ! 3D kernel free-surface
      integer, parameter :: nzwgt =  8    ! 3D vertical weights
      integer, parameter :: n2dps =  9    ! 2D engine Predictor Step
      integer, parameter :: n2dcs = 10    ! 2D engine Corrector Step
      integer, parameter :: n2dfx = 11    ! time-averaged 2D fluxes
      integer, parameter :: n3duv = 12    ! 3D momentum and fluxes
      integer, parameter :: n3dtv = 13    ! 3D tracer variables
!
!-----------------------------------------------------------------------
!  Nesting parameters.
!-----------------------------------------------------------------------
!
!  Nested grid connectivity switches. It is used to determine the
!  dimensions of the numerical kernel allocatable arrays. The arrays
!  have extra points due to the contact regions in any of sides
!  of the physical grid (1=iwest, 2=isouth, 3=ieast, 4=inorth).
!
      logical, allocatable :: contactregion(:,:)    ! [4,Ngrids]
!
!  Logical switch indicating which coarser grid is a donor to a
!  finer receiver grid (RefineScale(rg) > 0) external contact points.
!  This switch is in terms of the donor coarser grid.
!
      logical, allocatable :: donortofiner(:)       ! {Ngrids]
!$OMP THREADPRIVATE (DonorToFiner)
!
!  Switch indicating which refined grid(s), with RefineScale(ng) > 0,
!  include finer refined grids inside: telescoping refinement.
!
      logical, allocatable :: telescoping(:)        ! [Ngrids]
!$OMP THREADPRIVATE (Telescoping)
!
!  Switch to compute depth-dependent, vertical interpolation weights.
!  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 have 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 the
!  computations because of less distributed-memory communications.
!
      logical :: get_vweights
!$OMP THREADPRIVATE (get_Vweights)
!
!  If refinement, it contains the coarser donor grid number to finer
!  receiver grid external contact points.  The donor grid is always
!  coarser that receiver grid.  This variable is in terms of the
!  finer receiver grid.
!
      integer, allocatable :: coarserdonor(:)       ! [Ngrids]
!$OMP THREADPRIVATE (CoarserDonor)
!
!  If refinement and two-way exchange, it contains the donor finer
!  grid number to a coarser receiver grid.  This variable is in
!  terms of the coarser donor grid.
!
      integer, allocatable :: finerdonor(:)        ! [Ngrids]
!$OMP THREADPRIVATE (FinerDonor)
!
!  Number of refined time-steps. In most cases, the number of refined
!  time-step is the same as the refinement scale ratio for numerical
!  stability. However, the user is allowed to take larger divisible
!  time-step with respect to the donor grid. The variable below is
!  computed donor and receiver time-step ratio from standard input.
!  It is up to the user to determine the appropiate time-step for
!  stability.
!
      integer, allocatable :: refinesteps(:)        ! [Ngrids]
!$OMP THREADPRIVATE (RefineSteps)
!
!  Refined time-steps counter with respect the coarse grid (ng=1)
!  single time-step.
!
      integer, allocatable :: refinestepscounter(:) ! [Ngrids]
!
!  Interval used in the two-way exchange between fine and coarse
!  grids.
!
      real(r8), allocatable :: twowayinterval(:)    ! [Ngrids]
!
!  Donor and reciver grids for each contact region. These paremeters
!  are also duplicated in the T_NGC structure.
!
      integer, allocatable :: donor_grid(:)         ! [Ncontact]
      integer, allocatable :: receiver_grid(:)      ! [Ncontact]
!
!  Rolling index and time (seconds) used in the temporal interpolation
!  of contact point data.
!
      integer,  allocatable :: rollingindex(:)      ! [Ncontact]
      real(dp), allocatable :: rollingtime(:,:)     ! [Ncontact]
!$OMP THREADPRIVATE (RollingIndex, RollingTime)
!
!  If refinement, donor grid (I,J) indices at PSI points used to extract
!  refined grid. Values are set to -999 if not applicable.
!
!          +---------+    J_top
!          |         |
!          | Refined |
!          |    grid |
!          |         |
!          +---------+    J_bottom
!       I_left   I_right
!
      integer, allocatable :: i_left(:)             ! [Ngrids]
      integer, allocatable :: i_right(:)            ! [Ngrids]
      integer, allocatable :: j_bottom(:)           ! [Ngrids]
      integer, allocatable :: j_top(:)              ! [Ngrids]
!
!  Compact arrays used to unpack data from nested grids contact points
!  NetCDF file. They are allocated to the size "datum" dimension in
!  routine "set_contact".  The start and end indices for each C-type
!  variable are used to unpack from compact vector.
!
      integer :: ncdatum
      integer, allocatable :: ncpoints(:)           ! [Ncontact]
      integer, allocatable :: nstrr(:), nendr(:)    ! [Ncontact]
      integer, allocatable :: nstru(:), nendu(:)    ! [Ncontact]
      integer, allocatable :: nstrv(:), nendv(:)    ! [Ncontact]
!
      integer, allocatable :: contact_region(:)     ! [NCdatum]
      integer, allocatable :: on_boundary(:)        ! [NCdatum]
      integer, allocatable :: idg_cp(:)             ! [NCdatum]
      integer, allocatable :: jdg_cp(:)             ! [NCdatum]
      integer, allocatable :: irg_cp(:)             ! [NCdatum]
      integer, allocatable :: jrg_cp(:)             ! [NCdatum]
!
!-----------------------------------------------------------------------
!  Nested grid connectivity (NGC) structure.
!-----------------------------------------------------------------------
!
!  This structure is used to store all the connectivity information
!  between nested grids.  It will be used extensively when processing
!  contact region points between data donor and data receiver grids.
!  The nested grid contact region information is processed outside of
!  ROMS and read from a NetCDF file for functionality and efficiency.
!
!  In nested grids, the value in the contact region are interpolated
!  from the data donor grid cell using the following conventions at
!  the horizontal location in receiver grid (Irg,Jrg) and donor grid
!  cell (Idg,Jdg):
!
!  suffix 'dg' = donor grid
!         'rg' = receiver grid
!
!            4---------------3 (Idg+1,Jdg+1)   weight(1) = (1-p) * (1-q)
!            |        .      |                 weight(2) =    p  * (1-q)
!            |    1-q .      |                 weight(3) =    p  *  q
!            |        .      |                 weight(4) = (1-p) *  q
!            |        .   p  |
!        Jrg |....... x .....|  Linear interpolation:
!            |  1-p   .      |
!            |        . q    |  V(Irg,Jrg) = weight(1) * F(Idg  ,Jdg  )+
!            |        .      |               weight(2) * F(Idg+1,Jdg  )+
!  (Idg,Jdg) 1---------------2               weight(3) * F(Idg+1,Jdg+1)+
!                    Irg                     weight(4) * F(Idg  ,Jdg+1)
!
!  Notice that if p=0 and q=0 at all contact points, the donor and
!  receiver grids are coincident since weight(1)=1.0 and weight(2:3)=0.
!  Therefore, the above formula is generic for any nested grid
!  configuration.
!
!  If Land/Sea masking, the interpolation weights are rescaled in
!  "mask_weights" during initialization to account masked points in
!  the contact regions. If wetting and drying, the rescaling is done
!  at every time step since the land/sea masking is time dependent.
!
      integer :: ncontact          ! total number of contact regions
!
      TYPE t_ngc
        logical :: coincident      ! coincident donor/receiver, p=q=0
        logical :: interpolate     ! perform vertical interpolation
        integer :: donor_grid      ! data donor grid number
        integer :: receiver_grid   ! data receiver grid number
        integer :: npoints         ! number of points in contact region
        integer, pointer :: irg(:) ! receiver grid, I-contact point
        integer, pointer :: jrg(:) ! receiver grid, J-contact point
        integer, pointer :: idg(:) ! donor grid, cell I-left   index
        integer, pointer :: jdg(:) ! donor grid, cell J-bottom index
# ifdef SOLVE3D
        integer, pointer :: kdg(:,:)         ! donor grid, cell K-index
# endif
        real(r8), pointer :: lweight(:,:)          ! linear weights
# ifdef WET_DRY
        real(r8), pointer :: lweightunmasked(:,:)  ! Unmasked Lweight
# endif
# ifdef QUADRATIC_WEIGHTS
        real(r8), pointer :: qweight(:,:)          ! quadratic weights
#  ifdef WET_DRY
        real(r8), pointer :: qweightunmasked(:,:)  ! Unmasked Qweight
#  endif
# endif
# ifdef SOLVE3D
        real(r8), pointer :: vweight(:,:,:)        ! vertical weights
# endif
      END TYPE t_ngc
!
      TYPE (t_ngc), allocatable :: rcontact(:)  ! RHO-points, [Ncontact]
      TYPE (t_ngc), allocatable :: ucontact(:)  ! U-points,   [Ncontact]
      TYPE (t_ngc), allocatable :: vcontact(:)  ! V-points,   [Ncontact]
!
!-----------------------------------------------------------------------
!  Boundary Contact Points (BCP) structure, allocated as (4,Ncontact).
!  The first dimension is for domain edge (1=iwest,2=isouth, 3=ieast,
!  4=inorth).
!-----------------------------------------------------------------------
!
!  Currently, this structure is only used in refinement grids where the
!  coarser (donor) and finer (receiver) grids have coincident boundaries
!  but with different I- and J-indices.  However, it can be used in the
!  future for composite grids with coincient boundaries.
!
!  The variable "C2Bindex" is used to tell us which contact points in
!  the "Ucontact" and "Vcontact" structure are located at the physical
!  boundary of the relevant nested grid. For example at the boundary
!  edge of a grid with contact region "cr", we can get the mapping
!  between contact point "m" and grid physical boundary edge index
!  "i" or "j" as:
!
!     m = BRY_CONTACT(iwest, cr) % C2Bindex(j)
!     m = BRY_CONTACT(isouth,cr) % C2Bindex(i)
!     m = BRY_CONTACT(ieast, cr) % C2Bindex(j)
!     m = BRY_CONTACT(inorth,cr) % C2Bindex(i)
!
!  This mapping is set during intialization and facilitates efficient
!  processing of nesting contact data.
!
      TYPE t_bcp
        integer :: spv                     ! fill value, unwanted index
        integer :: ibmin                   ! viable minimum Ib
        integer :: ibmax                   ! viable maximum Ib
        integer :: jbmin                   ! viable minimum Jb
        integer :: jbmax                   ! viable maximum Jb
        integer, pointer :: ib(:)          ! I-boundary index
        integer, pointer :: jb(:)          ! J-boundary index
 
        integer, pointer :: c2bindex(:)    ! contact to boundary index
 
# ifdef NESTING_DEBUG
        real(r8), pointer :: mflux(:)      ! perimeter mass flux
# endif
# ifdef SOLVE3D
        real(r8), pointer :: tflux(:,:,:)  ! perimeter tracer flux
# endif
      END TYPE t_bcp
!
      TYPE (t_bcp), allocatable :: bry_contact(:,:)    ! [4,Ncontact]
!
!-----------------------------------------------------------------------
!  Nested Grid Metrics (NGM) structure for contact regions.  Usually,
!  there are contact points outside of the regular (physical) nested
!  grid domain. That is, such contact points are located in the
!  extended (numerical) regions.  These metrics values are computed
!  when designing and generating the application grids.
!
!  It is recommended to build an intermediary fine resolution grid
!  encompassing the study area first and extract/sample all the ROMS
!  application nested grids from it.  This would give a better handle
!  on volume conservation, bathymetry, land/sea masking and other
!  issues.
!-----------------------------------------------------------------------
!
!  These metrics are written the contact points NetCDF file and save
!  separated here. It is very tricky to load these values directly
!  to global grid metrics because of parallelization.
!
       TYPE t_ngm
         real(r8), pointer :: angler(:)    ! angle between XI and EAST
         real(r8), pointer :: dndx(:)      ! d(1/pn)/d(XI)
         real(r8), pointer :: dmde(:)      ! d(1/pm)/d(ETA)
         real(r8), pointer :: f(:)         ! Coriolis parameter
         real(r8), pointer :: h(:)         ! bathymetry
         real(r8), pointer :: rmask(:)     ! land/sea RHO-mask
         real(r8), pointer :: umask(:)     ! land/sea U-mask
         real(r8), pointer :: vmask(:)     ! land/sea V-mask
         real(r8), pointer :: pm(:)        ! XI-coordinate metric
         real(r8), pointer :: pn(:)        ! ETA-coordinate metric
         real(r8), pointer :: xr(:)        ! X RHO-coordinate (m or deg)
         real(r8), pointer :: yr(:)        ! Y RHO-coordinate (m or deg)
         real(r8), pointer :: xu(:)        ! X U-coordinate (m or deg)
         real(r8), pointer :: yu(:)        ! Y U-coordinate (m or deg)
         real(r8), pointer :: xv(:)        ! X V-coordinate (m or deg)
         real(r8), pointer :: yv(:)        ! Y V-coordinate (m or deg)
       END TYPE t_ngm
!
       TYPE (t_ngm), allocatable :: contact_metric(:)    ! [Ncontact]
!
!-----------------------------------------------------------------------
!  Composite grids structure.  It contains the donor grid data at the
!  receiver grid contact points. The donor grid data is extracted for
!  the cell containing the contact point: 4 horizontal values to
!  facilitate spatial interpolation.
!-----------------------------------------------------------------------
!
      TYPE t_composite
        real(r8), pointer :: bustr(:,:)       ! [4,Npoints]
        real(r8), pointer :: bvstr(:,:)       ! [4,Npoints)
 
        real(r8), pointer :: ubar(:,:,:)      ! [4,Npoints,2]
        real(r8), pointer :: vbar(:,:,:)      ! [4,Npoints,2]
        real(r8), pointer :: zeta(:,:,:)      ! [4,Npoints,2]
 
        real(r8), pointer :: rzeta(:,:)       ! [4,Npoints]
 
# ifdef SOLVE3D
        real(r8), pointer :: du_avg1(:,:)     ! [4,Npoints]
        real(r8), pointer :: dv_avg1(:,:)     ! [4,Npoints]
        real(r8), pointer :: zt_avg1(:,:)     ! [4,Npoints]
 
        real(r8), pointer :: u(:,:,:)         ! [4,k,Npoints]
        real(r8), pointer :: v(:,:,:)         ! [4,k,Npoints]
 
        real(r8), pointer :: huon(:,:,:)      ! [4,k,Npoints]
        real(r8), pointer :: hvom(:,:,:)      ! [4,k,Npoints]
 
        real(r8), pointer :: t(:,:,:,:)       ! [4,k,Npoints,itrc]
# endif
 
      END TYPE t_composite
!
      TYPE (t_composite), allocatable :: composite(:)  ! [Ncontact]
!
!-----------------------------------------------------------------------
!  Refinement grids structure: It contains the coarser grid data at the
!  finer grid contact points. The finer grid data is extracted for the
!  cell containing the contact point: 4 horizontal values and 2 time
!  records (t1:t2) to facilitate the space-time interpolation.
!-----------------------------------------------------------------------
!
      TYPE t_refined
        real(r8), pointer :: ubar(:,:,:)      ! [4,Npoints,t1:t2]
        real(r8), pointer :: vbar(:,:,:)      ! [4,Npoints,t1:t2]
        real(r8), pointer :: zeta(:,:,:)      ! [4,Npoints,t1:t2]
 
        real(r8), pointer :: u2d_flux(:,:,:)  ! [4,Npoints,t1:t2]
        real(r8), pointer :: v2d_flux(:,:,:)  ! [4,Npoints,t1:t2]
 
        real(r8), pointer :: on_u(:)          ! [Npoints]
        real(r8), pointer :: om_v(:)          ! [Npoints]
 
# ifdef SOLVE3D
        real(r8), pointer :: u(:,:,:,:)       ! [4,k,Npoints,t1:t2]
        real(r8), pointer :: v(:,:,:,:)       ! [4,k,Npoints,t1:t2]
 
        real(r8), pointer :: t(:,:,:,:,:)     ! [4,k,Npoints,t1:t2,itrc]
# endif
      END TYPE t_refined
!
      TYPE (t_refined), allocatable :: refined(:)      ! [Ncontact]
!
      CONTAINS
!
      SUBROUTINE allocate_nesting
!
!=======================================================================
!                                                                      !
!  This routine allocates and initializes nesting structure for 2D     !
!  state variables.                                                    !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_boundary
      USE mod_scalars
!
!  Local variable declarations.
!
      integer :: lbi, ubi, lbj, ubj
      integer :: imin, imax, jmin, jmax
      integer :: ccr, cr, dg, ng, rg
      integer :: i, ibry, ic, id, ir, j, jd, jr, m, my_tile
      integer :: ispval
 
      integer, allocatable :: ibmin(:,:), ibmax(:,:)
      integer, allocatable :: jbmin(:,:), jbmax(:,:)
!
!-----------------------------------------------------------------------
!  Unpack Boundary Contact Points structure (type T_BCP).
!-----------------------------------------------------------------------
!
!  Allocate boundary connectivity (type T_BCP) structure.
!
      allocate ( bry_contact(4,ncontact) )
!
!  Allocate arrays in boundary connectivity structure.
!
      my_tile=-1                           ! for global values
      DO cr=1,ncontact
        rg=receiver_grid(cr)
        lbi=bounds(rg)%LBi(my_tile)
        ubi=bounds(rg)%UBi(my_tile)
        lbj=bounds(rg)%LBj(my_tile)
        ubj=bounds(rg)%UBj(my_tile)
        DO ibry=1,4
          SELECT CASE (ibry)
            CASE (iwest, ieast)
              allocate ( bry_contact(ibry,cr) % Ib(lbj:ubj) )
              dmem(rg)=dmem(rg)+real(ubj-lbj,r8)
 
              allocate ( bry_contact(ibry,cr) % Jb(lbj:ubj) )
              dmem(rg)=dmem(rg)+real(ubj-lbj,r8)
 
              allocate ( bry_contact(ibry,cr) % C2Bindex(lbj:ubj) )
              dmem(rg)=dmem(rg)+real(ubj-lbj,r8)
 
# ifdef NESTING_DEBUG
              allocate ( bry_contact(ibry,cr) % Mflux(lbj:ubj) )
              dmem(rg)=dmem(rg)+real(ubj-lbj,r8)
# endif
 
# ifdef SOLVE3D
              allocate ( bry_contact(ibry,cr) % Tflux(lbj:ubj,          &
     &                                                n(rg),nt(rg)) )
              dmem(rg)=dmem(rg)+real((ubj-lbj)*n(rg)*nt(rg),r8)
# endif
            CASE (isouth, inorth)
              allocate ( bry_contact(ibry,cr) % Ib(lbi:ubi) )
              dmem(rg)=dmem(rg)+real(ubi-lbi,r8)
 
              allocate ( bry_contact(ibry,cr) % Jb(lbi:ubi) )
              dmem(rg)=dmem(rg)+real(ubi-lbi,r8)
 
              allocate ( bry_contact(ibry,cr) % C2Bindex(lbi:ubi) )
              dmem(rg)=dmem(rg)+real(ubi-lbi,r8)
 
# ifdef NESTING_DEBUG
              allocate ( bry_contact(ibry,cr) % Mflux(lbi:ubi) )
              dmem(rg)=dmem(rg)+real(ubi-lbi,r8)
# endif
 
# ifdef SOLVE3D
              allocate ( bry_contact(ibry,cr) % Tflux(lbi:ubi,          &
     &                                                n(rg),nt(rg)) )
              dmem(rg)=dmem(rg)+real((ubi-lbi)*n(rg)*nt(rg),r8)
# endif
          END SELECT
        END DO
      END DO
!
!  Initialize boundary connectivity structure: Boundary indices array
!  are initialized to its special value.
!
      ispval=-999
!
      IF (.not.allocated(ibmin)) THEN
        allocate ( ibmin(4,ncontact) )
        DO ng=1,ngrids
          dmem(ng)=dmem(ng)+4.0_r8*real(ncontact,r8)
        END DO
        ibmin = -ispval
      END IF
      IF (.not.allocated(ibmax)) THEN
        allocate ( ibmax(4,ncontact) )
        DO ng=1,ngrids
          dmem(ng)=dmem(ng)+4.0_r8*real(ncontact,r8)
        END DO
        ibmax = ispval
      END IF
      IF (.not.allocated(jbmin)) THEN
        allocate ( jbmin(4,ncontact) )
        DO ng=1,ngrids
          dmem(ng)=dmem(ng)+4.0_r8*real(ncontact,r8)
        END DO
        jbmin = -ispval
      END IF
      IF (.not.allocated(jbmax)) THEN
        allocate ( jbmax(4,ncontact) )
        DO ng=1,ngrids
          dmem(ng)=dmem(ng)+4.0_r8*real(ncontact,r8)
        END DO
        jbmax = ispval
      END IF
!
      DO cr=1,ncontact
        DO ibry=1,4
          bry_contact(ibry,cr) % spv = ispval
          bry_contact(ibry,cr) % Ib = ispval
          bry_contact(ibry,cr) % Jb = ispval
          bry_contact(ibry,cr) % C2Bindex = ispval
# ifdef NESTING_DEBUG
          bry_contact(ibry,cr) % Mflux = 0.0_r8
# endif
# ifdef SOLVE3D
          bry_contact(ibry,cr) % Tflux = 0.0_r8
# endif
        END DO
      END DO
!
!  Identify contact points located on the grid boundary.  Notice that
!  the conjugate contact region (CCR) is also processed but it is not
!  yet used. Also, the CCR indices (Ib,Jb) are in refinement overwriten
!  in the m-loop below because several finer grid contact points
!  (RefineScale) are contained in the coarser grid cell. The C2Bindex
!  in this case has the value for the last processed contact point
!  with contact region "cr".
!
      DO m=1,ncdatum
        cr=contact_region(m)
        dg=donor_grid(cr)
        rg=receiver_grid(cr)
        ibry=on_boundary(m)
        DO ic=1,ncontact
          IF ((dg.eq.receiver_grid(ic)).and.                            &
     &        (rg.eq.donor_grid(ic))) THEN
            ccr=ic                            ! conjugate contact region
            EXIT
          END IF
        END DO
        IF ((ibry.eq.iwest ).or.(ibry.eq.ieast )) THEN
          ir=irg_cp(m)
          jr=jrg_cp(m)
          ibmin(ibry,cr )=min(ir,ibmin(ibry,cr ))
          ibmax(ibry,cr )=max(ir,ibmax(ibry,cr ))
          jbmin(ibry,cr )=min(jr,jbmin(ibry,cr ))
          jbmax(ibry,cr )=max(jr,jbmax(ibry,cr ))
          bry_contact(ibry,cr ) % Ib(jr) = ir
          bry_contact(ibry,cr ) % Jb(jr) = jr
          bry_contact(ibry,cr ) % C2Bindex(jr) = m-nstru(cr)+1
!
          id=idg_cp(m)
          jd=jdg_cp(m)
          ibmin(ibry,ccr)=min(id,ibmin(ibry,ccr))
          ibmax(ibry,ccr)=max(id,ibmax(ibry,ccr))
          jbmin(ibry,ccr)=min(jd,jbmin(ibry,ccr))
          jbmax(ibry,ccr)=max(jd,jbmax(ibry,ccr))
          bry_contact(ibry,ccr) % Ib(jd) = id
          bry_contact(ibry,ccr) % Jb(jd) = jd
          bry_contact(ibry,ccr) % C2Bindex(jd) = m-nstru(cr)+1  ! same
        ELSE IF ((ibry.eq.isouth).or.(ibry.eq.inorth)) THEN
          ir=irg_cp(m)
          jr=jrg_cp(m)
          ibmin(ibry,cr )=min(ir,ibmin(ibry,cr))
          ibmax(ibry,cr )=max(ir,ibmax(ibry,cr))
          jbmin(ibry,cr )=min(jr,jbmin(ibry,cr))
          jbmax(ibry,cr )=max(jr,jbmax(ibry,cr))
          bry_contact(ibry,cr ) % Ib(ir) = ir
          bry_contact(ibry,cr ) % Jb(ir) = jr
          bry_contact(ibry,cr ) % C2Bindex(ir) = m-nstrv(cr)+1
!
          id=idg_cp(m)
          jd=jdg_cp(m)
          ibmin(ibry,ccr)=min(id,ibmin(ibry,ccr))
          ibmax(ibry,ccr)=max(id,ibmax(ibry,ccr))
          jbmin(ibry,ccr)=min(jd,jbmin(ibry,ccr))
          jbmax(ibry,ccr)=max(jd,jbmax(ibry,ccr))
          bry_contact(ibry,ccr) % Ib(id) = id
          bry_contact(ibry,ccr) % Jb(id) = jd
          bry_contact(ibry,ccr) % C2Bindex(id) = m-nstrv(cr)+1  ! same
        END IF
      END DO
!
!  Set minimum and maximum indices to process at each boundary.
!
      DO cr=1,ncontact
        DO ibry=1,4
          IF (abs(ibmin(ibry,cr)).eq.abs(ispval)) THEN
            bry_contact(ibry,cr) % Ibmin = ispval
          ELSE
            bry_contact(ibry,cr) % Ibmin = ibmin(ibry,cr)
          END IF
 
          IF (abs(ibmax(ibry,cr)).eq.abs(ispval)) THEN
            bry_contact(ibry,cr) % Ibmax = ispval
          ELSE
            bry_contact(ibry,cr) % Ibmax = ibmax(ibry,cr)
          END IF
 
          IF (abs(jbmin(ibry,cr)).eq.abs(ispval)) THEN
            bry_contact(ibry,cr) % Jbmin = ispval
          ELSE
            bry_contact(ibry,cr) % Jbmin = jbmin(ibry,cr)
          END IF
 
          IF (abs(jbmax(ibry,cr)).eq.abs(ispval)) THEN
            bry_contact(ibry,cr) % Jbmax = ispval
          ELSE
            bry_contact(ibry,cr) % Jbmax = jbmax(ibry,cr)
          END IF
        END DO
      END DO
!
!-----------------------------------------------------------------------
!  Deactivate boundary condition switches if contact point lay on the
!  physical nested grid boundary.
!-----------------------------------------------------------------------
!
      DO cr=1,ncontact
        rg=receiver_grid(cr)
        IF (refinedgrid(rg)) THEN
          IF (refinescale(rg).gt.0) THEN
            lbc_apply(rg) % west  = .false.    ! This is a refinement
            lbc_apply(rg) % south = .false.    ! grid, so we need to
            lbc_apply(rg) % east  = .false.    ! avoid applying lateral
            lbc_apply(rg) % north = .false.    ! boundary conditions
          END IF
        ELSE
          DO ibry=1,4
            imin=bry_contact(ibry,cr) % Ibmin  ! Deactivate full or
            imax=bry_contact(ibry,cr) % Ibmax  ! partial lateral
            jmin=bry_contact(ibry,cr) % Jbmin  ! boundary conditions
            jmax=bry_contact(ibry,cr) % Jbmax
            SELECT CASE (ibry)
              CASE (iwest)
                IF ((jmin.ne.ispval).and.(jmax.ne.ispval)) THEN
                  DO j=jmin,jmax
                    lbc_apply(rg) % west (j) = .false.
                  END DO
                END IF
              CASE (isouth)
                IF ((imin.ne.ispval).and.(imax.ne.ispval)) THEN
                  DO i=imin,imax
                    lbc_apply(rg) % south(i) = .false.
                  END DO
                END IF
              CASE (ieast)
                IF ((jmin.ne.ispval).and.(jmax.ne.ispval)) THEN
                  DO j=jmin,jmax
                    lbc_apply(rg) % east (j) = .false.
                  END DO
                END IF
              CASE (inorth)
                IF ((imin.ne.ispval).and.(imax.ne.ispval)) THEN
                  DO i=imin,imax
                    lbc_apply(rg) % north(i) = .false.
                  END DO
                END IF
            END SELECT
          END DO
        END IF
      END DO
!
      RETURN
      END SUBROUTINE allocate_nesting
!
      SUBROUTINE deallocate_nesting
!
!=======================================================================
!                                                                      !
!  This routine allocates and initializes nesting structure for 2D     !
!  state variables.                                                    !
!                                                                      !
!=======================================================================
 
# ifdef SUBOBJECT_DEALLOCATION
!
      USE destroy_mod, ONLY : destroy
# endif
!
!  Local variable declarations.
!
      integer ::  cr, ibry
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", deallocate_nesting"
 
# ifdef SUBOBJECT_DEALLOCATION
!
!-----------------------------------------------------------------------
!  Deallocate each variable in the derived-type T_BCP boundary
!  connectivity structure separately.
!-----------------------------------------------------------------------
!
      DO cr=1,ncontact
        DO ibry=1,4
          IF (.not.destroy(ng, bry_contact(ibry,cr)%Ib, myfile,         &
     &                     __line__, 'BRY_CONTACT(ibry,cr)%Ib')) RETURN
 
          IF (.not.destroy(ng, bry_contact(ibry,cr)%Jb, myfile,         &
     &                     __line__, 'BRY_CONTACT(ibry,cr)%Jb')) RETURN
 
          IF (.not.destroy(ng, bry_contact(ibry,cr)%C2Bindex, myfile,   &
     &                __line__, 'BRY_CONTACT(ibry,cr)%C2Bindex')) RETURN
 
#  ifdef NESTING_DEBUG
          IF (.not.destroy(ng, bry_contact(ibry,cr)%Mflux, myfile,      &
     &                __line__, 'BRY_CONTACT(ibry,cr)%Mflux')) RETURN
#  endif
 
#  ifdef SOLVE3D
          IF (.not.destroy(ng, bry_contact(ibry,cr)%Tflux, myfile,      &
     &                __line__, 'BRY_CONTACT(ibry,cr)%Tflux')) RETURN
#  endif
        END DO
      END DO
!
!-----------------------------------------------------------------------
!  Deallocate each variable in the derived-type T_NGC grid connectivity
!  structure separately.
!-----------------------------------------------------------------------
!
      DO cr=1,ncontact
        IF (.not.destroy(ng, rcontact(cr)%Irg, myfile,                  &
     &                   __line__, 'Rcontact(cr)%Irg')) RETURN
        IF (.not.destroy(ng, ucontact(cr)%Irg, myfile,                  &
     &                   __line__, 'Ucontact(cr)%Irg')) RETURN
        IF (.not.destroy(ng, vcontact(cr)%Irg, myfile,                  &
     &                   __line__, 'Vcontact(cr)%Irg')) RETURN
!
        IF (.not.destroy(ng, rcontact(cr)%Jrg, myfile,                  &
     &                   __line__, 'Rcontact(cr)%Jrg')) RETURN
        IF (.not.destroy(ng, ucontact(cr)%Jrg, myfile,                  &
     &                   __line__, 'Ucontact(cr)%Jrg')) RETURN
        IF (.not.destroy(ng, vcontact(cr)%Jrg, myfile,                  &
     &                   __line__, 'Vcontact(cr)%Jrg')) RETURN
!
        IF (.not.destroy(ng, rcontact(cr)%Idg, myfile,                  &
     &                   __line__, 'Rcontact(cr)%Idg')) RETURN
        IF (.not.destroy(ng, ucontact(cr)%Idg, myfile,                  &
     &                   __line__, 'Ucontact(cr)%Idg')) RETURN
        IF (.not.destroy(ng, vcontact(cr)%Idg, myfile,                  &
     &                   __line__, 'Vcontact(cr)%Idg')) RETURN
!
        IF (.not.destroy(ng, rcontact(cr)%Jdg, myfile,                  &
     &                   __line__, 'Rcontact(cr)%Jdg')) RETURN
        IF (.not.destroy(ng, ucontact(cr)%Jdg, myfile,                  &
     &                   __line__, 'Ucontact(cr)%Jdg')) RETURN
        IF (.not.destroy(ng, vcontact(cr)%Jdg, myfile,                  &
     &                   __line__, 'Vcontact(cr)%Jdg')) RETURN
 
#  ifdef SOLVE3D
!
        IF (.not.destroy(ng, rcontact(cr)%Kdg, myfile,                  &
     &                   __line__, 'Rcontact(cr)%Kdg')) RETURN
        IF (.not.destroy(ng, ucontact(cr)%Kdg, myfile,                  &
     &                   __line__, 'Ucontact(cr)%Kdg')) RETURN
        IF (.not.destroy(ng, vcontact(cr)%Kdg, myfile,                  &
     &                   __line__, 'Vcontact(cr)%Kdg')) RETURN
#  endif
!
        IF (.not.destroy(ng, rcontact(cr)%Lweight, myfile,              &
     &                   __line__, 'Rcontact(cr)%Lweight')) RETURN
        IF (.not.destroy(ng, ucontact(cr)%Lweight, myfile,              &
     &                   __line__, 'Ucontact(cr)%Lweight')) RETURN
        IF (.not.destroy(ng, vcontact(cr)%Lweight, myfile,              &
     &                   __line__, 'Vcontact(cr)%Lweight')) RETURN
!
#  ifdef WET_DRY
        IF (.not.destroy(ng, rcontact(cr)%LweightUnmasked, myfile,      &
     &              __line__, 'Rcontact(cr)%LweightUnmasked')) RETURN
        IF (.not.destroy(ng, ucontact(cr)%LweightUnmasked, myfile,      &
     &              __line__, 'Ucontact(cr)%LweightUnmasked')) RETURN
        IF (.not.destroy(ng, vcontact(cr)%LweightUnmasked, myfile,      &
     &              __line__, 'Vcontact(cr)%LweightUnmasked')) RETURN
#  endif
 
#  ifdef QUADRATIC_WEIGHTS
!
        IF (.not.destroy(ng, rcontact(cr)%Qweight, myfile,              &
     &                   __line__, 'Rcontact(cr)%Qweight')) RETURN
        IF (.not.destroy(ng, ucontact(cr)%Qweight, myfile,              &
     &                   __line__, 'Ucontact(cr)%Qweight')) RETURN
        IF (.not.destroy(ng, vcontact(cr)%Qweight, myfile,              &
     &                   __line__, 'Vcontact(cr)%Qweight')) RETURN
 
#   ifdef WET_DRY
!
        IF (.not.destroy(ng, rcontact(cr)%QweightUnmasked, myfile,      &
     &              __line__, 'Rcontact(cr)%QweightUnmasked')) RETURN
        IF (.not.destroy(ng, ucontact(cr)%QweightUnmasked, myfile,      &
     &              __line__, 'Ucontact(cr)%QweightUnmasked')) RETURN
        IF (.not.destroy(ng, vcontact(cr)%QweightUnmasked, myfile,      &
     &              __line__, 'Vcontact(cr)%QweightUnmasked')) RETURN
#   endif
#  endif
 
#  ifdef SOLVE3D
!
        IF (.not.destroy(ng, rcontact(cr)%Vweight, myfile,              &
     &                   __line__, 'Rcontact(cr)%Vweight')) RETURN
        IF (.not.destroy(ng, ucontact(cr)%Vweight, myfile,              &
     &                   __line__, 'Ucontact(cr)%Vweight')) RETURN
        IF (.not.destroy(ng, vcontact(cr)%Vweight, myfile,              &
     &                   __line__, 'Vcontact(cr)%Vweight')) RETURN
 
#   if defined TANGENT || defined TL_IOMS
!
        IF (.not.destroy(ng, rcontact(cr)%tl_Vweight, myfile,           &
     &                   __line__, 'Rcontact(cr)%tl_Vweight')) RETURN
        IF (.not.destroy(ng, ucontact(cr)%tl_Vweight, myfile,           &
     &                   __line__, 'Ucontact(cr)%tl_Vweight')) RETURN
        IF (.not.destroy(ng, vcontact(cr)%tl_Vweight, myfile,           &
     &                   __line__, 'Vcontact(cr)%tl_Vweight')) RETURN
#   endif
 
#   ifdef ADJOINT
!
        IF (.not.destroy(ng, rcontact(cr)%ad_Vweight, myfile,           &
     &                   __line__, 'Rcontact(cr)%ad_Vweight')) RETURN
        IF (.not.destroy(ng, ucontact(cr)%ad_Vweight, myfile,           &
     &                   __line__, 'Ucontact(cr)%ad_Vweight')) RETURN
        IF (.not.destroy(ng, vcontact(cr)%ad_Vweight, myfile,           &
     &                   __line__, 'Vcontact(cr)%ad_Vweight')) RETURN
#   endif
#  endif
      END DO
!
!-----------------------------------------------------------------------
!  Deallocate each variable in the derived-type T_NGC contact region
!  metrics structure separately.
!-----------------------------------------------------------------------
!
      DO cr=1,ncontact
        IF (.not.destroy(ng, contact_metric(cr)%angler, myfile,         &
     &                   __line__, 'CONTACT_METRIC(cr)%angler')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%dndx, myfile,           &
     &                   __line__, 'CONTACT_METRIC(cr)%dndx')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%dmde, myfile,           &
     &                   __line__, 'CONTACT_METRIC(cr)%dmde')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%f, myfile,              &
     &                   __line__, 'CONTACT_METRIC(cr)%f')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%h, myfile,              &
     &                   __line__, 'CONTACT_METRIC(cr)%h')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%rmask, myfile,          &
     &                   __line__, 'CONTACT_METRIC(cr)%rmask')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%umask, myfile,          &
     &                   __line__, 'CONTACT_METRIC(cr)%umask')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%vmask, myfile,          &
     &                   __line__, 'CONTACT_METRIC(cr)%vmask')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%pm, myfile,             &
     &                   __line__, 'CONTACT_METRIC(cr)%pm')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%pn, myfile,             &
     &                   __line__, 'CONTACT_METRIC(cr)%pn')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%Xr, myfile,             &
     &                   __line__, 'CONTACT_METRIC(cr)%Xr')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%Yr, myfile,             &
     &                   __line__, 'CONTACT_METRIC(cr)%Yr')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%Xu, myfile,             &
     &                   __line__, 'CONTACT_METRIC(cr)%Xu')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%Yu, myfile,             &
     &                   __line__, 'CONTACT_METRIC(cr)%Yu')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%Xv, myfile,             &
     &                   __line__, 'CONTACT_METRIC(cr)%Xv')) RETURN
 
        IF (.not.destroy(ng, contact_metric(cr)%Yv, myfile,             &
     &                   __line__, 'CONTACT_METRIC(cr)%Yv')) RETURN
      END DO
!
!-----------------------------------------------------------------------
!  Deallocate each variable in the derived-type T_COMPOSITE for
!  composite grids contact region structure separately.
!-----------------------------------------------------------------------
!
      DO cr=1,ncontact
        IF (.not.destroy(ng, composite(cr)%bustr, myfile,               &
     &                   __line__, 'COMPOSITE(cr)%bustr')) RETURN
        IF (.not.destroy(ng, composite(cr)%bvstr, myfile,               &
     &                   __line__, 'COMPOSITE(cr)%bvstr')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%ubar, myfile,                &
     &                   __line__, 'COMPOSITE(cr)%ubar')) RETURN
        IF (.not.destroy(ng, composite(cr)%vbar, myfile,                &
     &                   __line__, 'COMPOSITE(cr)%vbar')) RETURN
        IF (.not.destroy(ng, composite(cr)%zeta, myfile,                &
     &                   __line__, 'COMPOSITE(cr)%zeta')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%rzeta, myfile,               &
     &                   __line__, 'COMPOSITE(cr)%rzeta')) RETURN
 
#   if defined TANGENT || defined TL_IOMS
        IF (.not.destroy(ng, composite(cr)%tl_bustr, myfile,            &
     &                   __line__, 'COMPOSITE(cr)%tl_bustr')) RETURN
        IF (.not.destroy(ng, composite(cr)%tl_bvstr, myfile,            &
     &                   __line__, 'COMPOSITE(cr)%tl_bvstr')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%tl_ubar, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%tl_ubar')) RETURN
        IF (.not.destroy(ng, composite(cr)%tl_vbar, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%tl_vbar')) RETURN
        IF (.not.destroy(ng, composite(cr)%tl_zeta, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%tl_zeta')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%tl_rzeta, myfile,            &
     &                   __line__, 'COMPOSITE(cr)%tl_rzeta')) RETURN
#   endif
 
#   ifdef ADJOINT
        IF (.not.destroy(ng, composite(cr)%ad_bustr, myfile,            &
     &                   __line__, 'COMPOSITE(cr)%ad_bustr')) RETURN
        IF (.not.destroy(ng, composite(cr)%ad_bvstr, myfile,            &
     &                   __line__, 'COMPOSITE(cr)%ad_bvstr')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%ad_ubar, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%ad_ubar')) RETURN
        IF (.not.destroy(ng, composite(cr)%ad_vbar, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%ad_vbar')) RETURN
        IF (.not.destroy(ng, composite(cr)%ad_zeta, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%ad_zeta')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%ad_rzeta, myfile,            &
     &                   __line__, 'COMPOSITE(cr)%ad_rzeta')) RETURN
#   endif
 
#  ifdef SOLVE3D
        IF (.not.destroy(ng, composite(cr)%DU_avg1, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%DU_avg1')) RETURN
        IF (.not.destroy(ng, composite(cr)%DV_avg1, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%DV_avg1')) RETURN
        IF (.not.destroy(ng, composite(cr)%Zt_avg1, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%Zt_avg1')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%u, myfile,                   &
     &                   __line__, 'COMPOSITE(cr)%u')) RETURN
        IF (.not.destroy(ng, composite(cr)%v, myfile,                   &
     &                   __line__, 'COMPOSITE(cr)%v')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%Huon, myfile,                &
     &                   __line__, 'COMPOSITE(cr)%Huon')) RETURN
        IF (.not.destroy(ng, composite(cr)%Hvom, myfile,                &
     &                   __line__, 'COMPOSITE(cr)%Hvom')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%t, myfile,                   &
     &                   __line__, 'COMPOSITE(cr)%t')) RETURN
 
#   if defined TANGENT || defined TL_IOMS
        IF (.not.destroy(ng, composite(cr)%tl_DU_avg1, myfile,          &
     &                   __line__, 'COMPOSITE(cr)%tl_DU_avg1')) RETURN
        IF (.not.destroy(ng, composite(cr)%tl_DV_avg1, myfile,          &
     &                   __line__, 'COMPOSITE(cr)%tl_DV_avg1')) RETURN
        IF (.not.destroy(ng, composite(cr)%tl_Zt_avg1, myfile,          &
     &                   __line__, 'COMPOSITE(cr)%tl_Zt_avg1')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%tl_u, myfile,                &
     &                   __line__, 'COMPOSITE(cr)%tl_u')) RETURN
        IF (.not.destroy(ng, composite(cr)%tl_v, myfile,                &
     &                   __line__, 'COMPOSITE(cr)%tl_v')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%tl_Huon, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%tl_Huon')) RETURN
        IF (.not.destroy(ng, composite(cr)%tl_Hvom, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%tl_Hvom')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%tl_t, myfile,                &
     &                   __line__, 'COMPOSITE(cr)%tl_t')) RETURN
#   endif
 
#   ifdef ADJOINT
        IF (.not.destroy(ng, composite(cr)%ad_DU_avg1, myfile,          &
     &                   __line__, 'COMPOSITE(cr)%ad_DU_avg1')) RETURN
        IF (.not.destroy(ng, composite(cr)%ad_DV_avg1, myfile,          &
     &                   __line__, 'COMPOSITE(cr)%ad_DV_avg1')) RETURN
        IF (.not.destroy(ng, composite(cr)%ad_Zt_avg1, myfile,          &
     &                   __line__, 'COMPOSITE(cr)%ad_Zt_avg1')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%ad_u, myfile,                &
     &                   __line__, 'COMPOSITE(cr)%ad_u')) RETURN
        IF (.not.destroy(ng, composite(cr)%ad_v, myfile,                &
     &                   __line__, 'COMPOSITE(cr)%ad_v')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%ad_Huon, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%ad_Huon')) RETURN
        IF (.not.destroy(ng, composite(cr)%ad_Hvom, myfile,             &
     &                   __line__, 'COMPOSITE(cr)%ad_Hvom')) RETURN
 
        IF (.not.destroy(ng, composite(cr)%ad_t, myfile,                &
     &                   __line__, 'COMPOSITE(cr)%ad_t')) RETURN
#   endif
#  endif
      END DO
!
!-----------------------------------------------------------------------
!  Deallocate each variable in the derived-type T_REFINED for
!  refinement grids contact region structure separately.
!-----------------------------------------------------------------------
!
      DO cr=1,ncontact
        IF (.not.destroy(ng, refined(cr)%ubar, myfile,                  &
     &                   __line__, 'REFINED(cr)%ubar')) RETURN
        IF (.not.destroy(ng, refined(cr)%vbar, myfile,                  &
     &                   __line__, 'REFINED(cr)%vbar')) RETURN
        IF (.not.destroy(ng, refined(cr)%zeta, myfile,                  &
     &                   __line__, 'REFINED(cr)%zeta')) RETURN
 
        IF (.not.destroy(ng, refined(cr)%U2d_flux, myfile,              &
     &                   __line__, 'REFINED(cr)%U2d_flux')) RETURN
        IF (.not.destroy(ng, refined(cr)%V2d_flux, myfile,               &
     &                   __line__, 'REFINED(cr)%V2d_flux')) RETURN
 
        IF (.not.destroy(ng, refined(cr)%on_u, myfile,                  &
     &                   __line__, 'REFINED(cr)%on_u')) RETURN
        IF (.not.destroy(ng, refined(cr)%om_v, myfile,                  &
     &                   __line__, 'REFINED(cr)%om_v')) RETURN
 
#  if defined TANGENT || defined TL_IOMS
        IF (.not.destroy(ng, refined(cr)%tl_ubar, myfile,               &
     &                   __line__, 'REFINED(cr)%tl_ubar')) RETURN
        IF (.not.destroy(ng, refined(cr)%tl_vbar, myfile,               &
     &                   __line__, 'REFINED(cr)%tl_vbar')) RETURN
        IF (.not.destroy(ng, refined(cr)%tl_zeta, myfile,               &
     &                   __line__, 'REFINED(cr)%tl_zeta')) RETURN
 
        IF (.not.destroy(ng, refined(cr)%tl_U2d_flux, myfile,           &
     &                   __line__, 'REFINED(cr)%tl_U2d_flux')) RETURN
        IF (.not.destroy(ng, refined(cr)%tl_V2d_flux, myfile,           &
     &                   __line__, 'REFINED(cr)%tl_V2d_flux')) RETURN
#  endif
 
#  ifdef ADJOINT
        IF (.not.destroy(ng, refined(cr)%ad_ubar, myfile,               &
     &                   __line__, 'REFINED(cr)%ad_ubar')) RETURN
        IF (.not.destroy(ng, refined(cr)%ad_vbar, myfile,               &
     &                   __line__, 'REFINED(cr)%ad_vbar')) RETURN
        IF (.not.destroy(ng, refined(cr)%ad_zeta, myfile,               &
     &                   __line__, 'REFINED(cr)%ad_zeta')) RETURN
 
        IF (.not.destroy(ng, refined(cr)%ad_U2d_flux, myfile,           &
     &                   __line__, 'REFINED(cr)%ad_U2d_flux')) RETURN
        IF (.not.destroy(ng, refined(cr)%ad_V2d_flux, myfile,           &
     &                   __line__, 'REFINED(cr)%ad_V2d_flux')) RETURN
#  endif
 
#  ifdef SOLVE3D
        IF (.not.destroy(ng, refined(cr)%u, myfile,                     &
     &                   __line__, 'REFINED(cr)%u')) RETURN
        IF (.not.destroy(ng, refined(cr)%v, myfile,                     &
     &                   __line__, 'REFINED(cr)%v')) RETURN
 
        IF (.not.destroy(ng, refined(cr)%t, myfile,                     &
     &                   __line__, 'REFINED(cr)%t')) RETURN
 
#   if defined TANGENT || defined TL_IOMS
        IF (.not.destroy(ng, refined(cr)%tl_u, myfile,                  &
     &                   __line__, 'REFINED(cr)%tl_u')) RETURN
        IF (.not.destroy(ng, refined(cr)%tl_v, myfile,                  &
     &                   __line__, 'REFINED(cr)%tl_v')) RETURN
 
        IF (.not.destroy(ng, refined(cr)%tl_t, myfile,                  &
     &                   __line__, 'REFINED(cr)%tl_t')) RETURN
#   endif
 
#   ifdef ADJOINT
        IF (.not.destroy(ng, refined(cr)%ad_u, myfile,                  &
     &                   __line__, 'REFINED(cr)%ad_u')) RETURN
        IF (.not.destroy(ng, refined(cr)%ad_v, myfile,                  &
     &                   __line__, 'REFINED(cr)%ad_v')) RETURN
 
        IF (.not.destroy(ng, refined(cr)%ad_t, myfile,                  &
     &                   __line__, 'REFINED(cr)%ad_t')) RETURN
#   endif
#  endif
      END DO
# endif
!
!-----------------------------------------------------------------------
!  Deallocate derived-type structures:
!-----------------------------------------------------------------------
!
!  Boundary connectivity.
 
      IF (allocated(bry_contact)) deallocate ( bry_contact )
!
!  Grid connectivity.
!
      IF (allocated(rcontact)) deallocate ( rcontact )
      IF (allocated(ucontact)) deallocate ( ucontact )
      IF (allocated(vcontact)) deallocate ( vcontact )
!
!  Contact region metrics.
!
      IF (allocated(contact_metric)) deallocate ( contact_metric )
!
!  Composite grid contact regions.
!
      IF (allocated(composite)) deallocate ( composite )
!
!  Deallocate refinement grids contact regions structure.
!
      IF (allocated(refined)) deallocate ( refined )
!
!-----------------------------------------------------------------------
!  Deallocate other variables in module.
!-----------------------------------------------------------------------
!
      IF (allocated(contactregion)) THEN
        deallocate ( coarserdonor )
      END IF
 
      IF (allocated(finerdonor)) THEN
        deallocate ( finerdonor )
      END IF
 
      IF (allocated(donortofiner)) THEN
        deallocate ( donortofiner )
      END IF
 
      IF (allocated(refinesteps)) THEN
        deallocate ( refinesteps )
      END IF
 
      IF (allocated(refinestepscounter)) THEN
        deallocate ( refinestepscounter )
      END IF
 
      IF (allocated(twowayinterval)) THEN
        deallocate ( twowayinterval )
      END IF
 
      IF (allocated(telescoping)) THEN
        deallocate ( telescoping )
      END IF
 
      IF (allocated(rollingindex)) THEN
        deallocate ( rollingindex )
      END IF
 
      IF (allocated(rollingtime)) THEN
        deallocate ( rollingtime )
      END IF
!
      RETURN
      END SUBROUTINE deallocate_nesting
!
      SUBROUTINE initialize_nesting
!
!=======================================================================
!                                                                      !
!  This routine initializes time varying nesting structures.           !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_scalars
!
!  Local variable declarations.
!
      integer :: cr
 
      real(r8), parameter :: inival = 0.0_r8
!
!-----------------------------------------------------------------------
!  Initialize time-varying contact regions structures.  They are used
!  to process values from contact regions to global kernel arrays and
!  vice versa.
!-----------------------------------------------------------------------
!
!  Composite grids contact region structure.
!
      IF (any(compositegrid)) THEN
        DO cr=1,ncontact
          composite(cr) % bustr = inival
          composite(cr) % bvstr = inival
 
          composite(cr) % ubar = inival
          composite(cr) % vbar = inival
          composite(cr) % zeta = inival
 
          composite(cr) % rzeta = inival
 
# ifdef SOLVE3D
          composite(cr) % DU_avg1 = inival
          composite(cr) % DV_avg1 = inival
          composite(cr) % Zt_avg1 = inival
 
          composite(cr) % u = inival
          composite(cr) % v = inival
 
          composite(cr) % Huon = inival
          composite(cr) % Hvom = inival
 
          composite(cr) % t = inival
# endif
        END DO
      END IF
!
!  Refinement grids contact region structure.
!
      IF (any(refinedgrid(:))) THEN
        DO cr=1,ncontact
          refined(cr) % ubar = inival
          refined(cr) % vbar = inival
          refined(cr) % zeta = inival
 
          refined(cr) % U2d_flux = inival
          refined(cr) % V2d_flux = inival
 
          refined(cr) % on_u = inival
          refined(cr) % om_v = inival
 
# ifdef SOLVE3D
          refined(cr) % u = inival
          refined(cr) % v = inival
 
          refined(cr) % t = inival
# endif
        END DO
      END IF
 
      RETURN
      END SUBROUTINE initialize_nesting
#endif
      END MODULE mod_nesting