MOM_dyn_horgrid.F90

!> Contains a shareable dynamic type for describing horizontal grids and metric data
!! and utilty routines that work on this type.
module mom_dyn_horgrid
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_hor_index, only : hor_index_type
use mom_domains, only : mom_domain_type
use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
 
implicit none ; private
 
public create_dyn_horgrid, destroy_dyn_horgrid, set_derived_dyn_horgrid
public rescale_dyn_horgrid_bathymetry
 
! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.
 
!> Describes the horizontal ocean grid with only dynamic memory arrays
type, public :: dyn_horgrid_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.
  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, allocatable, dimension(:,:) :: &
    mask2dt, &   !< 0 for land points and 1 for ocean points on the h-grid [nondim].
    geolatt, &   !< The geographic latitude at q points [degrees of latitude] or [m].
    geolont, &   !< The geographic longitude at q points [degrees of longitude] or [m].
    dxt, &       !< dxT is delta x at h points [m].
    idxt, &      !< 1/dxT [m-1].
    dyt, &       !< dyT is delta y at h points [m].
    idyt, &      !< IdyT is 1/dyT [m-1].
    areat, &     !< The area of an h-cell [m2].
    iareat       !< 1/areaT [m-2].
  real, allocatable, dimension(:,:) :: sin_rot
                 !< The sine of the angular rotation between the local model grid's northward
                 !! and the true northward directions [nondim].
  real, allocatable, dimension(:,:) :: cos_rot
                 !< The cosine of the angular rotation between the local model grid's northward
                 !! and the true northward directions [nondim].
 
  real, allocatable, dimension(:,:) :: &
    mask2dcu, &  !< 0 for boundary points and 1 for ocean points on the u grid [nondim].
    geolatcu, &  !< The geographic latitude at u points [degrees of latitude] or [m].
    geoloncu, &  !< The geographic longitude at u points [degrees of longitude] or [m].
    dxcu, &      !< dxCu is delta x at u points [m].
    idxcu, &     !< 1/dxCu [m-1].
    dycu, &      !< dyCu is delta y at u points [m].
    idycu, &     !< 1/dyCu [m-1].
    dy_cu, &     !< The unblocked lengths of the u-faces of the h-cell [m].
    iareacu, &   !< The masked inverse areas of u-grid cells [m2].
    areacu       !< The areas of the u-grid cells [m2].
 
  real, allocatable, dimension(:,:) :: &
    mask2dcv, &  !< 0 for boundary points and 1 for ocean points on the v grid [nondim].
    geolatcv, &  !< The geographic latitude at v points [degrees of latitude] or [m].
    geoloncv, &  !< The geographic longitude at v points [degrees of longitude] or [m].
    dxcv, &      !< dxCv is delta x at v points [m].
    idxcv, &     !< 1/dxCv [m-1].
    dycv, &      !< dyCv is delta y at v points [m].
    idycv, &     !< 1/dyCv [m-1].
    dx_cv, &     !< The unblocked lengths of the v-faces of the h-cell [m].
    iareacv, &   !< The masked inverse areas of v-grid cells [m2].
    areacv       !< The areas of the v-grid cells [m2].
 
  real, allocatable, dimension(:,:) :: &
    mask2dbu, &  !< 0 for boundary points and 1 for ocean points on the q grid [nondim].
    geolatbu, &  !< The geographic latitude at q points [degrees of latitude] or [m].
    geolonbu, &  !< The geographic longitude at q points [degrees of longitude] or [m].
    dxbu, &      !< dxBu is delta x at q points [m].
    idxbu, &     !< 1/dxBu [m-1].
    dybu, &      !< dyBu is delta y at q points [m].
    idybu, &     !< 1/dyBu [m-1].
    areabu, &    !< areaBu is the area of a q-cell [m2]
    iareabu      !< IareaBu = 1/areaBu [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.
  real, pointer, dimension(:) :: 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.
  real, pointer, dimension(:) :: 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.
    ! Except on a Cartesian grid, these are usually  some variant of "degrees".
 
  real, allocatable, dimension(:,:) :: &
    bathyt        !< Ocean bottom depth at tracer points, in depth units [Z ~> 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, allocatable, dimension(:,:) :: &
    dblock_u, &   !< Topographic depths at u-points at which the flow is blocked [Z ~> m].
    dopen_u       !< Topographic depths at u-points at which the flow is open at width dy_Cu [Z ~> m].
  real, allocatable, dimension(:,:) :: &
    dblock_v, &   !< Topographic depths at v-points at which the flow is blocked [Z ~> m].
    dopen_v       !< Topographic depths at v-points at which the flow is open at width dx_Cv [Z ~> m].
  real, allocatable, dimension(:,:) :: &
    coriolisbu    !< The Coriolis parameter at corner points [T-1 ~> s-1].
  real, allocatable, dimension(:,:) :: &
    df_dx, &      !< Derivative d/dx f (Coriolis parameter) at h-points [T-1 m-1 ~> s-1 m-1].
    df_dy         !< Derivative d/dy f (Coriolis parameter) at h-points [T-1 m-1 ~> s-1 m-1].
 
  ! These variables are global sums that are useful for 1-d diagnostics
  real :: areat_global  !< Global sum of h-cell area [m2]
  real :: iareat_global !< Global sum of inverse h-cell area (1/areaT_global) [m-2]
 
  ! 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 [m].
  real :: max_depth     !< The maximum depth of the ocean [Z ~> m].
end type dyn_horgrid_type
 
contains
 
!---------------------------------------------------------------------
!> Allocate memory used by the dyn_horgrid_type and related structures.
subroutine create_dyn_horgrid(G, HI, bathymetry_at_vel)
  type(dyn_horgrid_type), pointer    :: g  !< A pointer to the dynamic horizontal grid type
  type(hor_index_type),   intent(in) :: hi !< A hor_index_type for array extents
  logical,        optional, intent(in) :: bathymetry_at_vel !< If true, there are
                             !! separate values for the basin depths at velocity
                             !! points.  Otherwise the effects of topography are
                             !! entirely determined from thickness points.
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, isg, ieg, jsg, jeg
 
  ! This subroutine allocates the lateral elements of the dyn_horgrid_type that
  ! are always used and zeros them out.
 
  if (associated(g)) then
    call mom_error(warning, "create_dyn_horgrid called with an associated horgrid_type.")
  else
    allocate(g)
  endif
 
  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
 
  g%bathymetry_at_vel = .false.
  if (present(bathymetry_at_vel)) g%bathymetry_at_vel = bathymetry_at_vel
 
  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
 
  allocate(g%dxT(isd:ied,jsd:jed))       ; g%dxT(:,:) = 0.0
  allocate(g%dxCu(isdb:iedb,jsd:jed))    ; g%dxCu(:,:) = 0.0
  allocate(g%dxCv(isd:ied,jsdb:jedb))    ; g%dxCv(:,:) = 0.0
  allocate(g%dxBu(isdb:iedb,jsdb:jedb))  ; g%dxBu(:,:) = 0.0
  allocate(g%IdxT(isd:ied,jsd:jed))      ; g%IdxT(:,:) = 0.0
  allocate(g%IdxCu(isdb:iedb,jsd:jed))   ; g%IdxCu(:,:) = 0.0
  allocate(g%IdxCv(isd:ied,jsdb:jedb))   ; g%IdxCv(:,:) = 0.0
  allocate(g%IdxBu(isdb:iedb,jsdb:jedb)) ; g%IdxBu(:,:) = 0.0
 
  allocate(g%dyT(isd:ied,jsd:jed))       ; g%dyT(:,:) = 0.0
  allocate(g%dyCu(isdb:iedb,jsd:jed))    ; g%dyCu(:,:) = 0.0
  allocate(g%dyCv(isd:ied,jsdb:jedb))    ; g%dyCv(:,:) = 0.0
  allocate(g%dyBu(isdb:iedb,jsdb:jedb))  ; g%dyBu(:,:) = 0.0
  allocate(g%IdyT(isd:ied,jsd:jed))      ; g%IdyT(:,:) = 0.0
  allocate(g%IdyCu(isdb:iedb,jsd:jed))   ; g%IdyCu(:,:) = 0.0
  allocate(g%IdyCv(isd:ied,jsdb:jedb))   ; g%IdyCv(:,:) = 0.0
  allocate(g%IdyBu(isdb:iedb,jsdb:jedb)) ; g%IdyBu(:,:) = 0.0
 
  allocate(g%areaT(isd:ied,jsd:jed))       ; g%areaT(:,:) = 0.0
  allocate(g%IareaT(isd:ied,jsd:jed))      ; g%IareaT(:,:) = 0.0
  allocate(g%areaBu(isdb:iedb,jsdb:jedb))  ; g%areaBu(:,:) = 0.0
  allocate(g%IareaBu(isdb:iedb,jsdb:jedb)) ; g%IareaBu(:,:) = 0.0
 
  allocate(g%mask2dT(isd:ied,jsd:jed))      ; g%mask2dT(:,:) = 0.0
  allocate(g%mask2dCu(isdb:iedb,jsd:jed))   ; g%mask2dCu(:,:) = 0.0
  allocate(g%mask2dCv(isd:ied,jsdb:jedb))   ; g%mask2dCv(:,:) = 0.0
  allocate(g%mask2dBu(isdb:iedb,jsdb:jedb)) ; g%mask2dBu(:,:) = 0.0
  allocate(g%geoLatT(isd:ied,jsd:jed))      ; g%geoLatT(:,:) = 0.0
  allocate(g%geoLatCu(isdb:iedb,jsd:jed))   ; g%geoLatCu(:,:) = 0.0
  allocate(g%geoLatCv(isd:ied,jsdb:jedb))   ; g%geoLatCv(:,:) = 0.0
  allocate(g%geoLatBu(isdb:iedb,jsdb:jedb)) ; g%geoLatBu(:,:) = 0.0
  allocate(g%geoLonT(isd:ied,jsd:jed))      ; g%geoLonT(:,:) = 0.0
  allocate(g%geoLonCu(isdb:iedb,jsd:jed))   ; g%geoLonCu(:,:) = 0.0
  allocate(g%geoLonCv(isd:ied,jsdb:jedb))   ; g%geoLonCv(:,:) = 0.0
  allocate(g%geoLonBu(isdb:iedb,jsdb:jedb)) ; g%geoLonBu(:,:) = 0.0
 
  allocate(g%dx_Cv(isd:ied,jsdb:jedb))     ; g%dx_Cv(:,:) = 0.0
  allocate(g%dy_Cu(isdb:iedb,jsd:jed))     ; g%dy_Cu(:,:) = 0.0
 
  allocate(g%areaCu(isdb:iedb,jsd:jed))  ; g%areaCu(:,:) = 0.0
  allocate(g%areaCv(isd:ied,jsdb:jedb))  ; g%areaCv(:,:) = 0.0
  allocate(g%IareaCu(isdb:iedb,jsd:jed)) ; g%IareaCu(:,:) = 0.0
  allocate(g%IareaCv(isd:ied,jsdb:jedb)) ; g%IareaCv(:,:) = 0.0
 
  allocate(g%bathyT(isd:ied, jsd:jed)) ; g%bathyT(:,:) = 0.0
  allocate(g%CoriolisBu(isdb:iedb, jsdb:jedb)) ; g%CoriolisBu(:,:) = 0.0
  allocate(g%dF_dx(isd:ied, jsd:jed)) ; g%dF_dx(:,:) = 0.0
  allocate(g%dF_dy(isd:ied, jsd:jed)) ; g%dF_dy(:,:) = 0.0
 
  allocate(g%sin_rot(isd:ied,jsd:jed)) ; g%sin_rot(:,:) = 0.0
  allocate(g%cos_rot(isd:ied,jsd:jed)) ; g%cos_rot(:,:) = 1.0
 
  if (g%bathymetry_at_vel) then
    allocate(g%Dblock_u(isdb:iedb, jsd:jed)) ; g%Dblock_u(:,:) = 0.0
    allocate(g%Dopen_u(isdb:iedb, jsd:jed))  ; g%Dopen_u(:,:) = 0.0
    allocate(g%Dblock_v(isd:ied, jsdb:jedb)) ; g%Dblock_v(:,:) = 0.0
    allocate(g%Dopen_v(isd:ied, jsdb:jedb))  ; g%Dopen_v(:,:) = 0.0
  endif
 
  ! gridLonB and gridLatB are used as edge values in some cases, so they
  ! always need to use symmetric memory allcoations.
  allocate(g%gridLonT(isg:ieg))   ; g%gridLonT(:) = 0.0
  allocate(g%gridLonB(isg-1:ieg)) ; g%gridLonB(:) = 0.0
  allocate(g%gridLatT(jsg:jeg))   ; g%gridLatT(:) = 0.0
  allocate(g%gridLatB(jsg-1:jeg)) ; g%gridLatB(:) = 0.0
 
end subroutine create_dyn_horgrid
 
!> rescale_dyn_horgrid_bathymetry permits a change in the internal units for the bathymetry on the
!! grid, both rescaling the depths and recording the new internal depth units.
subroutine rescale_dyn_horgrid_bathymetry(G, m_in_new_units)
  type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid type
  real,                   intent(in)    :: m_in_new_units !< The new internal representation of 1 m depth.
 
  ! Local variables
  real :: rescale
  integer :: i, j, isd, ied, jsd, jed, 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
 
  if (m_in_new_units == 1.0) return
  if (m_in_new_units < 0.0) &
    call mom_error(fatal, "rescale_grid_bathymetry: Negative depth units are not permitted.")
  if (m_in_new_units == 0.0) &
    call mom_error(fatal, "rescale_grid_bathymetry: Zero depth units are not permitted.")
 
  rescale = 1.0 / m_in_new_units
  do j=jsd,jed ; do i=isd,ied
    g%bathyT(i,j) = rescale*g%bathyT(i,j)
  enddo ; enddo
  if (g%bathymetry_at_vel) then ; do j=jsd,jed ; do i=isdb,iedb
    g%Dblock_u(i,j) = rescale*g%Dblock_u(i,j) ; g%Dopen_u(i,j) = rescale*g%Dopen_u(i,j)
  enddo ; enddo ; endif
  if (g%bathymetry_at_vel) then ; do j=jsdb,jedb ; do i=isd,ied
    g%Dblock_v(i,j) = rescale*g%Dblock_v(i,j) ; g%Dopen_v(i,j) = rescale*g%Dopen_v(i,j)
  enddo ; enddo ; endif
  g%max_depth = rescale*g%max_depth
 
end subroutine rescale_dyn_horgrid_bathymetry
 
!> set_derived_dyn_horgrid calculates metric terms that are derived from other metrics.
subroutine set_derived_dyn_horgrid(G)
  type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid type
!    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_dyn_horgrid
 
!> 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
 
!---------------------------------------------------------------------
!> Release memory used by the dyn_horgrid_type and related structures.
subroutine destroy_dyn_horgrid(G)
  type(dyn_horgrid_type), pointer :: g !< The dynamic horizontal grid type
 
  if (.not.associated(g)) then
    call mom_error(fatal, "destroy_dyn_horgrid called with an unassociated horgrid_type.")
  endif
 
  deallocate(g%dxT)  ; deallocate(g%dxCu)  ; deallocate(g%dxCv)  ; deallocate(g%dxBu)
  deallocate(g%IdxT) ; deallocate(g%IdxCu) ; deallocate(g%IdxCv) ; deallocate(g%IdxBu)
 
  deallocate(g%dyT)  ; deallocate(g%dyCu)  ; deallocate(g%dyCv)  ; deallocate(g%dyBu)
  deallocate(g%IdyT) ; deallocate(g%IdyCu) ; deallocate(g%IdyCv) ; deallocate(g%IdyBu)
 
  deallocate(g%areaT)  ; deallocate(g%IareaT)
  deallocate(g%areaBu) ; deallocate(g%IareaBu)
  deallocate(g%areaCu) ; deallocate(g%IareaCu)
  deallocate(g%areaCv)  ; deallocate(g%IareaCv)
 
  deallocate(g%mask2dT)  ; deallocate(g%mask2dCu)
  deallocate(g%mask2dCv) ; deallocate(g%mask2dBu)
 
  deallocate(g%geoLatT)  ; deallocate(g%geoLatCu)
  deallocate(g%geoLatCv) ; deallocate(g%geoLatBu)
  deallocate(g%geoLonT)  ; deallocate(g%geoLonCu)
  deallocate(g%geoLonCv) ; deallocate(g%geoLonBu)
 
  deallocate(g%dx_Cv) ; deallocate(g%dy_Cu)
 
  deallocate(g%bathyT)  ; deallocate(g%CoriolisBu)
  deallocate(g%dF_dx)  ; deallocate(g%dF_dy)
  deallocate(g%sin_rot) ; deallocate(g%cos_rot)
 
  if (allocated(g%Dblock_u)) deallocate(g%Dblock_u)
  if (allocated(g%Dopen_u)) deallocate(g%Dopen_u)
  if (allocated(g%Dblock_v)) deallocate(g%Dblock_v)
  if (allocated(g%Dopen_v)) deallocate(g%Dopen_v)
 
  deallocate(g%gridLonT) ; deallocate(g%gridLatT)
  deallocate(g%gridLonB) ; deallocate(g%gridLatB)
 
  deallocate(g%Domain%mpp_domain)
  deallocate(g%Domain)
 
  deallocate(g)
 
end subroutine destroy_dyn_horgrid
 
end module mom_dyn_horgrid