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)

subroutine, public	set_rho_params (CS, min_thickness, integrate_downward_for_e, interp_CS)

subroutine, public	build_rho_column (CS, remapCS, nz, depth, h, T, S, eqn_of_state, zInterface)

subroutine, public	old_inflate_layers_1d (minThickness, N, h)

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

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

Function/Subroutine Documentation

◆ build_rho_column()

subroutine, public coord_rho::build_rho_column	(	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
	)

Parameters

[in]	cs	Regridding control structure
[in]	remapcs	Remapping parameters and options
[in]	nz	Number of levels
[in]	depth	Depth of ocean bottom (positive in m)
[in]	h	Layer thicknesses, in m
[in]	s	T and S for column
	eqn_of_state	Equation of state structure
[in,out]	zinterface	Absolute positions of interfaces

Definition at line 86 of file coord_rho.F90.

References regrid_interp::build_and_interpolate_grid(), deviation_tolerance, nb_regridding_iterations, old_inflate_layers_1d(), and mom_remapping::remapping_core_h().

Referenced by mom_diag_remap::diag_remap_update().

   !< 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
 

Here is the call graph for this function:

Here is the caller graph for this function:

◆ end_coord_rho()

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

Definition at line 64 of file coord_rho.F90.

   type(rho_cs), pointer :: cs
 
   ! 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

Definition at line 47 of file coord_rho.F90.

References mom_error_handler::mom_error().

   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

Here is the call graph for this function:

◆ old_inflate_layers_1d()

subroutine, public coord_rho::old_inflate_layers_1d	(	real, intent(in)	minThickness,
		integer, intent(in)	N,
		real, dimension(:), intent(inout)	h
	)

Definition at line 248 of file coord_rho.F90.

Referenced by build_rho_column().

 
   ! 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
 

Here is the caller graph for this function:

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

Definition at line 73 of file coord_rho.F90.

References mom_error_handler::mom_error().

   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

Here is the call graph for this function:

Variable Documentation

◆ deviation_tolerance

real, parameter coord_rho::deviation_tolerance = 1e-10

private

Deviation tolerance between succesive grids in regridding iterations.

Definition at line 38 of file coord_rho.F90.

Referenced by build_rho_column().

38 real, parameter :: deviation_tolerance = 1e-10

◆ nb_regridding_iterations

integer, parameter coord_rho::nb_regridding_iterations = 1

Maximum number of regridding iterations.

Definition at line 36 of file coord_rho.F90.

Referenced by build_rho_column().

36 integer, parameter :: nb_regridding_iterations = 1

Detailed Description

Data Types

Functions/Subroutines

Variables

Function/Subroutine Documentation

◆ build_rho_column()

◆ end_coord_rho()

◆ init_coord_rho()

◆ old_inflate_layers_1d()

◆ set_rho_params()

Variable Documentation

◆ deviation_tolerance

◆ nb_regridding_iterations