MOM_dyn_horgrid.F90

Go to the documentation of this file.

module mom_dyn_horgrid

!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

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

type, public :: dyn_horgrid_type
  type(mom_domain_type), pointer :: domain => null()
  type(mom_domain_type), pointer :: domain_aux => null() ! A non-symmetric auxiliary domain type.

  ! These elements can be copied from a provided hor_index_type.
  type(hor_index_type)  :: hi   ! Make this a pointer?
  integer :: isc, iec, jsc, jec ! The range of the computational domain indices
  integer :: isd, ied, jsd, jed ! and data domain indices at tracer cell centers.
  integer :: isg, ieg, jsg, jeg ! The range of the global domain tracer cell indices.
  integer :: iscb, iecb, jscb, jecb ! The range of the computational domain indices
  integer :: isdb, iedb, jsdb, jedb ! and data domain indices at tracer cell vertices.
  integer :: isgb, iegb, jsgb, jegb ! The range of the global domain vertex indices.
  integer :: isd_global         ! The values of isd and jsd in the global
  integer :: jsd_global         ! (decomposition invariant) index space.
  integer :: idg_offset         ! The offset between the corresponding global
  integer :: jdg_offset         ! and local array 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. 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, idxt, & ! dxT is delta x at h points, in m, and IdxT is 1/dxT in m-1.
    dyt, idyt, & ! dyT is delta y at h points, in m, and IdyT is 1/dyT in m-1.
    areat, &     ! areaT is the area of an h-cell, in m2.
    iareat, &    ! IareaT = 1/areaT, in m-2.
    sin_rot, &   ! The sine and cosine of the angular rotation between the local
    cos_rot      ! model grid's northward and the true northward directions.

  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 in degrees of latitude or m.
    geoloncu, &  ! The geographic longitude at u points in degrees of longitude or m.
    dxcu, idxcu, & ! dxCu is delta x at u points, in m, and IdxCu is 1/dxCu in m-1.
    dycu, idycu, & ! dyCu is delta y at u points, in m, and IdyCu is 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.

  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 in degrees of latitude or m.
    geoloncv, &  !  The geographic longitude at v points in degrees of longitude or m.
    dxcv, idxcv, & ! dxCv is delta x at v points, in m, and IdxCv is 1/dxCv in m-1.
    dycv, idycv, & ! dyCv is delta y at v points, in m, and IdyCv is 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, allocatable, dimension(:,:) :: &
    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, idxbu, & ! dxBu is delta x at q points, in m, and IdxBu is 1/dxBu in m-1.
    dybu, idybu, & ! dyBu is delta y at q points, in m, and IdyBu is 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(), gridlatb => null() ! The latitude of T or B points for
                        ! the purpose of labeling the output axes.
                        ! On many grids these are the same as geoLatT & geoLatBu.
  real, pointer, dimension(:) :: &
    gridlont => null(), gridlonb => null() ! The longitude of T or B points for
                        ! the purpose of labeling the output axes.
                        ! On many grids these are the same as geoLonT & geoLonBu.
  character(len=40) :: &
    x_axis_units, &     !   The units that are used in labeling the coordinate
    y_axis_units        ! axes.  Except on a Cartesian grid, these are usually
                        ! some variant of "degrees".

  real, allocatable, dimension(:,:) :: &
    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, allocatable, dimension(:,:) :: &
    dblock_u, &   ! Topographic depths at u-points at which the flow is blocked
    dopen_u       ! (Dblock_u) and open at width dy_Cu (Dopen_u), both in m.
  real, allocatable, dimension(:,:) :: &
    dblock_v, &   ! Topographic depths at v-points at which the flow is blocked
    dopen_v       ! (Dblock_v) and open at width dx_Cv (Dopen_v), both in m.
  real, allocatable, dimension(:,:) :: &
    coriolisbu    ! The Coriolis parameter at corner points, in s-1.
  real, allocatable, dimension(:,:) :: &
    df_dx, df_dy  ! Derivatives of f (Coriolis parameter) at h-points, in 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 in m2
  real :: iareat_global ! Global sum of inverse h-cell area (1/areaT_global)
                        ! in m2

  ! 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 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%dx_Cv_obc(isd:ied,jsdb:jedb)) ; g%dx_Cv_obc(:,:) = 0.0
  allocate(g%dy_Cu_obc(isdb:iedb,jsd:jed)) ; g%dy_Cu_obc(:,:) = 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

!> 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%dx_Cv_obc) ; deallocate(g%dy_Cu_obc)

  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