MOM_interface_heights.F90

Go to the documentation of this file.

!> The module calculates interface heights, including free surface height.
module mom_interface_heights

!***********************************************************************
!*                   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                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, February 2001                                  *
!*                                                                     *
!*    The subroutines here calculate the heights of interfaces in a    *
!*  way that is consistent with the calculations in the pressure       *
!*  gradient acceleration calculation.  In a Boussinseq model this is  *
!*  pretty simple vertical sum, but in a non-Boussinesq model it uses  *
!*  integrals of the equation of state.                                *
!*                                                                     *
!*  Macros written all in capital letters are defined in MOM_memory.h. *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q, CoriolisBu                            *
!*    j+1  > o > o >   At ^:  v                                        *
!*    j    x ^ x ^ x   At >:  u                                        *
!*    j    > o > o >   At o:  h, bathyT                                *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use mom_error_handler, only : mom_error, fatal
use mom_file_parser, only : log_version
use mom_grid, only : ocean_grid_type
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : int_specific_vol_dp

implicit none ; private

#include <MOM_memory.h>

public find_eta

interface find_eta
  module procedure find_eta_2d, find_eta_3d
end interface find_eta

contains

subroutine find_eta_3d(h, tv, G_Earth, G, GV, eta, eta_bt, halo_size)
  type(ocean_grid_type),                      intent(in)  :: G         !< The ocean's grid
                                                                       !! structure.
  type(verticalgrid_type),                    intent(in)  :: GV        !< The ocean's vertical
                                                                       !! grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: h         !< Layer thicknesses, in H
                                                                       !! (usually m or kg m-2).
  type(thermo_var_ptrs),                      intent(in)  :: tv        !< A structure pointing to
                                                                       !! various thermodynamic
                                                                       !! variables.
  real,                                       intent(in)  :: G_Earth   !< Earth gravitational
                                                                       !! acceleration (m/s2).
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: eta       !< layer interface heights
                                                                       !! (meter).
  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in)  :: eta_bt    !< optional barotropic
             !! variable that gives the "correct" free surface height (Boussinesq) or total water
             !! column mass per unit aread (non-Boussinesq).  This is used to dilate the layer.
             !! thicknesses when calculating interfaceheights, in H (m or kg m-2).
  integer,                          optional, intent(in)  :: halo_size !< width of halo points on
                                                                       !! which to calculate eta.

!   This subroutine determines the heights of all interfaces between layers,
! using the appropriate form for consistency with the calculation of the
! pressure gradient forces.  Additionally, these height may be dilated for
! consistency with the corresponding time-average quantity from the barotropic
! calculation.

! Arguments:
!  (in)      h       - layer thickness H (meter or kg/m2)
!  (in)      tv      - structure pointing to thermodynamic variables
!  (in)      G_Earth - Earth gravitational acceleration (m/s2)
!  (in)      G       - ocean grid structure
!  (in)      GV      - The ocean's vertical grid structure.
!  (out)     eta     - layer interface heights (meter)
!  (in,opt)  eta_bt  - optional barotropic variable that gives the "correct"
!                      free surface height (Boussinesq) or total water column
!                      mass per unit aread (non-Boussinesq).  This is used to
!                      dilate the layer thicknesses when calculating interface
!                      heights, in H (m or kg m-2).
! (in,opt) halo_size - width of halo points on which to calculate eta

  real :: p(szi_(g),szj_(g),szk_(g)+1)
  real :: dz_geo(szi_(g),szj_(g),szk_(g)) ! The change in geopotential height
                                          ! across a layer, in m2 s-2.
  real :: dilate(szi_(g))                 ! non-dimensional dilation factor
  real :: htot(szi_(g))                   ! total thickness H
  real :: I_gEarth
  integer i, j, k, isv, iev, jsv, jev, nz, halo

  halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)

  isv = g%isc-halo ; iev = g%iec+halo ; jsv = g%jsc-halo ; jev = g%jec+halo
  nz = g%ke

  if ((isv<g%isd) .or. (iev>g%ied) .or. (jsv<g%jsd) .or. (jev>g%jed)) &
    call mom_error(fatal,"find_eta called with an overly large halo_size.")

  i_gearth = 1.0 / g_earth

!$OMP parallel default(none) shared(isv,iev,jsv,jev,nz,eta,G,GV,h,eta_bt,tv,p, &
!$OMP                               G_Earth,dz_geo,halo,I_gEarth) &
!$OMP                       private(dilate,htot)
!$OMP do
  do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -g%bathyT(i,j) ; enddo ; enddo

  if (gv%Boussinesq) then
!$OMP do
    do j=jsv,jev ; do k=nz,1,-1; do i=isv,iev
      eta(i,j,k) = eta(i,j,k+1) + h(i,j,k)*gv%H_to_m
    enddo ; enddo ; enddo
    if (present(eta_bt)) then
      ! Dilate the water column to agree with the free surface height
      ! that is used for the dynamics.
!$OMP do
      do j=jsv,jev
        do i=isv,iev
          dilate(i) = (eta_bt(i,j)*gv%H_to_m + g%bathyT(i,j)) / &
                      (eta(i,j,1) + g%bathyT(i,j))
        enddo
        do k=1,nz ; do i=isv,iev
          eta(i,j,k) = dilate(i) * (eta(i,j,k) + g%bathyT(i,j)) - g%bathyT(i,j)
        enddo ; enddo
      enddo
    endif
  else
    if (associated(tv%eqn_of_state)) then
      ! ### THIS SHOULD BE P_SURF, IF AVAILABLE.
!$OMP do
      do j=jsv,jev
        do i=isv,iev ; p(i,j,1) = 0.0 ; enddo
        do k=1,nz ; do i=isv,iev
          p(i,j,k+1) = p(i,j,k) + g_earth*gv%H_to_kg_m2*h(i,j,k)
        enddo ; enddo
      enddo
!$OMP do
      do k=1,nz
        call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), &
                                 0.0, g%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo)
      enddo
!$OMP do
      do j=jsv,jev
        do k=nz,1,-1 ; do i=isv,iev
          eta(i,j,k) = eta(i,j,k+1) + i_gearth * dz_geo(i,j,k)
        enddo ; enddo
      enddo
    else
!$OMP do
      do j=jsv,jev ;  do k=nz,1,-1; do i=isv,iev
        eta(i,j,k) = eta(i,j,k+1) + gv%H_to_kg_m2*h(i,j,k)/gv%Rlay(k)
      enddo ; enddo ; enddo
    endif
    if (present(eta_bt)) then
      ! Dilate the water column to agree with the free surface height
      ! from the time-averaged barotropic solution.
!$OMP do
      do j=jsv,jev
        do i=isv,iev ; htot(i) = gv%H_subroundoff ; enddo
        do k=1,nz ; do i=isv,iev ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
        do i=isv,iev ; dilate(i) = eta_bt(i,j) / htot(i) ; enddo
        do k=1,nz ; do i=isv,iev
          eta(i,j,k) = dilate(i) * (eta(i,j,k) + g%bathyT(i,j)) - g%bathyT(i,j)
        enddo ; enddo
      enddo
    endif
  endif
!$OMP end parallel

end subroutine find_eta_3d


subroutine find_eta_2d(h, tv, G_Earth, G, GV, eta, eta_bt, halo_size)
  type(ocean_grid_type),                      intent(in)  :: G        !< The ocean's grid structure.
  type(verticalgrid_type),                    intent(in)  :: GV       !< The ocean's vertical grid
                                                                      !! structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: h        !< Layer thicknesses, in H
                                                                      !! (usually m or kg m-2).
  type(thermo_var_ptrs),                      intent(in)  :: tv       !< A structure pointing to
                                                                      !! various thermodynamic
                                                                      !! variables.
  real,                                       intent(in)  :: G_Earth  !< Earth gravitational
                                                                      !! acceleration (m/s2).
  real, dimension(SZI_(G),SZJ_(G)),           intent(out) :: eta      !< free surface height
                                                                      !! relative to mean sea
                                                                      !! level (z=0) (m).
  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in)  :: eta_bt   !< optional barotropic
                   !! variable that gives the "correct" free surface height (Boussinesq) or total
                   !! water column mass per unit aread (non-Boussinesq), in H (m or kg m-2).
  integer,                          optional, intent(in)  :: halo_size !< width of halo points on
                                                                       !! which to calculate eta.

!   This subroutine determines the free surface height, using the appropriate
! form for consistency with the calculation of the pressure gradient forces.
! Additionally, the sea surface height may be adjusted for consistency with the
! corresponding time-average quantity from the barotropic calculation.

! Arguments:
!  (in)       h      - layer thickness (meter or kg/m2)
!  (in)      tv      - structure pointing to various thermodynamic variables
!  (in)      G_Earth - Earth gravitational acceleration (m/s2)
!  (in)      G       - ocean grid structure
!  (in)      GV      - The ocean's vertical grid structure.
!  (out)     eta     - free surface height relative to mean sea level (z=0) (m)
!  (in,opt)  eta_bt  - optional barotropic variable that gives the "correct"
!                      free surface height (Boussinesq) or total water column
!                      mass per unit aread (non-Boussinesq), in H (m or kg m-2)
!  (in,opt)  halo_size - width of halo points on which to calculate eta

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
    p     ! The pressure in Pa.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    dz_geo     ! The change in geopotential height across a layer, in m2 s-2.
  real :: htot(szi_(g))  ! The sum of all layers' thicknesses, in kg m-2 or m.
  real :: I_gEarth
  integer i, j, k, is, ie, js, je, nz, halo

  halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
  nz = g%ke

  i_gearth = 1.0 / g_earth

!$OMP parallel default(none) shared(is,ie,js,je,nz,eta,G,GV,eta_bt,h,tv,p, &
!$OMP                               G_Earth,dz_geo,halo,I_gEarth) &
!$OMP                       private(htot)
!$OMP do
  do j=js,je ; do i=is,ie ; eta(i,j) = -g%bathyT(i,j) ; enddo ; enddo

  if (gv%Boussinesq) then
    if (present(eta_bt)) then
!$OMP do
      do j=js,je ; do i=is,ie
        eta(i,j) = eta_bt(i,j)
      enddo ; enddo
    else
!$OMP do
      do j=js,je ; do k=1,nz ; do i=is,ie
        eta(i,j) = eta(i,j) + h(i,j,k)*gv%H_to_m
      enddo ; enddo ; enddo
    endif
  else
    if (associated(tv%eqn_of_state)) then
!$OMP do
      do j=js,je
        do i=is,ie ; p(i,j,1) = 0.0 ; enddo

        do k=1,nz ; do i=is,ie
          p(i,j,k+1) = p(i,j,k) + g_earth*gv%H_to_kg_m2*h(i,j,k)
        enddo ; enddo
      enddo
!$OMP do
      do k = 1, nz
        call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), 0.0, &
                                 g%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo)
      enddo
!$OMP do
      do j=js,je ; do k=1,nz ; do i=is,ie
          eta(i,j) = eta(i,j) + i_gearth * dz_geo(i,j,k)
      enddo ; enddo ; enddo
    else
!$OMP do
      do j=js,je ; do k=1,nz ; do i=is,ie
        eta(i,j) = eta(i,j) + gv%H_to_kg_m2*h(i,j,k)/gv%Rlay(k)
      enddo ; enddo ; enddo
    endif
    if (present(eta_bt)) then
      !   Dilate the water column to agree with the the time-averaged column
      ! mass from the barotropic solution.
!$OMP do
      do j=js,je
        do i=is,ie ; htot(i) = gv%H_subroundoff ; enddo
        do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
        do i=is,ie
          eta(i,j) = (eta_bt(i,j) / htot(i)) * (eta(i,j) + g%bathyT(i,j)) - &
                     g%bathyT(i,j)
        enddo
      enddo
    endif
  endif
!$OMP end parallel

end subroutine find_eta_2d

end module mom_interface_heights