coord_adapt Module Reference

Detailed Description

Regrid columns for the adaptive coordinate.

Data Types
type	adapt_cs
	Control structure for adaptive coordinates (coord_adapt). More...

Functions/Subroutines
subroutine, public	init_coord_adapt (CS, nk, coordinateResolution, m_to_H, kg_m3_to_R)
	Initialise an adapt_CS with parameters. More...
subroutine, public	end_coord_adapt (CS)
	Clean up the coordinate control structure. More...
subroutine, public	set_adapt_params (CS, adaptTimeRatio, adaptAlpha, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptDrho0, adaptDoMin)
	This subtroutine can be used to set the parameters for coord_adapt module. More...
subroutine, public	build_adapt_column (CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNext)

Function/Subroutine Documentation

build_adapt_column()

subroutine, public coord_adapt::build_adapt_column	(	type(adapt_cs), intent(in)	CS,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(thermo_var_ptrs), intent(in)	tv,
		integer, intent(in)	i,
		integer, intent(in)	j,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke+1), intent(in)	zInt,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke+1), intent(in)	tInt,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke+1), intent(in)	sInt,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)	h,
		real, dimension( gv %ke+1), intent(inout)	zNext
	)

Parameters

[in]	cs	The control structure for this module
[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	tv	A structure pointing to various thermodynamic variables
[in]	i	The i-index of the column to work on
[in]	j	The j-index of the column to work on
[in]	zint	Interface heights [H ~> m or kg m-2].
[in]	tint	Interface temperatures [degC]
[in]	sint	Interface salinities [ppt]
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in,out]	znext	updated interface positions

Definition at line 116 of file coord_adapt.F90.

  type(adapt_cs),                              intent(in)    :: cs   !< The control structure for this module
  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
  type(thermo_var_ptrs),                       intent(in)    :: tv   !< A structure pointing to various
                                                                     !! thermodynamic variables
  integer,                                     intent(in)    :: i    !< The i-index of the column to work on
  integer,                                     intent(in)    :: j    !< The j-index of the column to work on
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in)    :: zint !< Interface heights [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in)    :: tint !< Interface temperatures [degC]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in)    :: sint !< Interface salinities [ppt]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZK_(GV)+1),                 intent(inout) :: znext !< updated interface positions

  ! Local variables
  integer :: k, nz
  real :: h_up, b1, b_denom_1, d1, depth, nominal_z, stretching
  real :: drdz  ! The vertical density gradient [R H-1 ~> kg m-4 or m-1]
  real, dimension(SZK_(GV)+1) :: alpha ! drho/dT [R degC-1 ~> kg m-3 degC-1]
  real, dimension(SZK_(GV)+1) :: beta  ! drho/dS [R ppt-1 ~> kg m-3 ppt-1]
  real, dimension(SZK_(GV)+1) :: del2sigma ! Laplacian of in situ density times grid spacing [R ~> kg m-3]
  real, dimension(SZK_(GV)+1) :: dh_d2s ! Thickness change in response to del2sigma [H ~> m or kg m-2]
  real, dimension(SZK_(GV)) :: kgrid, c1 ! grid diffusivity on layers, and tridiagonal work array

  nz = cs%nk

  ! set bottom and surface of zNext
  znext(1) = 0.
  znext(nz+1) = zint(i,j,nz+1)

  ! local depth for scaling diffusivity
  depth = g%bathyT(i,j) * gv%Z_to_H

  ! initialize del2sigma and the thickness change response to it zero
  del2sigma(:) = 0.0 ; dh_d2s(:) = 0.0

  ! calculate del-squared of neutral density by a
  ! stencilled finite difference
  ! TODO: this needs to be adjusted to account for vanished layers near topography

  ! up (j-1)
  if (g%mask2dT(i,j-1) > 0.) then
    call calculate_density_derivs( &
         0.5 * (tint(i,j,2:nz) + tint(i,j-1,2:nz)), &
         0.5 * (sint(i,j,2:nz) + sint(i,j-1,2:nz)), &
         0.5 * (zint(i,j,2:nz) + zint(i,j-1,2:nz)) * (gv%H_to_RZ * gv%g_Earth), &
         alpha, beta, tv%eqn_of_state, (/2,nz/) )

    del2sigma(2:nz) = del2sigma(2:nz) + &
         (alpha(2:nz) * (tint(i,j-1,2:nz) - tint(i,j,2:nz)) + &
          beta(2:nz)  * (sint(i,j-1,2:nz) - sint(i,j,2:nz)))
  endif
  ! down (j+1)
  if (g%mask2dT(i,j+1) > 0.) then
    call calculate_density_derivs( &
         0.5 * (tint(i,j,2:nz) + tint(i,j+1,2:nz)), &
         0.5 * (sint(i,j,2:nz) + sint(i,j+1,2:nz)), &
         0.5 * (zint(i,j,2:nz) + zint(i,j+1,2:nz)) * (gv%H_to_RZ * gv%g_Earth), &
         alpha, beta, tv%eqn_of_state, (/2,nz/) )

    del2sigma(2:nz) = del2sigma(2:nz) + &
         (alpha(2:nz) * (tint(i,j+1,2:nz) - tint(i,j,2:nz)) + &
          beta(2:nz)  * (sint(i,j+1,2:nz) - sint(i,j,2:nz)))
  endif
  ! left (i-1)
  if (g%mask2dT(i-1,j) > 0.) then
    call calculate_density_derivs( &
         0.5 * (tint(i,j,2:nz) + tint(i-1,j,2:nz)), &
         0.5 * (sint(i,j,2:nz) + sint(i-1,j,2:nz)), &
         0.5 * (zint(i,j,2:nz) + zint(i-1,j,2:nz)) * (gv%H_to_RZ * gv%g_Earth), &
         alpha, beta, tv%eqn_of_state, (/2,nz/) )

    del2sigma(2:nz) = del2sigma(2:nz) + &
         (alpha(2:nz) * (tint(i-1,j,2:nz) - tint(i,j,2:nz)) + &
          beta(2:nz)  * (sint(i-1,j,2:nz) - sint(i,j,2:nz)))
  endif
  ! right (i+1)
  if (g%mask2dT(i+1,j) > 0.) then
    call calculate_density_derivs( &
         0.5 * (tint(i,j,2:nz) + tint(i+1,j,2:nz)), &
         0.5 * (sint(i,j,2:nz) + sint(i+1,j,2:nz)), &
         0.5 * (zint(i,j,2:nz) + zint(i+1,j,2:nz)) * (gv%H_to_RZ * gv%g_Earth), &
         alpha, beta, tv%eqn_of_state, (/2,nz/) )

    del2sigma(2:nz) = del2sigma(2:nz) + &
         (alpha(2:nz) * (tint(i+1,j,2:nz) - tint(i,j,2:nz)) + &
          beta(2:nz)  * (sint(i+1,j,2:nz) - sint(i,j,2:nz)))
  endif

  ! at this point, del2sigma contains the local neutral density curvature at
  ! h-points, on interfaces
  ! we need to divide by drho/dz to give an interfacial displacement
  !
  ! a positive curvature means we're too light relative to adjacent columns,
  ! so del2sigma needs to be positive too (push the interface deeper)
  call calculate_density_derivs(tint(i,j,:), sint(i,j,:), zint(i,j,:) * (gv%H_to_RZ * gv%g_Earth), &
       alpha, beta, tv%eqn_of_state, (/1,nz+1/) )
  do k = 2, nz
    ! TODO make lower bound here configurable
    dh_d2s(k) = del2sigma(k) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / &
         max(alpha(k) * (tv%T(i,j,k) - tv%T(i,j,k-1)) + &
             beta(k)  * (tv%S(i,j,k) - tv%S(i,j,k-1)), 1e-20*us%kg_m3_to_R)

    ! don't move the interface so far that it would tangle with another
    ! interface in the direction we're moving (or exceed a Nyquist limit
    ! that could cause oscillations of the interface)
    h_up = merge(h(i,j,k), h(i,j,k-1), dh_d2s(k) > 0.)
    dh_d2s(k) = 0.5 * cs%adaptAlpha * &
         sign(min(abs(del2sigma(k)), 0.5 * h_up), dh_d2s(k))

    ! update interface positions so we can diffuse them
    znext(k) = zint(i,j,k) + dh_d2s(k)
  enddo

  ! solve diffusivity equation to smooth grid
  ! upper diagonal coefficients: -kGrid(2:nz)
  ! lower diagonal coefficients: -kGrid(1:nz-1)
  ! diagonal coefficients:       1 + (kGrid(1:nz-1) + kGrid(2:nz))
  !
  ! first, calculate the diffusivities within layers
  do k = 1, nz
    ! calculate the dr bit of drdz
    drdz = 0.5 * (alpha(k) + alpha(k+1)) * (tint(i,j,k+1) - tint(i,j,k)) + &
           0.5 * (beta(k)  + beta(k+1))  * (sint(i,j,k+1) - sint(i,j,k))
    ! divide by dz from the new interface positions
    drdz = drdz / (znext(k) - znext(k+1) + gv%H_subroundoff)
    ! don't do weird stuff in unstably-stratified regions
    drdz = max(drdz, 0.)

    ! set vertical grid diffusivity
    kgrid(k) = (cs%adaptTimeRatio * nz**2 * depth) * &
         (cs%adaptZoomCoeff / (cs%adaptZoom + 0.5*(znext(k) + znext(k+1))) + &
         (cs%adaptBuoyCoeff * drdz / cs%adaptDrho0) + &
         max(1.0 - cs%adaptZoomCoeff - cs%adaptBuoyCoeff, 0.0) / depth)
  enddo

  ! initial denominator (first diagonal element)
  b1 = 1.0
  ! initial Q_1 = 1 - q_1 = 1 - 0/1
  d1 = 1.0
  ! work on all interior interfaces
  do k = 2, nz
    ! calculate numerator of Q_k
    b_denom_1 = 1. + d1 * kgrid(k-1)
    ! update denominator for k
    b1 = 1.0 / (b_denom_1 + kgrid(k))

    ! calculate q_k
    c1(k) = kgrid(k) * b1
    ! update Q_k = 1 - q_k
    d1 = b_denom_1 * b1

    ! update RHS
    znext(k) = b1 * (znext(k) + kgrid(k-1)*znext(k-1))
  enddo
  ! final substitution
  do k = nz, 2, -1
    znext(k) = znext(k) + c1(k)*znext(k+1)
  enddo

  if (cs%adaptDoMin) then
    nominal_z = 0.
    stretching = zint(i,j,nz+1) / depth

    do k = 2, nz+1
      nominal_z = nominal_z + cs%coordinateResolution(k-1) * stretching
      ! take the deeper of the calculated and nominal positions
      znext(k) = max(znext(k), nominal_z)
      ! interface can't go below topography
      znext(k) = min(znext(k), zint(i,j,nz+1))
    enddo
  endif

end_coord_adapt()

subroutine, public coord_adapt::end_coord_adapt

(

type(adapt_cs), pointer

CS

)

Clean up the coordinate control structure.

Parameters

cs	The control structure for this module

Definition at line 80 of file coord_adapt.F90.

  type(adapt_cs), pointer :: cs  !< The control structure for this module

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

init_coord_adapt()

subroutine, public coord_adapt::init_coord_adapt	(	type(adapt_cs), pointer	CS,
		integer, intent(in)	nk,
		real, dimension(:), intent(in)	coordinateResolution,
		real, intent(in)	m_to_H,
		real, intent(in)	kg_m3_to_R
	)

Initialise an adapt_CS with parameters.

Parameters

	cs	Unassociated pointer to hold the control structure
[in]	nk	Number of layers in the grid
[in]	coordinateresolution	Nominal near-surface resolution [m] or other units specified with m_to_H
[in]	m_to_h	A conversion factor from m to the units of thicknesses
[in]	kg_m3_to_r	A conversion factor from kg m-3 to the units of density

Definition at line 54 of file coord_adapt.F90.

  type(adapt_cs),     pointer    :: cs !< Unassociated pointer to hold the control structure
  integer,            intent(in) :: nk !< Number of layers in the grid
  real, dimension(:), intent(in) :: coordinateresolution !< Nominal near-surface resolution [m] or
                                       !! other units specified with m_to_H
  real,               intent(in) :: m_to_h !< A conversion factor from m to the units of thicknesses
  real,               intent(in) :: kg_m3_to_r !< A conversion factor from kg m-3 to the units of density

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

  cs%nk = nk
  cs%coordinateResolution(:) = coordinateresolution(:)

  ! Set real parameter default values
  cs%adaptTimeRatio = 1e-1 ! Nondim.
  cs%adaptAlpha     = 1.0  ! Nondim.
  cs%adaptZoom      = 200.0 * m_to_h    ! [H ~> m or kg m-2]
  cs%adaptZoomCoeff = 0.0  ! Nondim.
  cs%adaptBuoyCoeff = 0.0  ! Nondim.
  cs%adaptDrho0     = 0.5 * kg_m3_to_r  ! [R ~> kg m-3]

set_adapt_params()

subroutine, public coord_adapt::set_adapt_params	(	type(adapt_cs), pointer	CS,
		real, intent(in), optional	adaptTimeRatio,
		real, intent(in), optional	adaptAlpha,
		real, intent(in), optional	adaptZoom,
		real, intent(in), optional	adaptZoomCoeff,
		real, intent(in), optional	adaptBuoyCoeff,
		real, intent(in), optional	adaptDrho0,
		logical, intent(in), optional	adaptDoMin
	)

This subtroutine can be used to set the parameters for coord_adapt module.

Parameters

	cs	The control structure for this module
[in]	adapttimeratio	Ratio of optimisation and diffusion timescales
[in]	adaptalpha	Nondimensional coefficient determining how much optimisation to apply
[in]	adaptzoom	Near-surface zooming depth [H ~> m or kg m-2]
[in]	adaptzoomcoeff	Near-surface zooming coefficient
[in]	adaptbuoycoeff	Stratification-dependent diffusion coefficient
[in]	adaptdrho0	Reference density difference for stratification-dependent diffusion [R ~> kg m-3]
[in]	adaptdomin	If true, form a HYCOM1-like mixed layer by preventing interfaces from becoming shallower than the depths set by coordinateResolution

Definition at line 91 of file coord_adapt.F90.

  type(adapt_cs),    pointer    :: cs  !< The control structure for this module
  real,    optional, intent(in) :: adapttimeratio !< Ratio of optimisation and diffusion timescales
  real,    optional, intent(in) :: adaptalpha     !< Nondimensional coefficient determining
                                                  !! how much optimisation to apply
  real,    optional, intent(in) :: adaptzoom      !< Near-surface zooming depth [H ~> m or kg m-2]
  real,    optional, intent(in) :: adaptzoomcoeff !< Near-surface zooming coefficient
  real,    optional, intent(in) :: adaptbuoycoeff !< Stratification-dependent diffusion coefficient
  real,    optional, intent(in) :: adaptdrho0  !< Reference density difference for
                                               !! stratification-dependent diffusion [R ~> kg m-3]
  logical, optional, intent(in) :: adaptdomin  !< If true, form a HYCOM1-like mixed layer by
                                               !! preventing interfaces from becoming shallower than
                                               !! the depths set by coordinateResolution

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

  if (present(adapttimeratio)) cs%adaptTimeRatio = adapttimeratio
  if (present(adaptalpha)) cs%adaptAlpha = adaptalpha
  if (present(adaptzoom)) cs%adaptZoom = adaptzoom
  if (present(adaptzoomcoeff)) cs%adaptZoomCoeff = adaptzoomcoeff
  if (present(adaptbuoycoeff)) cs%adaptBuoyCoeff = adaptbuoycoeff
  if (present(adaptdrho0)) cs%adaptDrho0 = adaptdrho0
  if (present(adaptdomin)) cs%adaptDoMin = adaptdomin