Detailed Description

Regrid columns for the adaptive coordinate.

Data Types
type	adapt_cs

Functions/Subroutines
subroutine, public	init_coord_adapt (CS, nk, coordinateResolution)
	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)

subroutine, public	build_adapt_column (CS, G, GV, 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(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]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in]	tv	A structure pointing to various thermodynamic variables
[in]	h	Layer thicknesses, in H (usually m or kg m-2)

Definition at line 92 of file coord_adapt.F90.

References mom_eos::calculate_density_derivs().

   type(adapt_cs),                              intent(in)    :: cs
   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(thermo_var_ptrs),                       intent(in)    :: tv   !< A structure pointing to various thermodynamic variables
   integer,                                     intent(in)    :: i, j
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in)    :: zint, tint, sint
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in)    :: h    !< Layer thicknesses, in H (usually 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, drdz, nominal_z, stretching
   real, dimension(SZK_(GV)+1) :: alpha, beta, del2sigma ! drho/dT and drho/dS
   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%m_to_H
 
   ! initialize del2sigma to zero
   del2sigma(:) = 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_Pa, &
          alpha, beta, 2, nz - 1, tv%eqn_of_state)
 
     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_Pa, &
          alpha, beta, 2, nz - 1, tv%eqn_of_state)
 
     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_Pa, &
          alpha, beta, 2, nz - 1, tv%eqn_of_state)
 
     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_Pa, &
          alpha, beta, 2, nz - 1, tv%eqn_of_state)
 
     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_Pa, &
        alpha, beta, 1, nz + 1, tv%eqn_of_state)
   do k = 2, nz
     ! TODO make lower bound here configurable
     del2sigma(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)
 
     ! 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), del2sigma(k) > 0.)
     del2sigma(k) = 0.5 * cs%adaptAlpha * &
          sign(min(abs(del2sigma(k)), 0.5 * h_up), del2sigma(k))
 
     ! update interface positions so we can diffuse them
     znext(k) = zint(i,j,k) + del2sigma(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 * gv%m_to_H + 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

Here is the call graph for this function:

◆ end_coord_adapt()

subroutine, public coord_adapt::end_coord_adapt ( type(adapt_cs), pointer CS )

Clean up the coordinate control structure.

Definition at line 65 of file coord_adapt.F90.

Referenced by mom_regridding::end_regridding().

   type(adapt_cs), pointer :: cs
 
   ! nothing to do
   if (.not. associated(cs)) return
   deallocate(cs%coordinateResolution)
   deallocate(cs)

Here is the caller graph for this function:

◆ init_coord_adapt()

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

Initialise an adapt_CS with parameters.

Parameters

cs	Unassociated pointer to hold the control structure

Definition at line 51 of file coord_adapt.F90.

References mom_error_handler::mom_error().

   type(adapt_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_adapt: CS already associated")
   allocate(cs)
   allocate(cs%coordinateResolution(nk))
 
   cs%nk = nk
   cs%coordinateResolution(:) = coordinateresolution(:)

Here is the call graph for this function:

◆ 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
	)

Definition at line 75 of file coord_adapt.F90.

References mom_error_handler::mom_error().

   type(adapt_cs),    pointer    :: cs
   real, optional,    intent(in) :: adapttimeratio, adaptalpha, adaptzoom, adaptzoomcoeff
   real, optional,    intent(in) :: adaptbuoycoeff, adaptdrho0
   logical, optional, intent(in) :: adaptdomin
 
   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

Here is the call graph for this function:

Detailed Description

Data Types

Functions/Subroutines

Function/Subroutine Documentation

◆ build_adapt_column()

◆ end_coord_adapt()

◆ init_coord_adapt()

◆ set_adapt_params()