coord_zlike.F90

Go to the documentation of this file.

!> Regrid columns for a z-like coordinate (z-star, z-level)
module coord_zlike

use mom_error_handler, only : mom_error, fatal

implicit none ; private

!> Control structure containing required parameters for a z-like coordinate
type, public :: zlike_cs
  private

  !> Number of levels
  integer :: nk

  !> Minimum thickness allowed for layers
  real :: min_thickness

  !> Target coordinate resolution
  real, allocatable, dimension(:) :: coordinateresolution
end type zlike_cs

public init_coord_zlike, set_zlike_params, build_zstar_column, end_coord_zlike

contains

!> Initialise a zlike_CS with pointers to parameters
subroutine init_coord_zlike(CS, nk, coordinateResolution)
  type(zlike_cs),     pointer    :: CS !< Unassociated pointer to hold the control structure
  integer,            intent(in) :: nk
  real, dimension(:), intent(in) :: coordinateResolution

  if (associated(cs)) call mom_error(fatal, "init_coord_zlike: CS already associated!")
  allocate(cs)
  allocate(cs%coordinateResolution(nk))

  cs%nk                   = nk
  cs%coordinateResolution = coordinateresolution
end subroutine init_coord_zlike

subroutine end_coord_zlike(CS)
  type(zlike_cs), pointer :: CS

  ! nothing to do
  if (.not. associated(cs)) return
  deallocate(cs%coordinateResolution)
  deallocate(cs)
end subroutine end_coord_zlike

subroutine set_zlike_params(CS, min_thickness)
  type(zlike_cs), pointer    :: CS
  real, optional, intent(in) :: min_thickness

  if (.not. associated(cs)) call mom_error(fatal, "set_zlike_params: CS not associated")

  if (present(min_thickness)) cs%min_thickness = min_thickness
end subroutine set_zlike_params

!> Builds a z* coordinate with a minimum thickness
subroutine build_zstar_column(CS, nz, depth, total_thickness, zInterface, z_rigid_top, eta_orig)
  type(zlike_cs),        intent(in)    :: CS !< Coordinate control structure
  integer,               intent(in)    :: nz !< Number of levels
  real,                  intent(in)    :: depth !< Depth of ocean bottom (positive in m)
  real,                  intent(in)    :: total_thickness !< Column thickness (positive in m)
  real, dimension(nz+1), intent(inout) :: zInterface !< Absolute positions of interfaces
  real, optional,        intent(in)    :: z_rigid_top !< The height of a rigid top (negative in m)
  real, optional,        intent(in)    :: eta_orig !< The actual original height of the top (m)
  ! Local variables
  real :: eta, stretching, dh, min_thickness, z0_top, z_star
  integer :: k
  logical :: new_zstar_def

  new_zstar_def = .false.
  min_thickness = min( cs%min_thickness, total_thickness/real(nz) )
  z0_top = 0.
  if (present(z_rigid_top)) then
    z0_top = z_rigid_top
    new_zstar_def = .true.
  endif

  ! Position of free-surface (or the rigid top, for which eta ~ z0_top)
  eta = total_thickness - depth
  if (present(eta_orig)) eta = eta_orig

  ! Conventional z* coordinate:
  !   z* = (z-eta) / stretching   where stretching = (H+eta)/H
  !   z = eta + stretching * z*
  ! The above gives z*(z=eta) = 0, z*(z=-H) = -H.
  ! With a rigid top boundary at eta = z0_top then
  !   z* = z0 + (z-eta) / stretching   where stretching = (H+eta)/(H+z0)
  !   z = eta + stretching * (z*-z0) * stretching
  stretching = total_thickness / ( depth + z0_top )

  if (new_zstar_def) then
    ! z_star is the notional z* coordinate in absence of upper/lower topography
    z_star = 0. ! z*=0 at the free-surface
    zinterface(1) = eta ! The actual position of the top of the column
    do k = 2,nz
      z_star = z_star - cs%coordinateResolution(k-1)
      ! This ensures that z is below a rigid upper surface (ice shelf bottom)
      zinterface(k) = min( eta + stretching * ( z_star - z0_top ), z0_top )
      ! This ensures that the layer in inflated
      zinterface(k) = min( zinterface(k), zinterface(k-1) - min_thickness )
      ! This ensures that z is above or at the topography
      zinterface(k) = max( zinterface(k), -depth + real(nz+1-k) * min_thickness )
    enddo
    zinterface(nz+1) = -depth

  else
    ! Integrate down from the top for a notional new grid, ignoring topography
    ! The starting position is offset by z0_top which, if z0_top<0, will place
    ! interfaces above the rigid boundary.
    zinterface(1) = eta
    do k = 1,nz
      dh = stretching * cs%coordinateResolution(k) ! Notional grid spacing
      zinterface(k+1) = zinterface(k) - dh
    enddo

    ! Integrating up from the bottom adjusting interface position to accommodate
    ! inflating layers without disturbing the interface above
    zinterface(nz+1) = -depth
    do k = nz,1,-1
      if ( zinterface(k) < (zinterface(k+1) + min_thickness) ) then
        zinterface(k) = zinterface(k+1) + min_thickness
      endif
    enddo
  endif

end subroutine build_zstar_column

end module coord_zlike