Detailed Description

Regrid columns for the continuous isopycnal (rho) coordinate.

Data Types
type	rho_cs
	Control structure containing required parameters for the rho coordinate. More...

Functions/Subroutines
subroutine, public	init_coord_rho (CS, nk, ref_pressure, target_density, interp_CS)
	Initialise a rho_CS with pointers to parameters. More...

subroutine, public	end_coord_rho (CS)
	This subroutine deallocates memory in the control structure for the coord_rho module. More...

subroutine, public	set_rho_params (CS, min_thickness, integrate_downward_for_e, interp_CS)
	This subroutine can be used to set the parameters for the coord_rho module. More...

subroutine, public	build_rho_column (CS, nz, depth, h, T, S, eqn_of_state, z_interface, h_neglect, h_neglect_edge)
	Build a rho coordinate column. More...

subroutine	build_rho_column_iteratively (CS, remapCS, nz, depth, h, T, S, eqn_of_state, zInterface, h_neglect, h_neglect_edge)
	Iteratively build a rho coordinate column. More...

subroutine	copy_finite_thicknesses (nk, h_in, threshold, nout, h_out, mapping)
	Copy column thicknesses with vanished layers removed. More...

subroutine, public	old_inflate_layers_1d (min_thickness, nk, h)
	Inflate vanished layers to finite (nonzero) width. More...

Variables
integer, parameter	nb_regridding_iterations = 1
	Maximum number of regridding iterations.

real, parameter	deviation_tolerance = 1e-10
	Deviation tolerance between succesive grids in regridding iterations.

Function/Subroutine Documentation

◆ build_rho_column()

subroutine, public coord_rho::build_rho_column	(	type(rho_cs), intent(in)	CS,
		integer, intent(in)	nz,
		real, intent(in)	depth,
		real, dimension(nz), intent(in)	h,
		real, dimension(nz), intent(in)	T,
		real, dimension(nz), intent(in)	S,
		type(eos_type), pointer	eqn_of_state,
		real, dimension(cs%nk+1), intent(inout)	z_interface,
		real, intent(in), optional	h_neglect,
		real, intent(in), optional	h_neglect_edge
	)

Build a rho coordinate column.

Density profiles are calculated on the source grid.

Positions of target densities (for interfaces) are found by interpolation.

Parameters

[in]	cs	coord_rho control structure
[in]	nz	Number of levels on source grid (i.e. length of h, T, S)
[in]	depth	Depth of ocean bottom (positive in m)
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in]	t	T for source column
[in]	s	S for source column
	eqn_of_state	Equation of state structure
[in,out]	z_interface	Absolute positions of interfaces
[in]	h_neglect	A negligibly small width for the purpose of cell reconstructions in the same units as h0.
[in]	h_neglect_edge	A negligibly small width for the purpose of edge value calculations in the same units as h0.

Definition at line 96 of file coord_rho.F90.

   type(rho_CS),        intent(in)    :: CS !< coord_rho control structure
   integer,             intent(in)    :: nz !< Number of levels on source grid (i.e. length of  h, T, S)
   real,                intent(in)    :: depth !< Depth of ocean bottom (positive in m)
   real, dimension(nz), intent(in)    :: h  !< Layer thicknesses [H ~> m or kg m-2]
   real, dimension(nz), intent(in)    :: T  !< T for source column
   real, dimension(nz), intent(in)    :: S  !< S for source column
   type(EOS_type),      pointer       :: eqn_of_state !< Equation of state structure
   real, dimension(CS%nk+1), &
                        intent(inout) :: z_interface !< Absolute positions of interfaces
   real,      optional, intent(in)    :: h_neglect !< A negligibly small width for the
                                              !! purpose of cell reconstructions
                                              !! in the same units as h0.
   real,      optional, intent(in)    :: h_neglect_edge !< A negligibly small width
                                              !! for the purpose of edge value calculations
                                              !! in the same units as h0.
   ! Local variables
   integer :: k, count_nonzero_layers
   integer, dimension(nz) :: mapping
   real, dimension(nz) :: p, densities, h_nv
   real, dimension(nz+1) :: xTmp
   real, dimension(CS%nk) :: h_new ! New thicknesses
   real, dimension(CS%nk+1) :: x1
  
   ! Construct source column with vanished layers removed (stored in h_nv)
   call copy_finite_thicknesses(nz, h, cs%min_thickness, count_nonzero_layers, h_nv, mapping)
  
   if (count_nonzero_layers > 1) then
     xtmp(1) = 0.0
     do k = 1,count_nonzero_layers
       xtmp(k+1) = xtmp(k) + h_nv(k)
     enddo
  
     ! Compute densities on source column
     p(:) = cs%ref_pressure
     call calculate_density(t, s, p, densities, 1, nz, eqn_of_state)
     do k = 1,count_nonzero_layers
       densities(k) = densities(mapping(k))
     enddo
  
     ! Based on source column density profile, interpolate to generate a new grid
     call build_and_interpolate_grid(cs%interp_CS, densities, count_nonzero_layers, &
                                     h_nv, xtmp, cs%target_density, cs%nk, h_new, &
                                     x1, h_neglect, h_neglect_edge)
  
     ! Inflate vanished layers
     call old_inflate_layers_1d(cs%min_thickness, cs%nk, h_new)
  
     ! Comment: The following adjustment of h_new, and re-calculation of h_new via x1 needs to be removed
     x1(1) = 0.0 ; do k = 1,cs%nk ; x1(k+1) = x1(k) + h_new(k) ; enddo
     do k = 1,cs%nk
       h_new(k) = x1(k+1) - x1(k)
     enddo
  
   else ! count_nonzero_layers <= 1
     if (nz == cs%nk) then
       h_new(:) = h(:) ! This keeps old behavior
     else
       h_new(:) = 0.
       h_new(1) = h(1)
     endif
   endif
  
   ! Return interface positions
   if (cs%integrate_downward_for_e) then
     ! Remapping is defined integrating from zero
     z_interface(1) = 0.
     do k = 1,cs%nk
       z_interface(k+1) = z_interface(k) - h_new(k)
     enddo
   else
     ! The rest of the model defines grids integrating up from the bottom
     z_interface(cs%nk+1) = -depth
     do k = cs%nk,1,-1
       z_interface(k) = z_interface(k+1) + h_new(k)
     enddo
   endif
  

◆ build_rho_column_iteratively()

subroutine coord_rho::build_rho_column_iteratively	(	type(rho_cs), intent(in)	CS,
		type(remapping_cs), intent(in)	remapCS,
		integer, intent(in)	nz,
		real, intent(in)	depth,
		real, dimension(nz), intent(in)	h,
		real, dimension(nz), intent(in)	T,
		real, dimension(nz), intent(in)	S,
		type(eos_type), pointer	eqn_of_state,
		real, dimension(nz+1), intent(inout)	zInterface,
		real, intent(in), optional	h_neglect,
		real, intent(in), optional	h_neglect_edge
	)

private

Iteratively build a rho coordinate column.

The algorithm operates as follows within each column:

Given T & S within each layer, the layer densities are computed.
Based on these layer densities, a global density profile is reconstructed (this profile is monotonically increasing and may be discontinuous)
The new grid interfaces are determined based on the target interface densities.
T & S are remapped onto the new grid.

Return to step 1 until convergence or until the maximum number of iterations is reached, whichever comes first.

Parameters

[in]	cs	Regridding control structure
[in]	remapcs	Remapping parameters and options
[in]	nz	Number of levels
[in]	depth	Depth of ocean bottom [Z ~> m]
[in]	h	Layer thicknesses in Z coordinates [Z ~> m]
[in]	t	T for column [degC]
[in]	s	S for column [ppt]
	eqn_of_state	Equation of state structure
[in,out]	zinterface	Absolute positions of interfaces
[in]	h_neglect	A negligibly small width for the purpose of cell reconstructions in the same units as h [Z ~> m]
[in]	h_neglect_edge	A negligibly small width for the purpose of edge value calculations in the same units as h [Z ~> m]

Definition at line 189 of file coord_rho.F90.

   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 [Z ~> m]
   real, dimension(nz),   intent(in)    :: h  !< Layer thicknesses in Z coordinates [Z ~> m]
   real, dimension(nz),   intent(in)    :: T  !< T for column [degC]
   real, dimension(nz),   intent(in)    :: S  !< S for column [ppt]
   type(EOS_type),        pointer       :: eqn_of_state !< Equation of state structure
   real, dimension(nz+1), intent(inout) :: zInterface !< Absolute positions of interfaces
   real,        optional, intent(in)    :: h_neglect !< A negligibly small width for the
                                              !! purpose of cell reconstructions
                                              !! in the same units as h [Z ~> m]
   real,        optional, intent(in)    :: h_neglect_edge !< A negligibly small width
                                              !! for the purpose of edge value calculations
                                              !! in the same units as h [Z ~> m]
   ! Local variables
   integer   :: k, m
   integer   :: count_nonzero_layers
   real      :: deviation            ! When iterating to determine the final
                                     ! grid, this is the deviation between two
                                     ! successive grids.
   real      :: threshold
   real, dimension(nz) :: p, densities, T_tmp, S_tmp, Tmp
   integer, dimension(nz) :: mapping
   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 ) )
  
     ! Construct column with vanished layers removed
     call copy_finite_thicknesses(nz, h0, threshold, count_nonzero_layers, htmp, mapping)
     if ( count_nonzero_layers <= 1 ) then
       h1(:) = h0(:)
       exit  ! stop iterations here
     endif
  
     xtmp(1) = 0.0
     do k = 1,count_nonzero_layers
       xtmp(k+1) = xtmp(k) + htmp(k)
     enddo
  
     ! 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))
     enddo
  
     ! 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, h_neglect, h_neglect_edge)
  
     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) ; enddo
  
     ! Remap T and S from previous grid to new grid
     do k = 1,nz
       h1(k) = x1(k+1) - x1(k)
     enddo
  
     call remapping_core_h(remapcs, nz, h0, s, nz, h1, tmp, h_neglect, h_neglect_edge)
     s_tmp(:) = tmp(:)
  
     call remapping_core_h(remapcs, nz, h0, t, nz, h1, tmp, h_neglect, h_neglect_edge)
     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
     enddo
     deviation = sqrt( deviation / (nz-1) )
  
     m = m + 1
  
     ! Copy final grid onto start grid for next iteration
     h0(:) = h1(:)
  
   enddo ! 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 accommodate 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 accommodate inflating layers
       ! without disturbing the interface above
     enddo
   endif
  

◆ copy_finite_thicknesses()

subroutine coord_rho::copy_finite_thicknesses	(	integer, intent(in)	nk,
		real, dimension(nk), intent(in)	h_in,
		real, intent(in)	threshold,
		integer, intent(out)	nout,
		real, dimension(nk), intent(out)	h_out,
		integer, dimension(nk), intent(out)	mapping
	)

private

Copy column thicknesses with vanished layers removed.

Parameters

[in]	nk	Number of layer for h_in, T_in, S_in
[in]	h_in	Thickness of input column
[in]	threshold	Thickness threshold defining vanished layers
[out]	nout	Number of non-vanished layers
[out]	h_out	Thickness of output column
[out]	mapping	Index of k-out corresponding to k-in

Definition at line 306 of file coord_rho.F90.

   integer,                intent(in)  :: nk !< Number of layer for h_in, T_in, S_in
   real, dimension(nk),    intent(in)  :: h_in !< Thickness of input column
   real,                   intent(in)  :: threshold !< Thickness threshold defining vanished layers
   integer,                intent(out) :: nout !< Number of non-vanished layers
   real, dimension(nk),    intent(out) :: h_out !< Thickness of output column
   integer, dimension(nk), intent(out) :: mapping !< Index of k-out corresponding to k-in
   ! Local variables
   integer :: k, k_thickest
   real :: thickness_in_vanished, thickest_h_out
  
   ! Build up new grid
   nout = 0
   thickness_in_vanished = 0.0
   thickest_h_out = h_in(1)
   k_thickest = 1
   do k = 1, nk
     mapping(k) = nout ! Note k>=nout always
     h_out(k) = 0.  ! Make sure h_out is set everywhere
     if (h_in(k) > threshold) then
       ! For non-vanished layers
       nout = nout + 1
       mapping(nout) = k
       h_out(nout) = h_in(k)
       if (h_out(nout) > thickest_h_out) then
         thickest_h_out = h_out(nout)
         k_thickest = nout
       endif
     else
       ! Add up mass in vanished layers
       thickness_in_vanished = thickness_in_vanished + h_in(k)
     endif
   enddo
  
   ! No finite layers
   if (nout <= 1) return
  
   ! Adjust for any lost volume in vanished layers
   h_out(k_thickest) = h_out(k_thickest) + thickness_in_vanished
  

◆ end_coord_rho()

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

This subroutine deallocates memory in the control structure for the coord_rho module.

Parameters

cs	Coordinate control structure

Definition at line 65 of file coord_rho.F90.

   type(rho_CS), pointer :: CS !< Coordinate control structure
  
   ! nothing to do
   if (.not. associated(cs)) return
   deallocate(cs%target_density)
   deallocate(cs)

◆ init_coord_rho()

subroutine, public coord_rho::init_coord_rho	(	type(rho_cs), pointer	CS,
		integer, intent(in)	nk,
		real, intent(in)	ref_pressure,
		real, dimension(:), intent(in)	target_density,
		type(interp_cs_type), intent(in)	interp_CS
	)

Initialise a rho_CS with pointers to parameters.

Parameters

	cs	Unassociated pointer to hold the control structure
[in]	nk	Number of layers in the grid
[in]	ref_pressure	Coordinate reference pressure [Pa]
[in]	target_density	Nominal density of interfaces [kg m-3]
[in]	interp_cs	Controls for interpolation

Definition at line 47 of file coord_rho.F90.

   type(rho_CS),         pointer    :: CS !< Unassociated pointer to hold the control structure
   integer,              intent(in) :: nk !< Number of layers in the grid
   real,                 intent(in) :: ref_pressure !< Coordinate reference pressure [Pa]
   real, dimension(:),   intent(in) :: target_density !< Nominal density of interfaces [kg m-3]
   type(interp_CS_type), intent(in) :: interp_CS !< Controls for interpolation
  
   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

◆ old_inflate_layers_1d()

subroutine, public coord_rho::old_inflate_layers_1d	(	real, intent(in)	min_thickness,
		integer, intent(in)	nk,
		real, dimension(:), intent(inout)	h
	)

Inflate vanished layers to finite (nonzero) width.

Parameters

[in]	min_thickness	Minimum allowed thickness [H ~> m or kg m-2]
[in]	nk	Number of layers in the grid
[in,out]	h	Layer thicknesses [H ~> m or kg m-2]

Definition at line 350 of file coord_rho.F90.

  
   ! Argument
   real,               intent(in)    :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2]
   integer,            intent(in)    :: nk  !< Number of layers in the grid
   real, dimension(:), intent(inout) :: h   !< Layer thicknesses [H ~> m or kg m-2]
  
   ! 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,nk
     if ( h(k) > min_thickness ) then
       count_nonzero_layers = count_nonzero_layers + 1
     endif
   enddo
  
   ! If all layer thicknesses are greater than the threshold, exit routine
   if ( count_nonzero_layers == nk ) return
  
   ! If all thicknesses are zero, inflate them all and exit
   if ( count_nonzero_layers == 0 ) then
     do k = 1,nk
       h(k) = min_thickness
     enddo
     return
   endif
  
   ! Inflate zero layers
   correction = 0.0
   do k = 1,nk
     if ( h(k) <= min_thickness ) then
       delta = min_thickness - h(k)
       correction = correction + delta
       h(k) = h(k) + delta
     endif
   enddo
  
   ! Modify thicknesses of nonzero layers to ensure volume conservation
   maxthickness = h(1)
   k_found = 1
   do k = 1,nk
     if ( h(k) > maxthickness ) then
       maxthickness = h(k)
       k_found = k
     endif
   enddo
  
   h(k_found) = h(k_found) - correction
  

◆ set_rho_params()

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

This subroutine can be used to set the parameters for the coord_rho module.

Parameters

	cs	Coordinate control structure
[in]	min_thickness	Minimum allowed thickness [H ~> m or kg m-2]
[in]	integrate_downward_for_e	If true, integrate for interface positions from the top downward. If false, integrate from the bottom upward, as does the rest of the model.
[in]	interp_cs	Controls for interpolation

Definition at line 75 of file coord_rho.F90.

   type(rho_CS),      pointer    :: CS !< Coordinate control structure
   real,    optional, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2]
   logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface
                                       !! positions from the top downward.  If false, integrate
                                       !! from the bottom upward, as does the rest of the model.
  
   type(interp_CS_type), optional, intent(in) :: interp_CS !< Controls for interpolation
  
   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

Detailed Description

Data Types

Functions/Subroutines

Variables

Function/Subroutine Documentation

◆ build_rho_column()

◆ build_rho_column_iteratively()

◆ copy_finite_thicknesses()

◆ end_coord_rho()

◆ init_coord_rho()

◆ old_inflate_layers_1d()

◆ set_rho_params()