coord_rho.F90

Go to the documentation of this file.

!> Regrid columns for the continuous isopycnal (rho) coordinate
module coord_rho

use mom_error_handler, only : mom_error, fatal
use mom_remapping,     only : remapping_cs, remapping_core_h
use mom_eos,           only : eos_type, calculate_density
use regrid_interp,     only : interp_cs_type, build_and_interpolate_grid, degree_max

implicit none ; private

!> Control structure containing required parameters for the rho coordinate
type, public :: rho_cs
  private

  !> Number of layers
  integer :: nk

  !> Minimum thickness allowed for layers
  real :: min_thickness

  !> Reference pressure for density calculations
  real :: ref_pressure

  !> If true, integrate for interface positions from the top downward.
  !! If false, integrate from the bottom upward, as does the rest of the model.
  logical :: integrate_downward_for_e

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

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

!> Maximum number of regridding iterations
integer, parameter :: nb_regridding_iterations = 1
!> Deviation tolerance between succesive grids in regridding iterations
real, parameter    :: deviation_tolerance = 1e-10
! This CPP macro embeds some safety checks

public init_coord_rho, set_rho_params, build_rho_column, old_inflate_layers_1d, end_coord_rho

contains

!> Initialise a rho_CS with pointers to parameters
subroutine init_coord_rho(CS, nk, ref_pressure, target_density, interp_CS)
  type(rho_cs),         pointer    :: CS !< Unassociated pointer to hold the control structure
  integer,              intent(in) :: nk
  real,                 intent(in) :: ref_pressure
  real, dimension(:),   intent(in) :: target_density
  type(interp_cs_type), intent(in) :: interp_CS

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

  cs%nk                = nk
  cs%ref_pressure      = ref_pressure
  cs%target_density(:) = target_density(:)
  cs%interp_CS         = interp_cs
end subroutine init_coord_rho

subroutine end_coord_rho(CS)
  type(rho_cs), pointer :: CS

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

subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, interp_CS)
  type(rho_cs),                   pointer    :: CS
  real,                 optional, intent(in) :: min_thickness
  logical,              optional, intent(in) :: integrate_downward_for_e
  type(interp_cs_type), optional, intent(in) :: interp_CS

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

  if (present(min_thickness)) cs%min_thickness = min_thickness
  if (present(integrate_downward_for_e)) cs%integrate_downward_for_e = integrate_downward_for_e
  if (present(interp_cs)) cs%interp_CS = interp_cs
end subroutine set_rho_params

subroutine build_rho_column(CS, remapCS, nz, depth, h, T, S, eqn_of_state, zInterface)
  !< Build a rho coordinate column
  !!
  !! The algorithn operates as follows within each column:
  !!
  !! 1. Given T & S within each layer, the layer densities are computed.
  !! 2. Based on these layer densities, a global density profile is reconstructed
  !!    (this profile is monotonically increasing and may be discontinuous)
  !! 3. The new grid interfaces are determined based on the target interface
  !!    densities.
  !! 4. T & S are remapped onto the new grid.
  !! 5. Return to step 1 until convergence or until the maximum number of
  !!    iterations is reached, whichever comes first.

  type(rho_cs),          intent(in)    :: CS !< Regridding control structure
  type(remapping_cs),    intent(in)    :: remapCS !< Remapping parameters and options
  integer,               intent(in)    :: nz !< Number of levels
  real,                  intent(in)    :: depth !< Depth of ocean bottom (positive in m)
  real, dimension(nz),   intent(in)    :: h  !< Layer thicknesses, in m
  real, dimension(nz),   intent(in)    :: T, S !< T and S for column
  type(eos_type),        pointer       :: eqn_of_state !< Equation of state structure
  real, dimension(nz+1), intent(inout) :: zInterface !< Absolute positions of interfaces

  ! Local variables
  integer   :: k, m
  integer   :: map_index
  integer   :: k_found
  integer   :: count_nonzero_layers
  real      :: deviation            ! When iterating to determine the final
                                    ! grid, this is the deviation between two
                                    ! successive grids.
  real      :: threshold
  real      :: max_thickness
  real      :: correction
  real, dimension(nz) :: p, densities, T_tmp, S_tmp, Tmp
  integer, dimension(nz) :: mapping
  real :: dh
  real, dimension(nz) :: h0, h1, hTmp
  real, dimension(nz+1) :: x0, x1, xTmp

  threshold = cs%min_thickness
  p(:) = cs%ref_pressure
  t_tmp(:) = t(:)
  s_tmp(:) = s(:)
  h0(:) = h(:)

  ! Start iterations to build grid
  m = 1
  deviation = 1e10
  do while ( ( m <= nb_regridding_iterations ) .and. &
             ( deviation > deviation_tolerance ) )

    ! Count number of nonzero layers within current water column
    count_nonzero_layers = 0
    do k = 1,nz
      if ( h0(k) > threshold ) then
        count_nonzero_layers = count_nonzero_layers + 1
      end if
    end do

    ! If there is at most one nonzero layer, stop here (no regridding)
    if ( count_nonzero_layers <= 1 ) then
      h1(:) = h0(:)
      exit  ! stop iterations here
    end if

    ! Build new grid containing only nonzero layers
    map_index = 1
    correction = 0.0
    do k = 1,nz
      if ( h0(k) > threshold ) then
        mapping(map_index) = k
        htmp(map_index) = h0(k)
        map_index = map_index + 1
      else
        correction = correction + h0(k)
      end if
    end do

    max_thickness = htmp(1)
    k_found = 1
    do k = 1,count_nonzero_layers
      if ( htmp(k) > max_thickness ) then
        max_thickness = htmp(k)
        k_found = k
      end if
    end do

    htmp(k_found) = htmp(k_found) + correction

    xtmp(1) = 0.0
    do k = 1,count_nonzero_layers
      xtmp(k+1) = xtmp(k) + htmp(k)
    end do

    ! Compute densities within current water column
    call calculate_density( t_tmp, s_tmp, p, densities,&
                             1, nz, eqn_of_state )

    do k = 1,count_nonzero_layers
      densities(k) = densities(mapping(k))
    end do

    ! One regridding iteration
    ! Based on global density profile, interpolate to generate a new grid
    call build_and_interpolate_grid(cs%interp_CS, densities, count_nonzero_layers, &
         htmp, xtmp, cs%target_density, nz, h1, x1)

    call old_inflate_layers_1d( cs%min_thickness, nz, h1 )
    x1(1) = 0.0 ; do k = 1,nz ; x1(k+1) = x1(k) + h1(k) ; end do

    ! Remap T and S from previous grid to new grid
    do k = 1,nz
      h1(k) = x1(k+1) - x1(k)
    end do

    call remapping_core_h(remapcs, nz, h0, s, nz, h1, tmp)
    s_tmp(:) = tmp(:)

    call remapping_core_h(remapcs, nz, h0, t, nz, h1, tmp)
    t_tmp(:) = tmp(:)

    ! Compute the deviation between two successive grids
    deviation = 0.0
    x0(1) = 0.0
    x1(1) = 0.0
    do k = 2,nz
      x0(k) = x0(k-1) + h0(k-1)
      x1(k) = x1(k-1) + h1(k-1)
      deviation = deviation + (x0(k)-x1(k))**2
    end do
    deviation = sqrt( deviation / (nz-1) )

    m = m + 1

    ! Copy final grid onto start grid for next iteration
    h0(:) = h1(:)

  end do ! end regridding iterations

  if (cs%integrate_downward_for_e) then
    zinterface(1) = 0.
    do k = 1,nz
      zinterface(k+1) = zinterface(k) - h1(k)
      ! Adjust interface position to accomodate inflating layers
      ! without disturbing the interface above
    enddo
  else
    ! The rest of the model defines grids integrating up from the bottom
    zinterface(nz+1) = -depth
    do k = nz,1,-1
      zinterface(k) = zinterface(k+1) + h1(k)
      ! Adjust interface position to accomodate inflating layers
      ! without disturbing the interface above
    enddo
  endif

end subroutine build_rho_column

!------------------------------------------------------------------------------
! Inflate vanished layers to finite (nonzero) width
!------------------------------------------------------------------------------
subroutine old_inflate_layers_1d( minThickness, N, h )

  ! Argument
  real,                intent(in) :: minThickness
  integer,             intent(in) :: N
  real,                intent(inout) :: h(:)

  ! Local variable
  integer   :: k
  integer   :: k_found
  integer   :: count_nonzero_layers
  real      :: delta
  real      :: correction
  real      :: maxThickness

  ! Count number of nonzero layers
  count_nonzero_layers = 0
  do k = 1,n
    if ( h(k) > minthickness ) then
      count_nonzero_layers = count_nonzero_layers + 1
    end if
  end do

  ! If all layer thicknesses are greater than the threshold, exit routine
  if ( count_nonzero_layers == n ) return

  ! If all thicknesses are zero, inflate them all and exit
  if ( count_nonzero_layers == 0 ) then
    do k = 1,n
      h(k) = minthickness
    end do
    return
  end if

  ! Inflate zero layers
  correction = 0.0
  do k = 1,n
    if ( h(k) <= minthickness ) then
      delta = minthickness - h(k)
      correction = correction + delta
      h(k) = h(k) + delta
    end if
  end do

  ! Modify thicknesses of nonzero layers to ensure volume conservation
  maxthickness = h(1)
  k_found = 1
  do k = 1,n
    if ( h(k) > maxthickness ) then
      maxthickness = h(k)
      k_found = k
    end if
  end do

  h(k_found) = h(k_found) - correction

end subroutine old_inflate_layers_1d

end module coord_rho