coord_hycom.F90

Go to the documentation of this file.

!> Regrid columns for the HyCOM coordinate
module coord_hycom

use mom_error_handler, only : mom_error, fatal
use mom_eos,           only : eos_type, calculate_density
use regrid_interp,     only : interp_cs_type, build_and_interpolate_grid

implicit none ; private

!> Control structure containing required parameters for the HyCOM coordinate
type, public :: hycom_cs
  private

  !> Number of layers/levels
  integer :: nk

  !> Nominal near-surface resolution
  real, allocatable, dimension(:) :: coordinateresolution

  !> Nominal density of interfaces
  real, allocatable, dimension(:) :: target_density

  !> Maximum depths of interfaces
  real, allocatable, dimension(:) :: max_interface_depths

  !> Maximum thicknesses of layers
  real, allocatable, dimension(:) :: max_layer_thickness

  !> Interpolation control structure
  type(interp_cs_type) :: interp_cs
end type hycom_cs

public init_coord_hycom, set_hycom_params, build_hycom1_column, end_coord_hycom

contains

!> Initialise a hycom_CS with pointers to parameters
subroutine init_coord_hycom(CS, nk, coordinateResolution, target_density, interp_CS)
  type(hycom_cs),       pointer    :: CS !< Unassociated pointer to hold the control structure
  integer,              intent(in) :: nk !< Number of layers in generated grid
  real, dimension(:),   intent(in) :: coordinateResolution !< Z-space thicknesses (m)
  real, dimension(:),   intent(in) :: target_density !< Interface target densities (kg/m3)
  type(interp_cs_type), intent(in) :: interp_CS !< Controls for interpolation

  if (associated(cs)) call mom_error(fatal, "init_coord_hycom: CS already associated!")
  allocate(cs)
  allocate(cs%coordinateResolution(nk))
  allocate(cs%target_density(nk+1))

  cs%nk                      = nk
  cs%coordinateResolution(:) = coordinateresolution(:)
  cs%target_density(:)       = target_density(:)
  cs%interp_CS               = interp_cs
end subroutine init_coord_hycom

subroutine end_coord_hycom(CS)
  type(hycom_cs), pointer :: CS

  ! nothing to do
  if (.not. associated(cs)) return
  deallocate(cs%coordinateResolution)
  deallocate(cs%target_density)
  if (allocated(cs%max_interface_depths)) deallocate(cs%max_interface_depths)
  if (allocated(cs%max_layer_thickness)) deallocate(cs%max_layer_thickness)
  deallocate(cs)
end subroutine end_coord_hycom

subroutine set_hycom_params(CS, max_interface_depths, max_layer_thickness, interp_CS)
  type(hycom_cs),                 pointer    :: CS
  real, optional, dimension(:),   intent(in) :: max_interface_depths
  real, optional, dimension(:),   intent(in) :: max_layer_thickness
  type(interp_cs_type), optional, intent(in) :: interp_CS

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

  if (present(max_interface_depths)) then
    if (size(max_interface_depths) /= cs%nk+1) &
      call mom_error(fatal, "set_hycom_params: max_interface_depths inconsistent size")
    allocate(cs%max_interface_depths(cs%nk+1))
    cs%max_interface_depths(:) = max_interface_depths(:)
  endif

  if (present(max_layer_thickness)) then
    if (size(max_layer_thickness) /= cs%nk) &
      call mom_error(fatal, "set_hycom_params: max_layer_thickness inconsistent size")
    allocate(cs%max_layer_thickness(cs%nk))
    cs%max_layer_thickness(:) = max_layer_thickness(:)
  endif

  if (present(interp_cs)) cs%interp_CS = interp_cs
end subroutine set_hycom_params

!> Build a HyCOM coordinate column
subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, z_col, z_col_new)
  type(hycom_cs),        intent(in)    :: CS !< Coordinate control structure
  type(eos_type),        pointer       :: eqn_of_state !< Equation of state structure
  integer,               intent(in)    :: nz !< Number of levels
  real,                  intent(in)    :: depth !< Depth of ocean bottom (positive in H)
  real, dimension(nz),   intent(in)    :: T, S !< T and S for column
  real, dimension(nz),   intent(in)    :: h  !< Layer thicknesses, in m
  real, dimension(nz),   intent(in)    :: p_col !< Layer pressure in Pa
  real, dimension(nz+1), intent(in)    :: z_col ! Interface positions relative to the surface in H units (m or kg m-2)
  real, dimension(nz+1), intent(inout) :: z_col_new !< Absolute positions of interfaces

  ! Local variables
  integer   :: k
  real, dimension(nz) :: rho_col, h_col_new ! Layer quantities
  real :: stretching ! z* stretching, converts z* to z.
  real :: nominal_z ! Nominal depth of interface is using z* (m or Pa)
  real :: hNew
  logical :: maximum_depths_set ! If true, the maximum depths of interface have been set.
  logical :: maximum_h_set      ! If true, the maximum layer thicknesses have been set.

  maximum_depths_set = allocated(cs%max_interface_depths)
  maximum_h_set = allocated(cs%max_layer_thickness)

  ! Work bottom recording potential density
  call calculate_density(t, s, p_col, rho_col, 1, nz, eqn_of_state)
  ! This ensures the potential density profile is monotonic
  ! although not necessarily single valued.
  do k = nz-1, 1, -1
    rho_col(k) = min( rho_col(k), rho_col(k+1) )
  enddo

  ! Interpolates for the target interface position with the rho_col profile
  ! Based on global density profile, interpolate to generate a new grid
  call build_and_interpolate_grid(cs%interp_CS, rho_col, nz, h(:), z_col, &
       cs%target_density, nz, h_col_new, z_col_new)

  ! Sweep down the interfaces and make sure that the interface is at least
  ! as deep as a nominal target z* grid
  nominal_z = 0.
  stretching = z_col(nz+1) / depth ! Stretches z* to z
  do k = 2, nz+1
    nominal_z = nominal_z + cs%coordinateResolution(k-1) * stretching
    z_col_new(k) = max( z_col_new(k), nominal_z )
    z_col_new(k) = min( z_col_new(k), z_col(nz+1) )
  enddo

  if (maximum_depths_set .and. maximum_h_set) then ; do k=2,nz
    ! The loop bounds are 2 & nz so the top and bottom interfaces do not move.
    ! Recall that z_col_new is positive downward.
    z_col_new(k) = min(z_col_new(k), cs%max_interface_depths(k), &
                       z_col_new(k-1) + cs%max_layer_thickness(k-1))
  enddo ; elseif (maximum_depths_set) then ; do k=2,nz
    z_col_new(k) = min(z_col_new(k), cs%max_interface_depths(k))
  enddo ; elseif (maximum_h_set) then ; do k=2,nz
    z_col_new(k) = min(z_col_new(k), z_col_new(k-1) + cs%max_layer_thickness(k-1))
  enddo ; endif
end subroutine build_hycom1_column

end module coord_hycom