mom_regridding Module Reference

Detailed Description

Generates vertical grids as part of the ALE algorithm.

A vertical grid is defined solely by the cell thicknesses, \(h\). Most calculations in this module start with the coordinate at the bottom of the column set to -depth, and use a increasing value of coordinate with decreasing k. This is consistent with the rest of MOM6 that uses position, \(z\) which is a negative quantity for most of the ocean.

A change in grid is define through a change in position of the interfaces:

\[ z^n_{k+1/2} = z^{n-1}_{k+1/2} + \Delta z_{k+1/2} \]

with the positive upward coordinate convention

\[ z_{k-1/2} = z_{k+1/2} + h_k \]

so that

\[ h^n_k = h^{n-1}_k + ( \Delta z_{k-1/2} - \Delta z_{k+1/2} ) \]

Original date of creation: 2008.06.09 by L. White

Data Types
type	regridding_cs
	Regridding control structure. More...

Functions/Subroutines
subroutine, public	initialize_regridding (CS, GV, US, max_depth, param_file, mdl, coord_mode, param_prefix, param_suffix)
	Initialization and configures a regridding control structure based on customizable run-time parameters. More...
subroutine	check_grid_def (filename, varname, expected_units, msg, ierr)
	Do some basic checks on the vertical grid definition file, variable. More...
subroutine, public	end_regridding (CS)
	Deallocation of regridding memory. More...
subroutine, public	regridding_main (remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h, conv_adjust)
	Dispatching regridding routine for orchestrating regridding & remapping. More...
subroutine	calc_h_new_by_dz (CS, G, GV, h, dzInterface, h_new)
	Calculates h_new from h + delta_k dzInterface. More...
subroutine, public	check_remapping_grid (G, GV, h, dzInterface, msg)
	Check that the total thickness of two grids match. More...
subroutine, public	check_grid_column (nk, depth, h, dzInterface, msg)
	Check that the total thickness of new and old grids are consistent. More...
subroutine	filtered_grid_motion (CS, nk, z_old, z_new, dz_g)
	Returns the change in interface position motion after filtering and assuming the top and bottom interfaces do not move. The filtering is a function of depth, and is applied as the integrated average filtering over the trajectory of the interface. By design, this code can not give tangled interfaces provided that z_old and z_new are not already tangled. More...
subroutine	build_zstar_grid (CS, G, GV, h, dzInterface, frac_shelf_h)
	Builds a z*-ccordinate grid with partial steps (Adcroft and Campin, 2004). z* is defined as z* = (z-eta)/(H+eta)*H s.t. z*=0 when z=eta and z*=-H when z=-H . More...
subroutine	build_sigma_grid (CS, G, GV, h, dzInterface)
	This routine builds a grid based on terrain-following coordinates. More...
subroutine	build_rho_grid (G, GV, US, h, tv, dzInterface, remapCS, CS)
	This routine builds a new grid based on a given set of target interface densities. More...
subroutine	build_grid_hycom1 (G, GV, US, h, tv, h_new, dzInterface, CS)
	Builds a simple HyCOM-like grid with the deepest location of potential density interpolated from the column profile and a clipping of depth for each interface to a fixed z* or p* grid. This should probably be (optionally?) changed to find the nearest location of the target density. More...
subroutine	build_grid_adaptive (G, GV, US, h, tv, dzInterface, remapCS, CS)
	This subroutine builds an adaptive grid that follows density surfaces where possible, subject to constraints on the smoothness of interface heights. More...
subroutine	build_grid_slight (G, GV, US, h, tv, dzInterface, CS)
	Builds a grid that tracks density interfaces for water that is denser than the surface density plus an increment of some number of layers, and uses all lighter layers uniformly above this location. Note that this amounts to interpolating to find the depth of an arbitrary (non-integer) interface index which should make the results vary smoothly in space to the extent that the surface density and interior stratification vary smoothly in space. Over shallow topography, this will tend to give a uniform sigma-like coordinate. For sufficiently shallow water, a minimum grid spacing is used to avoid certain instabilities. More...
subroutine	adjust_interface_motion (CS, nk, h_old, dz_int)
	Adjust dz_Interface to ensure non-negative future thicknesses. More...
subroutine	build_grid_arbitrary (G, GV, h, dzInterface, h_new, CS)
subroutine, public	inflate_vanished_layers_old (CS, G, GV, h)
subroutine	convective_adjustment (G, GV, h, tv)
	Achieve convective adjustment by swapping layers. More...
real function, dimension(nk), public	uniformresolution (nk, coordMode, maxDepth, rhoLight, rhoHeavy)
	Return a uniform resolution vector in the units of the coordinate. More...
subroutine	initcoord (CS, GV, US, coord_mode)
	Initialize the coordinate resolutions by calling the appropriate initialization routine for the specified coordinate mode. More...
subroutine, public	setcoordinateresolution (dz, CS, scale)
	Set the fixed resolution data. More...
subroutine, public	set_target_densities_from_gv (GV, US, CS)
	Set target densities based on the old Rlay variable. More...
subroutine, public	set_target_densities (CS, rho_int)
	Set target densities based on vector of interface values. More...
subroutine, public	set_regrid_max_depths (CS, max_depths, units_to_H)
	Set maximum interface depths based on a vector of input values. More...
subroutine, public	set_regrid_max_thickness (CS, max_h, units_to_H)
	Set maximum layer thicknesses based on a vector of input values. More...
real function, dimension(cs%nk), public	getcoordinateresolution (CS, undo_scaling)
	Query the fixed resolution data. More...
real function, dimension(cs%nk+1), public	getcoordinateinterfaces (CS, undo_scaling)
	Query the target coordinate interface positions. More...
character(len=20) function, public	getcoordinateunits (CS)
	Query the target coordinate units. More...
character(len=20) function, public	getcoordinateshortname (CS)
	Query the short name of the coordinate. More...
subroutine, public	set_regrid_params (CS, boundary_extrapolation, min_thickness, old_grid_weight, interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, compress_fraction, ref_pressure, dz_min_surface, nz_fixed_surface, Rho_ML_avg_depth, nlay_ML_to_interior, fix_haloclines, halocline_filt_len, halocline_strat_tol, integrate_downward_for_e, remap_answers_2018, adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin, adaptDrho0)
	Can be used to set any of the parameters for MOM_regridding. More...
integer function, public	get_regrid_size (CS)
	Returns the number of levels/layers in the regridding control structure. More...
type(zlike_cs) function, public	get_zlike_cs (CS)
	This returns a copy of the zlike_CS stored in the regridding control structure. More...
type(sigma_cs) function, public	get_sigma_cs (CS)
	This returns a copy of the sigma_CS stored in the regridding control structure. More...
type(rho_cs) function, public	get_rho_cs (CS)
	This returns a copy of the rho_CS stored in the regridding control structure. More...
real function, dimension(cs%nk), public	getstaticthickness (CS, SSH, depth)
	Return coordinate-derived thicknesses for fixed coordinate systems. More...
subroutine	dz_function1 (string, dz)
	Parses a string and generates a dz(:) profile that goes like k**power. More...
integer function	rho_function1 (string, rho_target)
	Parses a string and generates a rho_target(:) profile with refined resolution downward and returns the number of levels. More...

Variables
character(len= *), parameter, public	regriddingcoordinatemodedoc = " LAYER - Isopycnal or stacked shallow water layers\n"// " ZSTAR, Z* - stretched geopotential z*\n"// " SIGMA_SHELF_ZSTAR - stretched geopotential z* ignoring shelf\n"// " SIGMA - terrain following coordinates\n"// " RHO - continuous isopycnal\n"// " HYCOM1 - HyCOM-like hybrid coordinate\n"// " SLIGHT - stretched coordinates above continuous isopycnal\n"// " ADAPTIVE - optimize for smooth neutral density surfaces"
	Documentation for coordinate options.
character(len= *), parameter, public	regriddinginterpschemedoc = " P1M_H2 (2nd-order accurate)\n"// " P1M_H4 (2nd-order accurate)\n"// " P1M_IH4 (2nd-order accurate)\n"// " PLM (2nd-order accurate)\n"// " PPM_H4 (3rd-order accurate)\n"// " PPM_IH4 (3rd-order accurate)\n"// " P3M_IH4IH3 (4th-order accurate)\n"// " P3M_IH6IH5 (4th-order accurate)\n"// " PQM_IH4IH3 (4th-order accurate)\n"// " PQM_IH6IH5 (5th-order accurate)"
	Documentation for regridding interpolation schemes.
character(len= *), parameter, public	regriddingdefaultinterpscheme = "P1M_H2"
	Default interpolation scheme.
logical, parameter, public	regriddingdefaultboundaryextrapolation = .false.
	Default mode for boundary extrapolation.
real, parameter, public	regriddingdefaultminthickness = 1.e-3
	Default minimum thickness for some coordinate generation modes.

Function/Subroutine Documentation

adjust_interface_motion()

subroutine mom_regridding::adjust_interface_motion ( type(regridding_cs), intent(in)  CS, integer, intent(in)  nk, real, dimension(nk), intent(in)  h_old, real, dimension(cs%nk+1), intent(inout)  dz_int  )

private

Adjust dz_Interface to ensure non-negative future thicknesses.

Parameters

[in]	cs	Regridding control structure
[in]	nk	Number of layers in h_old
[in]	h_old	Minium allowed thickness of h [H ~> m or kg m-2]
[in,out]	dz_int	Minium allowed thickness of h [H ~> m or kg m-2]

Definition at line 1665 of file MOM_regridding.F90.

  type(regridding_CS),      intent(in)    :: CS !< Regridding control structure
  integer,                  intent(in)    :: nk !< Number of layers in h_old
  real, dimension(nk),      intent(in)    :: h_old !< Minium allowed thickness of h [H ~> m or kg m-2]
  real, dimension(CS%nk+1), intent(inout) :: dz_int !< Minium allowed thickness of h [H ~> m or kg m-2]
  ! Local variables
  integer :: k
  real :: h_new, eps, h_total, h_err
 
  eps = 1. ; eps = epsilon(eps)
 
  h_total = 0. ; h_err = 0.
  do k = 1, min(cs%nk,nk)
    h_total = h_total + h_old(k)
    h_err = h_err + max( h_old(k), abs(dz_int(k)), abs(dz_int(k+1)) )*eps
    h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
    if (h_new < -3.0*h_err) then
      write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), &
        'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
        'h_new=',h_new,'h_err=',h_err
      call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
                     'implied h<0 is larger than roundoff!')
    endif
  enddo
  if (cs%nk>nk) then
    do k = nk+1, cs%nk
      h_err = h_err + max( abs(dz_int(k)), abs(dz_int(k+1)) )*eps
      h_new = ( dz_int(k) - dz_int(k+1) )
      if (h_new < -3.0*h_err) then
        write(0,*) 'h<0 at k=',k,'h_old was empty',&
          'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
          'h_new=',h_new,'h_err=',h_err
        call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
                       'implied h<0 is larger than roundoff!')
      endif
    enddo
  endif
  do k = min(cs%nk,nk),2,-1
    h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
    if (h_new<cs%min_thickness) &
      dz_int(k) = ( dz_int(k+1) - h_old(k) ) + cs%min_thickness ! Implies next h_new = min_thickness
    h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
    if (h_new<0.) &
      dz_int(k) = ( 1. - eps ) * ( dz_int(k+1) - h_old(k) ) ! Backup in case min_thickness==0
    h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
    if (h_new<0.) then
      write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), &
        'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
        'h_new=',h_new
      stop 'Still did not work!'
      call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
                     'Repeated adjustment for roundoff h<0 failed!')
    endif
  enddo
 !if (dz_int(1)/=0.) stop 'MOM_regridding: adjust_interface_motion() surface moved'
 

build_grid_adaptive()

subroutine mom_regridding::build_grid_adaptive ( type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)  h, type(thermo_var_ptrs), intent(in)  tv, real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke+1), intent(inout)  dzInterface, type(remapping_cs), intent(in)  remapCS, type(regridding_cs), intent(in)  CS  )

private

This subroutine builds an adaptive grid that follows density surfaces where possible, subject to constraints on the smoothness of interface heights.

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in]	tv	A structure pointing to various thermodynamic variables
[in,out]	dzinterface	The change in interface depth [H ~> m or kg m-2]
[in]	remapcs	The remapping control structure
[in]	cs	Regridding control structure

Definition at line 1532 of file MOM_regridding.F90.

  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
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs),                       intent(in)    :: tv   !< A structure pointing to various
                                                                     !! thermodynamic variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth
                                                                     !! [H ~> m or kg m-2]
  type(remapping_CS),                          intent(in)    :: remapCS !< The remapping control structure
  type(regridding_CS),                         intent(in)    :: CS   !< Regridding control structure
 
  ! local variables
  integer :: i, j, k, nz ! indices and dimension lengths
  ! temperature, salinity and pressure on interfaces
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt, sInt
  ! current interface positions and after tendency term is applied
  ! positive downward
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt ! Interface depths [H ~> m or kg m-2]
  real, dimension(SZK_(GV)+1) :: zNext  ! New interface depths [H ~> m or kg m-2]
 
  nz = gv%ke
 
  ! position surface at z = 0.
  zint(:,:,1) = 0.
 
  ! work on interior interfaces
  do k = 2, nz ; do j = g%jsc-2,g%jec+2 ; do i = g%isc-2,g%iec+2
    tint(i,j,k) = 0.5 * (tv%T(i,j,k-1) + tv%T(i,j,k))
    sint(i,j,k) = 0.5 * (tv%S(i,j,k-1) + tv%S(i,j,k))
    zint(i,j,k) = zint(i,j,k-1) + h(i,j,k-1) ! zInt in [H]
  enddo ; enddo ; enddo
 
  ! top and bottom temp/salt interfaces are just the layer
  ! average values
  tint(:,:,1) = tv%T(:,:,1) ; tint(:,:,nz+1) = tv%T(:,:,nz)
  sint(:,:,1) = tv%S(:,:,1) ; sint(:,:,nz+1) = tv%S(:,:,nz)
 
  ! set the bottom interface depth
  zint(:,:,nz+1)  = zint(:,:,nz) + h(:,:,nz)
 
  ! calculate horizontal density derivatives (alpha/beta)
  ! between cells in a 5-point stencil, columnwise
  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
    if (g%mask2dT(i,j) < 0.5) then
      dzinterface(i,j,:) = 0. ! land point, don't move interfaces, and skip
      cycle
    endif
 
    call build_adapt_column(cs%adapt_CS, g, gv, us, tv, i, j, zint, tint, sint, h, znext)
 
    call filtered_grid_motion(cs, nz, zint(i,j,:), znext, dzinterface(i,j,:))
    ! convert from depth to z
    do k = 1, nz+1 ; dzinterface(i,j,k) = -dzinterface(i,j,k) ; enddo
    call adjust_interface_motion(cs, nz, h(i,j,:), dzinterface(i,j,:))
  enddo ; enddo

build_grid_arbitrary()

subroutine mom_regridding::build_grid_arbitrary ( type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, real, dimension(szi_(g),szj_(g), szk_(gv)), intent(in)  h, real, dimension(szi_(g),szj_(g), szk_(gv)+1), intent(inout)  dzInterface, real, intent(inout)  h_new, type(regridding_cs), intent(in)  CS  )

private

Parameters

[in]	g	Ocean grid structure
[in]	gv	Ocean vertical grid structure
[in]	h	Original layer thicknesses [H ~> m or kg m-2]
[in,out]	dzinterface	The change in interface depth [H ~> m or kg m-2]
[in,out]	h_new	New layer thicknesses [H ~> m or kg m-2]
[in]	cs	Regridding control structure

Definition at line 1726 of file MOM_regridding.F90.

!------------------------------------------------------------------------------
! This routine builds a grid based on arbitrary rules
!------------------------------------------------------------------------------
 
  ! Arguments
  type(ocean_grid_type),                        intent(in)    :: G  !< Ocean grid structure
  type(verticalGrid_type),                      intent(in)    :: GV !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)),   intent(in)    :: h  !< Original layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface
                                                                    !! depth [H ~> m or kg m-2]
  real,                                         intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2]
  type(regridding_CS),                          intent(in)    :: CS !< Regridding control structure
 
  ! Local variables
  integer   :: i, j, k
  integer   :: nz
  real      :: z_inter(SZK_(GV)+1)
  real      :: total_height
  real      :: delta_h
  real      :: max_depth
  real      :: eta              ! local elevation
  real      :: local_depth
  real      :: x1, y1, x2, y2
  real      :: x, t
 
  nz = gv%ke
  max_depth = g%max_depth*gv%Z_to_H
 
  do j = g%jsc-1,g%jec+1
    do i = g%isc-1,g%iec+1
 
      ! Local depth
      local_depth = g%bathyT(i,j)*gv%Z_to_H
 
      ! Determine water column height
      total_height = 0.0
      do k = 1,nz
        total_height = total_height + h(i,j,k)
      enddo
 
      eta = total_height - local_depth
 
      ! Compute new thicknesses based on stretched water column
      delta_h = (max_depth + eta) / nz
 
      ! Define interfaces
      z_inter(1) = eta
      do k = 1,nz
        z_inter(k+1) = z_inter(k) - delta_h
      enddo
 
      ! Refine grid in the middle
      do k = 1,nz+1
        x1 = 0.35; y1 = 0.45; x2 = 0.65; y2 = 0.55
 
        x = - ( z_inter(k) - eta ) / max_depth
 
        if ( x <= x1 ) then
          t = y1*x/x1
        elseif ( (x > x1 ) .and. ( x < x2 )) then
          t = y1 + (y2-y1) * (x-x1) / (x2-x1)
        else
          t = y2 + (1.0-y2) * (x-x2) / (1.0-x2)
        endif
 
        z_inter(k) = -t * max_depth + eta
 
      enddo
 
      ! Modify interface heights to account for topography
      z_inter(nz+1) = - local_depth
 
      ! Modify interface heights to avoid layers of zero thicknesses
      do k = nz,1,-1
        if ( z_inter(k) < (z_inter(k+1) + cs%min_thickness) ) then
          z_inter(k) = z_inter(k+1) + cs%min_thickness
        endif
      enddo
 
      ! Chnage in interface position
      x = 0. ! Left boundary at x=0
      dzinterface(i,j,1) = 0.
      do k = 2,nz
        x = x + h(i,j,k)
        dzinterface(i,j,k) = z_inter(k) - x
      enddo
      dzinterface(i,j,nz+1) = 0.
 
    enddo
  enddo
 
stop 'OOOOOOPS' ! For some reason the gnu compiler will not let me delete this
                ! routine????
 

build_grid_hycom1()

subroutine mom_regridding::build_grid_hycom1 ( type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)  h, type(thermo_var_ptrs), intent(in)  tv, real, dimension( g %isd: g %ied, g %jsd: g %jed,cs%nk), intent(inout)  h_new, real, dimension( g %isd: g %ied, g %jsd: g %jed,cs%nk+1), intent(inout)  dzInterface, type(regridding_cs), intent(in)  CS  )

private

Builds a simple HyCOM-like grid with the deepest location of potential density interpolated from the column profile and a clipping of depth for each interface to a fixed z* or p* grid. This should probably be (optionally?) changed to find the nearest location of the target density.

Remarks

{ Based on Bleck, 2002: An oceanice general circulation model framed in hybrid isopycnic-Cartesian coordinates, Ocean Modelling 37, 55-88. http://dx.doi.org/10.1016/S1463-5003(01)00012-9 }

Parameters

[in]	g	Grid structure
[in]	gv	Ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	h	Existing model thickness [H ~> m or kg m-2]
[in]	tv	Thermodynamics structure
[in]	cs	Regridding control structure
[in,out]	h_new	New layer thicknesses [H ~> m or kg m-2]
[in,out]	dzinterface	Changes in interface position

Definition at line 1460 of file MOM_regridding.F90.

  type(ocean_grid_type),                     intent(in)    :: G  !< Grid structure
  type(verticalGrid_type),                   intent(in)    :: GV !< Ocean vertical grid structure
  type(unit_scale_type),                     intent(in)    :: US !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h  !< Existing model thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),                     intent(in)    :: tv !< Thermodynamics structure
  type(regridding_CS),                       intent(in)    :: CS !< Regridding control structure
  real, dimension(SZI_(G),SZJ_(G),CS%nk),    intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),CS%nk+1),  intent(inout) :: dzInterface !< Changes in interface position
 
  ! Local variables
  real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2]
  real, dimension(CS%nk+1) :: z_col_new ! New interface positions relative to the surface [H ~> m or kg m-2]
  real, dimension(SZK_(GV)+1) :: dz_col  ! The realized change in z_col [H ~> m or kg m-2]
  real, dimension(SZK_(GV))   :: p_col   ! Layer center pressure [R L2 T-2 ~> Pa]
  real :: ref_pres  ! The reference pressure [R L2 T-2 ~> Pa]
  integer   :: i, j, k, nki
  real :: depth
  real :: h_neglect, h_neglect_edge
 
  if (.not.cs%remap_answers_2018) then
    h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
  elseif (gv%Boussinesq) then
    h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
  else
    h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
  endif
 
  if (.not.cs%target_density_set) call mom_error(fatal, "build_grid_HyCOM1 : "//&
        "Target densities must be set before build_grid_HyCOM1 is called.")
 
  nki = min(gv%ke, cs%nk)
 
  ! Build grid based on target interface densities
  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
    if (g%mask2dT(i,j)>0.) then
 
      depth = g%bathyT(i,j) * gv%Z_to_H
 
      z_col(1) = 0. ! Work downward rather than bottom up
      do k = 1, gv%ke
        z_col(k+1) = z_col(k) + h(i,j,k)
        p_col(k) = tv%P_Ref + cs%compressibility_fraction * &
             ( 0.5 * ( z_col(k) + z_col(k+1) ) * (gv%H_to_RZ*gv%g_Earth) - tv%P_Ref )
      enddo
 
      call build_hycom1_column(cs%hycom_CS, tv%eqn_of_state, gv%ke, depth, &
                               h(i,j,:), tv%T(i,j,:), tv%S(i,j,:), p_col, &
                               z_col, z_col_new, zscale=gv%Z_to_H, &
                               h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
 
      ! Calculate the final change in grid position after blending new and old grids
      call filtered_grid_motion( cs, gv%ke, z_col, z_col_new, dz_col )
 
      ! This adjusts things robust to round-off errors
      dz_col(:) = -dz_col(:)
      call adjust_interface_motion( cs, gv%ke, h(i,j,:), dz_col(:) )
 
      dzinterface(i,j,1:nki+1) = dz_col(1:nki+1)
      if (nki<cs%nk) dzinterface(i,j,nki+2:cs%nk+1) = 0.
 
    else ! on land
      dzinterface(i,j,:) = 0.
    endif ! mask2dT
  enddo ; enddo ! i,j
 
  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
 

build_grid_slight()

subroutine mom_regridding::build_grid_slight ( type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)  h, type(thermo_var_ptrs), intent(in)  tv, real, dimension(szi_(g),szj_(g),szk_(gv)+1), intent(inout)  dzInterface, type(regridding_cs), intent(in)  CS  )

private

Builds a grid that tracks density interfaces for water that is denser than the surface density plus an increment of some number of layers, and uses all lighter layers uniformly above this location. Note that this amounts to interpolating to find the depth of an arbitrary (non-integer) interface index which should make the results vary smoothly in space to the extent that the surface density and interior stratification vary smoothly in space. Over shallow topography, this will tend to give a uniform sigma-like coordinate. For sufficiently shallow water, a minimum grid spacing is used to avoid certain instabilities.

Parameters

[in]	g	Grid structure
[in]	gv	Ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	h	Existing model thickness [H ~> m or kg m-2]
[in]	tv	Thermodynamics structure
[in,out]	dzinterface	Changes in interface position
[in]	cs	Regridding control structure

Definition at line 1599 of file MOM_regridding.F90.

  type(ocean_grid_type),                       intent(in)    :: G  !< Grid structure
  type(verticalGrid_type),                     intent(in)    :: GV !< Ocean vertical grid structure
  type(unit_scale_type),                       intent(in)    :: US !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in)    :: h  !< Existing model thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),                       intent(in)    :: tv !< Thermodynamics structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< Changes in interface position
  type(regridding_CS),                         intent(in)    :: CS !< Regridding control structure
 
  real, dimension(SZK_(GV)+1) :: z_col   ! Interface positions relative to the surface [H ~> m or kg m-2]
  real, dimension(SZK_(GV)+1) :: z_col_new ! Interface positions relative to the surface [H ~> m or kg m-2]
  real, dimension(SZK_(GV)+1) :: dz_col  ! The realized change in z_col [H ~> m or kg m-2]
  real, dimension(SZK_(GV))   :: p_col   ! Layer center pressure [R L2 T-2 ~> Pa]
  real :: depth
  integer :: i, j, k, nz
  real :: h_neglect, h_neglect_edge
 
  if (.not.cs%remap_answers_2018) then
    h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
  elseif (gv%Boussinesq) then
    h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
  else
    h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
  endif
 
  nz = gv%ke
 
  if (.not.cs%target_density_set) call mom_error(fatal, "build_grid_SLight : "//&
        "Target densities must be set before build_grid_SLight is called.")
 
  ! Build grid based on target interface densities
  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
    if (g%mask2dT(i,j)>0.) then
 
      depth = g%bathyT(i,j) * gv%Z_to_H
      z_col(1) = 0. ! Work downward rather than bottom up
      do k=1,nz
        z_col(k+1) = z_col(k) + h(i,j,k)
        p_col(k) = tv%P_Ref + cs%compressibility_fraction * &
                    ( 0.5 * ( z_col(k) + z_col(k+1) ) * (gv%H_to_RZ*gv%g_Earth) - tv%P_Ref )
      enddo
 
      call build_slight_column(cs%slight_CS, tv%eqn_of_state, gv%H_to_RZ*gv%g_Earth, &
                          gv%H_subroundoff, nz, depth, h(i, j, :), &
                          tv%T(i, j, :), tv%S(i, j, :), p_col, z_col, z_col_new, &
                          h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
 
      ! Calculate the final change in grid position after blending new and old grids
      call filtered_grid_motion( cs, nz, z_col, z_col_new, dz_col )
      do k=1,nz+1 ; dzinterface(i,j,k) = -dz_col(k) ; enddo
#ifdef __DO_SAFETY_CHECKS__
      if (dzinterface(i,j,1) /= 0.) stop 'build_grid_SLight: Surface moved?!'
      if (dzinterface(i,j,nz+1) /= 0.) stop 'build_grid_SLight: Bottom moved?!'
#endif
 
      ! This adjusts things robust to round-off errors
      call adjust_interface_motion( cs, nz, h(i,j,:), dzinterface(i,j,:) )
 
    else ! on land
      dzinterface(i,j,:) = 0.
    endif ! mask2dT
  enddo ; enddo ! i,j
 

build_rho_grid()

subroutine mom_regridding::build_rho_grid ( type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)  h, type(thermo_var_ptrs), intent(in)  tv, real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke+1), intent(inout)  dzInterface, type(remapping_cs), intent(in)  remapCS, type(regridding_cs), intent(in)  CS  )

private

This routine builds a new grid based on a given set of target interface densities.

Parameters

[in]	g	Ocean grid structure
[in]	gv	Ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in]	tv	Thermodynamics structure
[in,out]	dzinterface	The change in interface depth [H ~> m or kg m-2]
[in]	remapcs	The remapping control structure
[in]	cs	Regridding control structure

Definition at line 1327 of file MOM_regridding.F90.

!------------------------------------------------------------------------------
! This routine builds a new grid based on a given set of target interface
! densities (these target densities are computed by taking the mean value
! of given layer densities). 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.
!------------------------------------------------------------------------------
 
  ! Arguments
  type(ocean_grid_type),                        intent(in)    :: G  !< Ocean grid structure
  type(verticalGrid_type),                      intent(in)    :: GV !< Ocean vertical grid structure
  type(unit_scale_type),                        intent(in)    :: US !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),    intent(in)    :: h  !< Layer thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs),                        intent(in)    :: tv !< Thermodynamics structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth
                                                                    !! [H ~> m or kg m-2]
  type(remapping_CS),                           intent(in)    :: remapCS !< The remapping control structure
  type(regridding_CS),                          intent(in)    :: CS !< Regridding control structure
 
  ! Local variables
  integer :: nz
  integer :: i, j, k
  real    :: nominalDepth   ! Depth of the bottom of the ocean, positive downward [H ~> m or kg m-2]
  real, dimension(SZK_(GV)+1) :: zOld, zNew ! Old and new interface heights [H ~> m or kg m-2]
  real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2]
#ifdef __DO_SAFETY_CHECKS__
  real    :: totalThickness
  real    :: dh
#endif
 
  if (.not.cs%remap_answers_2018) then
    h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
  elseif (gv%Boussinesq) then
    h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
  else
    h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
  endif
 
  nz = gv%ke
 
  if (.not.cs%target_density_set) call mom_error(fatal, "build_rho_grid: "//&
        "Target densities must be set before build_rho_grid is called.")
 
  ! Build grid based on target interface densities
  do j = g%jsc-1,g%jec+1
    do i = g%isc-1,g%iec+1
 
      if (g%mask2dT(i,j)==0.) then
        dzinterface(i,j,:) = 0.
        cycle
      endif
 
 
      ! Local depth (G%bathyT is positive)
      nominaldepth = g%bathyT(i,j)*gv%Z_to_H
 
      call build_rho_column(cs%rho_CS, nz, nominaldepth, h(i, j, :), &
                            tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, znew, &
                            h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
 
      if (cs%integrate_downward_for_e) then
        zold(1) = 0.
        do k = 1,nz
          zold(k+1) = zold(k) - h(i,j,k)
        enddo
      else
        ! The rest of the model defines grids integrating up from the bottom
        zold(nz+1) = - nominaldepth
        do k = nz,1,-1
          zold(k) = zold(k+1) + h(i,j,k)
        enddo
      endif
 
      ! Calculate the final change in grid position after blending new and old grids
      call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )
 
#ifdef __DO_SAFETY_CHECKS__
      do k = 2,nz
        if (znew(k) > zold(1)) then
          write(0,*) 'zOld=',zold
          write(0,*) 'zNew=',znew
          call mom_error( fatal, 'MOM_regridding, build_rho_grid: '//&
               'interior interface above surface!' )
        endif
        if (znew(k) > znew(k-1)) then
          write(0,*) 'zOld=',zold
          write(0,*) 'zNew=',znew
          call mom_error( fatal, 'MOM_regridding, build_rho_grid: '//&
               'interior interfaces cross!' )
        endif
      enddo
 
      totalthickness = 0.0
      do k = 1,nz
        totalthickness = totalthickness + h(i,j,k)
      enddo
 
      dh=max(nominaldepth,totalthickness)
      if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh) then
        write(0,*) 'min_thickness=',cs%min_thickness
        write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
        write(0,*) 'zNew(1)-zOld(1) = ',znew(1)-zold(1),epsilon(dh),nz
        do k=1,nz+1
          write(0,*) k,zold(k),znew(k)
        enddo
        do k=1,nz
          write(0,*) k,h(i,j,k),znew(k)-znew(k+1)
        enddo
        call mom_error( fatal, &
               'MOM_regridding, build_rho_grid: top surface has moved!!!' )
      endif
#endif
 
    enddo  ! end loop on i
  enddo  ! end loop on j
 

build_sigma_grid()

subroutine mom_regridding::build_sigma_grid ( type(regridding_cs), intent(in)  CS, type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)  h, real, dimension( g %isd: g %ied, g %jsd: g %jed, cs%nk+1), intent(inout)  dzInterface  )

private

This routine builds a grid based on terrain-following coordinates.

Parameters

[in]	cs	Regridding control structure
[in]	g	Ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in,out]	dzinterface	The change in interface depth [H ~> m or kg m-2]

Definition at line 1248 of file MOM_regridding.F90.

!------------------------------------------------------------------------------
! This routine builds a grid based on terrain-following coordinates.
! The module parameter coordinateResolution(:) determines the resolution in
! sigma coordinate, dSigma(:). sigma-coordinates are defined by
!   sigma = (eta-z)/(H+eta)  s.t. sigma=0 at z=eta and sigma=1 at z=-H .
!------------------------------------------------------------------------------
 
  ! Arguments
  type(regridding_CS),                       intent(in)    :: CS !< Regridding control structure
  type(ocean_grid_type),                     intent(in)    :: G  !< Ocean grid structure
  type(verticalGrid_type),                   intent(in)    :: GV !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h  !< Layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
                                                                 !! [H ~> m or kg m-2]
 
  ! Local variables
  integer :: i, j, k
  integer :: nz
  real    :: nominalDepth, totalThickness, dh
  real, dimension(SZK_(GV)+1) :: zOld, zNew
 
  nz = gv%ke
 
  do i = g%isc-1,g%iec+1
    do j = g%jsc-1,g%jec+1
 
      if (g%mask2dT(i,j)==0.) then
        dzinterface(i,j,:) = 0.
        cycle
      endif
 
      ! The rest of the model defines grids integrating up from the bottom
      nominaldepth = g%bathyT(i,j)*gv%Z_to_H
 
      ! Determine water column height
      totalthickness = 0.0
      do k = 1,nz
        totalthickness = totalthickness + h(i,j,k)
      enddo
 
      call build_sigma_column(cs%sigma_CS, nominaldepth, totalthickness, znew)
 
      ! Calculate the final change in grid position after blending new and old grids
      zold(nz+1) =  -nominaldepth
      do k = nz,1,-1
        zold(k) = zold(k+1) + h(i, j, k)
      enddo
 
      call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )
 
#ifdef __DO_SAFETY_CHECKS__
      dh=max(nominaldepth,totalthickness)
      if (abs(znew(1)-zold(1))>(cs%nk-1)*0.5*epsilon(dh)*dh) then
        write(0,*) 'min_thickness=',cs%min_thickness
        write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
        write(0,*) 'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz,cs%nk
        do k=1,nz+1
          write(0,*) k,zold(k),znew(k)
        enddo
        do k=1,cs%nk
          write(0,*) k,h(i,j,k),znew(k)-znew(k+1),totalthickness*cs%coordinateResolution(k),cs%coordinateResolution(k)
        enddo
        call mom_error( fatal, &
               'MOM_regridding, build_sigma_grid: top surface has moved!!!' )
      endif
      dzinterface(i,j,1) = 0.
      dzinterface(i,j,cs%nk+1) = 0.
#endif
 
    enddo
  enddo
 

build_zstar_grid()

subroutine mom_regridding::build_zstar_grid ( type(regridding_cs), intent(in)  CS, type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)  h, real, dimension( g %isd: g %ied, g %jsd: g %jed, cs%nk+1), intent(inout)  dzInterface, real, dimension(:,:), optional, pointer  frac_shelf_h  )

private

Builds a z*-ccordinate grid with partial steps (Adcroft and Campin, 2004). z* is defined as z* = (z-eta)/(H+eta)*H s.t. z*=0 when z=eta and z*=-H when z=-H .

Parameters

[in]	cs	Regridding control structure
[in]	g	Ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in,out]	dzinterface	The change in interface depth [H ~> m or kg m-2].
	frac_shelf_h	Fractional ice shelf coverage [nondim].

Definition at line 1158 of file MOM_regridding.F90.

 
  ! Arguments
  type(regridding_CS),                       intent(in)    :: CS !< Regridding control structure
  type(ocean_grid_type),                     intent(in)    :: G  !< Ocean grid structure
  type(verticalGrid_type),                   intent(in)    :: GV !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h  !< Layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
                                                                 !! [H ~> m or kg m-2].
  real, dimension(:,:),            optional, pointer       :: frac_shelf_h !< Fractional ice shelf coverage [nondim].
  ! Local variables
  real    :: nominalDepth, totalThickness, dh  ! Depths and thicknesses [H ~> m or kg m-2]
  real, dimension(SZK_(GV)+1) :: zOld, zNew    ! Coordinate interface heights [H ~> m or kg m-2]
  integer :: i, j, k, nz
  logical :: ice_shelf
 
  nz = gv%ke
  ice_shelf = .false.
  if (present(frac_shelf_h)) then
    if (associated(frac_shelf_h)) ice_shelf = .true.
  endif
 
!$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h,ice_shelf) &
!$OMP                          private(nominalDepth,totalThickness,zNew,dh,zOld)
  do j = g%jsc-1,g%jec+1
    do i = g%isc-1,g%iec+1
 
      if (g%mask2dT(i,j)==0.) then
        dzinterface(i,j,:) = 0.
        cycle
      endif
 
      ! Local depth (G%bathyT is positive)
      nominaldepth = g%bathyT(i,j)*gv%Z_to_H
 
      ! Determine water column thickness
      totalthickness = 0.0
      do k = 1,nz
        totalthickness = totalthickness + h(i,j,k)
      enddo
 
      zold(nz+1) = - nominaldepth
      do k = nz,1,-1
        zold(k) = zold(k+1) + h(i,j,k)
      enddo
 
      if (ice_shelf) then
        if (frac_shelf_h(i,j) > 0.) then ! under ice shelf
          call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, znew, &
                                z_rigid_top = totalthickness-nominaldepth, &
                                eta_orig=zold(1), zscale=gv%Z_to_H)
        else
          call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, &
                                znew, zscale=gv%Z_to_H)
        endif
      else
        call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, &
                                znew, zscale=gv%Z_to_H)
      endif
 
      ! Calculate the final change in grid position after blending new and old grids
      call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )
 
#ifdef __DO_SAFETY_CHECKS__
      dh=max(nominaldepth,totalthickness)
      if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh) then
        write(0,*) 'min_thickness=',cs%min_thickness
        write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
        write(0,*) 'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz
        do k=1,nz+1
          write(0,*) k,zold(k),znew(k)
        enddo
        do k=1,nz
          write(0,*) k,h(i,j,k),znew(k)-znew(k+1),cs%coordinateResolution(k)
        enddo
        call mom_error( fatal, &
               'MOM_regridding, build_zstar_grid(): top surface has moved!!!' )
      endif
#endif
 
      call adjust_interface_motion( cs, nz, h(i,j,:), dzinterface(i,j,:) )
 
    enddo
  enddo
 

calc_h_new_by_dz()

subroutine mom_regridding::calc_h_new_by_dz ( type(regridding_cs), intent(in)  CS, type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)  h, real, dimension( g %isd: g %ied, g %jsd: g %jed,cs%nk+1), intent(in)  dzInterface, real, dimension( g %isd: g %ied, g %jsd: g %jed,cs%nk), intent(inout)  h_new  )

private

Calculates h_new from h + delta_k dzInterface.

Parameters

[in]	cs	Regridding control structure
[in]	g	Grid structure
[in]	gv	Ocean vertical grid structure
[in]	h	Old layer thicknesses (arbitrary units)
[in]	dzinterface	Change in interface positions (same as h)
[in,out]	h_new	New layer thicknesses (same as h)

Definition at line 906 of file MOM_regridding.F90.

  type(regridding_CS),                       intent(in)    :: CS !< Regridding control structure
  type(ocean_grid_type),                     intent(in)    :: G  !< Grid structure
  type(verticalGrid_type),                   intent(in)    :: GV !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h  !< Old layer thicknesses (arbitrary units)
  real, dimension(SZI_(G),SZJ_(G),CS%nk+1),  intent(in)    :: dzInterface !< Change in interface positions (same as h)
  real, dimension(SZI_(G),SZJ_(G),CS%nk),    intent(inout) :: h_new !< New layer thicknesses (same as h)
  ! Local variables
  integer :: i, j, k, nki
 
  nki = min(cs%nk, gv%ke)
 
  !$OMP parallel do default(shared)
  do j = g%jsc-1,g%jec+1
    do i = g%isc-1,g%iec+1
      if (g%mask2dT(i,j)>0.) then
        do k=1,nki
          h_new(i,j,k) = max( 0., h(i,j,k) + ( dzinterface(i,j,k) - dzinterface(i,j,k+1) ) )
        enddo
        if (cs%nk > gv%ke) then
          do k=nki+1, cs%nk
            h_new(i,j,k) = max( 0., dzinterface(i,j,k) - dzinterface(i,j,k+1) )
          enddo
        endif
      else
        h_new(i,j,1:nki) = h(i,j,1:nki)
        if (cs%nk > gv%ke) h_new(i,j,nki+1:cs%nk) = 0.
        ! On land points, why are we keeping the original h rather than setting to zero? -AJA
      endif
    enddo
  enddo
 

check_grid_column()

subroutine, public mom_regridding::check_grid_column	(	integer, intent(in)	nk,
		real, intent(in)	depth,
		real, dimension(nk), intent(in)	h,
		real, dimension(nk+1), intent(in)	dzInterface,
		character(len=*), intent(in)	msg
	)

Check that the total thickness of new and old grids are consistent.

Parameters

[in]	nk	Number of cells
[in]	depth	Depth of bottom [Z ~> m] or arbitrary units
[in]	h	Cell thicknesses [Z ~> m] or arbitrary units
[in]	dzinterface	Change in interface positions (same units as h)
[in]	msg	Message to append to errors

Definition at line 959 of file MOM_regridding.F90.

  integer,               intent(in) :: nk !< Number of cells
  real,                  intent(in) :: depth !< Depth of bottom [Z ~> m] or arbitrary units
  real, dimension(nk),   intent(in) :: h  !< Cell thicknesses [Z ~> m] or arbitrary units
  real, dimension(nk+1), intent(in) :: dzInterface !< Change in interface positions (same units as h)
  character(len=*),      intent(in) :: msg !< Message to append to errors
  ! Local variables
  integer :: k
  real    :: eps, total_h_old, total_h_new, h_new, z_old, z_new
 
  eps =1. ; eps = epsilon(eps)
 
  ! Total thickness of grid h
  total_h_old = 0.
  do k = 1,nk
    total_h_old = total_h_old + h(k)
  enddo
 
  ! Integrate upwards for the interfaces consistent with the rest of MOM6
  z_old = - depth
  if (depth == 0.) z_old = - total_h_old
  total_h_new = 0.
  do k = nk,1,-1
    z_old = z_old + h(k) ! Old interface position above layer k
    z_new = z_old + dzinterface(k) ! New interface position based on dzInterface
    h_new = h(k) + ( dzinterface(k) - dzinterface(k+1) ) ! New thickness
    if (h_new<0.) then
      write(0,*) 'k,h,hnew=',k,h(k),h_new
      write(0,*) 'dzI(k+1),dzI(k)=',dzinterface(k+1),dzinterface(k)
      call mom_error( fatal, 'MOM_regridding, check_grid_column: '//&
        'Negative layer thickness implied by re-gridding, '//trim(msg))
    endif
    total_h_new = total_h_new + h_new
 
  enddo
 
  ! Conservation by implied h_new
  if (abs(total_h_new-total_h_old)>real(nk-1)*0.5*(total_h_old+total_h_new)*eps) then
    write(0,*) 'nk=',nk
    do k = 1,nk
      write(0,*) 'k,h,hnew=',k,h(k),h(k)+(dzinterface(k)-dzinterface(k+1))
    enddo
    write(0,*) 'Hold,Hnew,Hnew-Hold=',total_h_old,total_h_new,total_h_new-total_h_old
    write(0,*) 'eps,(n)/2*eps*H=',eps,real(nk-1)*0.5*(total_h_old+total_h_new)*eps
    call mom_error( fatal, 'MOM_regridding, check_grid_column: '//&
      'Re-gridding did NOT conserve total thickness to within roundoff '//trim(msg))
  endif
 
  ! Check that the top and bottom are intentionally moving
  if (dzinterface(1) /= 0.) call mom_error( fatal, &
    'MOM_regridding, check_grid_column: Non-zero dzInterface at surface! '//trim(msg))
  if (dzinterface(nk+1) /= 0.) call mom_error( fatal, &
    'MOM_regridding, check_grid_column: Non-zero dzInterface at bottom! '//trim(msg))
 

check_grid_def()

subroutine mom_regridding::check_grid_def ( character(len=*), intent(in)  filename, character(len=*), intent(in)  varname, character(len=*), intent(in)  expected_units, character(len=*), intent(inout)  msg, logical, intent(out)  ierr  )

private

Do some basic checks on the vertical grid definition file, variable.

Parameters

[in]	filename	File name
[in]	varname	Variable name
[in]	expected_units	Expected units of variable
[in,out]	msg	Message to use for errors
[out]	ierr	True if an error occurs

Definition at line 737 of file MOM_regridding.F90.

  character(len=*), intent(in)    :: filename !< File name
  character(len=*), intent(in)    :: varname !< Variable name
  character(len=*), intent(in)    :: expected_units !< Expected units of variable
  character(len=*), intent(inout) :: msg !< Message to use for errors
  logical,          intent(out)   :: ierr !< True if an error occurs
  ! Local variables
  character (len=200) :: units, long_name
  integer :: ncid, status, intid, vid
  integer :: i
 
  ierr = .false.
  status = nf90_open(trim(filename), nf90_nowrite, ncid)
  if (status /= nf90_noerr) then
    ierr = .true.
    msg = 'File not found: '//trim(filename)
    return
  endif
 
  status = nf90_inq_varid(ncid, trim(varname), vid)
  if (status /= nf90_noerr) then
    ierr = .true.
    msg = 'Var not found: '//trim(varname)
    return
  endif
 
  status = nf90_get_att(ncid, vid, "units", units)
  if (status /= nf90_noerr) then
    ierr = .true.
    msg = 'Attribute not found: units'
    return
  endif
  ! NF90_GET_ATT can return attributes with null characters, which TRIM will not truncate.
  ! This loop replaces any null characters with a space so that the following check between
  ! the read units and the expected units will pass
  do i=1,len_trim(units)
    if (units(i:i) == char(0)) units(i:i) = " "
  enddo
 
  if (trim(units) /= trim(expected_units)) then
    if (trim(expected_units) == "meters") then
      if (trim(units) /= "m") then
        ierr = .true.
      endif
    else
      ierr = .true.
    endif
  endif
 
  if (ierr) then
    msg = 'Units incorrect: '//trim(units)//' /= '//trim(expected_units)
  endif
 

check_remapping_grid()

subroutine, public mom_regridding::check_remapping_grid	(	type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)	h,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke+1), intent(in)	dzInterface,
		character(len=*), intent(in)	msg
	)

Check that the total thickness of two grids match.

Parameters

[in]	g	Grid structure
[in]	gv	Ocean vertical grid structure
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in]	dzinterface	Change in interface positions [H ~> m or kg m-2]
[in]	msg	Message to append to errors

Definition at line 941 of file MOM_regridding.F90.

  type(ocean_grid_type),                       intent(in) :: G   !< Grid structure
  type(verticalGrid_type),                     intent(in) :: GV  !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in) :: h   !< Layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzInterface !< Change in interface positions
                                                                 !! [H ~> m or kg m-2]
  character(len=*),                            intent(in) :: msg !< Message to append to errors
  ! Local variables
  integer :: i, j
 
  !$OMP parallel do default(shared)
  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
    if (g%mask2dT(i,j)>0.) call check_grid_column( gv%ke, gv%Z_to_H*g%bathyT(i,j), h(i,j,:), dzinterface(i,j,:), msg )
  enddo ; enddo
 

convective_adjustment()

subroutine mom_regridding::convective_adjustment ( type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, real, dimension(szi_(g),szj_(g),szk_(gv)), intent(inout)  h, type(thermo_var_ptrs), intent(inout)  tv  )

private

Achieve convective adjustment by swapping layers.

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in,out]	h	Layer thicknesses [H ~> m or kg m-2]
[in,out]	tv	A structure pointing to various thermodynamic variables

Definition at line 1869 of file MOM_regridding.F90.

  type(ocean_grid_type),   intent(in)    :: G    !< The ocean's grid structure
  type(verticalGrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: h    !< Layer thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs),   intent(inout) :: tv   !< A structure pointing to various thermodynamic variables
!------------------------------------------------------------------------------
! Check each water column to see whether it is stratified. If not, sort the
! layers by successive swappings of water masses (bubble sort algorithm)
!------------------------------------------------------------------------------
 
  ! Local variables
  integer   :: i, j, k
  real      :: T0, T1       ! temperatures
  real      :: S0, S1       ! salinities
  real      :: r0, r1       ! densities
  real      :: h0, h1
  logical   :: stratified
  real, dimension(GV%ke) :: p_col, densities
 
  p_col(:) = 0.
 
  ! Loop on columns
  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
 
    ! Compute densities within current water column
    call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, densities, tv%eqn_of_state)
 
    ! Repeat restratification until complete
    do
      stratified = .true.
      do k = 1,gv%ke-1
        ! Gather information of current and next cells
        t0 = tv%T(i,j,k)  ; t1 = tv%T(i,j,k+1)
        s0 = tv%S(i,j,k)  ; s1 = tv%S(i,j,k+1)
        r0 = densities(k) ; r1 = densities(k+1)
        h0 = h(i,j,k) ; h1 = h(i,j,k+1)
        ! If the density of the current cell is larger than the density
        ! below it, we swap the cells and recalculate the densitiies
        ! within the swapped cells
        if ( r0 > r1 ) then
          tv%T(i,j,k) = t1 ; tv%T(i,j,k+1) = t0
          tv%S(i,j,k) = s1 ; tv%S(i,j,k+1) = s0
          h(i,j,k)    = h1 ; h(i,j,k+1)    = h0
          ! Recompute densities at levels k and k+1
          call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), densities(k), tv%eqn_of_state)
          call calculate_density( tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), &
                                  densities(k+1), tv%eqn_of_state )
          stratified = .false.
        endif
      enddo  ! k
 
      if ( stratified ) exit
    enddo
 
  enddo ; enddo  ! i & j
 

dz_function1()

subroutine mom_regridding::dz_function1 ( character(len=*), intent(in)  string, real, dimension(:), intent(inout)  dz  )

private

Parses a string and generates a dz(:) profile that goes like k**power.

Parameters

[in]	string	String with list of parameters in form dz_min, H_total, power, precision
[in,out]	dz	Profile of nominal thicknesses

Definition at line 2417 of file MOM_regridding.F90.

  character(len=*),   intent(in)    :: string !< String with list of parameters in form
                                              !! dz_min, H_total, power, precision
  real, dimension(:), intent(inout) :: dz     !< Profile of nominal thicknesses
  ! Local variables
  integer :: nk, k
  real    :: dz_min, power, prec, H_total
 
  nk = size(dz) ! Number of cells
  prec = -1024.
  read( string, *) dz_min, h_total, power, prec
  if (prec == -1024.) call mom_error(fatal,"dz_function1: "// &
          "Problem reading FNC1: string  ="//trim(string))
  ! Create profile of ( dz - dz_min )
  do k = 1, nk
    dz(k) = (real(k-1)/real(nk-1))**power
  enddo
  dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total
  dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
  dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total
  dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
  dz(nk) = dz(nk) + ( h_total - sum( dz(:) + dz_min ) ) ! Adjust bottommost layer
  dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
  dz(:) = dz(:) + dz_min ! Finally add in the constant dz_min
 

end_regridding()

subroutine, public mom_regridding::end_regridding

(

type(regridding_cs), intent(inout)

CS

)

Deallocation of regridding memory.

Parameters

[in,out]

cs

Regridding control structure

Definition at line 793 of file MOM_regridding.F90.

  type(regridding_CS), intent(inout) :: CS !< Regridding control structure
 
  if (associated(cs%zlike_CS))  call end_coord_zlike(cs%zlike_CS)
  if (associated(cs%sigma_CS))  call end_coord_sigma(cs%sigma_CS)
  if (associated(cs%rho_CS))    call end_coord_rho(cs%rho_CS)
  if (associated(cs%hycom_CS))  call end_coord_hycom(cs%hycom_CS)
  if (associated(cs%slight_CS)) call end_coord_slight(cs%slight_CS)
  if (associated(cs%adapt_CS))  call end_coord_adapt(cs%adapt_CS)
 
  deallocate( cs%coordinateResolution )
  if (allocated(cs%target_density)) deallocate( cs%target_density )
  if (allocated(cs%max_interface_depths) ) deallocate( cs%max_interface_depths )
  if (allocated(cs%max_layer_thickness) ) deallocate( cs%max_layer_thickness )
 

filtered_grid_motion()

subroutine mom_regridding::filtered_grid_motion ( type(regridding_cs), intent(in)  CS, integer, intent(in)  nk, real, dimension(nk+1), intent(in)  z_old, real, dimension(cs%nk+1), intent(in)  z_new, real, dimension(cs%nk+1), intent(inout)  dz_g  )

private

Returns the change in interface position motion after filtering and assuming the top and bottom interfaces do not move. The filtering is a function of depth, and is applied as the integrated average filtering over the trajectory of the interface. By design, this code can not give tangled interfaces provided that z_old and z_new are not already tangled.

Parameters

[in]	cs	Regridding control structure
[in]	nk	Number of cells in source grid
[in]	z_old	Old grid position [H ~> m or kg m-2]
[in]	z_new	New grid position [H ~> m or kg m-2]
[in,out]	dz_g	Change in interface positions [H ~> m or kg m-2]

Definition at line 1020 of file MOM_regridding.F90.

  type(regridding_CS),      intent(in)    :: CS !< Regridding control structure
  integer,                  intent(in)    :: nk !< Number of cells in source grid
  real, dimension(nk+1),    intent(in)    :: z_old !< Old grid position [H ~> m or kg m-2]
  real, dimension(CS%nk+1), intent(in)    :: z_new !< New grid position [H ~> m or kg m-2]
  real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions [H ~> m or kg m-2]
  ! Local variables
  real :: sgn  ! The sign convention for downward.
  real :: dz_tgt, zr1, z_old_k
  real :: Aq, Bq, dz0, z0, F0
  real :: zs, zd, dzwt, Idzwt
  real :: wtd, Iwtd
  real :: Int_zs, Int_zd, dInt_zs_zd
! For debugging:
  real, dimension(nk+1) :: z_act
!  real, dimension(nk+1) :: ddz_g_s, ddz_g_d
  logical :: debug = .false.
  integer :: k
 
  if ((z_old(nk+1) - z_old(1)) * (z_new(cs%nk+1) - z_new(1)) < 0.0) then
    call mom_error(fatal, "filtered_grid_motion: z_old and z_new use different sign conventions.")
  elseif ((z_old(nk+1) - z_old(1)) * (z_new(cs%nk+1) - z_new(1)) == 0.0) then
    ! This is a massless column, so do nothing and return.
    do k=1,cs%nk+1 ; dz_g(k) = 0.0 ; enddo ; return
  elseif ((z_old(nk+1) - z_old(1)) + (z_new(cs%nk+1) - z_new(1)) > 0.0) then
    sgn = 1.0
  else
    sgn = -1.0
  endif
 
  if (debug) then
    do k=2,cs%nk+1
      if (sgn*(z_new(k)-z_new(k-1)) < -5e-16*(abs(z_new(k))+abs(z_new(k-1))) ) &
        call mom_error(fatal, "filtered_grid_motion: z_new is tangled.")
    enddo
    do k=2,nk+1
      if (sgn*(z_old(k)-z_old(k-1)) < -5e-16*(abs(z_old(k))+abs(z_old(k-1))) ) &
        call mom_error(fatal, "filtered_grid_motion: z_old is tangled.")
    enddo
    ! ddz_g_s(:) = 0.0 ; ddz_g_d(:) = 0.0
  endif
 
  zs = cs%depth_of_time_filter_shallow
  zd = cs%depth_of_time_filter_deep
  wtd = 1.0 - cs%old_grid_weight
  iwtd = 1.0 / wtd
 
  dzwt = (zd - zs)
  idzwt = 0.0 ; if (abs(zd - zs) > 0.0) idzwt = 1.0 / (zd - zs)
  dint_zs_zd = 0.5*(1.0 + iwtd) * (zd - zs)
  aq = 0.5*(iwtd - 1.0)
 
  dz_g(1) = 0.0
  z_old_k = z_old(1)
  do k = 2,cs%nk+1
    if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model
    ! zr1 is positive and increases with depth, and dz_tgt is positive downward.
    dz_tgt = sgn*(z_new(k) - z_old_k)
    zr1 = sgn*(z_old_k - z_old(1))
 
    !   First, handle the two simple and common cases that do not pass through
    ! the adjustment rate transition zone.
    if ((zr1 > zd) .and. (zr1 + wtd * dz_tgt > zd)) then
      dz_g(k) = sgn * wtd * dz_tgt
    elseif ((zr1 < zs) .and. (zr1 + dz_tgt < zs)) then
      dz_g(k) = sgn * dz_tgt
    else
      ! Find the new value by inverting the equation
      !   integral(0 to dz_new) Iwt(z) dz = dz_tgt
      ! This is trivial where Iwt is a constant, and agrees with the two limits above.
 
      ! Take test values at the transition points to figure out which segment
      ! the new value will be found in.
      if (zr1 >= zd) then
        int_zd = iwtd*(zd - zr1)
        int_zs = int_zd - dint_zs_zd
      elseif (zr1 <= zs) then
        int_zs = (zs - zr1)
        int_zd = dint_zs_zd + (zs - zr1)
      else
!        Int_zd = (zd - zr1) * (Iwtd + 0.5*(1.0 - Iwtd) * (zd - zr1) / (zd - zs))
        int_zd = (zd - zr1) * (iwtd*(0.5*(zd+zr1) - zs) + 0.5*(zd - zr1)) * idzwt
        int_zs = (zs - zr1) * (0.5*iwtd * ((zr1 - zs)) + (zd - 0.5*(zr1+zs))) * idzwt
        ! It has been verified that  Int_zs = Int_zd - dInt_zs_zd to within roundoff.
      endif
 
      if (dz_tgt >= int_zd) then ! The new location is in the deep, slow region.
        dz_g(k) = sgn * ((zd-zr1) + wtd*(dz_tgt - int_zd))
      elseif (dz_tgt <= int_zs) then ! The new location is in the shallow region.
        dz_g(k) = sgn * ((zs-zr1) + (dz_tgt - int_zs))
      else  ! We need to solve a quadratic equation for z_new.
        ! For accuracy, do the integral from the starting depth or the nearest
        ! edge of the transition region.  The results with each choice are
        ! mathematically equivalent, but differ in roundoff, and this choice
        ! should minimize the likelihood of inadvertently overlapping interfaces.
        if (zr1 <= zs) then ; dz0 = zs-zr1 ; z0 = zs ; f0 = dz_tgt - int_zs
        elseif (zr1 >= zd) then ; dz0 = zd-zr1 ; z0 = zd ; f0 = dz_tgt - int_zd
        else ; dz0 = 0.0 ; z0 = zr1 ; f0 = dz_tgt ; endif
 
        bq = (dzwt + 2.0*aq*(z0-zs))
        ! Solve the quadratic: Aq*(zn-z0)**2 + Bq*(zn-z0) - F0*dzwt = 0
        ! Note that b>=0, and the two terms in the standard form cancel for the right root.
        dz_g(k) = sgn * (dz0 + 2.0*f0*dzwt / (bq + sqrt(bq**2 + 4.0*aq*f0*dzwt) ))
 
!       if (debug) then
!         dz0 = zs-zr1 ; z0 = zs ; F0 = dz_tgt - Int_zs ; Bq = (dzwt + 2.0*Aq*(z0-zs))
!         ddz_g_s(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k)
!         dz0 = zd-zr1 ; z0 = zd ; F0 = dz_tgt - Int_zd ; Bq = (dzwt + 2.0*Aq*(z0-zs))
!         ddz_g_d(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k)
!
!         if (abs(ddz_g_s(k)) > 1e-12*(abs(dz_g(k)) + abs(dz_g(k)+ddz_g_s(k)))) &
!           call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled (sc).")
!         if (abs(ddz_g_d(k) - ddz_g_s(k)) > 1e-12*(abs(dz_g(k)+ddz_g_d(k)) + abs(dz_g(k)+ddz_g_s(k)))) &
!           call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled.")
!       endif
      endif
 
    endif
  enddo
 !dz_g(CS%nk+1) = 0.0
 
  if (debug) then
    z_old_k = z_old(1)
    do k=1,cs%nk+1
      if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model
      z_act(k) = z_old_k + dz_g(k)
    enddo
    do k=2,cs%nk+1
      if (sgn*((z_act(k))-z_act(k-1)) < -1e-15*(abs(z_act(k))+abs(z_act(k-1))) ) &
        call mom_error(fatal, "filtered_grid_motion: z_output is tangled.")
    enddo
  endif
 

get_regrid_size()

integer function, public mom_regridding::get_regrid_size

(

type(regridding_cs), intent(inout)

CS

)

Returns the number of levels/layers in the regridding control structure.

Parameters

[in,out]

cs

Regridding control structure

Definition at line 2347 of file MOM_regridding.F90.

  type(regridding_CS), intent(inout) :: CS !< Regridding control structure
 
  get_regrid_size = cs%nk
 

get_rho_cs()

type(rho_cs) function, public mom_regridding::get_rho_cs

(

type(regridding_cs), intent(in)

CS

)

This returns a copy of the rho_CS stored in the regridding control structure.

Parameters

[in]

cs

Regridding control structure

Definition at line 2371 of file MOM_regridding.F90.

  type(regridding_CS), intent(in) :: CS !< Regridding control structure
  type(rho_CS) :: get_rho_CS
 
  get_rho_cs = cs%rho_CS

get_sigma_cs()

type(sigma_cs) function, public mom_regridding::get_sigma_cs

(

type(regridding_cs), intent(in)

CS

)

This returns a copy of the sigma_CS stored in the regridding control structure.

Parameters

[in]

cs

Regridding control structure

Definition at line 2363 of file MOM_regridding.F90.

  type(regridding_CS), intent(in) :: CS !< Regridding control structure
  type(sigma_CS) :: get_sigma_CS
 
  get_sigma_cs = cs%sigma_CS

get_zlike_cs()

type(zlike_cs) function, public mom_regridding::get_zlike_cs

(

type(regridding_cs), intent(in)

CS

)

This returns a copy of the zlike_CS stored in the regridding control structure.

Parameters

[in]

cs

Regridding control structure

Definition at line 2355 of file MOM_regridding.F90.

  type(regridding_CS), intent(in) :: CS !< Regridding control structure
  type(zlike_CS) :: get_zlike_CS
 
  get_zlike_cs = cs%zlike_CS

getcoordinateinterfaces()

real function, dimension(cs%nk+1), public mom_regridding::getcoordinateinterfaces	(	type(regridding_cs), intent(in)	CS,
		logical, intent(in), optional	undo_scaling
	)

Query the target coordinate interface positions.

Parameters

[in]	cs	Regridding control structure
[in]	undo_scaling	If present and true, undo any internal rescaling of the resolution data.

Returns

Interface positions in target coordinate

Definition at line 2140 of file MOM_regridding.F90.

  type(regridding_CS), intent(in) :: CS                      !< Regridding control structure
  logical,   optional, intent(in) :: undo_scaling            !< If present and true, undo any internal
                                                             !! rescaling of the resolution data.
  real, dimension(CS%nk+1)        :: getCoordinateInterfaces !< Interface positions in target coordinate
 
  integer :: k
  logical :: unscale
  unscale = .false. ; if (present(undo_scaling)) unscale = undo_scaling
 
  ! When using a coordinate with target densities, we need to get the actual
  ! densities, rather than computing the interfaces based on resolution
  if (cs%regridding_scheme == regridding_rho) then
    if (.not. cs%target_density_set) &
      call mom_error(fatal, 'MOM_regridding, getCoordinateInterfaces: '//&
                            'target densities not set!')
 
    if (unscale) then
      getcoordinateinterfaces(:) = cs%coord_scale * cs%target_density(:)
    else
      getcoordinateinterfaces(:) = cs%target_density(:)
    endif
  else
    if (unscale) then
      getcoordinateinterfaces(1) = 0.
      do k = 1, cs%nk
        getcoordinateinterfaces(k+1) = getcoordinateinterfaces(k) - &
                                       cs%coord_scale * cs%coordinateResolution(k)
      enddo
    else
      getcoordinateinterfaces(1) = 0.
      do k = 1, cs%nk
        getcoordinateinterfaces(k+1) = getcoordinateinterfaces(k) - &
                                       cs%coordinateResolution(k)
      enddo
    endif
    ! The following line has an "abs()" to allow ferret users to reference
    ! data by index. It is a temporary work around...  :(  -AJA
    getcoordinateinterfaces(:) = abs( getcoordinateinterfaces(:) )
  endif
 

getcoordinateresolution()

real function, dimension(cs%nk), public mom_regridding::getcoordinateresolution	(	type(regridding_cs), intent(in)	CS,
		logical, intent(in), optional	undo_scaling
	)

Query the fixed resolution data.

Parameters

[in]	cs	Regridding control structure
[in]	undo_scaling	If present and true, undo any internal rescaling of the resolution data.

Definition at line 2122 of file MOM_regridding.F90.

  type(regridding_CS), intent(in) :: CS !< Regridding control structure
  logical,   optional, intent(in) :: undo_scaling !< If present and true, undo any internal
                                        !! rescaling of the resolution data.
  real, dimension(CS%nk)          :: getCoordinateResolution
 
  logical :: unscale
  unscale = .false. ; if (present(undo_scaling)) unscale = undo_scaling
 
  if (unscale) then
    getcoordinateresolution(:) = cs%coord_scale * cs%coordinateResolution(:)
  else
    getcoordinateresolution(:) = cs%coordinateResolution(:)
  endif
 

getcoordinateshortname()

character(len=20) function, public mom_regridding::getcoordinateshortname

(

type(regridding_cs), intent(in)

CS

)

Query the short name of the coordinate.

Parameters

[in]

cs

Regridding control structure

Definition at line 2209 of file MOM_regridding.F90.

  type(regridding_CS), intent(in) :: CS !< Regridding control structure
  character(len=20)               :: getCoordinateShortName
 
  select case ( cs%regridding_scheme )
    case ( regridding_zstar )
      !getCoordinateShortName = 'z*'
      ! The following line is a temporary work around...  :(  -AJA
      getcoordinateshortname = 'pseudo-depth, -z*'
    case ( regridding_sigma_shelf_zstar )
      getcoordinateshortname = 'pseudo-depth, -z*/sigma'
    case ( regridding_sigma )
      getcoordinateshortname = 'sigma'
    case ( regridding_rho )
      getcoordinateshortname = 'rho'
    case ( regridding_arbitrary )
      getcoordinateshortname = 'coordinate'
    case ( regridding_hycom1 )
      getcoordinateshortname = 'z-rho'
    case ( regridding_slight )
      getcoordinateshortname = 's-rho'
    case ( regridding_adaptive )
      getcoordinateshortname = 'adaptive'
    case default
      call mom_error(fatal,'MOM_regridding, getCoordinateShortName: '//&
                     'Unknown regridding scheme selected!')
  end select ! type of grid
 

getcoordinateunits()

character(len=20) function, public mom_regridding::getcoordinateunits

(

type(regridding_cs), intent(in)

CS

)

Query the target coordinate units.

Parameters

[in]

cs

Regridding control structure

Definition at line 2185 of file MOM_regridding.F90.

  type(regridding_CS), intent(in) :: CS !< Regridding control structure
  character(len=20)               :: getCoordinateUnits
 
  select case ( cs%regridding_scheme )
    case ( regridding_zstar, regridding_hycom1, regridding_slight, regridding_adaptive )
      getcoordinateunits = 'meter'
    case ( regridding_sigma_shelf_zstar )
      getcoordinateunits = 'meter/fraction'
    case ( regridding_sigma )
      getcoordinateunits = 'fraction'
    case ( regridding_rho )
      getcoordinateunits = 'kg/m3'
    case ( regridding_arbitrary )
      getcoordinateunits = 'unknown'
    case default
      call mom_error(fatal,'MOM_regridding, getCoordinateUnits: '//&
                     'Unknown regridding scheme selected!')
  end select ! type of grid
 

getstaticthickness()

real function, dimension(cs%nk), public mom_regridding::getstaticthickness	(	type(regridding_cs), intent(in)	CS,
		real, intent(in)	SSH,
		real, intent(in)	depth
	)

Return coordinate-derived thicknesses for fixed coordinate systems.

Parameters

[in]	cs	Regridding control structure
[in]	ssh	The sea surface height, in the same units as depth
[in]	depth	The maximum depth of the grid, often [Z ~> m]

Returns

The returned thicknesses in the units of depth

Definition at line 2380 of file MOM_regridding.F90.

  type(regridding_CS), intent(in) :: CS !< Regridding control structure
  real,                intent(in) :: SSH   !< The sea surface height, in the same units as depth
  real,                intent(in) :: depth !< The maximum depth of the grid, often [Z ~> m]
  real, dimension(CS%nk)          :: getStaticThickness !< The returned thicknesses in the units of depth
  ! Local
  integer :: k
  real :: z, dz
 
  select case ( cs%regridding_scheme )
    case ( regridding_zstar, regridding_sigma_shelf_zstar, regridding_hycom1, regridding_slight, regridding_adaptive )
      if (depth>0.) then
        z = ssh
        do k = 1, cs%nk
          dz = cs%coordinateResolution(k) * ( 1. + ssh/depth ) ! Nominal dz*
          dz = max(dz, 0.)              ! Avoid negative incase ssh=-depth
          dz = min(dz, depth - z)       ! Clip if below topography
          z = z + dz                    ! Bottom of layer
          getstaticthickness(k) = dz
        enddo
      else
        getstaticthickness(:) = 0. ! On land ...
      endif
    case ( regridding_sigma )
      getstaticthickness(:) = cs%coordinateResolution(:) * ( depth + ssh )
    case ( regridding_rho )
      getstaticthickness(:) = 0. ! Not applicable
    case ( regridding_arbitrary )
      getstaticthickness(:) = 0.  ! Not applicable
    case default
      call mom_error(fatal,'MOM_regridding, getStaticThickness: '//&
                     'Unknown regridding scheme selected!')
  end select ! type of grid
 

inflate_vanished_layers_old()

subroutine, public mom_regridding::inflate_vanished_layers_old	(	type(regridding_cs), intent(in)	CS,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		real, dimension(szi_(g),szj_(g), szk_(gv)), intent(inout)	h
	)

Parameters

[in]	cs	Regridding control structure
[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in,out]	h	Layer thicknesses [H ~> m or kg m-2]

Definition at line 1828 of file MOM_regridding.F90.

!------------------------------------------------------------------------------
! This routine is called when initializing the regridding options. The
! objective is to make sure all layers are at least as thick as the minimum
! thickness allowed for regridding purposes (this parameter is set in the
! MOM_input file or defaulted to 1.0e-3). When layers are too thin, they
! are inflated up to the minmum thickness.
!------------------------------------------------------------------------------
 
  ! Arguments
  type(regridding_CS),                    intent(in)    :: CS   !< Regridding control structure
  type(ocean_grid_type),                  intent(in)    :: G    !< The ocean's grid structure
  type(verticalGrid_type),                intent(in)    :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h    !< Layer thicknesses [H ~> m or kg m-2]
 
  ! Local variables
  integer :: i, j, k
  real    :: hTmp(GV%ke)
 
  do i = g%isc-1,g%iec+1
    do j = g%jsc-1,g%jec+1
 
      ! Build grid for current column
      do k = 1,gv%ke
        htmp(k) = h(i,j,k)
      enddo
 
      call old_inflate_layers_1d( cs%min_thickness, gv%ke, htmp )
 
      ! Save modified grid
      do k = 1,gv%ke
        h(i,j,k) = htmp(k)
      enddo
 
    enddo
  enddo
 

initcoord()

subroutine mom_regridding::initcoord ( type(regridding_cs), intent(inout)  CS, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, character(len=*), intent(in)  coord_mode  )

private

Initialize the coordinate resolutions by calling the appropriate initialization routine for the specified coordinate mode.

Parameters

[in,out]	cs	Regridding control structure
[in]	coord_mode	A string indicating the coordinate mode. See the documentation for regrid_consts for the recognized values.
[in]	gv	Ocean vertical grid structure
[in]	us	A dimensional unit scaling type

Definition at line 1972 of file MOM_regridding.F90.

  type(regridding_CS),     intent(inout) :: CS !< Regridding control structure
  character(len=*),        intent(in)    :: coord_mode !< A string indicating the coordinate mode.
                                               !! See the documentation for regrid_consts
                                               !! for the recognized values.
  type(verticalGrid_type), intent(in)    :: GV !< Ocean vertical grid structure
  type(unit_scale_type),   intent(in)    :: US  !< A dimensional unit scaling type
 
  select case (coordinatemode(coord_mode))
  case (regridding_zstar)
    call init_coord_zlike(cs%zlike_CS, cs%nk, cs%coordinateResolution)
  case (regridding_sigma_shelf_zstar)
    call init_coord_zlike(cs%zlike_CS, cs%nk, cs%coordinateResolution)
  case (regridding_sigma)
    call init_coord_sigma(cs%sigma_CS, cs%nk, cs%coordinateResolution)
  case (regridding_rho)
    call init_coord_rho(cs%rho_CS, cs%nk, cs%ref_pressure, cs%target_density, cs%interp_CS)
  case (regridding_hycom1)
    call init_coord_hycom(cs%hycom_CS, cs%nk, cs%coordinateResolution, cs%target_density, &
                          cs%interp_CS)
  case (regridding_slight)
    call init_coord_slight(cs%slight_CS, cs%nk, cs%ref_pressure, cs%target_density, &
                           cs%interp_CS, gv%m_to_H)
  case (regridding_adaptive)
    call init_coord_adapt(cs%adapt_CS, cs%nk, cs%coordinateResolution, gv%m_to_H, us%kg_m3_to_R)
  end select

initialize_regridding()

subroutine, public mom_regridding::initialize_regridding	(	type(regridding_cs), intent(inout)	CS,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		real, intent(in)	max_depth,
		type(param_file_type), intent(in)	param_file,
		character(len=*), intent(in)	mdl,
		character(len=*), intent(in)	coord_mode,
		character(len=*), intent(in)	param_prefix,
		character(len=*), intent(in)	param_suffix
	)

Initialization and configures a regridding control structure based on customizable run-time parameters.

Parameters

[in,out]	cs	Regridding control structure
[in]	gv	Ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	max_depth	The maximum depth of the ocean [Z ~> m].
[in]	param_file	Parameter file
[in]	mdl	Name of calling module.
[in]	coord_mode	Coordinate mode
[in]	param_prefix	String to prefix to parameter names. If empty, causes main model parameters to be used.
[in]	param_suffix	String to append to parameter names.

Definition at line 180 of file MOM_regridding.F90.

  type(regridding_CS),        intent(inout) :: CS  !< Regridding control structure
  type(verticalGrid_type),    intent(in)    :: GV  !< Ocean vertical grid structure
  type(unit_scale_type),      intent(in)    :: US  !< A dimensional unit scaling type
  real,                       intent(in)    :: max_depth  !< The maximum depth of the ocean [Z ~> m].
  type(param_file_type),      intent(in)    :: param_file !< Parameter file
  character(len=*),           intent(in)    :: mdl        !< Name of calling module.
  character(len=*),           intent(in)    :: coord_mode !< Coordinate mode
  character(len=*),           intent(in)    :: param_prefix !< String to prefix to parameter names.
                                                            !! If empty, causes main model parameters to be used.
  character(len=*),           intent(in)    :: param_suffix !< String to append to parameter names.
 
  ! Local variables
  integer :: ke ! Number of levels
  character(len=80)  :: string, string2, varName ! Temporary strings
  character(len=40)  :: coord_units, param_name, coord_res_param ! Temporary strings
  character(len=200) :: inputdir, fileName
  character(len=320) :: message ! Temporary strings
  character(len=12) :: expected_units ! Temporary strings
  logical :: tmpLogical, fix_haloclines, set_max, do_sum, main_parameters
  logical :: coord_is_state_dependent, ierr
  logical :: default_2018_answers, remap_answers_2018
  real :: filt_len, strat_tol, index_scale, tmpReal, P_Ref
  real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z).
  real :: dz_fixed_sfc, Rho_avg_depth, nlay_sfc_int
  real :: adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha
  real :: adaptDrho0 ! Reference density difference for stratification-dependent diffusion. [R ~> kg m-3]
  integer :: nz_fixed_sfc, k, nzf(4)
  real, dimension(:), allocatable :: dz     ! Resolution (thickness) in units of coordinate, which may be [m]
                                            ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units.
  real, dimension(:), allocatable :: h_max  ! Maximum layer thicknesses [H ~> m or kg m-2]
  real, dimension(:), allocatable :: z_max  ! Maximum interface depths [H ~> m or kg m-2] or other
                                            ! units depending on the coordinate
  real, dimension(:), allocatable :: dz_max ! Thicknesses used to find maximum interface depths
                                            ! [H ~> m or kg m-2] or other units
  real, dimension(:), allocatable :: rho_target ! Target density used in HYBRID mode [kg m-3]
  ! Thicknesses [m] that give level centers corresponding to table 2 of WOA09
  real, dimension(40) :: woa09_dz = (/ 5.,  10.,  10.,  15.,  22.5, 25., 25.,  25.,  &
                                      37.5, 50.,  50.,  75., 100., 100., 100., 100., &
                                     100., 100., 100., 100., 100., 100., 100., 175., &
                                     250., 375., 500., 500., 500., 500., 500., 500., &
                                     500., 500., 500., 500., 500., 500., 500., 500. /)
 
  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
 
  main_parameters=.false.
  if (len_trim(param_prefix)==0) main_parameters=.true.
  if (main_parameters .and. len_trim(param_suffix)>0) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
              'Suffix provided without prefix for parameter names!')
 
  cs%nk = 0
  cs%regridding_scheme = coordinatemode(coord_mode)
  coord_is_state_dependent = state_dependent(coord_mode)
  maximum_depth = us%Z_to_m*max_depth
 
  if (main_parameters) then
    ! Read coordinate units parameter (main model = REGRIDDING_COORDINATE_UNITS)
    call get_param(param_file, mdl, "REGRIDDING_COORDINATE_UNITS", coord_units, &
                 "Units of the regridding coordinate.", default=coordinateunits(coord_mode))
  else
    coord_units=coordinateunits(coord_mode)
  endif
 
  if (coord_is_state_dependent) then
    if (main_parameters) then
      param_name = "INTERPOLATION_SCHEME"
      string2 = regriddingdefaultinterpscheme
    else
      param_name = trim(param_prefix)//"_INTERP_SCHEME_"//trim(param_suffix)
      string2 = 'PPM_H4' ! Default for diagnostics
    endif
    call get_param(param_file, mdl, "INTERPOLATION_SCHEME", string, &
                 "This sets the interpolation scheme to use to "//&
                 "determine the new grid. These parameters are "//&
                 "only relevant when REGRIDDING_COORDINATE_MODE is "//&
                 "set to a function of state. Otherwise, it is not "//&
                 "used. It can be one of the following schemes: "//&
                 trim(regriddinginterpschemedoc), default=trim(string2))
    call set_regrid_params(cs, interp_scheme=string)
 
    call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                 "This sets the default value for the various _2018_ANSWERS parameters.", &
                 default=.false.)
    call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", remap_answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the "//&
                 "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                 "forms of the same expressions.", default=default_2018_answers)
    call set_regrid_params(cs, remap_answers_2018=remap_answers_2018)
  endif
 
  if (main_parameters .and. coord_is_state_dependent) then
    call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", tmplogical, &
                 "When defined, a proper high-order reconstruction "//&
                 "scheme is used within boundary cells rather "//&
                 "than PCM. E.g., if PPM is used for remapping, a "//&
                 "PPM reconstruction will also be used within "//&
                 "boundary cells.", default=regriddingdefaultboundaryextrapolation)
    call set_regrid_params(cs, boundary_extrapolation=tmplogical)
  else
    call set_regrid_params(cs, boundary_extrapolation=.false.)
  endif
 
  ! Read coordinate configuration parameter (main model = ALE_COORDINATE_CONFIG)
  if (main_parameters) then
    param_name = "ALE_COORDINATE_CONFIG"
    coord_res_param = "ALE_RESOLUTION"
    string2 = 'UNIFORM'
  else
    param_name = trim(param_prefix)//"_DEF_"//trim(param_suffix)
    coord_res_param = trim(param_prefix)//"_RES_"//trim(param_suffix)
    string2 = 'UNIFORM'
    if (maximum_depth>3000.) string2='WOA09' ! For convenience
  endif
  call get_param(param_file, mdl, param_name, string, &
                 "Determines how to specify the coordinate "//&
                 "resolution. Valid options are:\n"//&
                 " PARAM       - use the vector-parameter "//trim(coord_res_param)//"\n"//&
                 " UNIFORM[:N] - uniformly distributed\n"//&
                 " FILE:string - read from a file. The string specifies\n"//&
                 "               the filename and variable name, separated\n"//&
                 "               by a comma or space, e.g. FILE:lev.nc,dz\n"//&
                 "               or FILE:lev.nc,interfaces=zw\n"//&
                 " WOA09[:N]   - the WOA09 vertical grid (approximately)\n"//&
                 " FNC1:string - FNC1:dz_min,H_total,power,precision\n"//&
                 " HYBRID:string - read from a file. The string specifies\n"//&
                 "               the filename and two variable names, separated\n"//&
                 "               by a comma or space, for sigma-2 and dz. e.g.\n"//&
                 "               HYBRID:vgrid.nc,sigma2,dz",&
                 default=trim(string2))
  message = "The distribution of vertical resolution for the target\n"//&
            "grid used for Eulerian-like coordinates. For example,\n"//&
            "in z-coordinate mode, the parameter is a list of level\n"//&
            "thicknesses (in m). In sigma-coordinate mode, the list\n"//&
            "is of non-dimensional fractions of the water column."
  if (index(trim(string),'UNIFORM')==1) then
    if (len_trim(string)==7) then
      ke = gv%ke ! Use model nk by default
      tmpreal = maximum_depth
    elseif (index(trim(string),'UNIFORM:')==1 .and. len_trim(string)>8) then
      ! Format is "UNIFORM:N" or "UNIFORM:N,dz"
      ke = extract_integer(string(9:len_trim(string)),'',1)
      tmpreal = extract_real(string(9:len_trim(string)),',',2,missing_value=maximum_depth)
    else
      call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
          'Unable to interpret "'//trim(string)//'".')
    endif
    allocate(dz(ke))
    dz(:) = uniformresolution(ke, coord_mode, tmpreal, &
                              us%R_to_kg_m3*(gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(min(2,ke)))), &
                              us%R_to_kg_m3*(gv%Rlay(ke) + 0.5*(gv%Rlay(ke)-gv%Rlay(max(ke-1,1)))) )
    if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
                   trim(message), units=trim(coord_units))
  elseif (trim(string)=='PARAM') then
    ! Read coordinate resolution (main model = ALE_RESOLUTION)
    ke = gv%ke ! Use model nk by default
    allocate(dz(ke))
    call get_param(param_file, mdl, coord_res_param, dz, &
                   trim(message), units=trim(coord_units), fail_if_missing=.true.)
  elseif (index(trim(string),'FILE:')==1) then
    ! FILE:filename,var_name is assumed to be reading level thickness variables
    ! FILE:filename,interfaces=var_name reads positions
    if (string(6:6)=='.' .or. string(6:6)=='/') then
      ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
      filename = trim( extractword(trim(string(6:80)), 1) )
    else
      ! Otherwise assume we should look for the file in INPUTDIR
      filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
    endif
    if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
            "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
 
    varname = trim( extractword(trim(string(6:)), 2) )
    if (len_trim(varname)==0) then
      if (field_exists(filename,'dz')) then; varname = 'dz'
      elseif (field_exists(filename,'dsigma')) then; varname = 'dsigma'
      elseif (field_exists(filename,'ztest')) then; varname = 'ztest'
      else ;  call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
                    "Coordinate variable not specified and none could be guessed.")
      endif
    endif
    ! This check fails when the variable is a dimension variable! -AJA
   !if (.not. field_exists(fileName,trim(varName))) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
   !             "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
    if (cs%regridding_scheme == regridding_sigma) then
      expected_units = 'nondim'
    elseif (cs%regridding_scheme == regridding_rho) then
      expected_units = 'kg m-3'
    else
      expected_units = 'meters'
    endif
    if (index(trim(varname),'interfaces=')==1) then
      varname=trim(varname(12:))
      call check_grid_def(filename, varname, expected_units, message, ierr)
      if (ierr) call mom_error(fatal,trim(mdl)//", initialize_regridding: "//&
                  "Unsupported format in grid definition '"//trim(filename)//"'. Error message "//trim(message))
      call field_size(trim(filename), trim(varname), nzf)
      ke = nzf(1)-1
      if (ke < 1) call mom_error(fatal, trim(mdl)//" initialize_regridding via Var "//trim(varname)//&
                                 "in FILE "//trim(filename)//" requires at least 2 target interface values.")
      if (cs%regridding_scheme == regridding_rho) then
        allocate(rho_target(ke+1))
        call mom_read_data(trim(filename), trim(varname), rho_target)
      else
        allocate(dz(ke))
        allocate(z_max(ke+1))
        call mom_read_data(trim(filename), trim(varname), z_max)
        dz(:) = abs(z_max(1:ke) - z_max(2:ke+1))
        deallocate(z_max)
      endif
    else
      ! Assume reading resolution
      call field_size(trim(filename), trim(varname), nzf)
      ke = nzf(1)
      allocate(dz(ke))
      call mom_read_data(trim(filename), trim(varname), dz)
    endif
    if (main_parameters .and. (ke/=gv%ke)) then
      call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
                 'Mismatch in number of model levels and "'//trim(string)//'".')
    endif
    if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
               trim(message), units=coordinateunits(coord_mode))
  elseif (index(trim(string),'FNC1:')==1) then
    ke = gv%ke; allocate(dz(ke))
    call dz_function1( trim(string(6:)), dz )
    if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
               trim(message), units=coordinateunits(coord_mode))
  elseif (index(trim(string),'RFNC1:')==1) then
    ! Function used for set target interface densities
    ke = rho_function1( trim(string(7:)), rho_target )
  elseif (index(trim(string),'HYBRID:')==1) then
    ke = gv%ke; allocate(dz(ke))
    ! The following assumes the FILE: syntax of above but without "FILE:" in the string
    allocate(rho_target(ke+1))
    filename = trim( extractword(trim(string(8:)), 1) )
    if (filename(1:1)/='.' .and. filename(1:1)/='/') filename = trim(inputdir) // trim( filename )
    if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: HYBRID "// &
      "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
    varname = trim( extractword(trim(string(8:)), 2) )
    if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)//", initialize_regridding: HYBRID "// &
      "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
    call mom_read_data(trim(filename), trim(varname), rho_target)
    varname = trim( extractword(trim(string(8:)), 3) )
    if (varname(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz
      call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz )
    else ! Read dz from file
      if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)//", initialize_regridding: HYBRID "// &
        "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
      call mom_read_data(trim(filename), trim(varname), dz)
    endif
    if (main_parameters) then
      call log_param(param_file, mdl, "!"//coord_res_param, dz, &
               trim(message), units=coordinateunits(coord_mode))
      call log_param(param_file, mdl, "!TARGET_DENSITIES", rho_target, &
               'HYBRID target densities for interfaces', units=coordinateunits(coord_mode))
    endif
  elseif (index(trim(string),'WOA09')==1) then
    if (len_trim(string)==5) then
      tmpreal = 0. ; ke = 0
      do while (tmpreal<maximum_depth)
        ke = ke + 1
        tmpreal = tmpreal + woa09_dz(ke)
      enddo
    elseif (index(trim(string),'WOA09:')==1) then
      if (len_trim(string)==6) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
                 'Expected string of form "WOA09:N" but got "'//trim(string)//'".')
      ke = extract_integer(string(7:len_trim(string)),'',1)
    endif
    if (ke>40 .or. ke<1) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
                 'For "WOA05:N" N must 0<N<41 but got "'//trim(string)//'".')
    allocate(dz(ke))
    dz(1:ke) = woa09_dz(1:ke)
  else
    call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
      "Unrecognized coordinate configuration"//trim(string))
  endif
 
  if (main_parameters) then
    ! This is a work around to apparently needed to work with the from_Z initialization...  ???
    if (coordinatemode(coord_mode) == regridding_zstar .or. &
        coordinatemode(coord_mode) == regridding_hycom1 .or. &
        coordinatemode(coord_mode) == regridding_slight .or. &
        coordinatemode(coord_mode) == regridding_adaptive) then
      ! Adjust target grid to be consistent with maximum_depth
      tmpreal = sum( dz(:) )
      if (tmpreal < maximum_depth) then
        dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
      elseif (tmpreal > maximum_depth) then
        if ( dz(ke) + ( maximum_depth - tmpreal ) > 0. ) then
          dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
        else
          call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
            "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string))
        endif
      endif
    endif
  endif
 
  cs%nk=ke
 
  ! Target resolution (for fixed coordinates)
  allocate( cs%coordinateResolution(cs%nk) ); cs%coordinateResolution(:) = -1.e30
  if (state_dependent(cs%regridding_scheme)) then
    ! Target values
    allocate( cs%target_density(cs%nk+1) ); cs%target_density(:) = -1.e30*us%kg_m3_to_R
  endif
 
  if (allocated(dz)) then
    if (coordinatemode(coord_mode) == regridding_sigma) then
      call setcoordinateresolution(dz, cs, scale=1.0)
    elseif (coordinatemode(coord_mode) == regridding_rho) then
      call setcoordinateresolution(dz, cs, scale=us%kg_m3_to_R)
      cs%coord_scale = us%R_to_kg_m3
    elseif (coordinatemode(coord_mode) == regridding_adaptive) then
      call setcoordinateresolution(dz, cs, scale=gv%m_to_H)
      cs%coord_scale = gv%H_to_m
    else
      call setcoordinateresolution(dz, cs, scale=us%m_to_Z)
      cs%coord_scale = us%Z_to_m
    endif
  endif
 
  if (allocated(rho_target)) then
    call set_target_densities(cs, us%kg_m3_to_R*rho_target)
    deallocate(rho_target)
 
  ! \todo This line looks like it would overwrite the target densities set just above?
  elseif (coordinatemode(coord_mode) == regridding_rho) then
    call set_target_densities_from_gv(gv, us, cs)
    call log_param(param_file, mdl, "!TARGET_DENSITIES", us%R_to_kg_m3*cs%target_density(:), &
             'RHO target densities for interfaces', units=coordinateunits(coord_mode))
  endif
 
  ! initialise coordinate-specific control structure
  call initcoord(cs, gv, us, coord_mode)
 
  if (main_parameters .and. coord_is_state_dependent) then
    call get_param(param_file, mdl, "P_REF", p_ref, &
                 "The pressure that is used for calculating the coordinate "//&
                 "density.  (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
                 "This is only used if USE_EOS and ENABLE_THERMODYNAMICS are true.", &
                 units="Pa", default=2.0e7, scale=us%kg_m3_to_R*us%m_s_to_L_T**2)
    call get_param(param_file, mdl, "REGRID_COMPRESSIBILITY_FRACTION", tmpreal, &
                 "When interpolating potential density profiles we can add "//&
                 "some artificial compressibility solely to make homogeneous "//&
                 "regions appear stratified.", units="nondim", default=0.)
    call set_regrid_params(cs, compress_fraction=tmpreal, ref_pressure=p_ref)
  endif
 
  if (main_parameters) then
    call get_param(param_file, mdl, "MIN_THICKNESS", tmpreal, &
                 "When regridding, this is the minimum layer "//&
                 "thickness allowed.", units="m", scale=gv%m_to_H, &
                 default=regriddingdefaultminthickness )
    call set_regrid_params(cs, min_thickness=tmpreal)
  else
    call set_regrid_params(cs, min_thickness=0.)
  endif
 
  if (coordinatemode(coord_mode) == regridding_slight) then
    ! Set SLight-specific regridding parameters.
    call get_param(param_file, mdl, "SLIGHT_DZ_SURFACE", dz_fixed_sfc, &
                 "The nominal thickness of fixed thickness near-surface "//&
                 "layers with the SLight coordinate.", units="m", default=1.0, scale=gv%m_to_H)
    call get_param(param_file, mdl, "SLIGHT_NZ_SURFACE_FIXED", nz_fixed_sfc, &
                 "The number of fixed-depth surface layers with the SLight "//&
                 "coordinate.", units="nondimensional", default=2)
    call get_param(param_file, mdl, "SLIGHT_SURFACE_AVG_DEPTH", rho_avg_depth, &
                 "The thickness of the surface region over which to average "//&
                 "when calculating the density to use to define the interior "//&
                 "with the SLight coordinate.", units="m", default=1.0, scale=gv%m_to_H)
    call get_param(param_file, mdl, "SLIGHT_NLAY_TO_INTERIOR", nlay_sfc_int, &
                 "The number of layers to offset the surface density when "//&
                 "defining where the interior ocean starts with SLight.", &
                 units="nondimensional", default=2.0)
    call get_param(param_file, mdl, "SLIGHT_FIX_HALOCLINES", fix_haloclines, &
                 "If true, identify regions above the reference pressure "//&
                 "where the reference pressure systematically underestimates "//&
                 "the stratification and use this in the definition of the "//&
                 "interior with the SLight coordinate.", default=.false.)
 
    call set_regrid_params(cs, dz_min_surface=dz_fixed_sfc, &
                nz_fixed_surface=nz_fixed_sfc, rho_ml_avg_depth=rho_avg_depth, &
                nlay_ml_to_interior=nlay_sfc_int, fix_haloclines=fix_haloclines)
    if (fix_haloclines) then
      ! Set additional parameters related to SLIGHT_FIX_HALOCLINES.
      call get_param(param_file, mdl, "HALOCLINE_FILTER_LENGTH", filt_len, &
                 "A length scale over which to smooth the temperature and "//&
                 "salinity before identifying erroneously unstable haloclines.", &
                 units="m", default=2.0, scale=gv%m_to_H)
      call get_param(param_file, mdl, "HALOCLINE_STRAT_TOL", strat_tol, &
                 "A tolerance for the ratio of the stratification of the "//&
                 "apparent coordinate stratification to the actual value "//&
                 "that is used to identify erroneously unstable haloclines. "//&
                 "This ratio is 1 when they are equal, and sensible values "//&
                 "are between 0 and 0.5.", units="nondimensional", default=0.2)
      call set_regrid_params(cs, halocline_filt_len=filt_len, &
                             halocline_strat_tol=strat_tol)
    endif
 
  endif
 
  if (coordinatemode(coord_mode) == regridding_adaptive) then
    call get_param(param_file, mdl, "ADAPT_TIME_RATIO", adapttimeratio, &
                 "Ratio of ALE timestep to grid timescale.", units="nondim", default=1.0e-1)
    call get_param(param_file, mdl, "ADAPT_ZOOM_DEPTH", adaptzoom, &
                 "Depth of near-surface zooming region.", units="m", default=200.0, scale=gv%m_to_H)
    call get_param(param_file, mdl, "ADAPT_ZOOM_COEFF", adaptzoomcoeff, &
                 "Coefficient of near-surface zooming diffusivity.", units="nondim", default=0.2)
    call get_param(param_file, mdl, "ADAPT_BUOY_COEFF", adaptbuoycoeff, &
                 "Coefficient of buoyancy diffusivity.", units="nondim", default=0.8)
    call get_param(param_file, mdl, "ADAPT_ALPHA", adaptalpha, &
                 "Scaling on optimization tendency.", units="nondim", default=1.0)
    call get_param(param_file, mdl, "ADAPT_DO_MIN_DEPTH", tmplogical, &
                 "If true, make a HyCOM-like mixed layer by preventing interfaces "//&
                 "from being shallower than the depths specified by the regridding coordinate.", &
                 default=.false.)
    call get_param(param_file, mdl, "ADAPT_DRHO0", adaptdrho0, &
                 "Reference density difference for stratification-dependent diffusion.", &
                 units="kg m-3", default=0.5, scale=us%kg_m3_to_R)
 
    call set_regrid_params(cs, adapttimeratio=adapttimeratio, adaptzoom=adaptzoom, &
         adaptzoomcoeff=adaptzoomcoeff, adaptbuoycoeff=adaptbuoycoeff, adaptalpha=adaptalpha, &
         adaptdomin=tmplogical, adaptdrho0=adaptdrho0)
  endif
 
  if (main_parameters .and. coord_is_state_dependent) then
    call get_param(param_file, mdl, "MAXIMUM_INT_DEPTH_CONFIG", string, &
                 "Determines how to specify the maximum interface depths.\n"//&
                 "Valid options are:\n"//&
                 " NONE        - there are no maximum interface depths\n"//&
                 " PARAM       - use the vector-parameter MAXIMUM_INTERFACE_DEPTHS\n"//&
                 " FILE:string - read from a file. The string specifies\n"//&
                 "               the filename and variable name, separated\n"//&
                 "               by a comma or space, e.g. FILE:lev.nc,Z\n"//&
                 " FNC1:string - FNC1:dz_min,H_total,power,precision",&
                 default='NONE')
    message = "The list of maximum depths for each interface."
    allocate(z_max(ke+1))
    allocate(dz_max(ke))
    if ( trim(string) == "NONE") then
      ! Do nothing.
    elseif ( trim(string) ==  "PARAM") then
      call get_param(param_file, mdl, "MAXIMUM_INTERFACE_DEPTHS", z_max, &
                   trim(message), units="m", scale=gv%m_to_H, fail_if_missing=.true.)
      call set_regrid_max_depths(cs, z_max)
    elseif (index(trim(string),'FILE:')==1) then
      if (string(6:6)=='.' .or. string(6:6)=='/') then
        ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
        filename = trim( extractword(trim(string(6:80)), 1) )
      else
        ! Otherwise assume we should look for the file in INPUTDIR
        filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
      endif
      if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
        "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
 
      do_sum = .false.
      varname = trim( extractword(trim(string(6:)), 2) )
      if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
        "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
      if (len_trim(varname)==0) then
        if (field_exists(filename,'z_max')) then; varname = 'z_max'
        elseif (field_exists(filename,'dz')) then; varname = 'dz' ; do_sum = .true.
        elseif (field_exists(filename,'dz_max')) then; varname = 'dz_max' ; do_sum = .true.
        else ; call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
          "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
        endif
      endif
      if (do_sum) then
        call mom_read_data(trim(filename), trim(varname), dz_max)
        z_max(1) = 0.0 ; do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ; enddo
      else
        call mom_read_data(trim(filename), trim(varname), z_max)
      endif
      call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, &
                 trim(message), units=coordinateunits(coord_mode))
      call set_regrid_max_depths(cs, z_max, gv%m_to_H)
    elseif (index(trim(string),'FNC1:')==1) then
      call dz_function1( trim(string(6:)), dz_max )
      if ((coordinatemode(coord_mode) == regridding_slight) .and. &
          (dz_fixed_sfc > 0.0)) then
        do k=1,nz_fixed_sfc ; dz_max(k) = dz_fixed_sfc ; enddo
      endif
      z_max(1) = 0.0 ; do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ; enddo
      call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, &
                 trim(message), units=coordinateunits(coord_mode))
      call set_regrid_max_depths(cs, z_max, gv%m_to_H)
    else
      call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
        "Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string))
    endif
    deallocate(z_max)
    deallocate(dz_max)
 
    ! Optionally specify maximum thicknesses for each layer, enforced by moving
    ! the interface below a layer downward.
    call get_param(param_file, mdl, "MAX_LAYER_THICKNESS_CONFIG", string, &
                   "Determines how to specify the maximum layer thicknesses.\n"//&
                   "Valid options are:\n"//&
                   " NONE        - there are no maximum layer thicknesses\n"//&
                   " PARAM       - use the vector-parameter MAX_LAYER_THICKNESS\n"//&
                   " FILE:string - read from a file. The string specifies\n"//&
                   "               the filename and variable name, separated\n"//&
                   "               by a comma or space, e.g. FILE:lev.nc,Z\n"//&
                   " FNC1:string - FNC1:dz_min,H_total,power,precision",&
                   default='NONE')
    message = "The list of maximum thickness for each layer."
    allocate(h_max(ke))
    if ( trim(string) == "NONE") then
      ! Do nothing.
    elseif ( trim(string) ==  "PARAM") then
      call get_param(param_file, mdl, "MAX_LAYER_THICKNESS", h_max, &
                   trim(message), units="m", fail_if_missing=.true., scale=gv%m_to_H)
      call set_regrid_max_thickness(cs, h_max)
    elseif (index(trim(string),'FILE:')==1) then
      if (string(6:6)=='.' .or. string(6:6)=='/') then
        ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
        filename = trim( extractword(trim(string(6:80)), 1) )
      else
        ! Otherwise assume we should look for the file in INPUTDIR
        filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
      endif
      if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
        "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
 
      varname = trim( extractword(trim(string(6:)), 2) )
      if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
        "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
      if (len_trim(varname)==0) then
        if (field_exists(filename,'h_max')) then; varname = 'h_max'
        elseif (field_exists(filename,'dz_max')) then; varname = 'dz_max'
        else ; call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
          "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
        endif
      endif
      call mom_read_data(trim(filename), trim(varname), h_max)
      call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
                 trim(message), units=coordinateunits(coord_mode))
      call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
    elseif (index(trim(string),'FNC1:')==1) then
      call dz_function1( trim(string(6:)), h_max )
      call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
                 trim(message), units=coordinateunits(coord_mode))
      call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
    else
      call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
        "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(string))
    endif
    deallocate(h_max)
  endif
 
  if (allocated(dz)) deallocate(dz)

regridding_main()

subroutine, public mom_regridding::regridding_main	(	type(remapping_cs), intent(in)	remapCS,
		type(regridding_cs), intent(in)	CS,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(inout)	h,
		type(thermo_var_ptrs), intent(inout)	tv,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, cs%nk), intent(inout)	h_new,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, cs%nk+1), intent(inout)	dzInterface,
		real, dimension(:,:), optional, pointer	frac_shelf_h,
		logical, intent(in), optional	conv_adjust
	)

Dispatching regridding routine for orchestrating regridding & remapping.

Parameters

[in]	remapcs	Remapping parameters and options
[in]	cs	Regridding control structure
[in]	g	Ocean grid structure
[in]	gv	Ocean vertical grid structure
[in,out]	h	Current 3D grid obtained after the last time step
[in,out]	tv	Thermodynamical variables (T, S, ...)
[in,out]	h_new	New 3D grid consistent with target coordinate
[in,out]	dzinterface	The change in position of each interface
	frac_shelf_h	Fractional ice shelf coverage
[in]	conv_adjust	If true, do convective adjustment

Definition at line 812 of file MOM_regridding.F90.

!------------------------------------------------------------------------------
! This routine takes care of (1) building a new grid and (2) remapping between
! the old grid and the new grid. The creation of the new grid can be based
! on z coordinates, target interface densities, sigma coordinates or any
! arbitrary coordinate system.
!   The MOM6 interface positions are always calculated from the bottom up by
! accumulating the layer thicknesses starting at z=-G%bathyT.  z increases
! upwards (decreasing k-index).
!   The new grid is defined by the change in position of those interfaces in z
!       dzInterface = zNew - zOld.
!   Thus, if the regridding inflates the top layer, hNew(1) > hOld(1), then the
! second interface moves downward, zNew(2) < zOld(2), and dzInterface(2) < 0.
!       hNew(k) = hOld(k) - dzInterface(k+1) + dzInterface(k)
! IMPORTANT NOTE:
!   This is the converse of the sign convention used in the remapping code!
!------------------------------------------------------------------------------
 
  ! Arguments
  type(remapping_CS),                         intent(in)    :: remapCS !< Remapping parameters and options
  type(regridding_CS),                        intent(in)    :: CS     !< Regridding control structure
  type(ocean_grid_type),                      intent(in)    :: G      !< Ocean grid structure
  type(verticalGrid_type),                    intent(in)    :: GV     !< Ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h      !< Current 3D grid obtained after
                                                                      !! the last time step
  type(thermo_var_ptrs),                      intent(inout) :: tv     !< Thermodynamical variables (T, S, ...)
  real, dimension(SZI_(G),SZJ_(G), CS%nk),    intent(inout) :: h_new  !< New 3D grid consistent with target coordinate
  real, dimension(SZI_(G),SZJ_(G), CS%nk+1),  intent(inout) :: dzInterface !< The change in position of each interface
  real, dimension(:,:),             optional, pointer       :: frac_shelf_h !< Fractional ice shelf coverage
  logical,                          optional, intent(in   ) :: conv_adjust !< If true, do convective adjustment
  ! Local variables
  real :: trickGnuCompiler
  logical :: use_ice_shelf
  logical :: do_convective_adjustment
 
  do_convective_adjustment = .true.
  if (present(conv_adjust)) do_convective_adjustment = conv_adjust
 
  use_ice_shelf = .false.
  if (present(frac_shelf_h)) then
    if (associated(frac_shelf_h)) use_ice_shelf = .true.
  endif
 
  select case ( cs%regridding_scheme )
 
    case ( regridding_zstar )
      if (use_ice_shelf) then
        call build_zstar_grid( cs, g, gv, h, dzinterface, frac_shelf_h )
      else
        call build_zstar_grid( cs, g, gv, h, dzinterface )
      endif
      call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
 
    case ( regridding_sigma_shelf_zstar)
      call build_zstar_grid( cs, g, gv, h, dzinterface )
      call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
 
    case ( regridding_sigma )
      call build_sigma_grid( cs, g, gv, h, dzinterface )
      call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
 
    case ( regridding_rho )
      if (do_convective_adjustment) call convective_adjustment(g, gv, h, tv)
      call build_rho_grid( g, gv, g%US, h, tv, dzinterface, remapcs, cs )
      call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
 
    case ( regridding_arbitrary )
      call build_grid_arbitrary( g, gv, h, dzinterface, trickgnucompiler, cs )
      call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
 
    case ( regridding_hycom1 )
      call build_grid_hycom1( g, gv, g%US, h, tv, h_new, dzinterface, cs )
 
    case ( regridding_slight )
      call build_grid_slight( g, gv, g%US, h, tv, dzinterface, cs )
      call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
 
    case ( regridding_adaptive )
      call build_grid_adaptive(g, gv, g%US, h, tv, dzinterface, remapcs, cs)
      call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
 
    case default
      call mom_error(fatal,'MOM_regridding, regridding_main: '//&
                     'Unknown regridding scheme selected!')
 
  end select ! type of grid
 
#ifdef __DO_SAFETY_CHECKS__
  call check_remapping_grid(g, gv, h, dzinterface,'in regridding_main')
#endif
 

rho_function1()

integer function mom_regridding::rho_function1 ( character(len=*), intent(in)  string, real, dimension(:), intent(inout), allocatable  rho_target  )

private

Parses a string and generates a rho_target(:) profile with refined resolution downward and returns the number of levels.

Parameters

[in]	string	String with list of parameters in form dz_min, H_total, power, precision
[in,out]	rho_target	Profile of interface densities [kg m-3]

Definition at line 2446 of file MOM_regridding.F90.

  character(len=*),   intent(in)    :: string !< String with list of parameters in form
                                              !! dz_min, H_total, power, precision
  real, dimension(:), allocatable, intent(inout) :: rho_target !< Profile of interface densities [kg m-3]
  ! Local variables
  integer :: nki, k, nk
  real    :: ddx, dx, rho_1, rho_2, rho_3, drho, rho_4, drho_min
 
  read( string, *) nk, rho_1, rho_2, rho_3, drho, rho_4, drho_min
  allocate(rho_target(nk+1))
  nki = nk + 1 - 4 ! Number of interfaces minus 4 specified values
  rho_target(1) = rho_1
  rho_target(2) = rho_2
  dx = 0.
  do k = 0, nki
    ddx = max( drho_min, real(nki-k)/real(nki*nki) )
    dx = dx + ddx
    rho_target(3+k) = rho_3 + (2. * drho) * dx
  enddo
  rho_target(nki+4) = rho_4
 
  rho_function1 = nk
 

set_regrid_max_depths()

subroutine, public mom_regridding::set_regrid_max_depths	(	type(regridding_cs), intent(inout)	CS,
		real, dimension(cs%nk+1), intent(in)	max_depths,
		real, intent(in), optional	units_to_H
	)

Set maximum interface depths based on a vector of input values.

Parameters

[in,out]	cs	Regridding control structure
[in]	max_depths	Maximum interface depths, in arbitrary units
[in]	units_to_h	A conversion factor for max_depths into H units

Definition at line 2056 of file MOM_regridding.F90.

  type(regridding_CS),      intent(inout) :: CS !< Regridding control structure
  real, dimension(CS%nk+1), intent(in)    :: max_depths !< Maximum interface depths, in arbitrary units
  real, optional,           intent(in)    :: units_to_H !< A conversion factor for max_depths into H units
  ! Local variables
  real :: val_to_H
  integer :: K
 
  if (.not.allocated(cs%max_interface_depths)) allocate(cs%max_interface_depths(1:cs%nk+1))
 
  val_to_h = 1.0 ; if (present(units_to_h)) val_to_h = units_to_h
  if (max_depths(cs%nk+1) < max_depths(1)) val_to_h = -1.0*val_to_h
 
  ! Check for sign reversals in the depths.
  if (max_depths(cs%nk+1) < max_depths(1)) then
    do k=1,cs%nk ; if (max_depths(k+1) > max_depths(k)) &
      call mom_error(fatal, "Unordered list of maximum depths sent to set_regrid_max_depths!")
    enddo
  else
    do k=1,cs%nk ; if (max_depths(k+1) < max_depths(k)) &
      call mom_error(fatal, "Unordered list of maximum depths sent to set_regrid_max_depths.")
    enddo
  endif
 
  do k=1,cs%nk+1
    cs%max_interface_depths(k) = val_to_h * max_depths(k)
  enddo
 
  ! set max depths for coordinate
  select case (cs%regridding_scheme)
  case (regridding_hycom1)
    call set_hycom_params(cs%hycom_CS, max_interface_depths=cs%max_interface_depths)
  case (regridding_slight)
    call set_slight_params(cs%slight_CS, max_interface_depths=cs%max_interface_depths)
  end select

set_regrid_max_thickness()

subroutine, public mom_regridding::set_regrid_max_thickness	(	type(regridding_cs), intent(inout)	CS,
		real, dimension(cs%nk+1), intent(in)	max_h,
		real, intent(in), optional	units_to_H
	)

Set maximum layer thicknesses based on a vector of input values.

Parameters

[in,out]	cs	Regridding control structure
[in]	max_h	Maximum interface depths, in arbitrary units
[in]	units_to_h	A conversion factor for max_h into H units

Definition at line 2094 of file MOM_regridding.F90.

  type(regridding_CS),      intent(inout) :: CS !< Regridding control structure
  real, dimension(CS%nk+1), intent(in)    :: max_h !< Maximum interface depths, in arbitrary units
  real, optional,           intent(in)    :: units_to_H !< A conversion factor for max_h into H units
  ! Local variables
  real :: val_to_H
  integer :: K
 
  if (.not.allocated(cs%max_layer_thickness)) allocate(cs%max_layer_thickness(1:cs%nk))
 
  val_to_h = 1.0 ; if (present( units_to_h)) val_to_h = units_to_h
 
  do k=1,cs%nk
    cs%max_layer_thickness(k) = val_to_h * max_h(k)
  enddo
 
  ! set max thickness for coordinate
  select case (cs%regridding_scheme)
  case (regridding_hycom1)
    call set_hycom_params(cs%hycom_CS, max_layer_thickness=cs%max_layer_thickness)
  case (regridding_slight)
    call set_slight_params(cs%slight_CS, max_layer_thickness=cs%max_layer_thickness)
  end select

set_regrid_params()

subroutine, public mom_regridding::set_regrid_params	(	type(regridding_cs), intent(inout)	CS,
		logical, intent(in), optional	boundary_extrapolation,
		real, intent(in), optional	min_thickness,
		real, intent(in), optional	old_grid_weight,
		character(len=*), intent(in), optional	interp_scheme,
		real, intent(in), optional	depth_of_time_filter_shallow,
		real, intent(in), optional	depth_of_time_filter_deep,
		real, intent(in), optional	compress_fraction,
		real, intent(in), optional	ref_pressure,
		real, intent(in), optional	dz_min_surface,
		integer, intent(in), optional	nz_fixed_surface,
		real, intent(in), optional	Rho_ML_avg_depth,
		real, intent(in), optional	nlay_ML_to_interior,
		logical, intent(in), optional	fix_haloclines,
		real, intent(in), optional	halocline_filt_len,
		real, intent(in), optional	halocline_strat_tol,
		logical, intent(in), optional	integrate_downward_for_e,
		logical, intent(in), optional	remap_answers_2018,
		real, intent(in), optional	adaptTimeRatio,
		real, intent(in), optional	adaptZoom,
		real, intent(in), optional	adaptZoomCoeff,
		real, intent(in), optional	adaptBuoyCoeff,
		real, intent(in), optional	adaptAlpha,
		logical, intent(in), optional	adaptDoMin,
		real, intent(in), optional	adaptDrho0
	)

Can be used to set any of the parameters for MOM_regridding.

Parameters

[in,out]	cs	Regridding control structure
[in]	boundary_extrapolation	Extrapolate in boundary cells
[in]	min_thickness	Minimum thickness allowed when building the new grid [H ~> m or kg m-2]
[in]	old_grid_weight	Weight given to old coordinate when time-filtering grid
[in]	interp_scheme	Interpolation method for state-dependent coordinates
[in]	depth_of_time_filter_shallow	Depth to start cubic [H ~> m or kg m-2]
[in]	depth_of_time_filter_deep	Depth to end cubic [H ~> m or kg m-2]
[in]	compress_fraction	Fraction of compressibility to add to potential density [nondim]
[in]	ref_pressure	The reference pressure for density-dependent coordinates [R L2 T-2 ~> Pa]
[in]	dz_min_surface	The fixed resolution in the topmost SLight_nkml_min layers [H ~> m or kg m-2]
[in]	nz_fixed_surface	The number of fixed-thickness layers at the top of the model
[in]	rho_ml_avg_depth	Averaging depth over which to determine mixed layer potential density [H ~> m or kg m-2]
[in]	nlay_ml_to_interior	Number of layers to offset the mixed layer density to find resolved stratification [nondim]
[in]	fix_haloclines	Detect regions with much weaker stratification in the coordinate
[in]	halocline_filt_len	Length scale over which to filter T & S when looking for spuriously unstable water mass profiles [H ~> m or kg m-2]
[in]	halocline_strat_tol	Value of the stratification ratio that defines a problematic halocline region.
[in]	integrate_downward_for_e	If true, integrate for interface positions downward from the top.
[in]	remap_answers_2018	If true, use the order of arithmetic and expressions that recover the remapping answers from 2018. Otherwise use more robust but mathematically equivalent expressions.
[in]	adapttimeratio	Ratio of the ALE timestep to the grid timescale [nondim].
[in]	adaptzoom	Depth of near-surface zooming region [H ~> m or kg m-2].
[in]	adaptzoomcoeff	Coefficient of near-surface zooming diffusivity [nondim].
[in]	adaptbuoycoeff	Coefficient of buoyancy diffusivity [nondim].
[in]	adaptalpha	Scaling factor on optimization tendency [nondim].
[in]	adaptdomin	If true, make a HyCOM-like mixed layer by preventing interfaces from being shallower than the depths specified by the regridding coordinate.
[in]	adaptdrho0	Reference density difference for stratification-dependent diffusion. [R ~> kg m-3]

Definition at line 2240 of file MOM_regridding.F90.

  type(regridding_CS), intent(inout) :: CS !< Regridding control structure
  logical, optional, intent(in) :: boundary_extrapolation !< Extrapolate in boundary cells
  real,    optional, intent(in) :: min_thickness    !< Minimum thickness allowed when building the
                                                    !! new grid [H ~> m or kg m-2]
  real,    optional, intent(in) :: old_grid_weight  !< Weight given to old coordinate when time-filtering grid
  character(len=*), optional, intent(in) :: interp_scheme !< Interpolation method for state-dependent coordinates
  real,    optional, intent(in) :: depth_of_time_filter_shallow !< Depth to start cubic [H ~> m or kg m-2]
  real,    optional, intent(in) :: depth_of_time_filter_deep !< Depth to end cubic [H ~> m or kg m-2]
  real,    optional, intent(in) :: compress_fraction !< Fraction of compressibility to add to potential density [nondim]
  real,    optional, intent(in) :: ref_pressure     !< The reference pressure for density-dependent
                                                    !! coordinates [R L2 T-2 ~> Pa]
  real,    optional, intent(in) :: dz_min_surface   !< The fixed resolution in the topmost
                                                    !! SLight_nkml_min layers [H ~> m or kg m-2]
  integer, optional, intent(in) :: nz_fixed_surface !< The number of fixed-thickness layers at the top of the model
  real,    optional, intent(in) :: Rho_ml_avg_depth !< Averaging depth over which to determine mixed layer potential
                                                    !! density [H ~> m or kg m-2]
  real,    optional, intent(in) :: nlay_ML_to_interior !< Number of layers to offset the mixed layer density to find
                                                    !! resolved stratification [nondim]
  logical, optional, intent(in) :: fix_haloclines   !< Detect regions with much weaker stratification in the coordinate
  real,    optional, intent(in) :: halocline_filt_len !< Length scale over which to filter T & S when looking for
                                                    !! spuriously unstable water mass profiles [H ~> m or kg m-2]
  real,    optional, intent(in) :: halocline_strat_tol !< Value of the stratification ratio that defines a problematic
                                                    !! halocline region.
  logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface positions downward
                                                    !! from the top.
  logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions
                                                    !! that recover the remapping answers from 2018.  Otherwise
                                                    !! use more robust but mathematically equivalent expressions.
  real,    optional, intent(in) :: adaptTimeRatio   !< Ratio of the ALE timestep to the grid timescale [nondim].
  real,    optional, intent(in) :: adaptZoom        !< Depth of near-surface zooming region [H ~> m or kg m-2].
  real,    optional, intent(in) :: adaptZoomCoeff   !< Coefficient of near-surface zooming diffusivity [nondim].
  real,    optional, intent(in) :: adaptBuoyCoeff   !< Coefficient of buoyancy diffusivity [nondim].
  real,    optional, intent(in) :: adaptAlpha       !< Scaling factor on optimization tendency [nondim].
  logical, optional, intent(in) :: adaptDoMin       !< If true, make a HyCOM-like mixed layer by
                                                    !! preventing interfaces from being shallower than
                                                    !! the depths specified by the regridding coordinate.
  real,    optional, intent(in) :: adaptDrho0       !< Reference density difference for stratification-dependent
                                                    !! diffusion. [R ~> kg m-3]
 
  if (present(interp_scheme)) call set_interp_scheme(cs%interp_CS, interp_scheme)
  if (present(boundary_extrapolation)) call set_interp_extrap(cs%interp_CS, boundary_extrapolation)
 
  if (present(old_grid_weight)) then
    if (old_grid_weight<0. .or. old_grid_weight>1.) &
      call mom_error(fatal,'MOM_regridding, set_regrid_params: Weight is out side the range 0..1!')
    cs%old_grid_weight = old_grid_weight
  endif
  if (present(depth_of_time_filter_shallow)) cs%depth_of_time_filter_shallow = depth_of_time_filter_shallow
  if (present(depth_of_time_filter_deep)) cs%depth_of_time_filter_deep = depth_of_time_filter_deep
  if (present(depth_of_time_filter_shallow) .or. present(depth_of_time_filter_deep)) then
    if (cs%depth_of_time_filter_deep<cs%depth_of_time_filter_shallow) call mom_error(fatal,'MOM_regridding, '//&
                     'set_regrid_params: depth_of_time_filter_deep<depth_of_time_filter_shallow!')
  endif
 
  if (present(min_thickness)) cs%min_thickness = min_thickness
  if (present(compress_fraction)) cs%compressibility_fraction = compress_fraction
  if (present(ref_pressure)) cs%ref_pressure = ref_pressure
  if (present(integrate_downward_for_e)) cs%integrate_downward_for_e = integrate_downward_for_e
  if (present(remap_answers_2018)) cs%remap_answers_2018 = remap_answers_2018
 
  select case (cs%regridding_scheme)
  case (regridding_zstar)
    if (present(min_thickness)) call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
  case (regridding_sigma_shelf_zstar)
    if (present(min_thickness)) call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
  case (regridding_sigma)
    if (present(min_thickness)) call set_sigma_params(cs%sigma_CS, min_thickness=min_thickness)
  case (regridding_rho)
    if (present(min_thickness)) call set_rho_params(cs%rho_CS, min_thickness=min_thickness)
    if (present(integrate_downward_for_e)) &
      call set_rho_params(cs%rho_CS, integrate_downward_for_e=integrate_downward_for_e)
    if (associated(cs%rho_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
      call set_rho_params(cs%rho_CS, interp_cs=cs%interp_CS)
  case (regridding_hycom1)
    if (associated(cs%hycom_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
      call set_hycom_params(cs%hycom_CS, interp_cs=cs%interp_CS)
  case (regridding_slight)
    if (present(min_thickness))       call set_slight_params(cs%slight_CS, min_thickness=min_thickness)
    if (present(dz_min_surface))      call set_slight_params(cs%slight_CS, dz_ml_min=dz_min_surface)
    if (present(nz_fixed_surface))    call set_slight_params(cs%slight_CS, nz_fixed_surface=nz_fixed_surface)
    if (present(rho_ml_avg_depth))    call set_slight_params(cs%slight_CS, rho_ml_avg_depth=rho_ml_avg_depth)
    if (present(nlay_ml_to_interior)) call set_slight_params(cs%slight_CS, nlay_ml_offset=nlay_ml_to_interior)
    if (present(fix_haloclines))      call set_slight_params(cs%slight_CS, fix_haloclines=fix_haloclines)
    if (present(halocline_filt_len))  call set_slight_params(cs%slight_CS, halocline_filter_length=halocline_filt_len)
    if (present(halocline_strat_tol)) call set_slight_params(cs%slight_CS, halocline_strat_tol=halocline_strat_tol)
    if (present(compress_fraction))   call set_slight_params(cs%slight_CS, compressibility_fraction=compress_fraction)
    if (associated(cs%slight_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
      call set_slight_params(cs%slight_CS, interp_cs=cs%interp_CS)
  case (regridding_adaptive)
    if (present(adapttimeratio)) call set_adapt_params(cs%adapt_CS, adapttimeratio=adapttimeratio)
    if (present(adaptzoom))      call set_adapt_params(cs%adapt_CS, adaptzoom=adaptzoom)
    if (present(adaptzoomcoeff)) call set_adapt_params(cs%adapt_CS, adaptzoomcoeff=adaptzoomcoeff)
    if (present(adaptbuoycoeff)) call set_adapt_params(cs%adapt_CS, adaptbuoycoeff=adaptbuoycoeff)
    if (present(adaptalpha))     call set_adapt_params(cs%adapt_CS, adaptalpha=adaptalpha)
    if (present(adaptdomin))     call set_adapt_params(cs%adapt_CS, adaptdomin=adaptdomin)
    if (present(adaptdrho0))     call set_adapt_params(cs%adapt_CS, adaptdrho0=adaptdrho0)
  end select
 

set_target_densities()

subroutine, public mom_regridding::set_target_densities	(	type(regridding_cs), intent(inout)	CS,
		real, dimension(cs%nk+1), intent(in)	rho_int
	)

Set target densities based on vector of interface values.

Parameters

[in,out]	cs	Regridding control structure
[in]	rho_int	Interface densities [R ~> kg m-3]

Definition at line 2042 of file MOM_regridding.F90.

  type(regridding_CS),      intent(inout) :: CS !< Regridding control structure
  real, dimension(CS%nk+1), intent(in)    :: rho_int !< Interface densities [R ~> kg m-3]
 
  if (size(cs%target_density)/=size(rho_int)) then
    call mom_error(fatal, "set_target_densities inconsistent args!")
  endif
 
  cs%target_density(:) = rho_int(:)
  cs%target_density_set = .true.
 

set_target_densities_from_gv()

subroutine, public mom_regridding::set_target_densities_from_gv	(	type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(regridding_cs), intent(inout)	CS
	)

Set target densities based on the old Rlay variable.

Parameters

[in]	gv	Ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in,out]	cs	Regridding control structure

Definition at line 2019 of file MOM_regridding.F90.

  type(verticalGrid_type), intent(in)    :: GV !< Ocean vertical grid structure
  type(unit_scale_type),   intent(in)    :: US !< A dimensional unit scaling type
  type(regridding_CS),     intent(inout) :: CS !< Regridding control structure
  ! Local variables
  integer :: k, nz
 
  nz = cs%nk
  if (nz == 1) then ! Set a broad range of bounds.  Regridding may not be meaningful in this case.
    cs%target_density(1)    = 0.0
    cs%target_density(2)    = 2.0*gv%Rlay(1)
  else
    cs%target_density(1)    = (gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(2)))
    cs%target_density(nz+1) = (gv%Rlay(nz) + 0.5*(gv%Rlay(nz)-gv%Rlay(nz-1)))
    do k=2,nz
      cs%target_density(k)  = cs%target_density(k-1) + cs%coordinateResolution(k)
    enddo
  endif
  cs%target_density_set = .true.
 

setcoordinateresolution()

subroutine, public mom_regridding::setcoordinateresolution	(	real, dimension(:), intent(in)	dz,
		type(regridding_cs), intent(inout)	CS,
		real, intent(in), optional	scale
	)

Set the fixed resolution data.

Parameters

[in]	dz	A vector of vertical grid spacings
[in,out]	cs	Regridding control structure
[in]	scale	A scaling factor converting dz to coordRes

Definition at line 2002 of file MOM_regridding.F90.

  real, dimension(:),  intent(in)    :: dz !< A vector of vertical grid spacings
  type(regridding_CS), intent(inout) :: CS !< Regridding control structure
  real,      optional, intent(in)    :: scale !< A scaling factor converting dz to coordRes
 
  if (size(dz)/=cs%nk) call mom_error( fatal, &
      'setCoordinateResolution: inconsistent number of levels' )
 
  if (present(scale)) then
    cs%coordinateResolution(:) = scale*dz(:)
  else
    cs%coordinateResolution(:) = dz(:)
  endif
 

uniformresolution()

real function, dimension(nk), public mom_regridding::uniformresolution	(	integer, intent(in)	nk,
		character(len=*), intent(in)	coordMode,
		real, intent(in)	maxDepth,
		real, intent(in)	rhoLight,
		real, intent(in)	rhoHeavy
	)

Return a uniform resolution vector in the units of the coordinate.

Parameters

[in]	nk	Number of cells in source grid
[in]	coordmode	A string indicating the coordinate mode. See the documenttion for regrid_consts for the recognized values.
[in]	maxdepth	The range of the grid values in some modes
[in]	rholight	The minimum value of the grid in RHO mode
[in]	rhoheavy	The maximum value of the grid in RHO mode

Returns

The returned uniform resolution grid.

Definition at line 1931 of file MOM_regridding.F90.

!------------------------------------------------------------------------------
! Calculate a vector of uniform resolution in the units of the coordinate
!------------------------------------------------------------------------------
  ! Arguments
  integer,          intent(in) :: nk !< Number of cells in source grid
  character(len=*), intent(in) :: coordMode !< A string indicating the coordinate mode.
                                            !! See the documenttion for regrid_consts
                                            !! for the recognized values.
  real,             intent(in) :: maxDepth  !< The range of the grid values in some modes
  real,             intent(in) :: rhoLight  !< The minimum value of the grid in RHO mode
  real,             intent(in) :: rhoHeavy  !< The maximum value of the grid in RHO mode
 
  real                         :: uniformResolution(nk) !< The returned uniform resolution grid.
 
  ! Local variables
  integer :: scheme
 
  scheme = coordinatemode(coordmode)
  select case ( scheme )
 
    case ( regridding_zstar, regridding_hycom1, regridding_slight, regridding_sigma_shelf_zstar, &
           regridding_adaptive )
      uniformresolution(:) = maxdepth / real(nk)
 
    case ( regridding_rho )
      uniformresolution(:) = (rhoheavy - rholight) / real(nk)
 
    case ( regridding_sigma )
      uniformresolution(:) = 1. / real(nk)
 
    case default
      call mom_error(fatal, "MOM_regridding, uniformResolution: "//&
       "Unrecognized choice for coordinate mode ("//trim(coordmode)//").")
 
  end select ! type of grid