MOM6
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, 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, 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, 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, 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, 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, adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin)
 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. More...
 
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. More...
 
character(len= *), parameter, public regriddingdefaultinterpscheme = "P1M_H2"
 Default interpolation scheme. More...
 
logical, parameter, public regriddingdefaultboundaryextrapolation = .false.
 Default mode for boundary extrapolation. More...
 
real, parameter, public regriddingdefaultminthickness = 1.e-3
 Default minimum thickness for some coordinate generation modes. More...
 

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]csRegridding control structure
[in]nkNumber of layers in h_old
[in]h_oldMinium allowed thickness of h [H ~> m or kg m-2]
[in,out]dz_intMinium allowed thickness of h [H ~> m or kg m-2]

Definition at line 1640 of file MOM_regridding.F90.

1640  type(regridding_CS), intent(in) :: CS !< Regridding control structure
1641  integer, intent(in) :: nk !< Number of layers in h_old
1642  real, dimension(nk), intent(in) :: h_old !< Minium allowed thickness of h [H ~> m or kg m-2]
1643  real, dimension(CS%nk+1), intent(inout) :: dz_int !< Minium allowed thickness of h [H ~> m or kg m-2]
1644  ! Local variables
1645  integer :: k
1646  real :: h_new, eps, h_total, h_err
1647 
1648  eps = 1. ; eps = epsilon(eps)
1649 
1650  h_total = 0. ; h_err = 0.
1651  do k = 1, min(cs%nk,nk)
1652  h_total = h_total + h_old(k)
1653  h_err = h_err + max( h_old(k), abs(dz_int(k)), abs(dz_int(k+1)) )*eps
1654  h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1655  if (h_new < -3.0*h_err) then
1656  write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), &
1657  'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
1658  'h_new=',h_new,'h_err=',h_err
1659  call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
1660  'implied h<0 is larger than roundoff!')
1661  endif
1662  enddo
1663  if (cs%nk>nk) then
1664  do k = nk+1, cs%nk
1665  h_err = h_err + max( abs(dz_int(k)), abs(dz_int(k+1)) )*eps
1666  h_new = ( dz_int(k) - dz_int(k+1) )
1667  if (h_new < -3.0*h_err) then
1668  write(0,*) 'h<0 at k=',k,'h_old was empty',&
1669  'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
1670  'h_new=',h_new,'h_err=',h_err
1671  call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
1672  'implied h<0 is larger than roundoff!')
1673  endif
1674  enddo
1675  endif
1676  do k = min(cs%nk,nk),2,-1
1677  h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1678  if (h_new<cs%min_thickness) &
1679  dz_int(k) = ( dz_int(k+1) - h_old(k) ) + cs%min_thickness ! Implies next h_new = min_thickness
1680  h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1681  if (h_new<0.) &
1682  dz_int(k) = ( 1. - eps ) * ( dz_int(k+1) - h_old(k) ) ! Backup in case min_thickness==0
1683  h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1684  if (h_new<0.) then
1685  write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), &
1686  'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
1687  'h_new=',h_new
1688  stop 'Still did not work!'
1689  call mom_error( fatal, 'MOM_regridding: adjust_interface_motion() - '//&
1690  'Repeated adjustment for roundoff h<0 failed!')
1691  endif
1692  enddo
1693  !if (dz_int(1)/=0.) stop 'MOM_regridding: adjust_interface_motion() surface moved'
1694 

References mom_error_handler::mom_error().

Referenced by build_grid_adaptive(), build_grid_hycom1(), build_grid_slight(), and build_zstar_grid().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ build_grid_adaptive()

subroutine mom_regridding::build_grid_adaptive ( 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,
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]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]tvA structure pointing to various thermodynamic variables
[in,out]dzinterfaceThe change in interface depth [H ~> m or kg m-2]
[in]remapcsThe remapping control structure
[in]csRegridding control structure

Definition at line 1510 of file MOM_regridding.F90.

1510  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1511  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1512  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1513  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
1514  !! thermodynamic variables
1515  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth
1516  !! [H ~> m or kg m-2]
1517  type(remapping_CS), intent(in) :: remapCS !< The remapping control structure
1518  type(regridding_CS), intent(in) :: CS !< Regridding control structure
1519 
1520  ! local variables
1521  integer :: i, j, k, nz ! indices and dimension lengths
1522  ! temperature, salinity and pressure on interfaces
1523  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt, sInt
1524  ! current interface positions and after tendency term is applied
1525  ! positive downward
1526  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt
1527  real, dimension(SZK_(GV)+1) :: zNext
1528 
1529  nz = gv%ke
1530 
1531  ! position surface at z = 0.
1532  zint(:,:,1) = 0.
1533 
1534  ! work on interior interfaces
1535  do k = 2, nz ; do j = g%jsc-2,g%jec+2 ; do i = g%isc-2,g%iec+2
1536  tint(i,j,k) = 0.5 * (tv%T(i,j,k-1) + tv%T(i,j,k))
1537  sint(i,j,k) = 0.5 * (tv%S(i,j,k-1) + tv%S(i,j,k))
1538  zint(i,j,k) = zint(i,j,k-1) + h(i,j,k-1) ! zInt in [H]
1539  enddo ; enddo ; enddo
1540 
1541  ! top and bottom temp/salt interfaces are just the layer
1542  ! average values
1543  tint(:,:,1) = tv%T(:,:,1) ; tint(:,:,nz+1) = tv%T(:,:,nz)
1544  sint(:,:,1) = tv%S(:,:,1) ; sint(:,:,nz+1) = tv%S(:,:,nz)
1545 
1546  ! set the bottom interface depth
1547  zint(:,:,nz+1) = zint(:,:,nz) + h(:,:,nz)
1548 
1549  ! calculate horizontal density derivatives (alpha/beta)
1550  ! between cells in a 5-point stencil, columnwise
1551  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
1552  if (g%mask2dT(i,j) < 0.5) then
1553  dzinterface(i,j,:) = 0. ! land point, don't move interfaces, and skip
1554  cycle
1555  endif
1556 
1557  call build_adapt_column(cs%adapt_CS, g, gv, tv, i, j, zint, tint, sint, h, znext)
1558 
1559  call filtered_grid_motion(cs, nz, zint(i,j,:), znext, dzinterface(i,j,:))
1560  ! convert from depth to z
1561  do k = 1, nz+1 ; dzinterface(i,j,k) = -dzinterface(i,j,k) ; enddo
1562  call adjust_interface_motion(cs, nz, h(i,j,:), dzinterface(i,j,:))
1563  enddo ; enddo

References adjust_interface_motion(), coord_adapt::build_adapt_column(), and filtered_grid_motion().

Referenced by regridding_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]gOcean grid structure
[in]gvOcean vertical grid structure
[in]hOriginal layer thicknesses [H ~> m or kg m-2]
[in,out]dzinterfaceThe change in interface depth [H ~> m or kg m-2]
[in,out]h_newNew layer thicknesses [H ~> m or kg m-2]
[in]csRegridding control structure

Definition at line 1701 of file MOM_regridding.F90.

1701 !------------------------------------------------------------------------------
1702 ! This routine builds a grid based on arbitrary rules
1703 !------------------------------------------------------------------------------
1704 
1705  ! Arguments
1706  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1707  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1708  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(in) :: h !< Original layer thicknesses [H ~> m or kg m-2]
1709  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface
1710  !! depth [H ~> m or kg m-2]
1711  real, intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2]
1712  type(regridding_CS), intent(in) :: CS !< Regridding control structure
1713 
1714  ! Local variables
1715  integer :: i, j, k
1716  integer :: nz
1717  real :: z_inter(SZK_(GV)+1)
1718  real :: total_height
1719  real :: delta_h
1720  real :: max_depth
1721  real :: eta ! local elevation
1722  real :: local_depth
1723  real :: x1, y1, x2, y2
1724  real :: x, t
1725 
1726  nz = gv%ke
1727  max_depth = g%max_depth*gv%Z_to_H
1728 
1729  do j = g%jsc-1,g%jec+1
1730  do i = g%isc-1,g%iec+1
1731 
1732  ! Local depth
1733  local_depth = g%bathyT(i,j)*gv%Z_to_H
1734 
1735  ! Determine water column height
1736  total_height = 0.0
1737  do k = 1,nz
1738  total_height = total_height + h(i,j,k)
1739  enddo
1740 
1741  eta = total_height - local_depth
1742 
1743  ! Compute new thicknesses based on stretched water column
1744  delta_h = (max_depth + eta) / nz
1745 
1746  ! Define interfaces
1747  z_inter(1) = eta
1748  do k = 1,nz
1749  z_inter(k+1) = z_inter(k) - delta_h
1750  enddo
1751 
1752  ! Refine grid in the middle
1753  do k = 1,nz+1
1754  x1 = 0.35; y1 = 0.45; x2 = 0.65; y2 = 0.55
1755 
1756  x = - ( z_inter(k) - eta ) / max_depth
1757 
1758  if ( x <= x1 ) then
1759  t = y1*x/x1
1760  elseif ( (x > x1 ) .and. ( x < x2 )) then
1761  t = y1 + (y2-y1) * (x-x1) / (x2-x1)
1762  else
1763  t = y2 + (1.0-y2) * (x-x2) / (1.0-x2)
1764  endif
1765 
1766  z_inter(k) = -t * max_depth + eta
1767 
1768  enddo
1769 
1770  ! Modify interface heights to account for topography
1771  z_inter(nz+1) = - local_depth
1772 
1773  ! Modify interface heights to avoid layers of zero thicknesses
1774  do k = nz,1,-1
1775  if ( z_inter(k) < (z_inter(k+1) + cs%min_thickness) ) then
1776  z_inter(k) = z_inter(k+1) + cs%min_thickness
1777  endif
1778  enddo
1779 
1780  ! Chnage in interface position
1781  x = 0. ! Left boundary at x=0
1782  dzinterface(i,j,1) = 0.
1783  do k = 2,nz
1784  x = x + h(i,j,k)
1785  dzinterface(i,j,k) = z_inter(k) - x
1786  enddo
1787  dzinterface(i,j,nz+1) = 0.
1788 
1789  enddo
1790  enddo
1791 
1792 stop 'OOOOOOPS' ! For some reason the gnu compiler will not let me delete this
1793  ! routine????
1794 

Referenced by regridding_main().

Here is the caller graph for this function:

◆ build_grid_hycom1()

subroutine mom_regridding::build_grid_hycom1 ( 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,
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]gGrid structure
[in]gvOcean vertical grid structure
[in]hExisting model thickness [H ~> m or kg m-2]
[in]tvThermodynamics structure
[in]csRegridding control structure
[in,out]h_newNew layer thicknesses [H ~> m or kg m-2]
[in,out]dzinterfaceChanges in interface position

Definition at line 1441 of file MOM_regridding.F90.

1441  type(ocean_grid_type), intent(in) :: G !< Grid structure
1442  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1443  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Existing model thickness [H ~> m or kg m-2]
1444  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1445  type(regridding_CS), intent(in) :: CS !< Regridding control structure
1446  real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2]
1447  real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< Changes in interface position
1448 
1449  ! Local variables
1450  real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2]
1451  real, dimension(CS%nk+1) :: z_col_new ! New interface positions relative to the surface [H ~> m or kg m-2]
1452  real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2]
1453  real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [Pa]
1454  integer :: i, j, k, nki
1455  real :: depth
1456  real :: h_neglect, h_neglect_edge
1457 
1458  !### Try replacing both of these with GV%H_subroundoff
1459  if (gv%Boussinesq) then
1460  h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1461  else
1462  h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1463  endif
1464 
1465  if (.not.cs%target_density_set) call mom_error(fatal, "build_grid_HyCOM1 : "//&
1466  "Target densities must be set before build_grid_HyCOM1 is called.")
1467 
1468  nki = min(gv%ke, cs%nk)
1469 
1470  ! Build grid based on target interface densities
1471  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
1472  if (g%mask2dT(i,j)>0.) then
1473 
1474  depth = g%bathyT(i,j) * gv%Z_to_H
1475 
1476  z_col(1) = 0. ! Work downward rather than bottom up
1477  do k = 1, gv%ke
1478  z_col(k+1) = z_col(k) + h(i,j,k)
1479  p_col(k) = cs%ref_pressure + cs%compressibility_fraction * &
1480  ( 0.5 * ( z_col(k) + z_col(k+1) ) * gv%H_to_Pa - cs%ref_pressure )
1481  enddo
1482 
1483  call build_hycom1_column(cs%hycom_CS, tv%eqn_of_state, gv%ke, depth, &
1484  h(i, j, :), tv%T(i, j, :), tv%S(i, j, :), p_col, &
1485  z_col, z_col_new, zscale=gv%Z_to_H, &
1486  h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1487 
1488  ! Calculate the final change in grid position after blending new and old grids
1489  call filtered_grid_motion( cs, gv%ke, z_col, z_col_new, dz_col )
1490 
1491  ! This adjusts things robust to round-off errors
1492  dz_col(:) = -dz_col(:)
1493  call adjust_interface_motion( cs, gv%ke, h(i,j,:), dz_col(:) )
1494 
1495  dzinterface(i,j,1:nki+1) = dz_col(1:nki+1)
1496  if (nki<cs%nk) dzinterface(i,j,nki+2:cs%nk+1) = 0.
1497 
1498  else ! on land
1499  dzinterface(i,j,:) = 0.
1500  endif ! mask2dT
1501  enddo ; enddo ! i,j
1502 
1503  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
1504 

References adjust_interface_motion(), coord_hycom::build_hycom1_column(), calc_h_new_by_dz(), filtered_grid_motion(), and mom_error_handler::mom_error().

Referenced by regridding_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ build_grid_slight()

subroutine mom_regridding::build_grid_slight ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
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]gGrid structure
[in]gvOcean vertical grid structure
[in]hExisting model thickness [H ~> m or kg m-2]
[in]tvThermodynamics structure
[in,out]dzinterfaceChanges in interface position
[in]csRegridding control structure

Definition at line 1576 of file MOM_regridding.F90.

1576  type(ocean_grid_type), intent(in) :: G !< Grid structure
1577  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1578  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Existing model thickness [H ~> m or kg m-2]
1579  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1580  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< Changes in interface position
1581  type(regridding_CS), intent(in) :: CS !< Regridding control structure
1582 
1583  real, dimension(SZK_(GV)+1) :: z_col ! Interface positions relative to the surface [H ~> m or kg m-2]
1584  real, dimension(SZK_(GV)+1) :: z_col_new ! Interface positions relative to the surface [H ~> m or kg m-2]
1585  real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2]
1586  real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [Pa]
1587  real :: depth
1588  integer :: i, j, k, nz
1589  real :: h_neglect, h_neglect_edge
1590 
1591  !### Try replacing both of these with GV%H_subroundoff
1592  if (gv%Boussinesq) then
1593  h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1594  else
1595  h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1596  endif
1597 
1598  nz = gv%ke
1599 
1600  if (.not.cs%target_density_set) call mom_error(fatal, "build_grid_SLight : "//&
1601  "Target densities must be set before build_grid_SLight is called.")
1602 
1603  ! Build grid based on target interface densities
1604  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
1605  if (g%mask2dT(i,j)>0.) then
1606 
1607  depth = g%bathyT(i,j) * gv%Z_to_H
1608  z_col(1) = 0. ! Work downward rather than bottom up
1609  do k=1,nz
1610  z_col(k+1) = z_col(k) + h(i,j,k)
1611  p_col(k) = cs%ref_pressure + cs%compressibility_fraction * &
1612  ( 0.5 * ( z_col(k) + z_col(k+1) ) * gv%H_to_Pa - cs%ref_pressure )
1613  enddo
1614 
1615  call build_slight_column(cs%slight_CS, tv%eqn_of_state, gv%H_to_Pa, &
1616  gv%H_subroundoff, nz, depth, h(i, j, :), &
1617  tv%T(i, j, :), tv%S(i, j, :), p_col, z_col, z_col_new, &
1618  h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1619 
1620  ! Calculate the final change in grid position after blending new and old grids
1621  call filtered_grid_motion( cs, nz, z_col, z_col_new, dz_col )
1622  do k=1,nz+1 ; dzinterface(i,j,k) = -dz_col(k) ; enddo
1623 #ifdef __DO_SAFETY_CHECKS__
1624  if (dzinterface(i,j,1) /= 0.) stop 'build_grid_SLight: Surface moved?!'
1625  if (dzinterface(i,j,nz+1) /= 0.) stop 'build_grid_SLight: Bottom moved?!'
1626 #endif
1627 
1628  ! This adjusts things robust to round-off errors
1629  call adjust_interface_motion( cs, nz, h(i,j,:), dzinterface(i,j,:) )
1630 
1631  else ! on land
1632  dzinterface(i,j,:) = 0.
1633  endif ! mask2dT
1634  enddo ; enddo ! i,j
1635 

References adjust_interface_motion(), coord_slight::build_slight_column(), filtered_grid_motion(), and mom_error_handler::mom_error().

Referenced by regridding_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ build_rho_grid()

subroutine mom_regridding::build_rho_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,
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]gOcean grid structure
[in]gvOcean vertical grid structure
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]tvThermodynamics structure
[in,out]dzinterfaceThe change in interface depth [H ~> m or kg m-2]
[in]remapcsThe remapping control structure
[in]csRegridding control structure

Definition at line 1311 of file MOM_regridding.F90.

1311 !------------------------------------------------------------------------------
1312 ! This routine builds a new grid based on a given set of target interface
1313 ! densities (these target densities are computed by taking the mean value
1314 ! of given layer densities). The algorithn operates as follows within each
1315 ! column:
1316 ! 1. Given T & S within each layer, the layer densities are computed.
1317 ! 2. Based on these layer densities, a global density profile is reconstructed
1318 ! (this profile is monotonically increasing and may be discontinuous)
1319 ! 3. The new grid interfaces are determined based on the target interface
1320 ! densities.
1321 ! 4. T & S are remapped onto the new grid.
1322 ! 5. Return to step 1 until convergence or until the maximum number of
1323 ! iterations is reached, whichever comes first.
1324 !------------------------------------------------------------------------------
1325 
1326  ! Arguments
1327  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1328  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1329  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1330  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1331  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth
1332  !! [H ~> m or kg m-2]
1333  type(remapping_CS), intent(in) :: remapCS !< The remapping control structure
1334  type(regridding_CS), intent(in) :: CS !< Regridding control structure
1335 
1336  ! Local variables
1337  integer :: nz
1338  integer :: i, j, k
1339  real :: nominalDepth, totalThickness
1340  real, dimension(SZK_(GV)+1) :: zOld, zNew
1341  real :: h_neglect, h_neglect_edge
1342 #ifdef __DO_SAFETY_CHECKS__
1343  real :: dh
1344 #endif
1345 
1346  !### Try replacing both of these with GV%H_subroundoff
1347  if (gv%Boussinesq) then
1348  h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1349  else
1350  h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1351  endif
1352 
1353  nz = gv%ke
1354 
1355  if (.not.cs%target_density_set) call mom_error(fatal, "build_rho_grid: "//&
1356  "Target densities must be set before build_rho_grid is called.")
1357 
1358  ! Build grid based on target interface densities
1359  do j = g%jsc-1,g%jec+1
1360  do i = g%isc-1,g%iec+1
1361 
1362  if (g%mask2dT(i,j)==0.) then
1363  dzinterface(i,j,:) = 0.
1364  cycle
1365  endif
1366 
1367 
1368  ! Local depth (G%bathyT is positive)
1369  nominaldepth = g%bathyT(i,j)*gv%Z_to_H
1370 
1371  call build_rho_column(cs%rho_CS, nz, nominaldepth, h(i, j, :), &
1372  tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, znew, &
1373  h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1374 
1375  if (cs%integrate_downward_for_e) then
1376  zold(1) = 0.
1377  do k = 1,nz
1378  zold(k+1) = zold(k) - h(i,j,k)
1379  enddo
1380  else
1381  ! The rest of the model defines grids integrating up from the bottom
1382  zold(nz+1) = - nominaldepth
1383  do k = nz,1,-1
1384  zold(k) = zold(k+1) + h(i,j,k)
1385  enddo
1386  endif
1387 
1388  ! Calculate the final change in grid position after blending new and old grids
1389  call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )
1390 
1391 #ifdef __DO_SAFETY_CHECKS__
1392  do k = 2,nz
1393  if (znew(k) > zold(1)) then
1394  write(0,*) 'zOld=',zold
1395  write(0,*) 'zNew=',znew
1396  call mom_error( fatal, 'MOM_regridding, build_rho_grid: '//&
1397  'interior interface above surface!' )
1398  endif
1399  if (znew(k) > znew(k-1)) then
1400  write(0,*) 'zOld=',zold
1401  write(0,*) 'zNew=',znew
1402  call mom_error( fatal, 'MOM_regridding, build_rho_grid: '//&
1403  'interior interfaces cross!' )
1404  endif
1405  enddo
1406 
1407  totalthickness = 0.0
1408  do k = 1,nz
1409  totalthickness = totalthickness + h(i,j,k)
1410  enddo
1411 
1412  dh=max(nominaldepth,totalthickness)
1413  if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh) then
1414  write(0,*) 'min_thickness=',cs%min_thickness
1415  write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
1416  write(0,*) 'zNew(1)-zOld(1) = ',znew(1)-zold(1),epsilon(dh),nz
1417  do k=1,nz+1
1418  write(0,*) k,zold(k),znew(k)
1419  enddo
1420  do k=1,nz
1421  write(0,*) k,h(i,j,k),znew(k)-znew(k+1)
1422  enddo
1423  call mom_error( fatal, &
1424  'MOM_regridding, build_rho_grid: top surface has moved!!!' )
1425  endif
1426 #endif
1427 
1428  enddo ! end loop on i
1429  enddo ! end loop on j
1430 

References filtered_grid_motion(), and mom_error_handler::mom_error().

Referenced by regridding_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]csRegridding control structure
[in]gOcean grid structure
[in]gvocean vertical grid structure
[in]hLayer thicknesses [H ~> m or kg m-2]
[in,out]dzinterfaceThe change in interface depth [H ~> m or kg m-2]

Definition at line 1232 of file MOM_regridding.F90.

1232 !------------------------------------------------------------------------------
1233 ! This routine builds a grid based on terrain-following coordinates.
1234 ! The module parameter coordinateResolution(:) determines the resolution in
1235 ! sigma coordinate, dSigma(:). sigma-coordinates are defined by
1236 ! sigma = (eta-z)/(H+eta) s.t. sigma=0 at z=eta and sigma=1 at z=-H .
1237 !------------------------------------------------------------------------------
1238 
1239  ! Arguments
1240  type(regridding_CS), intent(in) :: CS !< Regridding control structure
1241  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1242  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1243  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1244  real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
1245  !! [H ~> m or kg m-2]
1246 
1247  ! Local variables
1248  integer :: i, j, k
1249  integer :: nz
1250  real :: nominalDepth, totalThickness, dh
1251  real, dimension(SZK_(GV)+1) :: zOld, zNew
1252 
1253  nz = gv%ke
1254 
1255  do i = g%isc-1,g%iec+1
1256  do j = g%jsc-1,g%jec+1
1257 
1258  if (g%mask2dT(i,j)==0.) then
1259  dzinterface(i,j,:) = 0.
1260  cycle
1261  endif
1262 
1263  ! The rest of the model defines grids integrating up from the bottom
1264  nominaldepth = g%bathyT(i,j)*gv%Z_to_H
1265 
1266  ! Determine water column height
1267  totalthickness = 0.0
1268  do k = 1,nz
1269  totalthickness = totalthickness + h(i,j,k)
1270  enddo
1271 
1272  call build_sigma_column(cs%sigma_CS, nominaldepth, totalthickness, znew)
1273 
1274  ! Calculate the final change in grid position after blending new and old grids
1275  zold(nz+1) = -nominaldepth
1276  do k = nz,1,-1
1277  zold(k) = zold(k+1) + h(i, j, k)
1278  enddo
1279 
1280  call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )
1281 
1282 #ifdef __DO_SAFETY_CHECKS__
1283  dh=max(nominaldepth,totalthickness)
1284  if (abs(znew(1)-zold(1))>(cs%nk-1)*0.5*epsilon(dh)*dh) then
1285  write(0,*) 'min_thickness=',cs%min_thickness
1286  write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
1287  write(0,*) 'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz,cs%nk
1288  do k=1,nz+1
1289  write(0,*) k,zold(k),znew(k)
1290  enddo
1291  do k=1,cs%nk
1292  write(0,*) k,h(i,j,k),znew(k)-znew(k+1),totalthickness*cs%coordinateResolution(k),cs%coordinateResolution(k)
1293  enddo
1294  call mom_error( fatal, &
1295  'MOM_regridding, build_sigma_grid: top surface has moved!!!' )
1296  endif
1297  dzinterface(i,j,1) = 0.
1298  dzinterface(i,j,cs%nk+1) = 0.
1299 #endif
1300 
1301  enddo
1302  enddo
1303 

References coord_sigma::build_sigma_column(), filtered_grid_motion(), and mom_error_handler::mom_error().

Referenced by regridding_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]csRegridding control structure
[in]gOcean grid structure
[in]gvocean vertical grid structure
[in]hLayer thicknesses [H ~> m or kg m-2]
[in,out]dzinterfaceThe change in interface depth [H ~> m or kg m-2].
frac_shelf_hFractional ice shelf coverage.

Definition at line 1137 of file MOM_regridding.F90.

1137 
1138  ! Arguments
1139  type(regridding_CS), intent(in) :: CS !< Regridding control structure
1140  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1141  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1142  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1143  real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
1144  !! [H ~> m or kg m-2].
1145  real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage.
1146  ! Local variables
1147  integer :: i, j, k
1148  integer :: nz
1149  real :: nominalDepth, totalThickness, dh
1150  real, dimension(SZK_(GV)+1) :: zOld, zNew
1151  real :: minThickness
1152  logical :: ice_shelf
1153 
1154  nz = gv%ke
1155  minthickness = cs%min_thickness
1156  ice_shelf = .false.
1157  if (present(frac_shelf_h)) then
1158  if (associated(frac_shelf_h)) ice_shelf = .true.
1159  endif
1160 
1161 !$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h, &
1162 !$OMP ice_shelf,minThickness) &
1163 !$OMP private(nominalDepth,totalThickness, &
1164 !$OMP zNew,dh,zOld)
1165  do j = g%jsc-1,g%jec+1
1166  do i = g%isc-1,g%iec+1
1167 
1168  if (g%mask2dT(i,j)==0.) then
1169  dzinterface(i,j,:) = 0.
1170  cycle
1171  endif
1172 
1173  ! Local depth (G%bathyT is positive)
1174  nominaldepth = g%bathyT(i,j)*gv%Z_to_H
1175 
1176  ! Determine water column thickness
1177  totalthickness = 0.0
1178  do k = 1,nz
1179  totalthickness = totalthickness + h(i,j,k)
1180  enddo
1181 
1182  zold(nz+1) = - nominaldepth
1183  do k = nz,1,-1
1184  zold(k) = zold(k+1) + h(i,j,k)
1185  enddo
1186 
1187  if (ice_shelf) then
1188  if (frac_shelf_h(i,j) > 0.) then ! under ice shelf
1189  call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, znew, &
1190  z_rigid_top = totalthickness-nominaldepth, &
1191  eta_orig=zold(1), zscale=gv%Z_to_H)
1192  else
1193  call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, &
1194  znew, zscale=gv%Z_to_H)
1195  endif
1196  else
1197  call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, &
1198  znew, zscale=gv%Z_to_H)
1199  endif
1200 
1201  ! Calculate the final change in grid position after blending new and old grids
1202  call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )
1203 
1204 #ifdef __DO_SAFETY_CHECKS__
1205  dh=max(nominaldepth,totalthickness)
1206  if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh) then
1207  write(0,*) 'min_thickness=',minthickness
1208  write(0,*) 'nominalDepth=',nominaldepth,'totalThickness=',totalthickness
1209  write(0,*) 'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz
1210  do k=1,nz+1
1211  write(0,*) k,zold(k),znew(k)
1212  enddo
1213  do k=1,nz
1214  write(0,*) k,h(i,j,k),znew(k)-znew(k+1),cs%coordinateResolution(k)
1215  enddo
1216  call mom_error( fatal, &
1217  'MOM_regridding, build_zstar_grid(): top surface has moved!!!' )
1218  endif
1219 #endif
1220 
1221  call adjust_interface_motion( cs, nz, h(i,j,:), dzinterface(i,j,:) )
1222 
1223  enddo
1224  enddo
1225 

References adjust_interface_motion(), coord_zlike::build_zstar_column(), filtered_grid_motion(), and mom_error_handler::mom_error().

Referenced by regridding_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]csRegridding control structure
[in]gGrid structure
[in]gvOcean vertical grid structure
[in]hOld layer thicknesses (arbitrary units)
[in]dzinterfaceChange in interface positions (same as h)
[in,out]h_newNew layer thicknesses (same as h)

Definition at line 885 of file MOM_regridding.F90.

885  type(regridding_CS), intent(in) :: CS !< Regridding control structure
886  type(ocean_grid_type), intent(in) :: G !< Grid structure
887  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
888  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Old layer thicknesses (arbitrary units)
889  real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(in) :: dzInterface !< Change in interface positions (same as h)
890  real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses (same as h)
891  ! Local variables
892  integer :: i, j, k, nki
893 
894  nki = min(cs%nk, gv%ke)
895 
896  !$OMP parallel do default(shared)
897  do j = g%jsc-1,g%jec+1
898  do i = g%isc-1,g%iec+1
899  if (g%mask2dT(i,j)>0.) then
900  do k=1,nki
901  h_new(i,j,k) = max( 0., h(i,j,k) + ( dzinterface(i,j,k) - dzinterface(i,j,k+1) ) )
902  enddo
903  if (cs%nk > gv%ke) then
904  do k=nki+1, cs%nk
905  h_new(i,j,k) = max( 0., dzinterface(i,j,k) - dzinterface(i,j,k+1) )
906  enddo
907  endif
908  else
909  h_new(i,j,1:nki) = h(i,j,1:nki)
910  if (cs%nk > gv%ke) h_new(i,j,nki+1:cs%nk) = 0.
911  ! On land points, why are we keeping the original h rather than setting to zero? -AJA
912  endif
913  enddo
914  enddo
915 

Referenced by build_grid_hycom1(), and regridding_main().

Here is the caller graph for this function:

◆ 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]nkNumber of cells
[in]depthDepth of bottom [Z ~> m] or arbitrary units
[in]hCell thicknesses [Z ~> m] or arbitrary units
[in]dzinterfaceChange in interface positions (same units as h)
[in]msgMessage to append to errors

Definition at line 938 of file MOM_regridding.F90.

938  integer, intent(in) :: nk !< Number of cells
939  real, intent(in) :: depth !< Depth of bottom [Z ~> m] or arbitrary units
940  real, dimension(nk), intent(in) :: h !< Cell thicknesses [Z ~> m] or arbitrary units
941  real, dimension(nk+1), intent(in) :: dzInterface !< Change in interface positions (same units as h)
942  character(len=*), intent(in) :: msg !< Message to append to errors
943  ! Local variables
944  integer :: k
945  real :: eps, total_h_old, total_h_new, h_new, z_old, z_new
946 
947  eps =1. ; eps = epsilon(eps)
948 
949  ! Total thickness of grid h
950  total_h_old = 0.
951  do k = 1,nk
952  total_h_old = total_h_old + h(k)
953  enddo
954 
955  ! Integrate upwards for the interfaces consistent with the rest of MOM6
956  z_old = - depth
957  if (depth == 0.) z_old = - total_h_old
958  total_h_new = 0.
959  do k = nk,1,-1
960  z_old = z_old + h(k) ! Old interface position above layer k
961  z_new = z_old + dzinterface(k) ! New interface position based on dzInterface
962  h_new = h(k) + ( dzinterface(k) - dzinterface(k+1) ) ! New thickness
963  if (h_new<0.) then
964  write(0,*) 'k,h,hnew=',k,h(k),h_new
965  write(0,*) 'dzI(k+1),dzI(k)=',dzinterface(k+1),dzinterface(k)
966  call mom_error( fatal, 'MOM_regridding, check_grid_column: '//&
967  'Negative layer thickness implied by re-gridding, '//trim(msg))
968  endif
969  total_h_new = total_h_new + h_new
970 
971  enddo
972 
973  ! Conservation by implied h_new
974  if (abs(total_h_new-total_h_old)>real(nk-1)*0.5*(total_h_old+total_h_new)*eps) then
975  write(0,*) 'nk=',nk
976  do k = 1,nk
977  write(0,*) 'k,h,hnew=',k,h(k),h(k)+(dzinterface(k)-dzinterface(k+1))
978  enddo
979  write(0,*) 'Hold,Hnew,Hnew-Hold=',total_h_old,total_h_new,total_h_new-total_h_old
980  write(0,*) 'eps,(n)/2*eps*H=',eps,real(nk-1)*0.5*(total_h_old+total_h_new)*eps
981  call mom_error( fatal, 'MOM_regridding, check_grid_column: '//&
982  'Re-gridding did NOT conserve total thickness to within roundoff '//trim(msg))
983  endif
984 
985  ! Check that the top and bottom are intentionally moving
986  if (dzinterface(1) /= 0.) call mom_error( fatal, &
987  'MOM_regridding, check_grid_column: Non-zero dzInterface at surface! '//trim(msg))
988  if (dzinterface(nk+1) /= 0.) call mom_error( fatal, &
989  'MOM_regridding, check_grid_column: Non-zero dzInterface at bottom! '//trim(msg))
990 

References mom_error_handler::mom_error().

Referenced by check_remapping_grid().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]filenameFile name
[in]varnameVariable name
[in]expected_unitsExpected units of variable
[in,out]msgMessage to use for errors
[out]ierrTrue if an error occurs

Definition at line 716 of file MOM_regridding.F90.

716  character(len=*), intent(in) :: filename !< File name
717  character(len=*), intent(in) :: varname !< Variable name
718  character(len=*), intent(in) :: expected_units !< Expected units of variable
719  character(len=*), intent(inout) :: msg !< Message to use for errors
720  logical, intent(out) :: ierr !< True if an error occurs
721  ! Local variables
722  character (len=200) :: units, long_name
723  integer :: ncid, status, intid, vid
724  integer :: i
725 
726  ierr = .false.
727  status = nf90_open(trim(filename), nf90_nowrite, ncid)
728  if (status /= nf90_noerr) then
729  ierr = .true.
730  msg = 'File not found: '//trim(filename)
731  return
732  endif
733 
734  status = nf90_inq_varid(ncid, trim(varname), vid)
735  if (status /= nf90_noerr) then
736  ierr = .true.
737  msg = 'Var not found: '//trim(varname)
738  return
739  endif
740 
741  status = nf90_get_att(ncid, vid, "units", units)
742  if (status /= nf90_noerr) then
743  ierr = .true.
744  msg = 'Attribute not found: units'
745  return
746  endif
747  ! NF90_GET_ATT can return attributes with null characters, which TRIM will not truncate.
748  ! This loop replaces any null characters with a space so that the following check between
749  ! the read units and the expected units will pass
750  do i=1,len_trim(units)
751  if (units(i:i) == char(0)) units(i:i) = " "
752  enddo
753 
754  if (trim(units) /= trim(expected_units)) then
755  if (trim(expected_units) == "meters") then
756  if (trim(units) /= "m") then
757  ierr = .true.
758  endif
759  else
760  ierr = .true.
761  endif
762  endif
763 
764  if (ierr) then
765  msg = 'Units incorrect: '//trim(units)//' /= '//trim(expected_units)
766  endif
767 

Referenced by initialize_regridding().

Here is the caller graph for this function:

◆ 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]gGrid structure
[in]gvOcean vertical grid structure
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]dzinterfaceChange in interface positions [H ~> m or kg m-2]
[in]msgMessage to append to errors

Definition at line 920 of file MOM_regridding.F90.

920  type(ocean_grid_type), intent(in) :: G !< Grid structure
921  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
922  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
923  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzInterface !< Change in interface positions
924  !! [H ~> m or kg m-2]
925  character(len=*), intent(in) :: msg !< Message to append to errors
926  ! Local variables
927  integer :: i, j
928 
929  !$OMP parallel do default(shared)
930  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
931  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 )
932  enddo ; enddo
933 

References check_grid_column().

Referenced by regridding_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in,out]hLayer thicknesses [H ~> m or kg m-2]
[in,out]tvA structure pointing to various thermodynamic variables

Definition at line 1844 of file MOM_regridding.F90.

1844  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1845  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1846  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1847  intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
1848  type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
1849 !------------------------------------------------------------------------------
1850 ! Check each water column to see whether it is stratified. If not, sort the
1851 ! layers by successive swappings of water masses (bubble sort algorithm)
1852 !------------------------------------------------------------------------------
1853 
1854  ! Local variables
1855  integer :: i, j, k
1856  real :: T0, T1 ! temperatures
1857  real :: S0, S1 ! salinities
1858  real :: r0, r1 ! densities
1859  real :: h0, h1
1860  logical :: stratified
1861  real, dimension(GV%ke) :: p_col, densities
1862 
1863  p_col(:) = 0.
1864 
1865  ! Loop on columns
1866  do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1
1867 
1868  ! Compute densities within current water column
1869  call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, &
1870  densities, 1, gv%ke, tv%eqn_of_state )
1871 
1872  ! Repeat restratification until complete
1873  do
1874  stratified = .true.
1875  do k = 1,gv%ke-1
1876  ! Gather information of current and next cells
1877  t0 = tv%T(i,j,k) ; t1 = tv%T(i,j,k+1)
1878  s0 = tv%S(i,j,k) ; s1 = tv%S(i,j,k+1)
1879  r0 = densities(k) ; r1 = densities(k+1)
1880  h0 = h(i,j,k) ; h1 = h(i,j,k+1)
1881  ! If the density of the current cell is larger than the density
1882  ! below it, we swap the cells and recalculate the densitiies
1883  ! within the swapped cells
1884  if ( r0 > r1 ) then
1885  tv%T(i,j,k) = t1 ; tv%T(i,j,k+1) = t0
1886  tv%S(i,j,k) = s1 ; tv%S(i,j,k+1) = s0
1887  h(i,j,k) = h1 ; h(i,j,k+1) = h0
1888  ! Recompute densities at levels k and k+1
1889  call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), &
1890  densities(k), tv%eqn_of_state )
1891  call calculate_density( tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), &
1892  densities(k+1), tv%eqn_of_state )
1893  stratified = .false.
1894  endif
1895  enddo ! k
1896 
1897  if ( stratified ) exit
1898  enddo
1899 
1900  enddo ; enddo ! i & j
1901 

Referenced by regridding_main().

Here is the caller graph for this function:

◆ 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]stringString with list of parameters in form dz_min, H_total, power, precision
[in,out]dzProfile of nominal thicknesses

Definition at line 2380 of file MOM_regridding.F90.

2380  character(len=*), intent(in) :: string !< String with list of parameters in form
2381  !! dz_min, H_total, power, precision
2382  real, dimension(:), intent(inout) :: dz !< Profile of nominal thicknesses
2383  ! Local variables
2384  integer :: nk, k
2385  real :: dz_min, power, prec, H_total
2386 
2387  nk = size(dz) ! Number of cells
2388  prec = -1024.
2389  read( string, *) dz_min, h_total, power, prec
2390  if (prec == -1024.) call mom_error(fatal,"dz_function1: "// &
2391  "Problem reading FNC1: string ="//trim(string))
2392  ! Create profile of ( dz - dz_min )
2393  do k = 1, nk
2394  dz(k) = (real(k-1)/real(nk-1))**power
2395  enddo
2396  dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total
2397  dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
2398  dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total
2399  dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
2400  dz(nk) = dz(nk) + ( h_total - sum( dz(:) + dz_min ) ) ! Adjust bottommost layer
2401  dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
2402  dz(:) = dz(:) + dz_min ! Finally add in the constant dz_min
2403 

References mom_error_handler::mom_error().

Referenced by initialize_regridding().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ end_regridding()

subroutine, public mom_regridding::end_regridding ( type(regridding_cs), intent(inout)  CS)

Deallocation of regridding memory.

Parameters
[in,out]csRegridding control structure

Definition at line 772 of file MOM_regridding.F90.

772  type(regridding_CS), intent(inout) :: CS !< Regridding control structure
773 
774  if (associated(cs%zlike_CS)) call end_coord_zlike(cs%zlike_CS)
775  if (associated(cs%sigma_CS)) call end_coord_sigma(cs%sigma_CS)
776  if (associated(cs%rho_CS)) call end_coord_rho(cs%rho_CS)
777  if (associated(cs%hycom_CS)) call end_coord_hycom(cs%hycom_CS)
778  if (associated(cs%slight_CS)) call end_coord_slight(cs%slight_CS)
779  if (associated(cs%adapt_CS)) call end_coord_adapt(cs%adapt_CS)
780 
781  deallocate( cs%coordinateResolution )
782  if (allocated(cs%target_density)) deallocate( cs%target_density )
783  if (allocated(cs%max_interface_depths) ) deallocate( cs%max_interface_depths )
784  if (allocated(cs%max_layer_thickness) ) deallocate( cs%max_layer_thickness )
785 

References coord_adapt::end_coord_adapt(), coord_hycom::end_coord_hycom(), coord_sigma::end_coord_sigma(), coord_slight::end_coord_slight(), and coord_zlike::end_coord_zlike().

Here is the call graph for this function:

◆ 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]csRegridding control structure
[in]nkNumber of cells in source grid
[in]z_oldOld grid position [m]
[in]z_newNew grid position [m]
[in,out]dz_gChange in interface positions [m]

Definition at line 999 of file MOM_regridding.F90.

999  type(regridding_CS), intent(in) :: CS !< Regridding control structure
1000  integer, intent(in) :: nk !< Number of cells in source grid
1001  real, dimension(nk+1), intent(in) :: z_old !< Old grid position [m]
1002  real, dimension(CS%nk+1), intent(in) :: z_new !< New grid position [m]
1003  real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions [m]
1004  ! Local variables
1005  real :: sgn ! The sign convention for downward.
1006  real :: dz_tgt, zr1, z_old_k
1007  real :: Aq, Bq, dz0, z0, F0
1008  real :: zs, zd, dzwt, Idzwt
1009  real :: wtd, Iwtd
1010  real :: Int_zs, Int_zd, dInt_zs_zd
1011 ! For debugging:
1012  real, dimension(nk+1) :: z_act
1013 ! real, dimension(nk+1) :: ddz_g_s, ddz_g_d
1014  logical :: debug = .false.
1015  integer :: k
1016 
1017  if ((z_old(nk+1) - z_old(1)) * (z_new(cs%nk+1) - z_new(1)) < 0.0) then
1018  call mom_error(fatal, "filtered_grid_motion: z_old and z_new use different sign conventions.")
1019  elseif ((z_old(nk+1) - z_old(1)) * (z_new(cs%nk+1) - z_new(1)) == 0.0) then
1020  ! This is a massless column, so do nothing and return.
1021  do k=1,cs%nk+1 ; dz_g(k) = 0.0 ; enddo ; return
1022  elseif ((z_old(nk+1) - z_old(1)) + (z_new(cs%nk+1) - z_new(1)) > 0.0) then
1023  sgn = 1.0
1024  else
1025  sgn = -1.0
1026  endif
1027 
1028  if (debug) then
1029  do k=2,cs%nk+1
1030  if (sgn*(z_new(k)-z_new(k-1)) < -5e-16*(abs(z_new(k))+abs(z_new(k-1))) ) &
1031  call mom_error(fatal, "filtered_grid_motion: z_new is tangled.")
1032  enddo
1033  do k=2,nk+1
1034  if (sgn*(z_old(k)-z_old(k-1)) < -5e-16*(abs(z_old(k))+abs(z_old(k-1))) ) &
1035  call mom_error(fatal, "filtered_grid_motion: z_old is tangled.")
1036  enddo
1037  ! ddz_g_s(:) = 0.0 ; ddz_g_d(:) = 0.0
1038  endif
1039 
1040  zs = cs%depth_of_time_filter_shallow
1041  zd = cs%depth_of_time_filter_deep
1042  wtd = 1.0 - cs%old_grid_weight
1043  iwtd = 1.0 / wtd
1044 
1045  dzwt = (zd - zs)
1046  idzwt = 0.0 ; if (abs(zd - zs) > 0.0) idzwt = 1.0 / (zd - zs)
1047  dint_zs_zd = 0.5*(1.0 + iwtd) * (zd - zs)
1048  aq = 0.5*(iwtd - 1.0)
1049 
1050  dz_g(1) = 0.0
1051  z_old_k = z_old(1)
1052  do k = 2,cs%nk+1
1053  if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model
1054  ! zr1 is positive and increases with depth, and dz_tgt is positive downward.
1055  dz_tgt = sgn*(z_new(k) - z_old_k)
1056  zr1 = sgn*(z_old_k - z_old(1))
1057 
1058  ! First, handle the two simple and common cases that do not pass through
1059  ! the adjustment rate transition zone.
1060  if ((zr1 > zd) .and. (zr1 + wtd * dz_tgt > zd)) then
1061  dz_g(k) = sgn * wtd * dz_tgt
1062  elseif ((zr1 < zs) .and. (zr1 + dz_tgt < zs)) then
1063  dz_g(k) = sgn * dz_tgt
1064  else
1065  ! Find the new value by inverting the equation
1066  ! integral(0 to dz_new) Iwt(z) dz = dz_tgt
1067  ! This is trivial where Iwt is a constant, and agrees with the two limits above.
1068 
1069  ! Take test values at the transition points to figure out which segment
1070  ! the new value will be found in.
1071  if (zr1 >= zd) then
1072  int_zd = iwtd*(zd - zr1)
1073  int_zs = int_zd - dint_zs_zd
1074  elseif (zr1 <= zs) then
1075  int_zs = (zs - zr1)
1076  int_zd = dint_zs_zd + (zs - zr1)
1077  else
1078 ! Int_zd = (zd - zr1) * (Iwtd + 0.5*(1.0 - Iwtd) * (zd - zr1) / (zd - zs))
1079  int_zd = (zd - zr1) * (iwtd*(0.5*(zd+zr1) - zs) + 0.5*(zd - zr1)) * idzwt
1080  int_zs = (zs - zr1) * (0.5*iwtd * ((zr1 - zs)) + (zd - 0.5*(zr1+zs))) * idzwt
1081  ! It has been verified that Int_zs = Int_zd - dInt_zs_zd to within roundoff.
1082  endif
1083 
1084  if (dz_tgt >= int_zd) then ! The new location is in the deep, slow region.
1085  dz_g(k) = sgn * ((zd-zr1) + wtd*(dz_tgt - int_zd))
1086  elseif (dz_tgt <= int_zs) then ! The new location is in the shallow region.
1087  dz_g(k) = sgn * ((zs-zr1) + (dz_tgt - int_zs))
1088  else ! We need to solve a quadratic equation for z_new.
1089  ! For accuracy, do the integral from the starting depth or the nearest
1090  ! edge of the transition region. The results with each choice are
1091  ! mathematically equivalent, but differ in roundoff, and this choice
1092  ! should minimize the likelihood of inadvertently overlapping interfaces.
1093  if (zr1 <= zs) then ; dz0 = zs-zr1 ; z0 = zs ; f0 = dz_tgt - int_zs
1094  elseif (zr1 >= zd) then ; dz0 = zd-zr1 ; z0 = zd ; f0 = dz_tgt - int_zd
1095  else ; dz0 = 0.0 ; z0 = zr1 ; f0 = dz_tgt ; endif
1096 
1097  bq = (dzwt + 2.0*aq*(z0-zs))
1098  ! Solve the quadratic: Aq*(zn-z0)**2 + Bq*(zn-z0) - F0*dzwt = 0
1099  ! Note that b>=0, and the two terms in the standard form cancel for the right root.
1100  dz_g(k) = sgn * (dz0 + 2.0*f0*dzwt / (bq + sqrt(bq**2 + 4.0*aq*f0*dzwt) ))
1101 
1102 ! if (debug) then
1103 ! dz0 = zs-zr1 ; z0 = zs ; F0 = dz_tgt - Int_zs ; Bq = (dzwt + 2.0*Aq*(z0-zs))
1104 ! ddz_g_s(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k)
1105 ! dz0 = zd-zr1 ; z0 = zd ; F0 = dz_tgt - Int_zd ; Bq = (dzwt + 2.0*Aq*(z0-zs))
1106 ! ddz_g_d(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k)
1107 !
1108 ! if (abs(ddz_g_s(k)) > 1e-12*(abs(dz_g(k)) + abs(dz_g(k)+ddz_g_s(k)))) &
1109 ! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled (sc).")
1110 ! 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)))) &
1111 ! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled.")
1112 ! endif
1113  endif
1114 
1115  endif
1116  enddo
1117  !dz_g(CS%nk+1) = 0.0
1118 
1119  if (debug) then
1120  z_old_k = z_old(1)
1121  do k=1,cs%nk+1
1122  if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model
1123  z_act(k) = z_old_k + dz_g(k)
1124  enddo
1125  do k=2,cs%nk+1
1126  if (sgn*((z_act(k))-z_act(k-1)) < -1e-15*(abs(z_act(k))+abs(z_act(k-1))) ) &
1127  call mom_error(fatal, "filtered_grid_motion: z_output is tangled.")
1128  enddo
1129  endif
1130 

References mom_error_handler::mom_error().

Referenced by build_grid_adaptive(), build_grid_hycom1(), build_grid_slight(), build_rho_grid(), build_sigma_grid(), and build_zstar_grid().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]csRegridding control structure

Definition at line 2310 of file MOM_regridding.F90.

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

◆ 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]csRegridding control structure

Definition at line 2334 of file MOM_regridding.F90.

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

Referenced by mom_diag_remap::diag_remap_update().

Here is the caller graph for this function:

◆ 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]csRegridding control structure

Definition at line 2326 of file MOM_regridding.F90.

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

Referenced by mom_diag_remap::diag_remap_update().

Here is the caller graph for this function:

◆ 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]csRegridding control structure

Definition at line 2318 of file MOM_regridding.F90.

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

Referenced by mom_diag_remap::diag_remap_update().

Here is the caller graph for this function:

◆ 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]csRegridding control structure
[in]undo_scalingIf present and true, undo any internal rescaling of the resolution data.
Returns
Interface positions in target coordinate

Definition at line 2113 of file MOM_regridding.F90.

2113  type(regridding_CS), intent(in) :: CS !< Regridding control structure
2114  logical, optional, intent(in) :: undo_scaling !< If present and true, undo any internal
2115  !! rescaling of the resolution data.
2116  real, dimension(CS%nk+1) :: getCoordinateInterfaces !< Interface positions in target coordinate
2117 
2118  integer :: k
2119  logical :: unscale
2120  unscale = .false. ; if (present(undo_scaling)) unscale = undo_scaling
2121 
2122  ! When using a coordinate with target densities, we need to get the actual
2123  ! densities, rather than computing the interfaces based on resolution
2124  if (cs%regridding_scheme == regridding_rho) then
2125  if (.not. cs%target_density_set) &
2126  call mom_error(fatal, 'MOM_regridding, getCoordinateInterfaces: '//&
2127  'target densities not set!')
2128 
2129  if (unscale) then
2130  getcoordinateinterfaces(:) = cs%coord_scale * cs%target_density(:)
2131  else
2132  getcoordinateinterfaces(:) = cs%target_density(:)
2133  endif
2134  else
2135  if (unscale) then
2136  getcoordinateinterfaces(1) = 0.
2137  do k = 1, cs%nk
2138  getcoordinateinterfaces(k+1) = getcoordinateinterfaces(k) - &
2139  cs%coord_scale * cs%coordinateResolution(k)
2140  enddo
2141  else
2142  getcoordinateinterfaces(1) = 0.
2143  do k = 1, cs%nk
2144  getcoordinateinterfaces(k+1) = getcoordinateinterfaces(k) - &
2145  cs%coordinateResolution(k)
2146  enddo
2147  endif
2148  ! The following line has an "abs()" to allow ferret users to reference
2149  ! data by index. It is a temporary work around... :( -AJA
2150  getcoordinateinterfaces(:) = abs( getcoordinateinterfaces(:) )
2151  endif
2152 

References mom_error_handler::mom_error().

Here is the call graph for this function:

◆ 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]csRegridding control structure
[in]undo_scalingIf present and true, undo any internal rescaling of the resolution data.

Definition at line 2095 of file MOM_regridding.F90.

2095  type(regridding_CS), intent(in) :: CS !< Regridding control structure
2096  logical, optional, intent(in) :: undo_scaling !< If present and true, undo any internal
2097  !! rescaling of the resolution data.
2098  real, dimension(CS%nk) :: getCoordinateResolution
2099 
2100  logical :: unscale
2101  unscale = .false. ; if (present(undo_scaling)) unscale = undo_scaling
2102 
2103  if (unscale) then
2104  getcoordinateresolution(:) = cs%coord_scale * cs%coordinateResolution(:)
2105  else
2106  getcoordinateresolution(:) = cs%coordinateResolution(:)
2107  endif
2108 

◆ getcoordinateshortname()

character(len=20) function, public mom_regridding::getcoordinateshortname ( type(regridding_cs), intent(in)  CS)

Query the short name of the coordinate.

Parameters
[in]csRegridding control structure

Definition at line 2182 of file MOM_regridding.F90.

2182  type(regridding_CS), intent(in) :: CS !< Regridding control structure
2183  character(len=20) :: getCoordinateShortName
2184 
2185  select case ( cs%regridding_scheme )
2186  case ( regridding_zstar )
2187  !getCoordinateShortName = 'z*'
2188  ! The following line is a temporary work around... :( -AJA
2189  getcoordinateshortname = 'pseudo-depth, -z*'
2190  case ( regridding_sigma_shelf_zstar )
2191  getcoordinateshortname = 'pseudo-depth, -z*/sigma'
2192  case ( regridding_sigma )
2193  getcoordinateshortname = 'sigma'
2194  case ( regridding_rho )
2195  getcoordinateshortname = 'rho'
2196  case ( regridding_arbitrary )
2197  getcoordinateshortname = 'coordinate'
2198  case ( regridding_hycom1 )
2199  getcoordinateshortname = 'z-rho'
2200  case ( regridding_slight )
2201  getcoordinateshortname = 's-rho'
2202  case ( regridding_adaptive )
2203  getcoordinateshortname = 'adaptive'
2204  case default
2205  call mom_error(fatal,'MOM_regridding, getCoordinateShortName: '//&
2206  'Unknown regridding scheme selected!')
2207  end select ! type of grid
2208 

References mom_error_handler::mom_error(), regrid_consts::regridding_adaptive, regrid_consts::regridding_hycom1, and regrid_consts::regridding_slight.

Here is the call graph for this function:

◆ getcoordinateunits()

character(len=20) function, public mom_regridding::getcoordinateunits ( type(regridding_cs), intent(in)  CS)

Query the target coordinate units.

Parameters
[in]csRegridding control structure

Definition at line 2158 of file MOM_regridding.F90.

2158  type(regridding_CS), intent(in) :: CS !< Regridding control structure
2159  character(len=20) :: getCoordinateUnits
2160 
2161  select case ( cs%regridding_scheme )
2162  case ( regridding_zstar, regridding_hycom1, regridding_slight, regridding_adaptive )
2163  getcoordinateunits = 'meter'
2164  case ( regridding_sigma_shelf_zstar )
2165  getcoordinateunits = 'meter/fraction'
2166  case ( regridding_sigma )
2167  getcoordinateunits = 'fraction'
2168  case ( regridding_rho )
2169  getcoordinateunits = 'kg/m3'
2170  case ( regridding_arbitrary )
2171  getcoordinateunits = 'unknown'
2172  case default
2173  call mom_error(fatal,'MOM_regridding, getCoordinateUnits: '//&
2174  'Unknown regridding scheme selected!')
2175  end select ! type of grid
2176 

References mom_error_handler::mom_error(), regrid_consts::regridding_adaptive, regrid_consts::regridding_hycom1, and regrid_consts::regridding_slight.

Here is the call graph for this function:

◆ 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]csRegridding control structure
[in]sshThe sea surface height, in the same units as depth
[in]depthThe maximum depth of the grid, often [Z ~> m]
Returns
The returned thicknesses in the units of depth

Definition at line 2343 of file MOM_regridding.F90.

2343  type(regridding_CS), intent(in) :: CS !< Regridding control structure
2344  real, intent(in) :: SSH !< The sea surface height, in the same units as depth
2345  real, intent(in) :: depth !< The maximum depth of the grid, often [Z ~> m]
2346  real, dimension(CS%nk) :: getStaticThickness !< The returned thicknesses in the units of depth
2347  ! Local
2348  integer :: k
2349  real :: z, dz
2350 
2351  select case ( cs%regridding_scheme )
2352  case ( regridding_zstar, regridding_sigma_shelf_zstar, regridding_hycom1, regridding_slight, regridding_adaptive )
2353  if (depth>0.) then
2354  z = ssh
2355  do k = 1, cs%nk
2356  dz = cs%coordinateResolution(k) * ( 1. + ssh/depth ) ! Nominal dz*
2357  dz = max(dz, 0.) ! Avoid negative incase ssh=-depth
2358  dz = min(dz, depth - z) ! Clip if below topography
2359  z = z + dz ! Bottom of layer
2360  getstaticthickness(k) = dz
2361  enddo
2362  else
2363  getstaticthickness(:) = 0. ! On land ...
2364  endif
2365  case ( regridding_sigma )
2366  getstaticthickness(:) = cs%coordinateResolution(:) * ( depth + ssh )
2367  case ( regridding_rho )
2368  getstaticthickness(:) = 0. ! Not applicable
2369  case ( regridding_arbitrary )
2370  getstaticthickness(:) = 0. ! Not applicable
2371  case default
2372  call mom_error(fatal,'MOM_regridding, getStaticThickness: '//&
2373  'Unknown regridding scheme selected!')
2374  end select ! type of grid
2375 

References mom_error_handler::mom_error(), regrid_consts::regridding_adaptive, regrid_consts::regridding_hycom1, and regrid_consts::regridding_slight.

Referenced by mom_ale::ale_initthicknesstocoord().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]csRegridding control structure
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in,out]hLayer thicknesses [H ~> m or kg m-2]

Definition at line 1803 of file MOM_regridding.F90.

1803 !------------------------------------------------------------------------------
1804 ! This routine is called when initializing the regridding options. The
1805 ! objective is to make sure all layers are at least as thick as the minimum
1806 ! thickness allowed for regridding purposes (this parameter is set in the
1807 ! MOM_input file or defaulted to 1.0e-3). When layers are too thin, they
1808 ! are inflated up to the minmum thickness.
1809 !------------------------------------------------------------------------------
1810 
1811  ! Arguments
1812  type(regridding_CS), intent(in) :: CS !< Regridding control structure
1813  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1814  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1815  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
1816 
1817  ! Local variables
1818  integer :: i, j, k
1819  real :: hTmp(GV%ke)
1820 
1821  do i = g%isc-1,g%iec+1
1822  do j = g%jsc-1,g%jec+1
1823 
1824  ! Build grid for current column
1825  do k = 1,gv%ke
1826  htmp(k) = h(i,j,k)
1827  enddo
1828 
1829  call old_inflate_layers_1d( cs%min_thickness, gv%ke, htmp )
1830 
1831  ! Save modified grid
1832  do k = 1,gv%ke
1833  h(i,j,k) = htmp(k)
1834  enddo
1835 
1836  enddo
1837  enddo
1838 

References coord_rho::old_inflate_layers_1d().

Here is the call graph for this function:

◆ 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]csRegridding control structure
[in]coord_modeA string indicating the coordinate mode. See the documenttion for regrid_consts for the recognized values.
[in]gvOcean vertical grid structure
[in]usA dimensional unit scaling type

Definition at line 1949 of file MOM_regridding.F90.

1949  type(regridding_CS), intent(inout) :: CS !< Regridding control structure
1950  character(len=*), intent(in) :: coord_mode !< A string indicating the coordinate mode.
1951  !! See the documenttion for regrid_consts
1952  !! for the recognized values.
1953  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1954  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1955 
1956  select case (coordinatemode(coord_mode))
1957  case (regridding_zstar)
1958  call init_coord_zlike(cs%zlike_CS, cs%nk, cs%coordinateResolution)
1959  case (regridding_sigma_shelf_zstar)
1960  call init_coord_zlike(cs%zlike_CS, cs%nk, cs%coordinateResolution)
1961  case (regridding_sigma)
1962  call init_coord_sigma(cs%sigma_CS, cs%nk, cs%coordinateResolution)
1963  case (regridding_rho)
1964  call init_coord_rho(cs%rho_CS, cs%nk, cs%ref_pressure, cs%target_density, cs%interp_CS, &
1965  rho_scale=us%kg_m3_to_R)
1966  case (regridding_hycom1)
1967  call init_coord_hycom(cs%hycom_CS, cs%nk, cs%coordinateResolution, cs%target_density, &
1968  cs%interp_CS, rho_scale=us%kg_m3_to_R)
1969  case (regridding_slight)
1970  call init_coord_slight(cs%slight_CS, cs%nk, cs%ref_pressure, cs%target_density, &
1971  cs%interp_CS, gv%m_to_H, rho_scale=us%kg_m3_to_R)
1972  case (regridding_adaptive)
1973  call init_coord_adapt(cs%adapt_CS, cs%nk, cs%coordinateResolution, gv%m_to_H)
1974  end select

References coord_adapt::init_coord_adapt(), coord_hycom::init_coord_hycom(), coord_sigma::init_coord_sigma(), coord_slight::init_coord_slight(), coord_zlike::init_coord_zlike(), regrid_consts::regridding_adaptive, regrid_consts::regridding_hycom1, and regrid_consts::regridding_slight.

Referenced by initialize_regridding().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]csRegridding control structure
[in]gvOcean vertical grid structure
[in]usA dimensional unit scaling type
[in]max_depthThe maximum depth of the ocean [Z ~> m].
[in]param_fileParameter file
[in]mdlName of calling module.
[in]coord_modeCoordinate mode
[in]param_prefixString to prefix to parameter names. If empty, causes main model parameters to be used.
[in]param_suffixString to append to parameter names.

Definition at line 177 of file MOM_regridding.F90.

177  type(regridding_CS), intent(inout) :: CS !< Regridding control structure
178  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
179  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
180  real, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m].
181  type(param_file_type), intent(in) :: param_file !< Parameter file
182  character(len=*), intent(in) :: mdl !< Name of calling module.
183  character(len=*), intent(in) :: coord_mode !< Coordinate mode
184  character(len=*), intent(in) :: param_prefix !< String to prefix to parameter names.
185  !! If empty, causes main model parameters to be used.
186  character(len=*), intent(in) :: param_suffix !< String to append to parameter names.
187 
188  ! Local variables
189  integer :: ke ! Number of levels
190  character(len=80) :: string, string2, varName ! Temporary strings
191  character(len=40) :: coord_units, param_name, coord_res_param ! Temporary strings
192  character(len=200) :: inputdir, fileName
193  character(len=320) :: message ! Temporary strings
194  character(len=12) :: expected_units ! Temporary strings
195  logical :: tmpLogical, fix_haloclines, set_max, do_sum, main_parameters
196  logical :: coord_is_state_dependent, ierr
197  real :: filt_len, strat_tol, index_scale, tmpReal
198  real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z).
199  real :: dz_fixed_sfc, Rho_avg_depth, nlay_sfc_int
200  real :: adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha
201  integer :: nz_fixed_sfc, k, nzf(4)
202  real, dimension(:), allocatable :: dz ! Resolution (thickness) in units of coordinate, which may be [m]
203  ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units.
204  real, dimension(:), allocatable :: h_max ! Maximum layer thicknesses [H ~> m or kg m-2]
205  real, dimension(:), allocatable :: z_max ! Maximum interface depths [H ~> m or kg m-2] or other
206  ! units depending on the coordinate
207  real, dimension(:), allocatable :: dz_max ! Thicknesses used to find maximum interface depths
208  ! [H ~> m or kg m-2] or other units
209  real, dimension(:), allocatable :: rho_target ! Target density used in HYBRID mode [kg m-3]
210  ! Thicknesses [m] that give level centers corresponding to table 2 of WOA09
211  real, dimension(40) :: woa09_dz = (/ 5., 10., 10., 15., 22.5, 25., 25., 25., &
212  37.5, 50., 50., 75., 100., 100., 100., 100., &
213  100., 100., 100., 100., 100., 100., 100., 175., &
214  250., 375., 500., 500., 500., 500., 500., 500., &
215  500., 500., 500., 500., 500., 500., 500., 500. /)
216 
217  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
218  inputdir = slasher(inputdir)
219 
220  main_parameters=.false.
221  if (len_trim(param_prefix)==0) main_parameters=.true.
222  if (main_parameters .and. len_trim(param_suffix)>0) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
223  'Suffix provided without prefix for parameter names!')
224 
225  cs%nk = 0
226  cs%regridding_scheme = coordinatemode(coord_mode)
227  coord_is_state_dependent = state_dependent(coord_mode)
228  maximum_depth = us%Z_to_m*max_depth
229 
230  if (main_parameters) then
231  ! Read coordinate units parameter (main model = REGRIDDING_COORDINATE_UNITS)
232  call get_param(param_file, mdl, "REGRIDDING_COORDINATE_UNITS", coord_units, &
233  "Units of the regridding coordinate.", default=coordinateunits(coord_mode))
234  else
235  coord_units=coordinateunits(coord_mode)
236  endif
237 
238  if (coord_is_state_dependent) then
239  if (main_parameters) then
240  param_name = "INTERPOLATION_SCHEME"
241  string2 = regriddingdefaultinterpscheme
242  else
243  param_name = trim(param_prefix)//"_INTERP_SCHEME_"//trim(param_suffix)
244  string2 = 'PPM_H4' ! Default for diagnostics
245  endif
246  call get_param(param_file, mdl, "INTERPOLATION_SCHEME", string, &
247  "This sets the interpolation scheme to use to "//&
248  "determine the new grid. These parameters are "//&
249  "only relevant when REGRIDDING_COORDINATE_MODE is "//&
250  "set to a function of state. Otherwise, it is not "//&
251  "used. It can be one of the following schemes: "//&
252  trim(regriddinginterpschemedoc), default=trim(string2))
253  call set_regrid_params(cs, interp_scheme=string)
254  endif
255 
256  if (main_parameters .and. coord_is_state_dependent) then
257  call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", tmplogical, &
258  "When defined, a proper high-order reconstruction "//&
259  "scheme is used within boundary cells rather "//&
260  "than PCM. E.g., if PPM is used for remapping, a "//&
261  "PPM reconstruction will also be used within "//&
262  "boundary cells.", default=regriddingdefaultboundaryextrapolation)
263  call set_regrid_params(cs, boundary_extrapolation=tmplogical)
264  else
265  call set_regrid_params(cs, boundary_extrapolation=.false.)
266  endif
267 
268  ! Read coordinate configuration parameter (main model = ALE_COORDINATE_CONFIG)
269  if (main_parameters) then
270  param_name = "ALE_COORDINATE_CONFIG"
271  coord_res_param = "ALE_RESOLUTION"
272  string2 = 'UNIFORM'
273  else
274  param_name = trim(param_prefix)//"_DEF_"//trim(param_suffix)
275  coord_res_param = trim(param_prefix)//"_RES_"//trim(param_suffix)
276  string2 = 'UNIFORM'
277  if (maximum_depth>3000.) string2='WOA09' ! For convenience
278  endif
279  call get_param(param_file, mdl, param_name, string, &
280  "Determines how to specify the coordinate "//&
281  "resolution. Valid options are:\n"//&
282  " PARAM - use the vector-parameter "//trim(coord_res_param)//"\n"//&
283  " UNIFORM[:N] - uniformly distributed\n"//&
284  " FILE:string - read from a file. The string specifies\n"//&
285  " the filename and variable name, separated\n"//&
286  " by a comma or space, e.g. FILE:lev.nc,dz\n"//&
287  " or FILE:lev.nc,interfaces=zw\n"//&
288  " WOA09[:N] - the WOA09 vertical grid (approximately)\n"//&
289  " FNC1:string - FNC1:dz_min,H_total,power,precision\n"//&
290  " HYBRID:string - read from a file. The string specifies\n"//&
291  " the filename and two variable names, separated\n"//&
292  " by a comma or space, for sigma-2 and dz. e.g.\n"//&
293  " HYBRID:vgrid.nc,sigma2,dz",&
294  default=trim(string2))
295  message = "The distribution of vertical resolution for the target\n"//&
296  "grid used for Eulerian-like coordinates. For example,\n"//&
297  "in z-coordinate mode, the parameter is a list of level\n"//&
298  "thicknesses (in m). In sigma-coordinate mode, the list\n"//&
299  "is of non-dimensional fractions of the water column."
300  if (index(trim(string),'UNIFORM')==1) then
301  if (len_trim(string)==7) then
302  ke = gv%ke ! Use model nk by default
303  tmpreal = maximum_depth
304  elseif (index(trim(string),'UNIFORM:')==1 .and. len_trim(string)>8) then
305  ! Format is "UNIFORM:N" or "UNIFORM:N,dz"
306  ke = extract_integer(string(9:len_trim(string)),'',1)
307  tmpreal = extract_real(string(9:len_trim(string)),',',2,missing_value=maximum_depth)
308  else
309  call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
310  'Unable to interpret "'//trim(string)//'".')
311  endif
312  allocate(dz(ke))
313  dz(:) = uniformresolution(ke, coord_mode, tmpreal, &
314  us%R_to_kg_m3*(gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(min(2,ke)))), &
315  us%R_to_kg_m3*(gv%Rlay(ke) + 0.5*(gv%Rlay(ke)-gv%Rlay(max(ke-1,1)))) )
316  if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
317  trim(message), units=trim(coord_units))
318  elseif (trim(string)=='PARAM') then
319  ! Read coordinate resolution (main model = ALE_RESOLUTION)
320  ke = gv%ke ! Use model nk by default
321  allocate(dz(ke))
322  call get_param(param_file, mdl, coord_res_param, dz, &
323  trim(message), units=trim(coord_units), fail_if_missing=.true.)
324  elseif (index(trim(string),'FILE:')==1) then
325  ! FILE:filename,var_name is assumed to be reading level thickness variables
326  ! FILE:filename,interfaces=var_name reads positions
327  if (string(6:6)=='.' .or. string(6:6)=='/') then
328  ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
329  filename = trim( extractword(trim(string(6:80)), 1) )
330  else
331  ! Otherwise assume we should look for the file in INPUTDIR
332  filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
333  endif
334  if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
335  "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
336 
337  varname = trim( extractword(trim(string(6:)), 2) )
338  if (len_trim(varname)==0) then
339  if (field_exists(filename,'dz')) then; varname = 'dz'
340  elseif (field_exists(filename,'dsigma')) then; varname = 'dsigma'
341  elseif (field_exists(filename,'ztest')) then; varname = 'ztest'
342  else ; call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
343  "Coordinate variable not specified and none could be guessed.")
344  endif
345  endif
346  ! This check fails when the variable is a dimension variable! -AJA
347  !if (.not. field_exists(fileName,trim(varName))) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
348  ! "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
349  if (cs%regridding_scheme == regridding_sigma) then
350  expected_units = 'nondim'
351  elseif (cs%regridding_scheme == regridding_rho) then
352  expected_units = 'kg m-3'
353  else
354  expected_units = 'meters'
355  endif
356  if (index(trim(varname),'interfaces=')==1) then
357  varname=trim(varname(12:))
358  call check_grid_def(filename, varname, expected_units, message, ierr)
359  if (ierr) call mom_error(fatal,trim(mdl)//", initialize_regridding: "//&
360  "Unsupported format in grid definition '"//trim(filename)//"'. Error message "//trim(message))
361  call field_size(trim(filename), trim(varname), nzf)
362  ke = nzf(1)-1
363  if (cs%regridding_scheme == regridding_rho) then
364  allocate(rho_target(ke+1))
365  call mom_read_data(trim(filename), trim(varname), rho_target)
366  else
367  allocate(dz(ke))
368  allocate(z_max(ke+1))
369  call mom_read_data(trim(filename), trim(varname), z_max)
370  dz(:) = abs(z_max(1:ke) - z_max(2:ke+1))
371  deallocate(z_max)
372  endif
373  else
374  ! Assume reading resolution
375  call field_size(trim(filename), trim(varname), nzf)
376  ke = nzf(1)
377  allocate(dz(ke))
378  call mom_read_data(trim(filename), trim(varname), dz)
379  endif
380  if (main_parameters .and. ke/=gv%ke) then
381  call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
382  'Mismatch in number of model levels and "'//trim(string)//'".')
383  endif
384  if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
385  trim(message), units=coordinateunits(coord_mode))
386  elseif (index(trim(string),'FNC1:')==1) then
387  ke = gv%ke; allocate(dz(ke))
388  call dz_function1( trim(string(6:)), dz )
389  if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
390  trim(message), units=coordinateunits(coord_mode))
391  elseif (index(trim(string),'RFNC1:')==1) then
392  ! Function used for set target interface densities
393  ke = rho_function1( trim(string(7:)), rho_target )
394  elseif (index(trim(string),'HYBRID:')==1) then
395  ke = gv%ke; allocate(dz(ke))
396  ! The following assumes the FILE: syntax of above but without "FILE:" in the string
397  allocate(rho_target(ke+1))
398  filename = trim( extractword(trim(string(8:)), 1) )
399  if (filename(1:1)/='.' .and. filename(1:1)/='/') filename = trim(inputdir) // trim( filename )
400  if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: HYBRID "// &
401  "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
402  varname = trim( extractword(trim(string(8:)), 2) )
403  if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)//", initialize_regridding: HYBRID "// &
404  "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
405  call mom_read_data(trim(filename), trim(varname), rho_target)
406  varname = trim( extractword(trim(string(8:)), 3) )
407  if (varname(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz
408  call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz )
409  else ! Read dz from file
410  if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)//", initialize_regridding: HYBRID "// &
411  "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
412  call mom_read_data(trim(filename), trim(varname), dz)
413  endif
414  if (main_parameters) then
415  call log_param(param_file, mdl, "!"//coord_res_param, dz, &
416  trim(message), units=coordinateunits(coord_mode))
417  call log_param(param_file, mdl, "!TARGET_DENSITIES", rho_target, &
418  'HYBRID target densities for interfaces', units=coordinateunits(coord_mode))
419  endif
420  elseif (index(trim(string),'WOA09')==1) then
421  if (len_trim(string)==5) then
422  tmpreal = 0. ; ke = 0
423  do while (tmpreal<maximum_depth)
424  ke = ke + 1
425  tmpreal = tmpreal + woa09_dz(ke)
426  enddo
427  elseif (index(trim(string),'WOA09:')==1) then
428  if (len_trim(string)==6) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
429  'Expected string of form "WOA09:N" but got "'//trim(string)//'".')
430  ke = extract_integer(string(7:len_trim(string)),'',1)
431  endif
432  if (ke>40 .or. ke<1) call mom_error(fatal,trim(mdl)//', initialize_regridding: '// &
433  'For "WOA05:N" N must 0<N<41 but got "'//trim(string)//'".')
434  allocate(dz(ke))
435  dz(1:ke) = woa09_dz(1:ke)
436  else
437  call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
438  "Unrecognized coordinate configuration"//trim(string))
439  endif
440 
441  if (main_parameters) then
442  ! This is a work around to apparently needed to work with the from_Z initialization... ???
443  if (coordinatemode(coord_mode) == regridding_zstar .or. &
444  coordinatemode(coord_mode) == regridding_hycom1 .or. &
445  coordinatemode(coord_mode) == regridding_slight .or. &
446  coordinatemode(coord_mode) == regridding_adaptive) then
447  ! Adjust target grid to be consistent with maximum_depth
448  tmpreal = sum( dz(:) )
449  if (tmpreal < maximum_depth) then
450  dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
451  elseif (tmpreal > maximum_depth) then
452  if ( dz(ke) + ( maximum_depth - tmpreal ) > 0. ) then
453  dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
454  else
455  call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
456  "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string))
457  endif
458  endif
459  endif
460  endif
461 
462  cs%nk=ke
463 
464  ! Target resolution (for fixed coordinates)
465  allocate( cs%coordinateResolution(cs%nk) ); cs%coordinateResolution(:) = -1.e30
466  if (state_dependent(cs%regridding_scheme)) then
467  ! Target values
468  allocate( cs%target_density(cs%nk+1) ); cs%target_density(:) = -1.e30*us%kg_m3_to_R
469  endif
470 
471  if (allocated(dz)) then
472  if (coordinatemode(coord_mode) == regridding_sigma) then
473  call setcoordinateresolution(dz, cs, scale=1.0)
474  elseif (coordinatemode(coord_mode) == regridding_rho) then
475  call setcoordinateresolution(dz, cs, scale=us%kg_m3_to_R)
476  cs%coord_scale = us%R_to_kg_m3
477  elseif (coordinatemode(coord_mode) == regridding_adaptive) then
478  call setcoordinateresolution(dz, cs, scale=gv%m_to_H)
479  cs%coord_scale = gv%H_to_m
480  else
481  call setcoordinateresolution(dz, cs, scale=us%m_to_Z)
482  cs%coord_scale = us%Z_to_m
483  endif
484  endif
485 
486  if (allocated(rho_target)) then
487  call set_target_densities(cs, us%kg_m3_to_R*rho_target)
488  deallocate(rho_target)
489 
490  ! \todo This line looks like it would overwrite the target densities set just above?
491  elseif (coordinatemode(coord_mode) == regridding_rho) then
492  call set_target_densities_from_gv(gv, us, cs)
493  call log_param(param_file, mdl, "!TARGET_DENSITIES", us%R_to_kg_m3*cs%target_density(:), &
494  'RHO target densities for interfaces', units=coordinateunits(coord_mode))
495  endif
496 
497  ! initialise coordinate-specific control structure
498  call initcoord(cs, gv, us, coord_mode)
499 
500  if (main_parameters .and. coord_is_state_dependent) then
501  call get_param(param_file, mdl, "REGRID_COMPRESSIBILITY_FRACTION", tmpreal, &
502  "When interpolating potential density profiles we can add "//&
503  "some artificial compressibility solely to make homogeneous "//&
504  "regions appear stratified.", default=0.)
505  call set_regrid_params(cs, compress_fraction=tmpreal)
506  endif
507 
508  if (main_parameters) then
509  call get_param(param_file, mdl, "MIN_THICKNESS", tmpreal, &
510  "When regridding, this is the minimum layer "//&
511  "thickness allowed.", units="m", scale=gv%m_to_H, &
512  default=regriddingdefaultminthickness )
513  call set_regrid_params(cs, min_thickness=tmpreal)
514  else
515  call set_regrid_params(cs, min_thickness=0.)
516  endif
517 
518  if (coordinatemode(coord_mode) == regridding_slight) then
519  ! Set SLight-specific regridding parameters.
520  call get_param(param_file, mdl, "SLIGHT_DZ_SURFACE", dz_fixed_sfc, &
521  "The nominal thickness of fixed thickness near-surface "//&
522  "layers with the SLight coordinate.", units="m", default=1.0, scale=gv%m_to_H)
523  call get_param(param_file, mdl, "SLIGHT_NZ_SURFACE_FIXED", nz_fixed_sfc, &
524  "The number of fixed-depth surface layers with the SLight "//&
525  "coordinate.", units="nondimensional", default=2)
526  call get_param(param_file, mdl, "SLIGHT_SURFACE_AVG_DEPTH", rho_avg_depth, &
527  "The thickness of the surface region over which to average "//&
528  "when calculating the density to use to define the interior "//&
529  "with the SLight coordinate.", units="m", default=1.0, scale=gv%m_to_H)
530  call get_param(param_file, mdl, "SLIGHT_NLAY_TO_INTERIOR", nlay_sfc_int, &
531  "The number of layers to offset the surface density when "//&
532  "defining where the interior ocean starts with SLight.", &
533  units="nondimensional", default=2.0)
534  call get_param(param_file, mdl, "SLIGHT_FIX_HALOCLINES", fix_haloclines, &
535  "If true, identify regions above the reference pressure "//&
536  "where the reference pressure systematically underestimates "//&
537  "the stratification and use this in the definition of the "//&
538  "interior with the SLight coordinate.", default=.false.)
539 
540  call set_regrid_params(cs, dz_min_surface=dz_fixed_sfc, &
541  nz_fixed_surface=nz_fixed_sfc, rho_ml_avg_depth=rho_avg_depth, &
542  nlay_ml_to_interior=nlay_sfc_int, fix_haloclines=fix_haloclines)
543  if (fix_haloclines) then
544  ! Set additional parameters related to SLIGHT_FIX_HALOCLINES.
545  call get_param(param_file, mdl, "HALOCLINE_FILTER_LENGTH", filt_len, &
546  "A length scale over which to smooth the temperature and "//&
547  "salinity before identifying erroneously unstable haloclines.", &
548  units="m", default=2.0)
549  call get_param(param_file, mdl, "HALOCLINE_STRAT_TOL", strat_tol, &
550  "A tolerance for the ratio of the stratification of the "//&
551  "apparent coordinate stratification to the actual value "//&
552  "that is used to identify erroneously unstable haloclines. "//&
553  "This ratio is 1 when they are equal, and sensible values "//&
554  "are between 0 and 0.5.", units="nondimensional", default=0.2)
555  call set_regrid_params(cs, halocline_filt_len=filt_len, &
556  halocline_strat_tol=strat_tol)
557  endif
558 
559  endif
560 
561  if (coordinatemode(coord_mode) == regridding_adaptive) then
562  call get_param(param_file, mdl, "ADAPT_TIME_RATIO", adapttimeratio, &
563  "Ratio of ALE timestep to grid timescale.", units="s", default=1e-1) !### Should the units be "nondim"?
564  call get_param(param_file, mdl, "ADAPT_ZOOM_DEPTH", adaptzoom, &
565  "Depth of near-surface zooming region.", units="m", default=200.0, scale=gv%m_to_H)
566  call get_param(param_file, mdl, "ADAPT_ZOOM_COEFF", adaptzoomcoeff, &
567  "Coefficient of near-surface zooming diffusivity.", &
568  units="nondim", default=0.2)
569  call get_param(param_file, mdl, "ADAPT_BUOY_COEFF", adaptbuoycoeff, &
570  "Coefficient of buoyancy diffusivity.", &
571  units="nondim", default=0.8)
572  call get_param(param_file, mdl, "ADAPT_ALPHA", adaptalpha, &
573  "Scaling on optimization tendency.", &
574  units="nondim", default=1.0)
575  call get_param(param_file, mdl, "ADAPT_DO_MIN_DEPTH", tmplogical, &
576  "If true, make a HyCOM-like mixed layer by preventing interfaces "//&
577  "from being shallower than the depths specified by the regridding coordinate.", &
578  default=.false.)
579 
580  call set_regrid_params(cs, adapttimeratio=adapttimeratio, adaptzoom=adaptzoom, &
581  adaptzoomcoeff=adaptzoomcoeff, adaptbuoycoeff=adaptbuoycoeff, adaptalpha=adaptalpha, &
582  adaptdomin=tmplogical)
583  endif
584 
585  if (main_parameters .and. coord_is_state_dependent) then
586  call get_param(param_file, mdl, "MAXIMUM_INT_DEPTH_CONFIG", string, &
587  "Determines how to specify the maximum interface depths.\n"//&
588  "Valid options are:\n"//&
589  " NONE - there are no maximum interface depths\n"//&
590  " PARAM - use the vector-parameter MAXIMUM_INTERFACE_DEPTHS\n"//&
591  " FILE:string - read from a file. The string specifies\n"//&
592  " the filename and variable name, separated\n"//&
593  " by a comma or space, e.g. FILE:lev.nc,Z\n"//&
594  " FNC1:string - FNC1:dz_min,H_total,power,precision",&
595  default='NONE')
596  message = "The list of maximum depths for each interface."
597  allocate(z_max(ke+1))
598  allocate(dz_max(ke))
599  if ( trim(string) == "NONE") then
600  ! Do nothing.
601  elseif ( trim(string) == "PARAM") then
602  call get_param(param_file, mdl, "MAXIMUM_INTERFACE_DEPTHS", z_max, &
603  trim(message), units="m", scale=gv%m_to_H, fail_if_missing=.true.)
604  call set_regrid_max_depths(cs, z_max)
605  elseif (index(trim(string),'FILE:')==1) then
606  if (string(6:6)=='.' .or. string(6:6)=='/') then
607  ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
608  filename = trim( extractword(trim(string(6:80)), 1) )
609  else
610  ! Otherwise assume we should look for the file in INPUTDIR
611  filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
612  endif
613  if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
614  "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
615 
616  do_sum = .false.
617  varname = trim( extractword(trim(string(6:)), 2) )
618  if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
619  "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
620  if (len_trim(varname)==0) then
621  if (field_exists(filename,'z_max')) then; varname = 'z_max'
622  elseif (field_exists(filename,'dz')) then; varname = 'dz' ; do_sum = .true.
623  elseif (field_exists(filename,'dz_max')) then; varname = 'dz_max' ; do_sum = .true.
624  else ; call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
625  "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
626  endif
627  endif
628  if (do_sum) then
629  call mom_read_data(trim(filename), trim(varname), dz_max)
630  z_max(1) = 0.0 ; do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ; enddo
631  else
632  call mom_read_data(trim(filename), trim(varname), z_max)
633  endif
634  call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, &
635  trim(message), units=coordinateunits(coord_mode))
636  call set_regrid_max_depths(cs, z_max, gv%m_to_H)
637  elseif (index(trim(string),'FNC1:')==1) then
638  call dz_function1( trim(string(6:)), dz_max )
639  if ((coordinatemode(coord_mode) == regridding_slight) .and. &
640  (dz_fixed_sfc > 0.0)) then
641  do k=1,nz_fixed_sfc ; dz_max(k) = dz_fixed_sfc ; enddo
642  endif
643  z_max(1) = 0.0 ; do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ; enddo
644  call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, &
645  trim(message), units=coordinateunits(coord_mode))
646  call set_regrid_max_depths(cs, z_max, gv%m_to_H)
647  else
648  call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
649  "Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string))
650  endif
651  deallocate(z_max)
652  deallocate(dz_max)
653 
654  ! Optionally specify maximum thicknesses for each layer, enforced by moving
655  ! the interface below a layer downward.
656  call get_param(param_file, mdl, "MAX_LAYER_THICKNESS_CONFIG", string, &
657  "Determines how to specify the maximum layer thicknesses.\n"//&
658  "Valid options are:\n"//&
659  " NONE - there are no maximum layer thicknesses\n"//&
660  " PARAM - use the vector-parameter MAX_LAYER_THICKNESS\n"//&
661  " FILE:string - read from a file. The string specifies\n"//&
662  " the filename and variable name, separated\n"//&
663  " by a comma or space, e.g. FILE:lev.nc,Z\n"//&
664  " FNC1:string - FNC1:dz_min,H_total,power,precision",&
665  default='NONE')
666  message = "The list of maximum thickness for each layer."
667  allocate(h_max(ke))
668  if ( trim(string) == "NONE") then
669  ! Do nothing.
670  elseif ( trim(string) == "PARAM") then
671  call get_param(param_file, mdl, "MAX_LAYER_THICKNESS", h_max, &
672  trim(message), units="m", fail_if_missing=.true., scale=gv%m_to_H)
673  call set_regrid_max_thickness(cs, h_max)
674  elseif (index(trim(string),'FILE:')==1) then
675  if (string(6:6)=='.' .or. string(6:6)=='/') then
676  ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
677  filename = trim( extractword(trim(string(6:80)), 1) )
678  else
679  ! Otherwise assume we should look for the file in INPUTDIR
680  filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
681  endif
682  if (.not. file_exists(filename)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
683  "Specified file not found: Looking for '"//trim(filename)//"' ("//trim(string)//")")
684 
685  varname = trim( extractword(trim(string(6:)), 2) )
686  if (.not. field_exists(filename,varname)) call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
687  "Specified field not found: Looking for '"//trim(varname)//"' ("//trim(string)//")")
688  if (len_trim(varname)==0) then
689  if (field_exists(filename,'h_max')) then; varname = 'h_max'
690  elseif (field_exists(filename,'dz_max')) then; varname = 'dz_max'
691  else ; call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
692  "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
693  endif
694  endif
695  call mom_read_data(trim(filename), trim(varname), h_max)
696  call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
697  trim(message), units=coordinateunits(coord_mode))
698  call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
699  elseif (index(trim(string),'FNC1:')==1) then
700  call dz_function1( trim(string(6:)), h_max )
701  call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
702  trim(message), units=coordinateunits(coord_mode))
703  call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
704  else
705  call mom_error(fatal,trim(mdl)//", initialize_regridding: "// &
706  "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(string))
707  endif
708  deallocate(h_max)
709  endif
710 
711  if (allocated(dz)) deallocate(dz)

References check_grid_def(), dz_function1(), mom_string_functions::extract_integer(), mom_string_functions::extract_real(), mom_string_functions::extractword(), initcoord(), mom_error_handler::mom_error(), regrid_consts::regridding_adaptive, regrid_consts::regridding_hycom1, regrid_consts::regridding_slight, regriddingdefaultboundaryextrapolation, regriddingdefaultinterpscheme, regriddingdefaultminthickness, regriddinginterpschemedoc, rho_function1(), set_regrid_max_depths(), set_regrid_max_thickness(), set_regrid_params(), set_target_densities(), set_target_densities_from_gv(), setcoordinateresolution(), and uniformresolution().

Here is the call graph for this function:

◆ 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]remapcsRemapping parameters and options
[in]csRegridding control structure
[in]gOcean grid structure
[in]gvOcean vertical grid structure
[in,out]hCurrent 3D grid obtained after the last time step
[in,out]tvThermodynamical variables (T, S, ...)
[in,out]h_newNew 3D grid consistent with target coordinate
[in,out]dzinterfaceThe change in position of each interface
frac_shelf_hFractional ice shelf coverage
[in]conv_adjustIf true, do convective adjustment

Definition at line 791 of file MOM_regridding.F90.

791 !------------------------------------------------------------------------------
792 ! This routine takes care of (1) building a new grid and (2) remapping between
793 ! the old grid and the new grid. The creation of the new grid can be based
794 ! on z coordinates, target interface densities, sigma coordinates or any
795 ! arbitrary coordinate system.
796 ! The MOM6 interface positions are always calculated from the bottom up by
797 ! accumulating the layer thicknesses starting at z=-G%bathyT. z increases
798 ! upwards (decreasing k-index).
799 ! The new grid is defined by the change in position of those interfaces in z
800 ! dzInterface = zNew - zOld.
801 ! Thus, if the regridding inflates the top layer, hNew(1) > hOld(1), then the
802 ! second interface moves downward, zNew(2) < zOld(2), and dzInterface(2) < 0.
803 ! hNew(k) = hOld(k) - dzInterface(k+1) + dzInterface(k)
804 ! IMPORTANT NOTE:
805 ! This is the converse of the sign convention used in the remapping code!
806 !------------------------------------------------------------------------------
807 
808  ! Arguments
809  type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options
810  type(regridding_CS), intent(in) :: CS !< Regridding control structure
811  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
812  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
813  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after
814  !! the last time step
815  type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamical variables (T, S, ...)
816  real, dimension(SZI_(G),SZJ_(G), CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target coordinate
817  real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in position of each interface
818  real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage
819  logical, optional, intent(in ) :: conv_adjust !< If true, do convective adjustment
820  ! Local variables
821  real :: trickGnuCompiler
822  logical :: use_ice_shelf
823  logical :: do_convective_adjustment
824 
825  do_convective_adjustment = .true.
826  if (present(conv_adjust)) do_convective_adjustment = conv_adjust
827 
828  use_ice_shelf = .false.
829  if (present(frac_shelf_h)) then
830  if (associated(frac_shelf_h)) use_ice_shelf = .true.
831  endif
832 
833  select case ( cs%regridding_scheme )
834 
835  case ( regridding_zstar )
836  if (use_ice_shelf) then
837  call build_zstar_grid( cs, g, gv, h, dzinterface, frac_shelf_h )
838  else
839  call build_zstar_grid( cs, g, gv, h, dzinterface )
840  endif
841  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
842 
843  case ( regridding_sigma_shelf_zstar)
844  call build_zstar_grid( cs, g, gv, h, dzinterface )
845  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
846 
847  case ( regridding_sigma )
848  call build_sigma_grid( cs, g, gv, h, dzinterface )
849  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
850 
851  case ( regridding_rho )
852  if (do_convective_adjustment) call convective_adjustment(g, gv, h, tv)
853  call build_rho_grid( g, gv, h, tv, dzinterface, remapcs, cs )
854  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
855 
856  case ( regridding_arbitrary )
857  call build_grid_arbitrary( g, gv, h, dzinterface, trickgnucompiler, cs )
858  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
859 
860  case ( regridding_hycom1 )
861  call build_grid_hycom1( g, gv, h, tv, h_new, dzinterface, cs )
862 
863  case ( regridding_slight )
864  call build_grid_slight( g, gv, h, tv, dzinterface, cs )
865  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
866 
867  case ( regridding_adaptive )
868  call build_grid_adaptive(g, gv, h, tv, dzinterface, remapcs, cs)
869  call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
870 
871  case default
872  call mom_error(fatal,'MOM_regridding, regridding_main: '//&
873  'Unknown regridding scheme selected!')
874 
875  end select ! type of grid
876 
877 #ifdef __DO_SAFETY_CHECKS__
878  call check_remapping_grid(g, gv, h, dzinterface,'in regridding_main')
879 #endif
880 

References build_grid_adaptive(), build_grid_arbitrary(), build_grid_hycom1(), build_grid_slight(), build_rho_grid(), build_sigma_grid(), build_zstar_grid(), calc_h_new_by_dz(), check_remapping_grid(), convective_adjustment(), mom_error_handler::mom_error(), regrid_consts::regridding_adaptive, regrid_consts::regridding_hycom1, and regrid_consts::regridding_slight.

Here is the call graph for this function:

◆ 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]stringString with list of parameters in form dz_min, H_total, power, precision
[in,out]rho_targetProfile of interface densities [kg m-3]

Definition at line 2409 of file MOM_regridding.F90.

2409  character(len=*), intent(in) :: string !< String with list of parameters in form
2410  !! dz_min, H_total, power, precision
2411  real, dimension(:), allocatable, intent(inout) :: rho_target !< Profile of interface densities [kg m-3]
2412  ! Local variables
2413  integer :: nki, k, nk
2414  real :: ddx, dx, rho_1, rho_2, rho_3, drho, rho_4, drho_min
2415 
2416  read( string, *) nk, rho_1, rho_2, rho_3, drho, rho_4, drho_min
2417  allocate(rho_target(nk+1))
2418  nki = nk + 1 - 4 ! Number of interfaces minus 4 specified values
2419  rho_target(1) = rho_1
2420  rho_target(2) = rho_2
2421  dx = 0.
2422  do k = 0, nki
2423  ddx = max( drho_min, real(nki-k)/real(nki*nki) )
2424  dx = dx + ddx
2425  rho_target(3+k) = rho_3 + (2. * drho) * dx
2426  enddo
2427  rho_target(nki+4) = rho_4
2428 
2429  rho_function1 = nk
2430 

Referenced by initialize_regridding().

Here is the caller graph for this function:

◆ 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]csRegridding control structure
[in]max_depthsMaximum interface depths, in arbitrary units
[in]units_to_hA conversion factor for max_depths into H units

Definition at line 2029 of file MOM_regridding.F90.

2029  type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2030  real, dimension(CS%nk+1), intent(in) :: max_depths !< Maximum interface depths, in arbitrary units
2031  real, optional, intent(in) :: units_to_H !< A conversion factor for max_depths into H units
2032  ! Local variables
2033  real :: val_to_H
2034  integer :: K
2035 
2036  if (.not.allocated(cs%max_interface_depths)) allocate(cs%max_interface_depths(1:cs%nk+1))
2037 
2038  val_to_h = 1.0 ; if (present(units_to_h)) val_to_h = units_to_h
2039  if (max_depths(cs%nk+1) < max_depths(1)) val_to_h = -1.0*val_to_h
2040 
2041  ! Check for sign reversals in the depths.
2042  if (max_depths(cs%nk+1) < max_depths(1)) then
2043  do k=1,cs%nk ; if (max_depths(k+1) > max_depths(k)) &
2044  call mom_error(fatal, "Unordered list of maximum depths sent to set_regrid_max_depths!")
2045  enddo
2046  else
2047  do k=1,cs%nk ; if (max_depths(k+1) < max_depths(k)) &
2048  call mom_error(fatal, "Unordered list of maximum depths sent to set_regrid_max_depths.")
2049  enddo
2050  endif
2051 
2052  do k=1,cs%nk+1
2053  cs%max_interface_depths(k) = val_to_h * max_depths(k)
2054  enddo
2055 
2056  ! set max depths for coordinate
2057  select case (cs%regridding_scheme)
2058  case (regridding_hycom1)
2059  call set_hycom_params(cs%hycom_CS, max_interface_depths=cs%max_interface_depths)
2060  case (regridding_slight)
2061  call set_slight_params(cs%slight_CS, max_interface_depths=cs%max_interface_depths)
2062  end select

References mom_error_handler::mom_error(), regrid_consts::regridding_hycom1, regrid_consts::regridding_slight, coord_hycom::set_hycom_params(), and coord_slight::set_slight_params().

Referenced by initialize_regridding().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]csRegridding control structure
[in]max_hMaximum interface depths, in arbitrary units
[in]units_to_hA conversion factor for max_h into H units

Definition at line 2067 of file MOM_regridding.F90.

2067  type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2068  real, dimension(CS%nk+1), intent(in) :: max_h !< Maximum interface depths, in arbitrary units
2069  real, optional, intent(in) :: units_to_H !< A conversion factor for max_h into H units
2070  ! Local variables
2071  real :: val_to_H
2072  integer :: K
2073 
2074  if (.not.allocated(cs%max_layer_thickness)) allocate(cs%max_layer_thickness(1:cs%nk))
2075 
2076  val_to_h = 1.0 ; if (present( units_to_h)) val_to_h = units_to_h
2077 
2078  do k=1,cs%nk
2079  cs%max_layer_thickness(k) = val_to_h * max_h(k)
2080  enddo
2081 
2082  ! set max thickness for coordinate
2083  select case (cs%regridding_scheme)
2084  case (regridding_hycom1)
2085  call set_hycom_params(cs%hycom_CS, max_layer_thickness=cs%max_layer_thickness)
2086  case (regridding_slight)
2087  call set_slight_params(cs%slight_CS, max_layer_thickness=cs%max_layer_thickness)
2088  end select

References regrid_consts::regridding_hycom1, regrid_consts::regridding_slight, coord_hycom::set_hycom_params(), and coord_slight::set_slight_params().

Referenced by initialize_regridding().

Here is the call graph for this function:
Here is the caller graph for this function:

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

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

Parameters
[in,out]csRegridding control structure
[in]boundary_extrapolationExtrapolate in boundary cells
[in]min_thicknessMinimum thickness allowed when building the new grid [H ~> m or kg m-2]
[in]old_grid_weightWeight given to old coordinate when time-filtering grid
[in]interp_schemeInterpolation method for state-dependent coordinates
[in]depth_of_time_filter_shallowDepth to start cubic [H ~> m or kg m-2]
[in]depth_of_time_filter_deepDepth to end cubic [H ~> m or kg m-2]
[in]compress_fractionFraction of compressibility to add to potential density
[in]dz_min_surfaceThe fixed resolution in the topmost SLight_nkml_min layers [H ~> m or kg m-2]
[in]nz_fixed_surfaceThe number of fixed-thickness layers at the top of the model
[in]rho_ml_avg_depthAveraging depth over which to determine mixed layer potential density [H ~> m or kg m-2]
[in]nlay_ml_to_interiorNumber of layers to offset the mixed layer density to find resolved stratification [nondim]
[in]fix_haloclinesDetect regions with much weaker stratification in the coordinate
[in]halocline_filt_lenLength scale over which to filter T & S when looking for spuriously unstable water mass profiles [m]
[in]halocline_strat_tolValue of the stratification ratio that defines a problematic halocline region.
[in]integrate_downward_for_eIf true, integrate for interface positions downward from the top.
[in]adapttimeratioRatio of the ALE timestep to the grid timescale [nondim].
[in]adaptzoomDepth of near-surface zooming region [H ~> m or kg m-2].
[in]adaptzoomcoeffCoefficient of near-surface zooming diffusivity [nondim].
[in]adaptbuoycoeffCoefficient of buoyancy diffusivity [nondim].
[in]adaptalphaScaling factor on optimization tendency [nondim].
[in]adaptdominIf true, make a HyCOM-like mixed layer by preventing interfaces from being shallower than the depths specified by the regridding coordinate.

Definition at line 2218 of file MOM_regridding.F90.

2218  type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2219  logical, optional, intent(in) :: boundary_extrapolation !< Extrapolate in boundary cells
2220  real, optional, intent(in) :: min_thickness !< Minimum thickness allowed when building the
2221  !! new grid [H ~> m or kg m-2]
2222  real, optional, intent(in) :: old_grid_weight !< Weight given to old coordinate when time-filtering grid
2223  character(len=*), optional, intent(in) :: interp_scheme !< Interpolation method for state-dependent coordinates
2224  real, optional, intent(in) :: depth_of_time_filter_shallow !< Depth to start cubic [H ~> m or kg m-2]
2225  real, optional, intent(in) :: depth_of_time_filter_deep !< Depth to end cubic [H ~> m or kg m-2]
2226  real, optional, intent(in) :: compress_fraction !< Fraction of compressibility to add to potential density
2227  real, optional, intent(in) :: dz_min_surface !< The fixed resolution in the topmost
2228  !! SLight_nkml_min layers [H ~> m or kg m-2]
2229  integer, optional, intent(in) :: nz_fixed_surface !< The number of fixed-thickness layers at the top of the model
2230  real, optional, intent(in) :: Rho_ml_avg_depth !< Averaging depth over which to determine mixed layer potential
2231  !! density [H ~> m or kg m-2]
2232  real, optional, intent(in) :: nlay_ML_to_interior !< Number of layers to offset the mixed layer density to find
2233  !! resolved stratification [nondim]
2234  logical, optional, intent(in) :: fix_haloclines !< Detect regions with much weaker stratification in the coordinate
2235  real, optional, intent(in) :: halocline_filt_len !< Length scale over which to filter T & S when looking for
2236  !! spuriously unstable water mass profiles [m]
2237  real, optional, intent(in) :: halocline_strat_tol !< Value of the stratification ratio that defines a problematic
2238  !! halocline region.
2239  logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface positions downward
2240  !! from the top.
2241  real, optional, intent(in) :: adaptTimeRatio !< Ratio of the ALE timestep to the grid timescale [nondim].
2242  real, optional, intent(in) :: adaptZoom !< Depth of near-surface zooming region [H ~> m or kg m-2].
2243  real, optional, intent(in) :: adaptZoomCoeff !< Coefficient of near-surface zooming diffusivity [nondim].
2244  real, optional, intent(in) :: adaptBuoyCoeff !< Coefficient of buoyancy diffusivity [nondim].
2245  real, optional, intent(in) :: adaptAlpha !< Scaling factor on optimization tendency [nondim].
2246  logical, optional, intent(in) :: adaptDoMin !< If true, make a HyCOM-like mixed layer by
2247  !! preventing interfaces from being shallower than
2248  !! the depths specified by the regridding coordinate.
2249 
2250  if (present(interp_scheme)) call set_interp_scheme(cs%interp_CS, interp_scheme)
2251  if (present(boundary_extrapolation)) call set_interp_extrap(cs%interp_CS, boundary_extrapolation)
2252 
2253  if (present(old_grid_weight)) then
2254  if (old_grid_weight<0. .or. old_grid_weight>1.) &
2255  call mom_error(fatal,'MOM_regridding, set_regrid_params: Weight is out side the range 0..1!')
2256  cs%old_grid_weight = old_grid_weight
2257  endif
2258  if (present(depth_of_time_filter_shallow)) cs%depth_of_time_filter_shallow = depth_of_time_filter_shallow
2259  if (present(depth_of_time_filter_deep)) cs%depth_of_time_filter_deep = depth_of_time_filter_deep
2260  if (present(depth_of_time_filter_shallow) .or. present(depth_of_time_filter_deep)) then
2261  if (cs%depth_of_time_filter_deep<cs%depth_of_time_filter_shallow) call mom_error(fatal,'MOM_regridding, '//&
2262  'set_regrid_params: depth_of_time_filter_deep<depth_of_time_filter_shallow!')
2263  endif
2264 
2265  if (present(min_thickness)) cs%min_thickness = min_thickness
2266  if (present(compress_fraction)) cs%compressibility_fraction = compress_fraction
2267  if (present(integrate_downward_for_e)) cs%integrate_downward_for_e = integrate_downward_for_e
2268 
2269  select case (cs%regridding_scheme)
2270  case (regridding_zstar)
2271  if (present(min_thickness)) call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
2272  case (regridding_sigma_shelf_zstar)
2273  if (present(min_thickness)) call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
2274  case (regridding_sigma)
2275  if (present(min_thickness)) call set_sigma_params(cs%sigma_CS, min_thickness=min_thickness)
2276  case (regridding_rho)
2277  if (present(min_thickness)) call set_rho_params(cs%rho_CS, min_thickness=min_thickness)
2278  if (present(integrate_downward_for_e)) &
2279  call set_rho_params(cs%rho_CS, integrate_downward_for_e=integrate_downward_for_e)
2280  if (associated(cs%rho_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
2281  call set_rho_params(cs%rho_CS, interp_cs=cs%interp_CS)
2282  case (regridding_hycom1)
2283  if (associated(cs%hycom_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
2284  call set_hycom_params(cs%hycom_CS, interp_cs=cs%interp_CS)
2285  case (regridding_slight)
2286  if (present(min_thickness)) call set_slight_params(cs%slight_CS, min_thickness=min_thickness)
2287  if (present(dz_min_surface)) call set_slight_params(cs%slight_CS, dz_ml_min=dz_min_surface)
2288  if (present(nz_fixed_surface)) call set_slight_params(cs%slight_CS, nz_fixed_surface=nz_fixed_surface)
2289  if (present(rho_ml_avg_depth)) call set_slight_params(cs%slight_CS, rho_ml_avg_depth=rho_ml_avg_depth)
2290  if (present(nlay_ml_to_interior)) call set_slight_params(cs%slight_CS, nlay_ml_offset=nlay_ml_to_interior)
2291  if (present(fix_haloclines)) call set_slight_params(cs%slight_CS, fix_haloclines=fix_haloclines)
2292  if (present(halocline_filt_len)) call set_slight_params(cs%slight_CS, halocline_filter_length=halocline_filt_len)
2293  if (present(halocline_strat_tol)) call set_slight_params(cs%slight_CS, halocline_strat_tol=halocline_strat_tol)
2294  if (present(compress_fraction)) call set_slight_params(cs%slight_CS, compressibility_fraction=compress_fraction)
2295  if (associated(cs%slight_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
2296  call set_slight_params(cs%slight_CS, interp_cs=cs%interp_CS)
2297  case (regridding_adaptive)
2298  if (present(adapttimeratio)) call set_adapt_params(cs%adapt_CS, adapttimeratio=adapttimeratio)
2299  if (present(adaptzoom)) call set_adapt_params(cs%adapt_CS, adaptzoom=adaptzoom)
2300  if (present(adaptzoomcoeff)) call set_adapt_params(cs%adapt_CS, adaptzoomcoeff=adaptzoomcoeff)
2301  if (present(adaptbuoycoeff)) call set_adapt_params(cs%adapt_CS, adaptbuoycoeff=adaptbuoycoeff)
2302  if (present(adaptalpha)) call set_adapt_params(cs%adapt_CS, adaptalpha=adaptalpha)
2303  if (present(adaptdomin)) call set_adapt_params(cs%adapt_CS, adaptdomin=adaptdomin)
2304  end select
2305 

References mom_error_handler::mom_error(), regrid_consts::regridding_adaptive, regrid_consts::regridding_hycom1, regrid_consts::regridding_slight, coord_adapt::set_adapt_params(), coord_hycom::set_hycom_params(), regrid_interp::set_interp_extrap(), regrid_interp::set_interp_scheme(), coord_sigma::set_sigma_params(), coord_slight::set_slight_params(), and coord_zlike::set_zlike_params().

Referenced by mom_oda_driver_mod::init_oda(), and initialize_regridding().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]csRegridding control structure
[in]rho_intInterface densities [R ~> kg m-3]

Definition at line 2015 of file MOM_regridding.F90.

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

References mom_error_handler::mom_error().

Referenced by initialize_regridding().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]gvOcean vertical grid structure
[in]usA dimensional unit scaling type
[in,out]csRegridding control structure

Definition at line 1997 of file MOM_regridding.F90.

1997  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1998  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1999  type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2000  ! Local variables
2001  integer :: k, nz
2002 
2003  nz = cs%nk
2004  cs%target_density(1) = (gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(2)))
2005  cs%target_density(nz+1) = (gv%Rlay(nz) + 0.5*(gv%Rlay(nz)-gv%Rlay(nz-1)))
2006  do k = 2,nz
2007  cs%target_density(k) = cs%target_density(k-1) + cs%coordinateResolution(k)
2008  enddo
2009  cs%target_density_set = .true.
2010 

Referenced by initialize_regridding().

Here is the caller graph for this function:

◆ 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]dzA vector of vertical grid spacings
[in,out]csRegridding control structure
[in]scaleA scaling factor converting dz to coordRes

Definition at line 1980 of file MOM_regridding.F90.

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

References mom_error_handler::mom_error().

Referenced by initialize_regridding().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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]nkNumber of cells in source grid
[in]coordmodeA string indicating the coordinate mode. See the documenttion for regrid_consts for the recognized values.
[in]maxdepthThe range of the grid values in some modes
[in]rholightThe minimum value of the grid in RHO mode
[in]rhoheavyThe maximum value of the grid in RHO mode
Returns
The returned uniform resolution grid.

Definition at line 1908 of file MOM_regridding.F90.

1908 !------------------------------------------------------------------------------
1909 ! Calculate a vector of uniform resolution in the units of the coordinate
1910 !------------------------------------------------------------------------------
1911  ! Arguments
1912  integer, intent(in) :: nk !< Number of cells in source grid
1913  character(len=*), intent(in) :: coordMode !< A string indicating the coordinate mode.
1914  !! See the documenttion for regrid_consts
1915  !! for the recognized values.
1916  real, intent(in) :: maxDepth !< The range of the grid values in some modes
1917  real, intent(in) :: rhoLight !< The minimum value of the grid in RHO mode
1918  real, intent(in) :: rhoHeavy !< The maximum value of the grid in RHO mode
1919 
1920  real :: uniformResolution(nk) !< The returned uniform resolution grid.
1921 
1922  ! Local variables
1923  integer :: scheme
1924 
1925  scheme = coordinatemode(coordmode)
1926  select case ( scheme )
1927 
1928  case ( regridding_zstar, regridding_hycom1, regridding_slight, regridding_sigma_shelf_zstar, &
1929  regridding_adaptive )
1930  uniformresolution(:) = maxdepth / real(nk)
1931 
1932  case ( regridding_rho )
1933  uniformresolution(:) = (rhoheavy - rholight) / real(nk)
1934 
1935  case ( regridding_sigma )
1936  uniformresolution(:) = 1. / real(nk)
1937 
1938  case default
1939  call mom_error(fatal, "MOM_regridding, uniformResolution: "//&
1940  "Unrecognized choice for coordinate mode ("//trim(coordmode)//").")
1941 
1942  end select ! type of grid
1943 

References mom_error_handler::mom_error(), regrid_consts::regridding_adaptive, regrid_consts::regridding_hycom1, and regrid_consts::regridding_slight.

Referenced by initialize_regridding().

Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ regriddingcoordinatemodedoc

character(len=*), parameter, public mom_regridding::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.

Definition at line 141 of file MOM_regridding.F90.

141 character(len=*), parameter, public :: regriddingCoordinateModeDoc = &
142  " LAYER - Isopycnal or stacked shallow water layers\n"//&
143  " ZSTAR, Z* - stretched geopotential z*\n"//&
144  " SIGMA_SHELF_ZSTAR - stretched geopotential z* ignoring shelf\n"//&
145  " SIGMA - terrain following coordinates\n"//&
146  " RHO - continuous isopycnal\n"//&
147  " HYCOM1 - HyCOM-like hybrid coordinate\n"//&
148  " SLIGHT - stretched coordinates above continuous isopycnal\n"//&
149  " ADAPTIVE - optimize for smooth neutral density surfaces"

◆ regriddingdefaultboundaryextrapolation

logical, parameter, public mom_regridding::regriddingdefaultboundaryextrapolation = .false.

Default mode for boundary extrapolation.

Definition at line 167 of file MOM_regridding.F90.

167 logical, parameter, public :: regriddingDefaultBoundaryExtrapolation = .false.

Referenced by initialize_regridding().

◆ regriddingdefaultinterpscheme

character(len=*), parameter, public mom_regridding::regriddingdefaultinterpscheme = "P1M_H2"

Default interpolation scheme.

Definition at line 165 of file MOM_regridding.F90.

165 character(len=*), parameter, public :: regriddingDefaultInterpScheme = "P1M_H2"

Referenced by initialize_regridding().

◆ regriddingdefaultminthickness

real, parameter, public mom_regridding::regriddingdefaultminthickness = 1.e-3

Default minimum thickness for some coordinate generation modes.

Definition at line 169 of file MOM_regridding.F90.

169 real, parameter, public :: regriddingDefaultMinThickness = 1.e-3

Referenced by initialize_regridding().

◆ regriddinginterpschemedoc

character(len=*), parameter, public mom_regridding::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.

Definition at line 152 of file MOM_regridding.F90.

152 character(len=*), parameter, public :: regriddingInterpSchemeDoc = &
153  " P1M_H2 (2nd-order accurate)\n"//&
154  " P1M_H4 (2nd-order accurate)\n"//&
155  " P1M_IH4 (2nd-order accurate)\n"//&
156  " PLM (2nd-order accurate)\n"//&
157  " PPM_H4 (3rd-order accurate)\n"//&
158  " PPM_IH4 (3rd-order accurate)\n"//&
159  " P3M_IH4IH3 (4th-order accurate)\n"//&
160  " P3M_IH6IH5 (4th-order accurate)\n"//&
161  " PQM_IH4IH3 (4th-order accurate)\n"//&
162  " PQM_IH6IH5 (5th-order accurate)"

Referenced by initialize_regridding().