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, dev_tol)
	Iteratively build a rho coordinate column. More...

subroutine	copy_finite_thicknesses (nk, h_in, thresh, 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...

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 downward) [H ~> m or kg m-2]
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in]	t	Temperature for source column [degC]
[in]	s	Salinity for source column [ppt]
	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 [H ~> m or kg m-2]
[in]	h_neglect_edge	A negligibly small width for the purpose of edge value calculations [H ~> m or kg m-2]

Definition at line 92 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 downward) [H ~> m or kg m-2]
   real, dimension(nz), intent(in)    :: h  !< Layer thicknesses [H ~> m or kg m-2]
   real, dimension(nz), intent(in)    :: T  !< Temperature for source column [degC]
   real, dimension(nz), intent(in)    :: S  !< Salinity for source column [ppt]
   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 [H ~> m or kg m-2]
   real,      optional, intent(in)    :: h_neglect_edge !< A negligibly small width for the purpose
                                              !! of edge value calculations [H ~> m or kg m-2]
 
   ! Local variables
   integer :: k, count_nonzero_layers
   integer, dimension(nz) :: mapping
   real, dimension(nz) :: pres     ! Pressures used to calculate density [R L2 T-2 ~> Pa]
   real, dimension(nz) :: h_nv     ! Thicknesses of non-vanishing layers [H ~> m or kg m-2]
   real, dimension(nz) :: densities ! Layer density [R ~> kg m-3]
   real, dimension(nz+1) :: xTmp   ! Temporary positions [H ~> m or kg m-2]
   real, dimension(CS%nk) :: h_new ! New thicknesses [H ~> m or kg m-2]
   real, dimension(CS%nk+1) :: x1  ! Interface heights [H ~> m or kg m-2]
 
   ! 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
     pres(:) = cs%ref_pressure
     call calculate_density(t, s, pres, densities, 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,
		real, intent(in), optional	dev_tol
	)

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]
[in]	dev_tol	The tolerance for the deviation between successive grids for determining when the iterative solver has converged [Z ~> m]

Definition at line 188 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]
   real,        optional, intent(in)    :: dev_tol !< The tolerance for the deviation between
                                              !! successive grids for determining when the
                                              !! iterative solver has converged [Z ~> m]
 
   ! Local variables
   real, dimension(nz+1) :: x0, x1, xTmp ! Temporary interface heights [Z ~> m]
   real, dimension(nz) :: pres       ! The pressure used in the equation of state [R L2 T-2 ~> Pa].
   real, dimension(nz) :: densities  ! Layer densities [R ~> kg m-3]
   real, dimension(nz) :: T_tmp, S_tmp ! A temporary profile of temperature [degC] and salinity [ppt].
   real, dimension(nz) :: Tmp        ! A temporary variable holding a remapped variable.
   real, dimension(nz) :: h0, h1, hTmp ! Temporary thicknesses [Z ~> m]
   real :: deviation            ! When iterating to determine the final grid, this is the
                                ! deviation between two successive grids [Z ~> m].
   real :: deviation_tol        ! Deviation tolerance between succesive grids in
                                ! regridding iterations [Z ~> m]
   real :: threshold            ! The minimum thickness for a layer to be considered to exist [Z ~> m]
   integer, dimension(nz) :: mapping ! The indices of the massive layers in the initial column.
   integer :: k, m, count_nonzero_layers
 
   !  Maximum number of regridding iterations
   integer, parameter :: NB_REGRIDDING_ITERATIONS = 1
 
   threshold = cs%min_thickness
   pres(:) = cs%ref_pressure
   t_tmp(:) = t(:)
   s_tmp(:) = s(:)
   h0(:) = h(:)
 
   ! Start iterations to build grid
   m = 1
   deviation_tol = 1.0e-15*depth ; if (present(dev_tol)) deviation_tol = dev_tol
 
   do m=1,nb_regridding_iterations
 
     ! 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, pres, densities, 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) )
 
     if ( deviation <= deviation_tol ) exit
 
     ! 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)	thresh,
		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 [H ~> m or kg m-2] or [Z ~> m]
[in]	thresh	Thickness threshold defining vanished layers [H ~> m or kg m-2] or [Z ~> m]
[out]	nout	Number of non-vanished layers
[out]	h_out	Thickness of output column [H ~> m or kg m-2] or [Z ~> m]
[out]	mapping	Index of k-out corresponding to k-in

Definition at line 313 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 [H ~> m or kg m-2] or [Z ~> m]
   real,                   intent(in)  :: thresh  !< Thickness threshold defining vanished
                                                  !! layers [H ~> m or kg m-2] or [Z ~> m]
   integer,                intent(out) :: nout    !< Number of non-vanished layers
   real, dimension(nk),    intent(out) :: h_out   !< Thickness of output column [H ~> m or kg m-2] or [Z ~> m]
   integer, dimension(nk), intent(out) :: mapping !< Index of k-out corresponding to k-in
   ! Local variables
   integer :: k, k_thickest
   real :: thickness_in_vanished ! Summed thicknesses in discarded layers [H ~> m or kg m-2] or [Z ~> m]
   real :: thickest_h_out        ! Thickness of the thickest layer [H ~> m or kg m-2] or [Z ~> m]
 
   ! 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) > thresh) 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 61 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 [R L2 T-2 ~> Pa]
[in]	target_density	Nominal density of interfaces [R ~> kg m-3]
[in]	interp_cs	Controls for interpolation

Definition at line 42 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 [R L2 T-2 ~> Pa]
   real, dimension(:),   intent(in) :: target_density !< Nominal density of interfaces [R ~> 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 359 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 71 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

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