coord_zlike.F90

!> Regrid columns for a z-like coordinate (z-star, z-level)
module coord_zlike
 
! This file is part of MOM6. See LICENSE.md for the license.
 
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 to be generated
  integer :: nk
 
  !> Minimum thickness allowed for layers, in the same thickness units (perhaps [H ~> m or kg m-2])
  !! that will be used in all subsequent calls to build_zstar_column with this structure.
  real :: min_thickness
 
  !> Target coordinate resolution, usually in [Z ~> m]
  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 !< Number of levels in the grid
  real, dimension(:), intent(in) :: coordinateresolution !< Target coordinate resolution [Z ~> m]
 
  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
 
!> Deallocates the zlike control structure
subroutine end_coord_zlike(CS)
  type(zlike_cs), pointer :: cs !< Coordinate control structure
 
  ! Nothing to do
  if (.not. associated(cs)) return
  deallocate(cs%coordinateResolution)
  deallocate(cs)
end subroutine end_coord_zlike
 
!> Set parameters in the zlike structure
subroutine set_zlike_params(CS, min_thickness)
  type(zlike_cs), pointer    :: cs !< Coordinate control structure
  real, optional, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2]
 
  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, depth, total_thickness, zInterface, &
                              z_rigid_top, eta_orig, zScale)
  type(zlike_cs),           intent(in)    :: cs !< Coordinate control structure
  real,                     intent(in)    :: depth !< Depth of ocean bottom (positive in the output units)
  real,                     intent(in)    :: total_thickness !< Column thickness (positive in the same units as depth)
  real, dimension(CS%nk+1), intent(inout) :: zinterface !< Absolute positions of interfaces
  real, optional,           intent(in)    :: z_rigid_top !< The height of a rigid top (negative in the
                                                         !! same units as depth)
  real, optional,           intent(in)    :: eta_orig !< The actual original height of the top in the
                                                      !! same units as depth
  real, optional,           intent(in)    :: zscale !< Scaling factor from the target coordinate resolution
                                                    !! in Z to desired units for zInterface, perhaps Z_to_H
  ! Local variables
  real :: eta, stretching, dh, min_thickness, z0_top, z_star, z_scale
  integer :: k
  logical :: new_zstar_def
 
  z_scale = 1.0 ; if (present(zscale)) z_scale = zscale
 
  new_zstar_def = .false.
  min_thickness = min( cs%min_thickness, total_thickness/real(cs%nk) )
  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,cs%nk
      z_star = z_star - cs%coordinateResolution(k-1)*z_scale
      ! 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(cs%nk+1-k) * min_thickness )
    enddo
    zinterface(cs%nk+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,cs%nk
      dh = stretching * cs%coordinateResolution(k)*z_scale ! 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(cs%nk+1) = -depth
    do k = cs%nk,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