mom_interface_heights Module Reference

Detailed Description

Functions for calculating interface heights, including free surface height.

Data Types
interface	find_eta
	Calculates the heights of sruface or all interfaces from layer thicknesses. More...

Functions/Subroutines
subroutine	find_eta_3d (h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
	Calculates 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. More...
subroutine	find_eta_2d (h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
	Calculates 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. More...

Function/Subroutine Documentation

find_eta_2d()

subroutine mom_interface_heights::find_eta_2d ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h, type(thermo_var_ptrs), intent(in)  tv, type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)  eta, real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in), optional  eta_bt, integer, intent(in), optional  halo_size, real, intent(in), optional  eta_to_m  )

private

Calculates 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.

Parameters

[in]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	us	A dimensional unit scaling type
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in]	tv	A structure pointing to various thermodynamic variables.
[out]	eta	free surface height relative to mean sea level (z=0) often [Z ~> m].
[in]	eta_bt	optional barotropic variable that gives the "correct" free surface height (Boussinesq) or total water column mass per unit area (non-Boussinesq) [H ~> m or kg m-2].
[in]	halo_size	width of halo points on which to calculate eta.
[in]	eta_to_m	The conversion factor from the units of eta to m; by default this is USZ_to_m.

Definition at line 149 of file MOM_interface_heights.F90.

  type(ocean_grid_type),                      intent(in)  :: g   !< The ocean's grid structure.
  type(verticalgrid_type),                    intent(in)  :: gv  !< The ocean's vertical grid structure.
  type(unit_scale_type),                      intent(in)  :: us  !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: h   !< Layer thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs),                      intent(in)  :: tv  !< A structure pointing to various
                                                                 !! thermodynamic variables.
  real, dimension(SZI_(G),SZJ_(G)),           intent(out) :: eta !< free surface height relative to
                                                                 !! mean sea level (z=0) often [Z ~> 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 area (non-Boussinesq) [H ~> m or kg m-2].
  integer,                          optional, intent(in)  :: halo_size !< width of halo points on
                                                                 !! which to calculate eta.
  real,                             optional, intent(in)  :: eta_to_m  !< The conversion factor from
             !! the units of eta to m; by default this is US%Z_to_m.
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
    p          ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    dz_geo     ! The change in geopotential height across a layer [L2 T-2 ~> m2 s-2].
  real :: htot(szi_(g))  ! The sum of all layers' thicknesses [H ~> m or kg m-2].
  real :: i_gearth          ! The inverse of the gravitational acceleration times the
                            ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1]
  real :: z_to_eta, h_to_eta, h_to_rho_eta ! Unit conversion factors with obvious names.
  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

  z_to_eta = 1.0 ; if (present(eta_to_m)) z_to_eta = us%Z_to_m / eta_to_m
  h_to_eta = gv%H_to_Z * z_to_eta
  h_to_rho_eta =  gv%H_to_RZ * z_to_eta
  i_gearth = z_to_eta / gv%g_Earth

!$OMP parallel default(shared) private(htot)
!$OMP do
  do j=js,je ; do i=is,ie ; eta(i,j) = -z_to_eta*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) = h_to_eta*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)*h_to_eta
      enddo ; enddo ; enddo
    endif
  else
    if (associated(tv%eqn_of_state)) then
!$OMP do
      do j=js,je
        if (associated(tv%p_surf)) then
          do i=is,ie ; p(i,j,1) = tv%p_surf(i,j) ; enddo
        else
          do i=is,ie ; p(i,j,1) = 0.0 ; enddo
        endif

        do k=1,nz ; do i=is,ie
          p(i,j,k+1) = p(i,j,k) + gv%g_Earth*gv%H_to_RZ*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, us, 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) + h_to_rho_eta*h(i,j,k) / gv%Rlay(k)
      enddo ; enddo ; enddo
    endif
    if (present(eta_bt)) then
      !   Dilate the water column to agree with 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) + z_to_eta*g%bathyT(i,j)) - &
                     z_to_eta*g%bathyT(i,j)
        enddo
      enddo
    endif
  endif
!$OMP end parallel

find_eta_3d()

subroutine mom_interface_heights::find_eta_3d ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h, type(thermo_var_ptrs), intent(in)  tv, type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(out)  eta, real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in), optional  eta_bt, integer, intent(in), optional  halo_size, real, intent(in), optional  eta_to_m  )

private

Calculates 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.

Parameters

[in]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	us	A dimensional unit scaling type
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in]	tv	A structure pointing to various thermodynamic variables.
[out]	eta	layer interface heights [Z ~> m] or 1/eta_to_m m).
[in]	eta_bt	optional barotropic variable that gives the "correct" free surface height (Boussinesq) or total water column mass per unit area (non-Boussinesq). This is used to dilate the layer. thicknesses when calculating interfaceheights [H ~> m or kg m-2].
[in]	halo_size	width of halo points on which to calculate eta.
[in]	eta_to_m	The conversion factor from the units of eta to m; by default this is USZ_to_m.

Definition at line 32 of file MOM_interface_heights.F90.

  type(ocean_grid_type),                      intent(in)  :: g   !< The ocean's grid structure.
  type(verticalgrid_type),                    intent(in)  :: gv  !< The ocean's vertical grid structure.
  type(unit_scale_type),                      intent(in)  :: us  !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: h   !< Layer thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs),                      intent(in)  :: tv  !< A structure pointing to various
                                                                 !! thermodynamic variables.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: eta !< layer interface heights
                                                                 !! [Z ~> m] or 1/eta_to_m 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 area (non-Boussinesq).  This is used to dilate the layer.
             !! thicknesses when calculating interfaceheights [H ~> m or kg m-2].
  integer,                          optional, intent(in)  :: halo_size !< width of halo points on
                                                                 !! which to calculate eta.
  real,                             optional, intent(in)  :: eta_to_m  !< The conversion factor from
             !! the units of eta to m; by default this is US%Z_to_m.

  ! Local variables
  real :: p(szi_(g),szj_(g),szk_(g)+1)    ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa]
  real :: dz_geo(szi_(g),szj_(g),szk_(g)) ! The change in geopotential height
                                          ! across a layer [L2 T-2 ~> m2 s-2].
  real :: dilate(szi_(g))                 ! non-dimensional dilation factor
  real :: htot(szi_(g))                   ! total thickness [H ~> m or kg m-2]
  real :: i_gearth          ! The inverse of the gravitational acceleration times the
                            ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1]
  real :: z_to_eta, h_to_eta, h_to_rho_eta ! Unit conversion factors with obvious names.
  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.")

  z_to_eta = 1.0 ; if (present(eta_to_m)) z_to_eta = us%Z_to_m / eta_to_m
  h_to_eta = gv%H_to_Z * z_to_eta
  h_to_rho_eta =  gv%H_to_RZ * z_to_eta
  i_gearth = z_to_eta / gv%g_Earth

!$OMP parallel default(shared) private(dilate,htot)
!$OMP do
  do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -z_to_eta*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)*h_to_eta
    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)*h_to_eta + z_to_eta*g%bathyT(i,j)) / &
                      (eta(i,j,1) + z_to_eta*g%bathyT(i,j))
        enddo
        do k=1,nz ; do i=isv,iev
          eta(i,j,k) = dilate(i) * (eta(i,j,k) + z_to_eta*g%bathyT(i,j)) - z_to_eta*g%bathyT(i,j)
        enddo ; enddo
      enddo
    endif
  else
    if (associated(tv%eqn_of_state)) then
!$OMP do
      do j=jsv,jev
        if (associated(tv%p_surf)) then
          do i=isv,iev ; p(i,j,1) = tv%p_surf(i,j) ; enddo
        else
          do i=isv,iev ; p(i,j,1) = 0.0 ; enddo
        endif
        do k=1,nz ; do i=isv,iev
          p(i,j,k+1) = p(i,j,k) + gv%g_Earth*gv%H_to_RZ*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, us, 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) + h_to_rho_eta*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) + z_to_eta*g%bathyT(i,j)) - z_to_eta*g%bathyT(i,j)
        enddo ; enddo
      enddo
    endif
  endif
!$OMP end parallel