MOM_grid.F90

Go to the documentation of this file.

!> Provides the ocean grid type
module mom_grid

! This file is part of MOM6. See LICENSE.md for the license.

use mom_hor_index, only : hor_index_type, hor_index_init
use mom_domains, only : mom_domain_type, get_domain_extent, compute_block_extent
use mom_error_handler, only : mom_error, mom_mesg, fatal
use mom_file_parser, only : get_param, log_param, log_version, param_file_type

implicit none ; private

#include <MOM_memory.h>

public mom_grid_init, mom_grid_end, set_derived_metrics, set_first_direction
public ispointincell, hor_index_type

!> Ocean grid type. See mom_grid for details.
type, public :: ocean_grid_type
  type(mom_domain_type), pointer :: domain => null() !< Ocean model domain
  type(mom_domain_type), pointer :: domain_aux => null() !< A non-symmetric auxiliary domain type.
  type(hor_index_type) :: hi !< Horizontal index ranges

  integer :: isc !< The start i-index of cell centers within the computational domain
  integer :: iec !< The end i-index of cell centers within the computational domain
  integer :: jsc !< The start j-index of cell centers within the computational domain
  integer :: jec !< The end j-index of cell centers within the computational domain

  integer :: isd !< The start i-index of cell centers within the data domain
  integer :: ied !< The end i-index of cell centers within the data domain
  integer :: jsd !< The start j-index of cell centers within the data domain
  integer :: jed !< The end j-index of cell centers within the data domain

  integer :: isg !< The start i-index of cell centers within the global domain
  integer :: ieg !< The end i-index of cell centers within the global domain
  integer :: jsg !< The start j-index of cell centers within the global domain
  integer :: jeg !< The end j-index of cell centers within the global domain

  integer :: iscb !< The start i-index of cell vertices within the computational domain
  integer :: iecb !< The end i-index of cell vertices within the computational domain
  integer :: jscb !< The start j-index of cell vertices within the computational domain
  integer :: jecb !< The end j-index of cell vertices within the computational domain

  integer :: isdb !< The start i-index of cell vertices within the data domain
  integer :: iedb !< The end i-index of cell vertices within the data domain
  integer :: jsdb !< The start j-index of cell vertices within the data domain
  integer :: jedb !< The end j-index of cell vertices within the data domain

  integer :: isgb !< The start i-index of cell vertices within the global domain
  integer :: iegb !< The end i-index of cell vertices within the global domain
  integer :: jsgb !< The start j-index of cell vertices within the global domain
  integer :: jegb !< The end j-index of cell vertices within the global domain

  integer :: isd_global !< The value of isd in the global index space (decompoistion invariant).
  integer :: jsd_global !< The value of isd in the global index space (decompoistion invariant).
  integer :: idg_offset !< The offset between the corresponding global and local i-indices.
  integer :: jdg_offset !< The offset between the corresponding global and local j-indices.
  integer :: ke         !< The number of layers in the vertical.
  logical :: symmetric  !< True if symmetric memory is used.
  logical :: nonblocking_updates  !< If true, non-blocking halo updates are
                                  !! allowed.  The default is .false. (for now).
  integer :: first_direction !< An integer that indicates which direction is
                             !! to be updated first in directionally split
                             !! parts of the calculation.  This can be altered
                             !! during the course of the run via calls to
                             !! set_first_direction.

  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    mask2dt, &   !< 0 for land points and 1 for ocean points on the h-grid. Nd.
    geolatt, &   !< The geographic latitude at q points in degrees of latitude or m.
    geolont, &   !< The geographic longitude at q points in degrees of longitude or m.
    dxt, &       !< dxT is delta x at h points, in m.
    idxt, &      !< 1/dxT in m-1.
    dyt, &       !< dyT is delta y at h points, in m, and IdyT is 1/dyT in m-1.
    idyt, &      !< dyT is delta y at h points, in m, and IdyT is 1/dyT in m-1.
    areat, &     !< The area of an h-cell, in m2.
    iareat, &    !< 1/areaT, in m-2.
    sin_rot, &   !< The sine of the angular rotation between the local model grid's northward
                 !! and the true northward directions.
    cos_rot      !< The cosine of the angular rotation between the local model grid's northward
                 !! and the true northward directions.

  real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: &
    mask2dcu, &  !< 0 for boundary points and 1 for ocean points on the u grid.  Nondim.
    geolatcu, &  !< The geographic latitude at u points in degrees of latitude or m.
    geoloncu, &  !< The geographic longitude at u points in degrees of longitude or m.
    dxcu, &      !< dxCu is delta x at u points, in m.
    idxcu, &     !< 1/dxCu in m-1.
    dycu, &      !< dyCu is delta y at u points, in m.
    idycu, &     !< 1/dyCu in m-1.
    dy_cu, &     !< The unblocked lengths of the u-faces of the h-cell in m.
    dy_cu_obc, & !< The unblocked lengths of the u-faces of the h-cell in m for OBC.
    iareacu, &   !< The masked inverse areas of u-grid cells in m2.
    areacu       !< The areas of the u-grid cells in m2.
  !> \todo dy_Cu_obc is not used?

  real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: &
    mask2dcv, &  !< 0 for boundary points and 1 for ocean points on the v grid.  Nondim.
    geolatcv, &  !< The geographic latitude at v points in degrees of latitude or m.
    geoloncv, &  !< The geographic longitude at v points in degrees of longitude or m.
    dxcv, &      !< dxCv is delta x at v points, in m.
    idxcv, &     !< 1/dxCv in m-1.
    dycv, &      !< dyCv is delta y at v points, in m.
    idycv, &     !< 1/dyCv in m-1.
    dx_cv, &     !< The unblocked lengths of the v-faces of the h-cell in m.
    dx_cv_obc, & !< The unblocked lengths of the v-faces of the h-cell in m for OBC.
    iareacv, &   !< The masked inverse areas of v-grid cells in m2.
    areacv       !< The areas of the v-grid cells in m2.

  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
    mask2dbu, &  !< 0 for boundary points and 1 for ocean points on the q grid.  Nondim.
    geolatbu, &  !< The geographic latitude at q points in degrees of latitude or m.
    geolonbu, &  !< The geographic longitude at q points in degrees of longitude or m.
    dxbu, &      !< dxBu is delta x at q points, in m.
    idxbu, &     !< 1/dxBu in m-1.
    dybu, &      !< dyBu is delta y at q points, in m.
    idybu, &     !< 1/dyBu in m-1.
    areabu, &    !< areaBu is the area of a q-cell, in m2
    iareabu      !< IareaBu = 1/areaBu in m-2.

  real, pointer, dimension(:) :: &
    gridlatt => null(), & !< The latitude of T points for the purpose of labeling the output axes.
                          !! On many grids this is the same as geoLatT.
    gridlatb => null()    !< The latitude of B points for the purpose of labeling the output axes.
                          !! On many grids this is the same as geoLatBu.
  real, pointer, dimension(:) :: &
    gridlont => null(), & !< The longitude of T points for the purpose of labeling the output axes.
                          !! On many grids this is the same as geoLonT.
    gridlonb => null()    !< The longitude of B points for the purpose of labeling the output axes.
                          !! On many grids this is the same as geoLonBu.
  character(len=40) :: &
    x_axis_units, &     !< The units that are used in labeling the x coordinate axes.
    y_axis_units        !< The units that are used in labeling the y coordinate axes.

  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    bathyt        !< Ocean bottom depth at tracer points, in m.

  logical :: bathymetry_at_vel  !< If true, there are separate values for the
                  !! basin depths at velocity points.  Otherwise the effects of
                  !! of topography are entirely determined from thickness points.
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: &
    dblock_u, &   !< Topographic depths at u-points at which the flow is blocked, in m.
    dopen_u       !< Topographic depths at u-points at which the flow is open at width dy_Cu, in m.
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: &
    dblock_v, &   !< Topographic depths at v-points at which the flow is blocked, in m.
    dopen_v       !< Topographic depths at v-points at which the flow is open at width dx_Cv, in m.
  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
    coriolisbu    !< The Coriolis parameter at corner points, in s-1.
  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    df_dx, &      !< Derivative d/dx f (Coriolis parameter) at h-points, in s-1 m-1.
    df_dy         !< Derivative d/dy f (Coriolis parameter) at h-points, in s-1 m-1.
  real :: g_earth !< The gravitational acceleration in m s-2.

  ! These variables are global sums that are useful for 1-d diagnostics
  real :: areat_global  !< Global sum of h-cell area in m2
  real :: iareat_global !< Global sum of inverse h-cell area (1/areaT_global) in m2.

  ! These variables are for block structures.
  integer :: nblocks
  type(hor_index_type), pointer :: block(:) => null() ! store indices for each block

  ! These parameters are run-time parameters that are used during some
  ! initialization routines (but not all)
  real :: south_lat     !< The latitude (or y-coordinate) of the first v-line
  real :: west_lon      !< The longitude (or x-coordinate) of the first u-line
  real :: len_lat = 0.  !< The latitudinal (or y-coord) extent of physical domain
  real :: len_lon = 0.  !< The longitudinal (or x-coord) extent of physical domain
  real :: rad_earth = 6.378e6 !< The radius of the planet in meters.
  real :: max_depth     !< The maximum depth of the ocean in meters.
end type ocean_grid_type

contains

!> MOM_grid_init initializes the ocean grid array sizes and grid memory.
subroutine mom_grid_init(G, param_file, HI, global_indexing, bathymetry_at_vel)
  type(ocean_grid_type), intent(inout) :: G          !< The horizontal grid type
  type(param_file_type), intent(in)    :: param_file !< Parameter file handle
  type(hor_index_type), &
                  optional, intent(in) :: HI !< A hor_index_type for array extents
  logical,        optional, intent(in) :: global_indexing !< If true use global index
                             !! values instead of having the data domain on each
                             !! processor start at 1.
  logical,        optional, intent(in) :: bathymetry_at_vel !< If true, there are
                             !! separate values for the ocean bottom depths at
                             !! velocity points.  Otherwise the effects of topography
                             !! are entirely determined from thickness points.

! This include declares and sets the variable "version".
#include "version_variable.h"
  integer :: isd, ied, jsd, jed, nk
  integer :: IsdB, IedB, JsdB, JedB
  integer :: ied_max, jed_max
  integer :: niblock, njblock, nihalo, njhalo, nblocks, n, i, j
  logical :: local_indexing  ! If false use global index values instead of having
                             ! the data domain on each processor start at 1.

  integer, allocatable, dimension(:) :: ibegin, iend, jbegin, jend
  character(len=40)  :: mod_nm  = "MOM_grid" ! This module's name.


  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mod_nm, version, &
                   "Parameters providing information about the lateral grid.")


  call get_param(param_file, mod_nm, "NIBLOCK", niblock, "The number of blocks "// &
                 "in the x-direction on each processor (for openmp).", default=1, &
                 layoutparam=.true.)
  call get_param(param_file, mod_nm, "NJBLOCK", njblock, "The number of blocks "// &
                 "in the y-direction on each processor (for openmp).", default=1, &
                 layoutparam=.true.)

  if (present(hi)) then
    g%HI = hi

    g%isc = hi%isc ; g%iec = hi%iec ; g%jsc = hi%jsc ; g%jec = hi%jec
    g%isd = hi%isd ; g%ied = hi%ied ; g%jsd = hi%jsd ; g%jed = hi%jed
    g%isg = hi%isg ; g%ieg = hi%ieg ; g%jsg = hi%jsg ; g%jeg = hi%jeg

    g%IscB = hi%IscB ; g%IecB = hi%IecB ; g%JscB = hi%JscB ; g%JecB = hi%JecB
    g%IsdB = hi%IsdB ; g%IedB = hi%IedB ; g%JsdB = hi%JsdB ; g%JedB = hi%JedB
    g%IsgB = hi%IsgB ; g%IegB = hi%IegB ; g%JsgB = hi%JsgB ; g%JegB = hi%JegB

    g%idg_offset = hi%idg_offset ; g%jdg_offset = hi%jdg_offset
    g%isd_global = g%isd + hi%idg_offset ; g%jsd_global = g%jsd + hi%jdg_offset
    g%symmetric = hi%symmetric
  else
    local_indexing = .true.
    if (present(global_indexing)) local_indexing = .not.global_indexing
    call hor_index_init(g%Domain, g%HI, param_file, &
                        local_indexing=local_indexing)

    ! get_domain_extent ensures that domains start at 1 for compatibility between
    ! static and dynamically allocated arrays, unless global_indexing is true.
    call get_domain_extent(g%Domain, g%isc, g%iec, g%jsc, g%jec, &
                           g%isd, g%ied, g%jsd, g%jed, &
                           g%isg, g%ieg, g%jsg, g%jeg, &
                           g%idg_offset, g%jdg_offset, g%symmetric, &
                           local_indexing=local_indexing)
    g%isd_global = g%isd+g%idg_offset ; g%jsd_global = g%jsd+g%jdg_offset
  endif

  g%nonblocking_updates = g%Domain%nonblocking_updates

  ! Set array sizes for fields that are discretized at tracer cell boundaries.
  g%IscB = g%isc ; g%JscB = g%jsc
  g%IsdB = g%isd ; g%JsdB = g%jsd
  g%IsgB = g%isg ; g%JsgB = g%jsg
  if (g%symmetric) then
    g%IscB = g%isc-1 ; g%JscB = g%jsc-1
    g%IsdB = g%isd-1 ; g%JsdB = g%jsd-1
    g%IsgB = g%isg-1 ; g%JsgB = g%jsg-1
  endif
  g%IecB = g%iec ; g%JecB = g%jec
  g%IedB = g%ied ; g%JedB = g%jed
  g%IegB = g%ieg ; g%JegB = g%jeg

  call mom_mesg("  MOM_grid.F90, MOM_grid_init: allocating metrics", 5)

  call allocate_metrics(g)

  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  g%bathymetry_at_vel = .false.
  if (present(bathymetry_at_vel)) g%bathymetry_at_vel = bathymetry_at_vel
  if (g%bathymetry_at_vel) then
    alloc_(g%Dblock_u(isdb:iedb, jsd:jed)) ; g%Dblock_u(:,:) = 0.0
    alloc_(g%Dopen_u(isdb:iedb, jsd:jed))  ; g%Dopen_u(:,:) = 0.0
    alloc_(g%Dblock_v(isd:ied, jsdb:jedb)) ; g%Dblock_v(:,:) = 0.0
    alloc_(g%Dopen_v(isd:ied, jsdb:jedb))  ; g%Dopen_v(:,:) = 0.0
  endif

! setup block indices.
  nihalo = g%Domain%nihalo
  njhalo = g%Domain%njhalo
  nblocks = niblock * njblock
  if (nblocks < 1) call mom_error(fatal, "MOM_grid_init: " // &
       "nblocks(=NI_BLOCK*NJ_BLOCK) must be no less than 1")

  allocate(ibegin(niblock), iend(niblock), jbegin(njblock), jend(njblock))
  call compute_block_extent(g%HI%isc,g%HI%iec,niblock,ibegin,iend)
  call compute_block_extent(g%HI%jsc,g%HI%jec,njblock,jbegin,jend)
  !-- make sure the last block is the largest.
  do i = 1, niblock-1
    if (iend(i)-ibegin(i) > iend(niblock)-ibegin(niblock) ) call mom_error(fatal, &
       "MOM_grid_init: the last block size in x-direction is not the largest")
  enddo
  do j = 1, njblock-1
    if (jend(j)-jbegin(j) > jend(njblock)-jbegin(njblock) ) call mom_error(fatal, &
       "MOM_grid_init: the last block size in y-direction is not the largest")
  enddo

  g%nblocks = nblocks
  allocate(g%Block(nblocks))
  ied_max = 1 ; jed_max = 1
  do n = 1,nblocks
    ! Copy all information from the array index type describing the local grid.
    g%Block(n) = g%HI

    i = mod((n-1), niblock) + 1
    j = (n-1)/niblock + 1
    !--- isd and jsd are always 1 for each block to permit array reuse.
    g%Block(n)%isd = 1 ; g%Block(n)%jsd = 1
    g%Block(n)%isc = g%Block(n)%isd+nihalo
    g%Block(n)%jsc = g%Block(n)%jsd+njhalo
    g%Block(n)%iec = g%Block(n)%isc + iend(i) - ibegin(i)
    g%Block(n)%jec = g%Block(n)%jsc + jend(j) - jbegin(j)
    g%Block(n)%ied = g%Block(n)%iec + nihalo
    g%Block(n)%jed = g%Block(n)%jec + njhalo
    g%Block(n)%IscB = g%Block(n)%isc; g%Block(n)%IecB = g%Block(n)%iec
    g%Block(n)%JscB = g%Block(n)%jsc; g%Block(n)%JecB = g%Block(n)%jec
    !   For symmetric memory domains, the first block will have the extra point
    ! at the lower boundary of its computational domain.
    if (g%symmetric) then
      if (i==1) g%Block(n)%IscB = g%Block(n)%IscB-1
      if (j==1) g%Block(n)%JscB = g%Block(n)%JscB-1
    endif
    g%Block(n)%IsdB = g%Block(n)%isd; g%Block(n)%IedB = g%Block(n)%ied
    g%Block(n)%JsdB = g%Block(n)%jsd; g%Block(n)%JedB = g%Block(n)%jed
    !--- For symmetric memory domain, every block will have an extra point
    !--- at the lower boundary of its data domain.
    if (g%symmetric) then
      g%Block(n)%IsdB = g%Block(n)%IsdB-1
      g%Block(n)%JsdB = g%Block(n)%JsdB-1
    endif
    g%Block(n)%idg_offset = (ibegin(i) - g%Block(n)%isc) + g%HI%idg_offset
    g%Block(n)%jdg_offset = (jbegin(j) - g%Block(n)%jsc) + g%HI%jdg_offset
    ! Find the largest values of ied and jed so that all blocks will have the
    ! same size in memory.
    ied_max = max(ied_max, g%Block(n)%ied)
    jed_max = max(jed_max, g%Block(n)%jed)
  enddo

  ! Reset all of the data domain sizes to match the largest for array reuse,
  ! recalling that all block have isd=jed=1 for array reuse.
  do n = 1,nblocks
    g%Block(n)%ied = ied_max ; g%Block(n)%IedB = ied_max
    g%Block(n)%jed = jed_max ; g%Block(n)%JedB = jed_max
  enddo

  !-- do some bounds error checking
  if ( g%block(nblocks)%ied+g%block(nblocks)%idg_offset > g%HI%ied + g%HI%idg_offset ) &
        call mom_error(fatal, "MOM_grid_init: G%ied_bk > G%ied")
  if ( g%block(nblocks)%jed+g%block(nblocks)%jdg_offset > g%HI%jed + g%HI%jdg_offset ) &
        call mom_error(fatal, "MOM_grid_init: G%jed_bk > G%jed")

end subroutine mom_grid_init


!> set_derived_metrics calculates metric terms that are derived from other metrics.
subroutine set_derived_metrics(G)
  type(ocean_grid_type), intent(inout) :: G    !< The horizontal grid structure
!    Various inverse grid spacings and derived areas are calculated within this
!  subroutine.
  integer :: i, j, isd, ied, jsd, jed
  integer :: IsdB, IedB, JsdB, JedB

  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  do j=jsd,jed ; do i=isd,ied
    if (g%dxT(i,j) < 0.0) g%dxT(i,j) = 0.0
    if (g%dyT(i,j) < 0.0) g%dyT(i,j) = 0.0
    g%IdxT(i,j) = adcroft_reciprocal(g%dxT(i,j))
    g%IdyT(i,j) = adcroft_reciprocal(g%dyT(i,j))
    g%IareaT(i,j) = adcroft_reciprocal(g%areaT(i,j))
  enddo ; enddo

  do j=jsd,jed ; do i=isdb,iedb
    if (g%dxCu(i,j) < 0.0) g%dxCu(i,j) = 0.0
    if (g%dyCu(i,j) < 0.0) g%dyCu(i,j) = 0.0
    g%IdxCu(i,j) = adcroft_reciprocal(g%dxCu(i,j))
    g%IdyCu(i,j) = adcroft_reciprocal(g%dyCu(i,j))
  enddo ; enddo

  do j=jsdb,jedb ; do i=isd,ied
    if (g%dxCv(i,j) < 0.0) g%dxCv(i,j) = 0.0
    if (g%dyCv(i,j) < 0.0) g%dyCv(i,j) = 0.0
    g%IdxCv(i,j) = adcroft_reciprocal(g%dxCv(i,j))
    g%IdyCv(i,j) = adcroft_reciprocal(g%dyCv(i,j))
  enddo ; enddo

  do j=jsdb,jedb ; do i=isdb,iedb
    if (g%dxBu(i,j) < 0.0) g%dxBu(i,j) = 0.0
    if (g%dyBu(i,j) < 0.0) g%dyBu(i,j) = 0.0

    g%IdxBu(i,j) = adcroft_reciprocal(g%dxBu(i,j))
    g%IdyBu(i,j) = adcroft_reciprocal(g%dyBu(i,j))
    ! areaBu has usually been set to a positive area elsewhere.
    if (g%areaBu(i,j) <= 0.0) g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
    g%IareaBu(i,j) =  adcroft_reciprocal(g%areaBu(i,j))
  enddo ; enddo
end subroutine set_derived_metrics

!> Adcroft_reciprocal(x) = 1/x for |x|>0 or 0 for x=0.
function adcroft_reciprocal(val) result(I_val)
  real, intent(in) :: val  !< The value being inverted.
  real :: I_val            !< The Adcroft reciprocal of val.

  i_val = 0.0 ; if (val /= 0.0) i_val = 1.0/val
end function adcroft_reciprocal

!> Returns true if the coordinates (x,y) are within the h-cell (i,j)
logical function ispointincell(G, i, j, x, y)
  type(ocean_grid_type), intent(in) :: G !< Grid type
  integer,               intent(in) :: i !< i index of cell to test
  integer,               intent(in) :: j !< j index of cell to test
  real,                  intent(in) :: x !< x coordinate of point
  real,                  intent(in) :: y !< y coordinate of point
  ! Local variables
  real :: xNE, xNW, xSE, xSW, yNE, yNW, ySE, ySW
  real :: p0, p1, p2, p3, l0, l1, l2, l3
  ispointincell = .false.
  xne = g%geoLonBu(i  ,j  ) ; yne = g%geoLatBu(i  ,j  )
  xnw = g%geoLonBu(i-1,j  ) ; ynw = g%geoLatBu(i-1,j  )
  xse = g%geoLonBu(i  ,j-1) ; yse = g%geoLatBu(i  ,j-1)
  xsw = g%geoLonBu(i-1,j-1) ; ysw = g%geoLatBu(i-1,j-1)
  ! This is a crude calculation that assume a geographic coordinate system
  if (x<min(xne,xnw,xse,xsw) .or. x>max(xne,xnw,xse,xsw) .or. &
      y<min(yne,ynw,yse,ysw) .or. y>max(yne,ynw,yse,ysw) ) then
    return ! Avoid the more complicated calculation
  endif
  l0 = (x-xsw)*(yse-ysw) - (y-ysw)*(xse-xsw)
  l1 = (x-xse)*(yne-yse) - (y-yse)*(xne-xse)
  l2 = (x-xne)*(ynw-yne) - (y-yne)*(xnw-xne)
  l3 = (x-xnw)*(ysw-ynw) - (y-ynw)*(xsw-xnw)

  p0 = sign(1., l0) ; if (l0 == 0.) p0=0.
  p1 = sign(1., l1) ; if (l1 == 0.) p1=0.
  p2 = sign(1., l2) ; if (l2 == 0.) p2=0.
  p3 = sign(1., l3) ; if (l3 == 0.) p3=0.

  if ( (abs(p0)+abs(p2)) + (abs(p1)+abs(p3)) == abs((p0+p2) + (p1+p3)) ) then
    ispointincell=.true.
  endif
end function ispointincell

subroutine set_first_direction(G, y_first)
  type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure
  integer,               intent(in) :: y_first

  g%first_direction = y_first
end subroutine set_first_direction

!> Allocate memory used by the ocean_grid_type and related structures.
subroutine allocate_metrics(G)
  type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isg, ieg, jsg, jeg

  ! This subroutine allocates the lateral elements of the ocean_grid_type that
  ! are always used and zeros them out.

  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg

  alloc_(g%dxT(isd:ied,jsd:jed))       ; g%dxT(:,:) = 0.0
  alloc_(g%dxCu(isdb:iedb,jsd:jed))    ; g%dxCu(:,:) = 0.0
  alloc_(g%dxCv(isd:ied,jsdb:jedb))    ; g%dxCv(:,:) = 0.0
  alloc_(g%dxBu(isdb:iedb,jsdb:jedb))  ; g%dxBu(:,:) = 0.0
  alloc_(g%IdxT(isd:ied,jsd:jed))      ; g%IdxT(:,:) = 0.0
  alloc_(g%IdxCu(isdb:iedb,jsd:jed))   ; g%IdxCu(:,:) = 0.0
  alloc_(g%IdxCv(isd:ied,jsdb:jedb))   ; g%IdxCv(:,:) = 0.0
  alloc_(g%IdxBu(isdb:iedb,jsdb:jedb)) ; g%IdxBu(:,:) = 0.0

  alloc_(g%dyT(isd:ied,jsd:jed))       ; g%dyT(:,:) = 0.0
  alloc_(g%dyCu(isdb:iedb,jsd:jed))    ; g%dyCu(:,:) = 0.0
  alloc_(g%dyCv(isd:ied,jsdb:jedb))    ; g%dyCv(:,:) = 0.0
  alloc_(g%dyBu(isdb:iedb,jsdb:jedb))  ; g%dyBu(:,:) = 0.0
  alloc_(g%IdyT(isd:ied,jsd:jed))      ; g%IdyT(:,:) = 0.0
  alloc_(g%IdyCu(isdb:iedb,jsd:jed))   ; g%IdyCu(:,:) = 0.0
  alloc_(g%IdyCv(isd:ied,jsdb:jedb))   ; g%IdyCv(:,:) = 0.0
  alloc_(g%IdyBu(isdb:iedb,jsdb:jedb)) ; g%IdyBu(:,:) = 0.0

  alloc_(g%areaT(isd:ied,jsd:jed))       ; g%areaT(:,:) = 0.0
  alloc_(g%IareaT(isd:ied,jsd:jed))      ; g%IareaT(:,:) = 0.0
  alloc_(g%areaBu(isdb:iedb,jsdb:jedb))  ; g%areaBu(:,:) = 0.0
  alloc_(g%IareaBu(isdb:iedb,jsdb:jedb)) ; g%IareaBu(:,:) = 0.0

  alloc_(g%mask2dT(isd:ied,jsd:jed))      ; g%mask2dT(:,:) = 0.0
  alloc_(g%mask2dCu(isdb:iedb,jsd:jed))   ; g%mask2dCu(:,:) = 0.0
  alloc_(g%mask2dCv(isd:ied,jsdb:jedb))   ; g%mask2dCv(:,:) = 0.0
  alloc_(g%mask2dBu(isdb:iedb,jsdb:jedb)) ; g%mask2dBu(:,:) = 0.0
  alloc_(g%geoLatT(isd:ied,jsd:jed))      ; g%geoLatT(:,:) = 0.0
  alloc_(g%geoLatCu(isdb:iedb,jsd:jed))   ; g%geoLatCu(:,:) = 0.0
  alloc_(g%geoLatCv(isd:ied,jsdb:jedb))   ; g%geoLatCv(:,:) = 0.0
  alloc_(g%geoLatBu(isdb:iedb,jsdb:jedb)) ; g%geoLatBu(:,:) = 0.0
  alloc_(g%geoLonT(isd:ied,jsd:jed))      ; g%geoLonT(:,:) = 0.0
  alloc_(g%geoLonCu(isdb:iedb,jsd:jed))   ; g%geoLonCu(:,:) = 0.0
  alloc_(g%geoLonCv(isd:ied,jsdb:jedb))   ; g%geoLonCv(:,:) = 0.0
  alloc_(g%geoLonBu(isdb:iedb,jsdb:jedb)) ; g%geoLonBu(:,:) = 0.0

  alloc_(g%dx_Cv(isd:ied,jsdb:jedb))     ; g%dx_Cv(:,:) = 0.0
  alloc_(g%dy_Cu(isdb:iedb,jsd:jed))     ; g%dy_Cu(:,:) = 0.0
  alloc_(g%dx_Cv_obc(isd:ied,jsdb:jedb)) ; g%dx_Cv_obc(:,:) = 0.0
  alloc_(g%dy_Cu_obc(isdb:iedb,jsd:jed)) ; g%dy_Cu_obc(:,:) = 0.0

  alloc_(g%areaCu(isdb:iedb,jsd:jed))  ; g%areaCu(:,:) = 0.0
  alloc_(g%areaCv(isd:ied,jsdb:jedb))  ; g%areaCv(:,:) = 0.0
  alloc_(g%IareaCu(isdb:iedb,jsd:jed)) ; g%IareaCu(:,:) = 0.0
  alloc_(g%IareaCv(isd:ied,jsdb:jedb)) ; g%IareaCv(:,:) = 0.0

  alloc_(g%bathyT(isd:ied, jsd:jed)) ; g%bathyT(:,:) = 0.0
  alloc_(g%CoriolisBu(isdb:iedb, jsdb:jedb)) ; g%CoriolisBu(:,:) = 0.0
  alloc_(g%dF_dx(isd:ied, jsd:jed)) ; g%dF_dx(:,:) = 0.0
  alloc_(g%dF_dy(isd:ied, jsd:jed)) ; g%dF_dy(:,:) = 0.0

  alloc_(g%sin_rot(isd:ied,jsd:jed)) ; g%sin_rot(:,:) = 0.0
  alloc_(g%cos_rot(isd:ied,jsd:jed)) ; g%cos_rot(:,:) = 1.0

  allocate(g%gridLonT(isg:ieg))   ; g%gridLonT(:) = 0.0
  allocate(g%gridLonB(g%IsgB:g%IegB)) ; g%gridLonB(:) = 0.0
  allocate(g%gridLatT(jsg:jeg))   ; g%gridLatT(:) = 0.0
  allocate(g%gridLatB(g%JsgB:g%JegB)) ; g%gridLatB(:) = 0.0

end subroutine allocate_metrics

!> Release memory used by the ocean_grid_type and related structures.
subroutine mom_grid_end(G)
  type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type

  dealloc_(g%dxT)  ; dealloc_(g%dxCu)  ; dealloc_(g%dxCv)  ; dealloc_(g%dxBu)
  dealloc_(g%IdxT) ; dealloc_(g%IdxCu) ; dealloc_(g%IdxCv) ; dealloc_(g%IdxBu)

  dealloc_(g%dyT)  ; dealloc_(g%dyCu)  ; dealloc_(g%dyCv)  ; dealloc_(g%dyBu)
  dealloc_(g%IdyT) ; dealloc_(g%IdyCu) ; dealloc_(g%IdyCv) ; dealloc_(g%IdyBu)

  dealloc_(g%areaT)  ; dealloc_(g%IareaT)
  dealloc_(g%areaBu) ; dealloc_(g%IareaBu)
  dealloc_(g%areaCu) ; dealloc_(g%IareaCu)
  dealloc_(g%areaCv)  ; dealloc_(g%IareaCv)

  dealloc_(g%mask2dT)  ; dealloc_(g%mask2dCu)
  dealloc_(g%mask2dCv) ; dealloc_(g%mask2dBu)

  dealloc_(g%geoLatT)  ; dealloc_(g%geoLatCu)
  dealloc_(g%geoLatCv) ; dealloc_(g%geoLatBu)
  dealloc_(g%geoLonT)  ; dealloc_(g%geoLonCu)
  dealloc_(g%geoLonCv) ; dealloc_(g%geoLonBu)

  dealloc_(g%dx_Cv) ; dealloc_(g%dy_Cu)
  dealloc_(g%dx_Cv_obc) ; dealloc_(g%dy_Cu_obc)

  dealloc_(g%bathyT)  ; dealloc_(g%CoriolisBu)
  dealloc_(g%dF_dx)  ; dealloc_(g%dF_dy)
  dealloc_(g%sin_rot) ; dealloc_(g%cos_rot)

  if (g%bathymetry_at_vel) then
    dealloc_(g%Dblock_u) ; dealloc_(g%Dopen_u)
    dealloc_(g%Dblock_v) ; dealloc_(g%Dopen_v)
  endif

  deallocate(g%gridLonT) ; deallocate(g%gridLatT)
  deallocate(g%gridLonB) ; deallocate(g%gridLatB)

  deallocate(g%Domain%mpp_domain)
  deallocate(g%Domain)

end subroutine mom_grid_end

!> \namespace mom_grid
!!
!! Grid metrics and their inverses are labelled according to their staggered location on a Arakawa C (or B) grid.
!! - Metrics centered on h- or T-points are labelled T, e.g. dxT is the distance across the cell in the x-direction.
!! - Metrics centered on u-points are labelled Cu (C-grid u location). e.g. dyCu is the y-distance between two corners of a T-cell.
!! - Metrics centered on v-points are labelled Cv (C-grid v location). e.g. dyCv is the y-distance between two -points.
!! - Metrics centered on q-points are labelled Bu (B-grid u,v location). e.g. areaBu is the area centered on a q-point.
!!
!! \image html Grid_metrics.png "The labelling of distances (grid metrics) at various staggered location on an T-cell and around a q-point.
!!
!! Areas centered at T-, u-, v- and q- points are `areaT`, `areaCu`, `areaCv` and `areaBu` respectively.
!!
!! The reciprocal of metrics are pre-calculated and also stored in the ocean_grid_type with a I prepended to the name.
!! For example, `1./areaT` is called `IareaT`, and `1./dyCv` is `IdyCv`.
!!
!! Geographic latitude and longitude (or model coordinates if not on a sphere) are stored in `geoLatT`, `geoLonT` for T-points.
!! u-, v- and q- point coordinates are follow same pattern of replacing T with Cu, Cv and Bu respectively.
!!
!! Each location also has a 2D mask indicating whether the entire column is land or ocean.
!! `mask2dT` is 1 if the column is wet or 0 if the T-cell is land.
!! `mask2dCu` is 1 if both neighboring column are ocean, and 0 if either is land.

end module mom_grid