From 37f84b771bd7e4b5c4e9c4be596a12e5e50f002c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 13 Mar 2023 09:31:41 -0400 Subject: [PATCH] +(*)Remove coord_SLight Removed the coord_slight module and all calls to it, and obsoleted all run-time parameters that are exclusively related to it. This code was an attempt from 2015 to define an appropriate hybrid vertical coordinate for global climate modeling, but it never worked very well (usually falling apart in the second year), and it has not been used in any publication or active model for many years. The test case that exercised this coordinate in the MOM6-examples test suite is also being removed via MOM6-examples PR #388. The coord_SLight code is being eliminated altogether now to simplify the MOM6 code base and reduce the volume of untested and unused code. All answers in all known MOM6 configurations in active use are bitwise identical, although there is a remote chance that someone somewhere might be using the SLIGHT coordinate. --- src/ALE/MOM_regridding.F90 | 188 +----- src/ALE/coord_slight.F90 | 733 ------------------------ src/ALE/regrid_consts.F90 | 6 - src/diagnostics/MOM_obsolete_params.F90 | 8 + src/framework/MOM_diag_remap.F90 | 4 - 5 files changed, 18 insertions(+), 921 deletions(-) delete mode 100644 src/ALE/coord_slight.F90 diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 53072909a5..2fc65d0298 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -22,7 +22,7 @@ module MOM_regridding use regrid_consts, only : REGRIDDING_LAYER, REGRIDDING_ZSTAR use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SIGMA use regrid_consts, only : REGRIDDING_ARBITRARY, REGRIDDING_SIGMA_SHELF_ZSTAR -use regrid_consts, only : REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, REGRIDDING_SLIGHT, REGRIDDING_ADAPTIVE +use regrid_consts, only : REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, REGRIDDING_ADAPTIVE use regrid_interp, only : interp_CS_type, set_interp_scheme, set_interp_extrap use coord_zlike, only : init_coord_zlike, zlike_CS, set_zlike_params, build_zstar_column, end_coord_zlike @@ -30,7 +30,6 @@ module MOM_regridding use coord_rho, only : init_coord_rho, rho_CS, set_rho_params, build_rho_column, end_coord_rho use coord_rho, only : old_inflate_layers_1d use coord_hycom, only : init_coord_hycom, hycom_CS, set_hycom_params, build_hycom1_column, end_coord_hycom -use coord_slight, only : init_coord_slight, slight_CS, set_slight_params, build_slight_column, end_coord_slight use coord_adapt, only : init_coord_adapt, adapt_CS, set_adapt_params, build_adapt_column, end_coord_adapt use MOM_hybgen_regrid, only : hybgen_regrid, hybgen_regrid_CS, init_hybgen_regrid, end_hybgen_regrid use MOM_hybgen_regrid, only : write_Hybgen_coord_file @@ -61,7 +60,7 @@ module MOM_regridding !! This array is the nominal coordinate of interfaces and is the !! running sum of coordinateResolution, in [R ~> kg m-3]. i.e. !! target_density(k+1) = coordinateResolution(k) + coordinateResolution(k) - !! It is only used in "rho", "SLight" or "Hycom" mode. + !! It is only used in "rho" or "Hycom" mode. real, dimension(:), allocatable :: target_density !> A flag to indicate that the target_density arrays has been filled with data. @@ -129,7 +128,6 @@ module MOM_regridding type(sigma_CS), pointer :: sigma_CS => null() !< Control structure for sigma coordinate generator type(rho_CS), pointer :: rho_CS => null() !< Control structure for rho coordinate generator type(hycom_CS), pointer :: hycom_CS => null() !< Control structure for hybrid coordinate generator - type(slight_CS), pointer :: slight_CS => null() !< Control structure for Slight-coordinate generator type(adapt_CS), pointer :: adapt_CS => null() !< Control structure for adaptive coordinate generator type(hybgen_regrid_CS), pointer :: hybgen_CS => NULL() !< Control structure for hybgen regridding @@ -157,7 +155,6 @@ module MOM_regridding " RHO - continuous isopycnal\n"//& " HYCOM1 - HyCOM-like hybrid coordinate\n"//& " HYBGEN - Hybrid coordinate from the Hycom hybgen code\n"//& - " SLIGHT - stretched coordinates above continuous isopycnal\n"//& " ADAPTIVE - optimize for smooth neutral density surfaces" !> Documentation for regridding interpolation schemes @@ -208,18 +205,17 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m character(len=200) :: inputdir, fileName character(len=320) :: message ! Temporary strings character(len=12) :: expected_units, alt_units ! Temporary strings - logical :: tmpLogical, fix_haloclines, do_sum, main_parameters + logical :: tmpLogical, do_sum, main_parameters logical :: coord_is_state_dependent, ierr integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. logical :: remap_answers_2018 integer :: remap_answer_date ! The vintage of the remapping expressions to use. - real :: filt_len, strat_tol, tmpReal, P_Ref + real :: tmpReal, P_Ref real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z). - real :: dz_fixed_sfc, Rho_avg_depth, nlay_sfc_int real :: adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha real :: adaptDrho0 ! Reference density difference for stratification-dependent diffusion. [R ~> kg m-3] - integer :: nz_fixed_sfc, k, nzf(4) + integer :: k, nzf(4) real, dimension(:), allocatable :: dz ! Resolution (thickness) in units of coordinate, which may be [m] ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units. real, dimension(:), allocatable :: h_max ! Maximum layer thicknesses [H ~> m or kg m-2] @@ -488,7 +484,6 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m if (coordinateMode(coord_mode) == REGRIDDING_ZSTAR .or. & coordinateMode(coord_mode) == REGRIDDING_HYCOM1 .or. & coordinateMode(coord_mode) == REGRIDDING_HYBGEN .or. & - coordinateMode(coord_mode) == REGRIDDING_SLIGHT .or. & coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then ! Adjust target grid to be consistent with maximum_depth tmpReal = sum( dz(:) ) @@ -589,49 +584,6 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m default=.false.) endif - if (coordinateMode(coord_mode) == REGRIDDING_SLIGHT) then - ! Set SLight-specific regridding parameters. - call get_param(param_file, mdl, "SLIGHT_DZ_SURFACE", dz_fixed_sfc, & - "The nominal thickness of fixed thickness near-surface "//& - "layers with the SLight coordinate.", units="m", default=1.0, scale=GV%m_to_H) - call get_param(param_file, mdl, "SLIGHT_NZ_SURFACE_FIXED", nz_fixed_sfc, & - "The number of fixed-depth surface layers with the SLight "//& - "coordinate.", units="nondimensional", default=2) - call get_param(param_file, mdl, "SLIGHT_SURFACE_AVG_DEPTH", Rho_avg_depth, & - "The thickness of the surface region over which to average "//& - "when calculating the density to use to define the interior "//& - "with the SLight coordinate.", units="m", default=1.0, scale=GV%m_to_H) - call get_param(param_file, mdl, "SLIGHT_NLAY_TO_INTERIOR", nlay_sfc_int, & - "The number of layers to offset the surface density when "//& - "defining where the interior ocean starts with SLight.", & - units="nondimensional", default=2.0) - call get_param(param_file, mdl, "SLIGHT_FIX_HALOCLINES", fix_haloclines, & - "If true, identify regions above the reference pressure "//& - "where the reference pressure systematically underestimates "//& - "the stratification and use this in the definition of the "//& - "interior with the SLight coordinate.", default=.false.) - - call set_regrid_params(CS, dz_min_surface=dz_fixed_sfc, & - nz_fixed_surface=nz_fixed_sfc, Rho_ML_avg_depth=Rho_avg_depth, & - nlay_ML_to_interior=nlay_sfc_int, fix_haloclines=fix_haloclines) - if (fix_haloclines) then - ! Set additional parameters related to SLIGHT_FIX_HALOCLINES. - call get_param(param_file, mdl, "HALOCLINE_FILTER_LENGTH", filt_len, & - "A length scale over which to smooth the temperature and "//& - "salinity before identifying erroneously unstable haloclines.", & - units="m", default=2.0, scale=GV%m_to_H) - call get_param(param_file, mdl, "HALOCLINE_STRAT_TOL", strat_tol, & - "A tolerance for the ratio of the stratification of the "//& - "apparent coordinate stratification to the actual value "//& - "that is used to identify erroneously unstable haloclines. "//& - "This ratio is 1 when they are equal, and sensible values "//& - "are between 0 and 0.5.", units="nondimensional", default=0.2) - call set_regrid_params(CS, halocline_filt_len=filt_len, & - halocline_strat_tol=strat_tol) - endif - - endif - if (coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then call get_param(param_file, mdl, "ADAPT_TIME_RATIO", adaptTimeRatio, & "Ratio of ALE timestep to grid timescale.", units="nondim", default=1.0e-1) @@ -710,10 +662,6 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m call set_regrid_max_depths(CS, z_max, GV%m_to_H) elseif (index(trim(string),'FNC1:')==1) then call dz_function1( trim(string(6:)), dz_max ) - if ((coordinateMode(coord_mode) == REGRIDDING_SLIGHT) .and. & - (dz_fixed_sfc > 0.0)) then - do k=1,nz_fixed_sfc ; dz_max(k) = dz_fixed_sfc ; enddo - endif z_max(1) = 0.0 ; do K=1,ke ; z_max(K+1) = z_max(K) + dz_max(K) ; enddo call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, & trim(message), units=coordinateUnits(coord_mode)) @@ -795,7 +743,6 @@ subroutine end_regridding(CS) if (associated(CS%sigma_CS)) call end_coord_sigma(CS%sigma_CS) if (associated(CS%rho_CS)) call end_coord_rho(CS%rho_CS) if (associated(CS%hycom_CS)) call end_coord_hycom(CS%hycom_CS) - if (associated(CS%slight_CS)) call end_coord_slight(CS%slight_CS) if (associated(CS%adapt_CS)) call end_coord_adapt(CS%adapt_CS) if (associated(CS%hybgen_CS)) call end_hybgen_regrid(CS%hybgen_CS) @@ -869,9 +816,6 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, & case ( REGRIDDING_HYBGEN ) call hybgen_regrid(G, GV, G%US, h, tv, CS%hybgen_CS, dzInterface, PCM_cell) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) - case ( REGRIDDING_SLIGHT ) - call build_grid_SLight( G, GV, G%US, h, tv, dzInterface, CS ) - call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_ADAPTIVE ) call build_grid_adaptive(G, GV, G%US, h, tv, dzInterface, remapCS, CS) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) @@ -915,7 +859,7 @@ subroutine regridding_preadjust_reqs(CS, do_conv_adj, do_hybgen_unmix, hybgen_CS select case ( CS%regridding_scheme ) case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_SIGMA, REGRIDDING_ARBITRARY, & - REGRIDDING_HYCOM1, REGRIDDING_SLIGHT, REGRIDDING_ADAPTIVE ) + REGRIDDING_HYCOM1, REGRIDDING_ADAPTIVE ) do_conv_adj = .false. ; do_hybgen_unmix = .false. case ( REGRIDDING_RHO ) do_conv_adj = .true. ; do_hybgen_unmix = .false. @@ -1662,84 +1606,6 @@ subroutine build_grid_adaptive(G, GV, US, h, tv, dzInterface, remapCS, CS) enddo ; enddo end subroutine build_grid_adaptive -!> 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. -subroutine build_grid_SLight(G, GV, US, h, tv, dzInterface, CS) - type(ocean_grid_type), intent(in) :: G !< Grid structure - type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Existing model thickness [H ~> m or kg m-2] - type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure - type(regridding_CS), intent(in) :: CS !< Regridding control structure - real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< Changes in interface position - - real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2] - real, dimension(SZK_(GV)) :: p_col ! Layer center pressure in the input column [R L2 T-2 ~> Pa] - real, dimension(CS%nk+1) :: z_col_new ! New interface positions relative to the surface [H ~> m or kg m-2] - real, dimension(CS%nk+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2] - - ! Local variables - real :: depth ! Depth of the ocean relative to the mean sea surface height in thickness units [H ~> m or kg m-2] - integer :: i, j, k, nz - real :: h_neglect, h_neglect_edge - - if (CS%remap_answer_date >= 20190101) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then - h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 - else - h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 - endif - - nz = GV%ke - - call assert((GV%ke == CS%nk), "build_grid_SLight is only written to work "//& - "with the same number of input and target layers.") - call assert(CS%target_density_set, "build_grid_SLight : "//& - "Target densities must be set before build_grid_SLight is called.") - - ! Build grid based on target interface densities - do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 - if (G%mask2dT(i,j)>0.) then - - depth = (G%bathyT(i,j)+G%Z_ref) * GV%Z_to_H - z_col(1) = 0. ! Work downward rather than bottom up - do K=1,nz - z_col(K+1) = z_col(K) + h(i,j,k) - p_col(k) = tv%P_Ref + CS%compressibility_fraction * & - ( 0.5 * ( z_col(K) + z_col(K+1) ) * (GV%H_to_RZ*GV%g_Earth) - tv%P_Ref ) - enddo - - call build_slight_column(CS%slight_CS, tv%eqn_of_state, GV%H_to_RZ*GV%g_Earth, & - GV%H_subroundoff, nz, depth, h(i, j, :), & - tv%T(i, j, :), tv%S(i, j, :), p_col, z_col, z_col_new, & - h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) - - ! Calculate the final change in grid position after blending new and old grids - call filtered_grid_motion( CS, nz, z_col, z_col_new, dz_col ) - do K=1,nz+1 ; dzInterface(i,j,K) = -dz_col(K) ; enddo -#ifdef __DO_SAFETY_CHECKS__ - if (dzInterface(i,j,1) /= 0.) stop 'build_grid_SLight: Surface moved?!' - if (dzInterface(i,j,nz+1) /= 0.) stop 'build_grid_SLight: Bottom moved?!' -#endif - - ! This adjusts things robust to round-off errors - call adjust_interface_motion( CS, nz, h(i,j,:), dzInterface(i,j,:) ) - - else ! on land - dzInterface(i,j,:) = 0. - endif ! mask2dT - enddo ; enddo ! i,j - -end subroutine build_grid_SLight - !> Adjust dz_Interface to ensure non-negative future thicknesses subroutine adjust_interface_motion( CS, nk, h_old, dz_int ) type(regridding_CS), intent(in) :: CS !< Regridding control structure @@ -2033,7 +1899,7 @@ function uniformResolution(nk,coordMode,maxDepth,rhoLight,rhoHeavy) scheme = coordinateMode(coordMode) select case ( scheme ) - case ( REGRIDDING_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, REGRIDDING_SLIGHT, & + case ( REGRIDDING_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, & REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_ADAPTIVE ) uniformResolution(:) = maxDepth / real(nk) @@ -2076,9 +1942,6 @@ subroutine initCoord(CS, GV, US, coord_mode, param_file) CS%interp_CS) case (REGRIDDING_HYBGEN) call init_hybgen_regrid(CS%hybgen_CS, GV, US, param_file) - case (REGRIDDING_SLIGHT) - call init_coord_slight(CS%slight_CS, CS%nk, CS%ref_pressure, CS%target_density, & - CS%interp_CS, GV%m_to_H) case (REGRIDDING_ADAPTIVE) call init_coord_adapt(CS%adapt_CS, CS%nk, CS%coordinateResolution, GV%m_to_H, US%kg_m3_to_R) end select @@ -2172,8 +2035,6 @@ subroutine set_regrid_max_depths( CS, max_depths, units_to_H ) select case (CS%regridding_scheme) case (REGRIDDING_HYCOM1) call set_hycom_params(CS%hycom_CS, max_interface_depths=CS%max_interface_depths) - case (REGRIDDING_SLIGHT) - call set_slight_params(CS%slight_CS, max_interface_depths=CS%max_interface_depths) end select end subroutine set_regrid_max_depths @@ -2198,8 +2059,6 @@ subroutine set_regrid_max_thickness( CS, max_h, units_to_H ) select case (CS%regridding_scheme) case (REGRIDDING_HYCOM1) call set_hycom_params(CS%hycom_CS, max_layer_thickness=CS%max_layer_thickness) - case (REGRIDDING_SLIGHT) - call set_slight_params(CS%slight_CS, max_layer_thickness=CS%max_layer_thickness) end select end subroutine set_regrid_max_thickness @@ -2311,7 +2170,7 @@ function getCoordinateUnits( CS ) character(len=20) :: getCoordinateUnits select case ( CS%regridding_scheme ) - case ( REGRIDDING_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, REGRIDDING_SLIGHT, & + case ( REGRIDDING_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, & REGRIDDING_ADAPTIVE ) getCoordinateUnits = 'meter' case ( REGRIDDING_SIGMA_SHELF_ZSTAR ) @@ -2352,8 +2211,6 @@ function getCoordinateShortName( CS ) getCoordinateShortName = 'z-rho' case ( REGRIDDING_HYBGEN ) getCoordinateShortName = 'hybrid' - case ( REGRIDDING_SLIGHT ) - getCoordinateShortName = 's-rho' case ( REGRIDDING_ADAPTIVE ) getCoordinateShortName = 'adaptive' case default @@ -2366,8 +2223,7 @@ end function getCoordinateShortName !> Can be used to set any of the parameters for MOM_regridding. subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_grid_weight, & interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, & - compress_fraction, ref_pressure, dz_min_surface, nz_fixed_surface, Rho_ML_avg_depth, & - nlay_ML_to_interior, fix_haloclines, halocline_filt_len, halocline_strat_tol, & + compress_fraction, ref_pressure, & integrate_downward_for_e, remap_answers_2018, remap_answer_date, & adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin, adaptDrho0) type(regridding_CS), intent(inout) :: CS !< Regridding control structure @@ -2381,18 +2237,6 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri real, optional, intent(in) :: compress_fraction !< Fraction of compressibility to add to potential density [nondim] real, optional, intent(in) :: ref_pressure !< The reference pressure for density-dependent !! coordinates [R L2 T-2 ~> Pa] - real, optional, intent(in) :: dz_min_surface !< The fixed resolution in the topmost - !! SLight_nkml_min layers [H ~> m or kg m-2] - integer, optional, intent(in) :: nz_fixed_surface !< The number of fixed-thickness layers at the top of the model - real, optional, intent(in) :: Rho_ml_avg_depth !< Averaging depth over which to determine mixed layer potential - !! density [H ~> m or kg m-2] - real, optional, intent(in) :: nlay_ML_to_interior !< Number of layers to offset the mixed layer density to find - !! resolved stratification [nondim] - logical, optional, intent(in) :: fix_haloclines !< Detect regions with much weaker stratification in the coordinate - real, optional, intent(in) :: halocline_filt_len !< Length scale over which to filter T & S when looking for - !! spuriously unstable water mass profiles [H ~> m or kg m-2] - real, optional, intent(in) :: halocline_strat_tol !< Value of the stratification ratio that defines a problematic - !! halocline region. logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface positions downward !! from the top. logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions @@ -2457,18 +2301,6 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri call set_hycom_params(CS%hycom_CS, interp_CS=CS%interp_CS) case (REGRIDDING_HYBGEN) ! Do nothing for now. - case (REGRIDDING_SLIGHT) - if (present(min_thickness)) call set_slight_params(CS%slight_CS, min_thickness=min_thickness) - if (present(dz_min_surface)) call set_slight_params(CS%slight_CS, dz_ml_min=dz_min_surface) - if (present(nz_fixed_surface)) call set_slight_params(CS%slight_CS, nz_fixed_surface=nz_fixed_surface) - if (present(Rho_ML_avg_depth)) call set_slight_params(CS%slight_CS, Rho_ML_avg_depth=Rho_ML_avg_depth) - if (present(nlay_ML_to_interior)) call set_slight_params(CS%slight_CS, nlay_ML_offset=nlay_ML_to_interior) - if (present(fix_haloclines)) call set_slight_params(CS%slight_CS, fix_haloclines=fix_haloclines) - if (present(halocline_filt_len)) call set_slight_params(CS%slight_CS, halocline_filter_length=halocline_filt_len) - if (present(halocline_strat_tol)) call set_slight_params(CS%slight_CS, halocline_strat_tol=halocline_strat_tol) - if (present(compress_fraction)) call set_slight_params(CS%slight_CS, compressibility_fraction=compress_fraction) - if (associated(CS%slight_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) & - call set_slight_params(CS%slight_CS, interp_CS=CS%interp_CS) case (REGRIDDING_ADAPTIVE) if (present(adaptTimeRatio)) call set_adapt_params(CS%adapt_CS, adaptTimeRatio=adaptTimeRatio) if (present(adaptZoom)) call set_adapt_params(CS%adapt_CS, adaptZoom=adaptZoom) @@ -2526,7 +2358,7 @@ function getStaticThickness( CS, SSH, depth ) select case ( CS%regridding_scheme ) case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, & - REGRIDDING_SLIGHT, REGRIDDING_ADAPTIVE ) + REGRIDDING_ADAPTIVE ) if (depth>0.) then z = ssh do k = 1, CS%nk diff --git a/src/ALE/coord_slight.F90 b/src/ALE/coord_slight.F90 deleted file mode 100644 index 4b4ac8a153..0000000000 --- a/src/ALE/coord_slight.F90 +++ /dev/null @@ -1,733 +0,0 @@ -!> Regrid columns for the SLight coordinate -module coord_slight - -! This file is part of MOM6. See LICENSE.md for the license. - -use MOM_error_handler, only : MOM_error, FATAL -use MOM_EOS, only : EOS_type, calculate_compress -use MOM_EOS, only : calculate_density, calculate_density_derivs -use regrid_interp, only : interp_CS_type, regridding_set_ppolys -use regrid_interp, only : NR_ITERATIONS, NR_TOLERANCE, DEGREE_MAX - -implicit none ; private - -!> Control structure containing required parameters for the SLight coordinate -type, public :: slight_CS ; private - - !> Number of layers/levels - integer :: nk - - !> Minimum thickness allowed when building the new grid through regridding [H ~> m or kg m-2] - real :: min_thickness - - !> Reference pressure for potential density calculations [R L2 T-2 ~> Pa] - real :: ref_pressure - - !> Fraction (between 0 and 1) of compressibility to add to potential density - !! profiles when interpolating for target grid positions. [nondim] - real :: compressibility_fraction - - ! The following 4 parameters were introduced for use with the SLight coordinate: - !> Depth over which to average to determine the mixed layer potential density [H ~> m or kg m-2] - real :: Rho_ML_avg_depth - - !> Number of layers to offset the mixed layer density to find resolved stratification [nondim] - real :: nlay_ml_offset - - !> The number of fixed-thickness layers at the top of the model - integer :: nz_fixed_surface = 2 - - !> The fixed resolution in the topmost SLight_nkml_min layers [H ~> m or kg m-2] - real :: dz_ml_min - - !> If true, detect regions with much weaker stratification in the coordinate - !! than based on in-situ density, and use a stretched coordinate there. - logical :: fix_haloclines = .false. - - !> A length scale over which to filter T & S when looking for spuriously - !! unstable water mass profiles [H ~> m or kg m-2]. - real :: halocline_filter_length - - !> A value of the stratification ratio that defines a problematic halocline region [nondim]. - real :: halocline_strat_tol - - !> Nominal density of interfaces [R ~> kg m-3]. - real, allocatable, dimension(:) :: target_density - - !> Maximum depths of interfaces [H ~> m or kg m-2]. - real, allocatable, dimension(:) :: max_interface_depths - - !> Maximum thicknesses of layers [H ~> m or kg m-2]. - real, allocatable, dimension(:) :: max_layer_thickness - - !> Interpolation control structure - type(interp_CS_type) :: interp_CS -end type slight_CS - -public init_coord_slight, set_slight_params, build_slight_column, end_coord_slight - -contains - -!> Initialise a slight_CS with pointers to parameters -subroutine init_coord_slight(CS, nk, ref_pressure, target_density, interp_CS, m_to_H) - type(slight_CS), pointer :: CS !< Unassociated pointer to hold the control structure - integer, intent(in) :: nk !< Number of layers in the grid - real, intent(in) :: ref_pressure !< Coordinate reference pressure [R L2 T-2 ~> Pa] - real, dimension(:), intent(in) :: target_density !< Nominal density of interfaces [R ~> kg m-3] - type(interp_CS_type), intent(in) :: interp_CS !< Controls for interpolation - real, optional, intent(in) :: m_to_H !< A conversion factor from m to the units of thicknesses - - real :: m_to_H_rescale ! A unit conversion factor. - - if (associated(CS)) call MOM_error(FATAL, "init_coord_slight: CS already associated!") - allocate(CS) - allocate(CS%target_density(nk+1)) - - m_to_H_rescale = 1.0 ; if (present(m_to_H)) m_to_H_rescale = m_to_H - - CS%nk = nk - CS%ref_pressure = ref_pressure - CS%target_density(:) = target_density(:) - CS%interp_CS = interp_CS - - ! Set real parameter default values - CS%compressibility_fraction = 0. ! Nondim. - CS%Rho_ML_avg_depth = 1.0 * m_to_H_rescale - CS%nlay_ml_offset = 2.0 ! Nondim. - CS%dz_ml_min = 1.0 * m_to_H_rescale - CS%halocline_filter_length = 2.0 * m_to_H_rescale - CS%halocline_strat_tol = 0.25 ! Nondim. - -end subroutine init_coord_slight - -!> This subroutine deallocates memory in the control structure for the coord_slight module -subroutine end_coord_slight(CS) - type(slight_CS), pointer :: CS !< Coordinate control structure - - ! nothing to do - if (.not. associated(CS)) return - deallocate(CS%target_density) - deallocate(CS) -end subroutine end_coord_slight - -!> This subroutine can be used to set the parameters for the coord_slight module -subroutine set_slight_params(CS, max_interface_depths, max_layer_thickness, & - min_thickness, compressibility_fraction, dz_ml_min, & - nz_fixed_surface, Rho_ML_avg_depth, nlay_ML_offset, fix_haloclines, & - halocline_filter_length, halocline_strat_tol, interp_CS) - type(slight_CS), pointer :: CS !< Coordinate control structure - real, dimension(:), & - optional, intent(in) :: max_interface_depths !< Maximum depths of interfaces [H ~> m or kg m-2] - real, dimension(:), & - optional, intent(in) :: max_layer_thickness !< Maximum thicknesses of layers [H ~> m or kg m-2] - real, optional, intent(in) :: min_thickness !< Minimum thickness allowed when building the - !! new grid through regridding [H ~> m or kg m-2] - real, optional, intent(in) :: compressibility_fraction !< Fraction (between 0 and 1) of - !! compressibility to add to potential density profiles when - !! interpolating for target grid positions. [nondim] - real, optional, intent(in) :: dz_ml_min !< The fixed resolution in the topmost - !! SLight_nkml_min layers [H ~> m or kg m-2] - integer, optional, intent(in) :: nz_fixed_surface !< The number of fixed-thickness layers at the - !! top of the model - real, optional, intent(in) :: Rho_ML_avg_depth !< Depth over which to average to determine - !! the mixed layer potential density [H ~> m or kg m-2] - real, optional, intent(in) :: nlay_ML_offset !< Number of layers to offset the mixed layer - !! density to find resolved stratification [nondim] - logical, optional, intent(in) :: fix_haloclines !< If true, detect regions with much weaker than - !! based on in-situ density, and use a stretched coordinate there. - real, optional, intent(in) :: halocline_filter_length !< A length scale over which to filter T & S - !! when looking for spuriously unstable water mass profiles [H ~> m or kg m-2]. - real, optional, intent(in) :: halocline_strat_tol !< A value of the stratification ratio that - !! defines a problematic halocline region [nondim]. - type(interp_CS_type), & - optional, intent(in) :: interp_CS !< Controls for interpolation - - if (.not. associated(CS)) call MOM_error(FATAL, "set_slight_params: CS not associated") - - if (present(max_interface_depths)) then - if (size(max_interface_depths) /= CS%nk+1) & - call MOM_error(FATAL, "set_slight_params: max_interface_depths inconsistent size") - allocate(CS%max_interface_depths(CS%nk+1)) - CS%max_interface_depths(:) = max_interface_depths(:) - endif - - if (present(max_layer_thickness)) then - if (size(max_layer_thickness) /= CS%nk) & - call MOM_error(FATAL, "set_slight_params: max_layer_thickness inconsistent size") - allocate(CS%max_layer_thickness(CS%nk)) - CS%max_layer_thickness(:) = max_layer_thickness(:) - endif - - if (present(min_thickness)) CS%min_thickness = min_thickness - if (present(compressibility_fraction)) CS%compressibility_fraction = compressibility_fraction - - if (present(dz_ml_min)) CS%dz_ml_min = dz_ml_min - if (present(nz_fixed_surface)) CS%nz_fixed_surface = nz_fixed_surface - if (present(Rho_ML_avg_depth)) CS%Rho_ML_avg_depth = Rho_ML_avg_depth - if (present(nlay_ML_offset)) CS%nlay_ML_offset = nlay_ML_offset - if (present(fix_haloclines)) CS%fix_haloclines = fix_haloclines - if (present(halocline_filter_length)) CS%halocline_filter_length = halocline_filter_length - if (present(halocline_strat_tol)) then - if (halocline_strat_tol > 1.0) call MOM_error(FATAL, "set_slight_params: "//& - "HALOCLINE_STRAT_TOL must not exceed 1.0.") - CS%halocline_strat_tol = halocline_strat_tol - endif - - if (present(interp_CS)) CS%interp_CS = interp_CS -end subroutine set_slight_params - -!> Build a SLight coordinate column -subroutine build_slight_column(CS, eqn_of_state, H_to_pres, H_subroundoff, & - nz, depth, h_col, T_col, S_col, p_col, z_col, z_col_new, & - h_neglect, h_neglect_edge) - type(slight_CS), intent(in) :: CS !< Coordinate control structure - type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure - real, intent(in) :: H_to_pres !< A conversion factor from thicknesses to - !! scaled pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1] - real, intent(in) :: H_subroundoff !< GV%H_subroundoff - integer, intent(in) :: nz !< Number of levels - real, intent(in) :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2]) - real, dimension(nz), intent(in) :: T_col !< T for column [C ~> degC] - real, dimension(nz), intent(in) :: S_col !< S for column [S ~> ppt] - real, dimension(nz), intent(in) :: h_col !< Layer thicknesses [H ~> m or kg m-2] - real, dimension(nz), intent(in) :: p_col !< Layer center pressure [R L2 T-2 ~> Pa] - real, dimension(nz+1), intent(in) :: z_col !< Interface positions relative to the surface [H ~> m or kg m-2] - real, dimension(nz+1), intent(inout) :: z_col_new !< Absolute positions of interfaces [H ~> m or kg m-2] - real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose of - !! cell reconstructions [H ~> m or kg m-2]. - real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose - !! of edge value calculations [H ~> m or kg m-2]. - ! Local variables - real, dimension(nz) :: rho_col ! Layer densities [R ~> kg m-3] - real, dimension(nz) :: T_f, S_f ! Filtered layer temperature [C ~> degC] and salinity [S ~> ppt] - logical, dimension(nz+1) :: reliable ! If true, this interface is in a reliable position. - real, dimension(nz+1) :: T_int, S_int ! Temperature [C ~> degC] and salinity [S ~> ppt] interpolated to interfaces. - real, dimension(nz+1) :: rho_tmp ! A temporary density [R ~> kg m-3] - real, dimension(nz+1) :: drho_dp ! The partial derivative of density with pressure [T2 L-2 ~> kg m-3 Pa-1] - real, dimension(nz+1) :: p_IS, p_R ! Pressures [R L2 T-2 ~> Pa] - real, dimension(nz+1) :: drhoIS_dT ! The partial derivative of in situ density with temperature - ! in [R C-1 ~> kg m-3 degC-1] - real, dimension(nz+1) :: drhoIS_dS ! The partial derivative of in situ density with salinity - ! in [R S-1 ~> kg m-3 ppt-1] - real, dimension(nz+1) :: drhoR_dT ! The partial derivative of reference density with temperature - ! in [R C-1 ~> kg m-3 degC-1] - real, dimension(nz+1) :: drhoR_dS ! The partial derivative of reference density with salinity - ! in [R S-1 ~> kg m-3 ppt-1] - real, dimension(nz+1) :: strat_rat - real :: H_to_cPa ! A conversion factor from thicknesses to the compressibility fraction times - ! the units of pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1] - real :: drIS, drR ! In situ and reference density differences [R ~> kg m-3] - real :: Fn_now, I_HStol, Fn_zero_val ! Nondimensional variables [nondim] - real :: z_int_unst ! The depth where the stratification allows the interior grid to start [H ~> m or kg m-2] - real :: dz ! A uniform layer thickness in very shallow water [H ~> m or kg m-2]. - real :: dz_ur ! The total thickness of an unstable region [H ~> m or kg m-2]. - real :: wgt, cowgt ! A weight and its complement [nondim]. - real :: rho_ml_av ! The average potential density in a near-surface region [R ~> kg m-3]. - real :: H_ml_av ! A thickness to try to use in taking the near-surface average [H ~> m or kg m-2]. - real :: rho_x_z ! A cumulative integral of a density [R H ~> kg m-2 or kg2 m-5]. - real :: z_wt ! The thickness actually used in taking the near-surface average [H ~> m or kg m-2]. - real :: k_interior ! The (real) value of k where the interior grid starts [nondim]. - real :: k_int2 ! The (real) value of k where the interior grid starts [nondim]. - real :: z_interior ! The depth where the interior grid starts [H ~> m or kg m-2]. - real :: z_ml_fix ! The depth at which the fixed-thickness near-surface layers end [H ~> m or kg m-2]. - real :: dz_dk ! The thickness of layers between the fixed-thickness - ! near-surface layars and the interior [H ~> m or kg m-2]. - real :: Lfilt ! A filtering lengthscale [H ~> m or kg m-2]. - logical :: maximum_depths_set ! If true, the maximum depths of interface have been set. - logical :: maximum_h_set ! If true, the maximum layer thicknesses have been set. - real :: h_tr, b_denom_1, b1, d1 ! Temporary variables used by the tridiagonal solver. - real, dimension(nz) :: c1 ! Temporary variables used by the tridiagonal solver. - integer :: kur1, kur2 ! The indicies at the top and bottom of an unreliable region. - integer :: kur_ss ! The index to start with in the search for the next unstable region. - integer :: k, nkml - - maximum_depths_set = allocated(CS%max_interface_depths) - maximum_h_set = allocated(CS%max_layer_thickness) - - if (z_col(nz+1) - z_col(1) < nz*CS%min_thickness) then - ! This is a nearly massless total depth, so distribute the water evenly. - dz = (z_col(nz+1) - z_col(1)) / real(nz) - do K=2,nz ; z_col_new(K) = z_col(1) + dz*real(K-1) ; enddo - else - call calculate_density(T_col, S_col, p_col, rho_col, eqn_of_state) - - ! Find the locations of the target potential densities, flagging - ! locations in apparently unstable regions as not reliable. - call rho_interfaces_col(rho_col, h_col, z_col, CS%target_density, nz, & - z_col_new, CS, reliable, debug=.true., & - h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) - - ! Ensure that the interfaces are at least CS%min_thickness apart. - if (CS%min_thickness > 0.0) then - ! Move down interfaces below overly thin layers. - do K=2,nz ; if (z_col_new(K) < z_col_new(K-1) + CS%min_thickness) then - z_col_new(K) = z_col_new(K-1) + CS%min_thickness - endif ; enddo - ! Now move up any interfaces that are too close to the bottom. - do K=nz,2,-1 ; if (z_col_new(K) > z_col_new(K+1) - CS%min_thickness) then - z_col_new(K) = z_col_new(K+1) - CS%min_thickness - else - exit ! No more interfaces can be too close to the bottom. - endif ; enddo - endif - - ! Fix up the unreliable regions. - kur_ss = 2 ! reliable(1) and reliable(nz+1) must always be true. - do - ! Search for the uppermost unreliable interface postion. - kur1 = nz+2 - do K=kur_ss,nz ; if (.not.reliable(K)) then - kur1 = K ; exit - endif ; enddo - if (kur1 > nz) exit ! Everything is now reliable. - - kur2 = kur1-1 ! For error checking. - do K=kur1+1,nz+1 ; if (reliable(K)) then - kur2 = K-1 ; kur_ss = K ; exit - endif ; enddo - if (kur2 < kur1) call MOM_error(FATAL, "Bad unreliable range.") - - dz_ur = z_col_new(kur2+1) - z_col_new(kur1-1) - ! drho = CS%target_density(kur2+1) - CS%target_density(kur1-1) - ! Perhaps reset the wgt and cowgt depending on how bad the old interface - ! locations were. - wgt = 1.0 ; cowgt = 0.0 ! = 1.0-wgt - do K=kur1,kur2 - z_col_new(K) = cowgt*z_col_new(K) + & - wgt * (z_col_new(kur1-1) + dz_ur*(K - (kur1-1)) / ((kur2 - kur1) + 2)) - enddo - enddo - - ! Determine which interfaces are in the s-space region and the depth extent - ! of this region. - z_wt = 0.0 ; rho_x_z = 0.0 - H_ml_av = CS%Rho_ml_avg_depth - do k=1,nz - if (z_wt + h_col(k) >= H_ml_av) then - rho_x_z = rho_x_z + rho_col(k) * (H_ml_av - z_wt) - z_wt = H_ml_av - exit - else - rho_x_z = rho_x_z + rho_col(k) * h_col(k) - z_wt = z_wt + h_col(k) - endif - enddo - if (z_wt > 0.0) rho_ml_av = rho_x_z / z_wt - - nkml = CS%nz_fixed_surface - ! Find the interface that matches rho_ml_av. - if (rho_ml_av <= CS%target_density(nkml)) then - k_interior = CS%nlay_ml_offset + real(nkml) - elseif (rho_ml_av > CS%target_density(nz+1)) then - k_interior = real(nz+1) - else ; do K=nkml,nz - if ((rho_ml_av >= CS%target_density(K)) .and. & - (rho_ml_av < CS%target_density(K+1))) then - k_interior = (CS%nlay_ml_offset + K) + & - (rho_ml_av - CS%target_density(K)) / & - (CS%target_density(K+1) - CS%target_density(K)) - exit - endif - enddo ; endif - if (k_interior > real(nz+1)) k_interior = real(nz+1) - - ! Linearly interpolate to find z_interior. This could be made more sophisticated. - K = int(ceiling(k_interior)) - z_interior = (K-k_interior)*z_col_new(K-1) + (1.0+(k_interior-K))*z_col_new(K) - - if (CS%fix_haloclines) then - ! ! Identify regions above the reference pressure where the chosen - ! ! potential density significantly underestimates the actual - ! ! stratification, and use these to find a second estimate of - ! ! z_int_unst and k_interior. - - if (CS%halocline_filter_length > 0.0) then - Lfilt = CS%halocline_filter_length - - ! Filter the temperature and salnity with a fixed lengthscale. - h_tr = h_col(1) + H_subroundoff - b1 = 1.0 / (h_tr + Lfilt) ; d1 = h_tr * b1 - T_f(1) = (b1*h_tr)*T_col(1) ; S_f(1) = (b1*h_tr)*S_col(1) - do k=2,nz - c1(k) = Lfilt * b1 - h_tr = h_col(k) + H_subroundoff ; b_denom_1 = h_tr + d1*Lfilt - b1 = 1.0 / (b_denom_1 + Lfilt) ; d1 = b_denom_1 * b1 - T_f(k) = b1 * (h_tr*T_col(k) + Lfilt*T_f(k-1)) - S_f(k) = b1 * (h_tr*S_col(k) + Lfilt*S_f(k-1)) - enddo - do k=nz-1,1,-1 - T_f(k) = T_f(k) + c1(k+1)*T_f(k+1) ; S_f(k) = S_f(k) + c1(k+1)*S_f(k+1) - enddo - else - do k=1,nz ; T_f(k) = T_col(k) ; S_f(k) = S_col(k) ; enddo - endif - - T_int(1) = T_f(1) ; S_int(1) = S_f(1) - do K=2,nz - T_int(K) = 0.5*(T_f(k-1) + T_f(k)) ; S_int(K) = 0.5*(S_f(k-1) + S_f(k)) - p_IS(K) = z_col(K) * H_to_pres - p_R(K) = CS%ref_pressure + CS%compressibility_fraction * ( p_IS(K) - CS%ref_pressure ) - enddo - T_int(nz+1) = T_f(nz) ; S_int(nz+1) = S_f(nz) - p_IS(nz+1) = z_col(nz+1) * H_to_pres - call calculate_density_derivs(T_int, S_int, p_IS, drhoIS_dT, drhoIS_dS, & - eqn_of_state, (/2,nz/) ) - call calculate_density_derivs(T_int, S_int, p_R, drhoR_dT, drhoR_dS, & - eqn_of_state, (/2,nz/) ) - if (CS%compressibility_fraction > 0.0) then - call calculate_compress(T_int, S_int, p_R(:), rho_tmp, drho_dp, eqn_of_state, (/2,nz/)) - else - do K=2,nz ; drho_dp(K) = 0.0 ; enddo - endif - - H_to_cPa = CS%compressibility_fraction * H_to_pres - strat_rat(1) = 1.0 - do K=2,nz - drIS = drhoIS_dT(K) * (T_f(k) - T_f(k-1)) + & - drhoIS_dS(K) * (S_f(k) - S_f(k-1)) - drR = (drhoR_dT(K) * (T_f(k) - T_f(k-1)) + & - drhoR_dS(K) * (S_f(k) - S_f(k-1))) + & - drho_dp(K) * (H_to_cPa*0.5*(h_col(k) + h_col(k-1))) - - if (drIS <= 0.0) then - strat_rat(K) = 2.0 ! Maybe do this? => ; if (drR < 0.0) strat_rat(K) = -2.0 - else - strat_rat(K) = 2.0*max(drR,0.0) / (drIS + abs(drR)) - endif - enddo - strat_rat(nz+1) = 1.0 - - z_int_unst = 0.0 ; Fn_now = 0.0 - Fn_zero_val = min(2.0*CS%halocline_strat_tol, & - 0.5*(1.0 + CS%halocline_strat_tol)) - if (CS%halocline_strat_tol > 0.0) then - ! Use Adcroft's reciprocal rule. - I_HStol = 0.0 ; if (Fn_zero_val - CS%halocline_strat_tol > 0.0) & - I_HStol = 1.0 / (Fn_zero_val - CS%halocline_strat_tol) - do k=nz,1,-1 ; if (CS%ref_pressure > p_IS(k+1)) then - z_int_unst = z_int_unst + Fn_now * h_col(k) - if (strat_rat(K) <= Fn_zero_val) then - if (strat_rat(K) <= CS%halocline_strat_tol) then ; Fn_now = 1.0 - else - Fn_now = max(Fn_now, (Fn_zero_val - strat_rat(K)) * I_HStol) - endif - endif - endif ; enddo - else - do k=nz,1,-1 ; if (CS%ref_pressure > p_IS(k+1)) then - z_int_unst = z_int_unst + Fn_now * h_col(k) - if (strat_rat(K) <= CS%halocline_strat_tol) Fn_now = 1.0 - endif ; enddo - endif - - if (z_interior < z_int_unst) then - ! Find a second estimate of the extent of the s-coordinate region. - kur1 = max(int(ceiling(k_interior)),2) - if (z_col_new(kur1-1) < z_interior) then - k_int2 = kur1 - do K = kur1,nz+1 ; if (z_col_new(K) >= z_int_unst) then - ! This is linear interpolation again. - if (z_col_new(K-1) >= z_int_unst) & - call MOM_error(FATAL,"build_grid_SLight, bad halocline structure.") - k_int2 = real(K-1) + (z_int_unst - z_col_new(K-1)) / & - (z_col_new(K) - z_col_new(K-1)) - exit - endif ; enddo - if (z_col_new(nz+1) < z_int_unst) then - ! This should be unnecessary. - z_int_unst = z_col_new(nz+1) ; k_int2 = real(nz+1) - endif - - ! Now take the larger values. - if (k_int2 > k_interior) then - k_interior = k_int2 ; z_interior = z_int_unst - endif - endif - endif - endif ! fix_haloclines - - z_col_new(1) = 0.0 - do K=2,nkml+1 - z_col_new(K) = min((K-1)*CS%dz_ml_min, & - z_col_new(nz+1) - CS%min_thickness*(nz+1-K)) - enddo - z_ml_fix = z_col_new(nkml+1) - if (z_interior > z_ml_fix) then - dz_dk = (z_interior - z_ml_fix) / (k_interior - (nkml+1)) - do K=nkml+2,int(floor(k_interior)) - z_col_new(K) = z_ml_fix + dz_dk * (K - (nkml+1)) - enddo - else ! The fixed-thickness z-region penetrates into the interior. - do K=nkml+2,nz - if (z_col_new(K) <= z_col_new(CS%nz_fixed_surface+1)) then - z_col_new(K) = z_col_new(CS%nz_fixed_surface+1) - else ; exit ; endif - enddo - endif - - if (maximum_depths_set .and. maximum_h_set) then ; do k=2,nz - ! The loop bounds are 2 & nz so the top and bottom interfaces do not move. - ! Recall that z_col_new is positive downward. - z_col_new(K) = min(z_col_new(K), CS%max_interface_depths(K), & - z_col_new(K-1) + CS%max_layer_thickness(k-1)) - enddo ; elseif (maximum_depths_set) then ; do K=2,nz - z_col_new(K) = min(z_col_new(K), CS%max_interface_depths(K)) - enddo ; elseif (maximum_h_set) then ; do k=2,nz - z_col_new(K) = min(z_col_new(K), z_col_new(K-1) + CS%max_layer_thickness(k-1)) - enddo ; endif - - endif ! Total thickness exceeds nz*CS%min_thickness. - -end subroutine build_slight_column - -!> Finds the new interface locations in a column of water that match the -!! prescribed target densities. -subroutine rho_interfaces_col(rho_col, h_col, z_col, rho_tgt, nz, z_col_new, & - CS, reliable, debug, h_neglect, h_neglect_edge) - integer, intent(in) :: nz !< Number of layers - real, dimension(nz), intent(in) :: rho_col !< Initial layer reference densities [R ~> kg m-3]. - real, dimension(nz), intent(in) :: h_col !< Initial layer thicknesses [H ~> m or kg m-2]. - real, dimension(nz+1), intent(in) :: z_col !< Initial interface heights [H ~> m or kg m-2]. - real, dimension(nz+1), intent(in) :: rho_tgt !< Interface target densities. - real, dimension(nz+1), intent(inout) :: z_col_new !< New interface heights [H ~> m or kg m-2]. - type(slight_CS), intent(in) :: CS !< Coordinate control structure - logical, dimension(nz+1), intent(inout) :: reliable !< If true, the interface positions - !! are well defined from a stable region. - logical, optional, intent(in) :: debug !< If present and true, do debugging checks. - real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose of - !! cell reconstructions [H ~> m or kg m-2] - real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose - !! of edge value calculations [H ~> m or kg m-2] - - real, dimension(nz+1) :: ru_max_int ! The maximum and minimum densities in - real, dimension(nz+1) :: ru_min_int ! an unstable region around an interface [R ~> kg m-3]. - real, dimension(nz) :: ru_max_lay ! The maximum and minimum densities in - real, dimension(nz) :: ru_min_lay ! an unstable region containing a layer [R ~> kg m-3]. - real, dimension(nz,2) :: ppoly_i_E ! Edge value of polynomial [R ~> kg m-3] - real, dimension(nz,2) :: ppoly_i_S ! Edge slope of polynomial [R H-1 ~> kg m-4 or m-1] - real, dimension(nz,DEGREE_MAX+1) :: ppoly_i_coefficients ! Coefficients of polynomial [R ~> kg m-3] - logical, dimension(nz) :: unstable_lay ! If true, this layer is in an unstable region. - logical, dimension(nz+1) :: unstable_int ! If true, this interface is in an unstable region. - real :: rt ! The current target density [R ~> kg m-3]. - real :: zf ! The fractional z-position within a layer of the target density [nondim]. - real :: rfn ! The target density relative to the interpolated density [R ~> kg m-3] - real :: a(5) ! Coefficients of a local polynomial minus the target density [R ~> kg m-3]. - real :: zf1, zf2 ! Two previous estimates of zf [nondim] - real :: rfn1, rfn2 ! Values of rfn at zf1 and zf2 [R ~> kg m-3] - real :: drfn_dzf ! The partial derivative of rfn with zf [R ~> kg m-3] - real :: sgn, delta_zf, zf_prev ! [nondim] - real :: tol ! The tolerance for convergence of zf [nondim] - logical :: k_found ! If true, the position has been found. - integer :: k_layer ! The index of the stable layer containing an interface. - integer :: ppoly_degree - integer :: k, k1, k1_min, itt, max_itt, m - - real :: z_sgn ! 1 or -1, depending on whether z increases with increasing K. - logical :: debugging - - debugging = .false. ; if (present(debug)) debugging = debug - max_itt = NR_ITERATIONS - tol = NR_TOLERANCE - - z_sgn = 1.0 ; if ( z_col(1) > z_col(nz+1) ) z_sgn = -1.0 - if (debugging) then - do K=1,nz - if (abs((z_col(K+1) - z_col(K)) - z_sgn*h_col(k)) > & - 1.0e-14*(abs(z_col(K+1)) + abs(z_col(K)) + abs(h_col(k))) ) & - call MOM_error(FATAL, "rho_interfaces_col: Inconsistent z_col and h_col") - enddo - endif - - if ( z_col(1) == z_col(nz+1) ) then - ! This is a massless column! - do K=1,nz+1 ; z_col_new(K) = z_col(1) ; reliable(K) = .true. ; enddo - return - endif - - ! This sets up the piecewise polynomials based on the rho_col profile. - call regridding_set_ppolys(CS%interp_CS, rho_col, nz, h_col, ppoly_i_E, ppoly_i_S, & - ppoly_i_coefficients, ppoly_degree, h_neglect, h_neglect_edge) - - ! Determine the density ranges of unstably stratified segments. - ! Interfaces that start out in an unstably stratified segment can - ! only escape if they are outside of the bounds of that segment, and no - ! interfaces are ever mapped into an unstable segment. - unstable_int(1) = .false. - ru_max_int(1) = ppoly_i_E(1,1) - - unstable_lay(1) = (ppoly_i_E(1,1) > ppoly_i_E(1,2)) - ru_max_lay(1) = max(ppoly_i_E(1,1), ppoly_i_E(1,2)) - - do K=2,nz - unstable_int(K) = (ppoly_i_E(k-1,2) > ppoly_i_E(k,1)) - ru_max_int(K) = max(ppoly_i_E(k-1,2), ppoly_i_E(k,1)) - ru_min_int(K) = min(ppoly_i_E(k-1,2), ppoly_i_E(k,1)) - if (unstable_int(K) .and. unstable_lay(k-1)) & - ru_max_int(K) = max(ru_max_lay(k-1), ru_max_int(K)) - - unstable_lay(k) = (ppoly_i_E(k,1) > ppoly_i_E(k,2)) - ru_max_lay(k) = max(ppoly_i_E(k,1), ppoly_i_E(k,2)) - ru_min_lay(k) = min(ppoly_i_E(k,1), ppoly_i_E(k,2)) - if (unstable_lay(k) .and. unstable_int(K)) & - ru_max_lay(k) = max(ru_max_int(K), ru_max_lay(k)) - enddo - unstable_int(nz+1) = .false. - ru_min_int(nz+1) = ppoly_i_E(nz,2) - - do K=nz,1,-1 - if (unstable_lay(k) .and. unstable_int(K+1)) & - ru_min_lay(k) = min(ru_min_int(K+1), ru_min_lay(k)) - - if (unstable_int(K) .and. unstable_lay(k)) & - ru_min_int(K) = min(ru_min_lay(k), ru_min_int(K)) - enddo - - z_col_new(1) = z_col(1) ; reliable(1) = .true. - k1_min = 1 - do K=2,nz ! Find the locations of the various target densities for the interfaces. - rt = rho_tgt(K) - k_layer = -1 - k_found = .false. - - ! Many light layers are found at the top, so start there. - if (rt <= ppoly_i_E(k1_min,1)) then - z_col_new(K) = z_col(k1_min) - k_found = .true. - ! Do not change k1_min for the next layer. - elseif (k1_min == nz+1) then - z_col_new(K) = z_col(nz+1) - else - ! Start with the previous location and search outward. - if (unstable_int(K) .and. (rt >= ru_min_int(K)) .and. (rt <= ru_max_int(K))) then - ! This interface started in an unstable region and should not move due to remapping. - z_col_new(K) = z_col(K) ; reliable(K) = .false. - k1_min = K ; k_found = .true. - elseif ((rt >= ppoly_i_E(k-1,2)) .and. (rt <= ppoly_i_E(k,1))) then - ! This interface is already in the right place and does not move. - z_col_new(K) = z_col(K) ; reliable(K) = .true. - k1_min = K ; k_found = .true. - elseif (rt < ppoly_i_E(k-1,2)) then ! Search upward - do k1=K-1,k1_min,-1 - ! Check whether rt is in layer k. - if ((rt < ppoly_i_E(k1,2)) .and. (rt > ppoly_i_E(k1,1))) then - ! rt is in layer k. - k_layer = k1 - k1_min = k1 ; k_found = .true. ; exit - elseif (unstable_lay(k1) .and. (rt >= ru_min_lay(k1)) .and. (rt <= ru_max_lay(K1))) then - ! rt would be found at unstable layer that it can not penetrate. - ! It is possible that this can never happen? - z_col_new(K) = z_col(K1+1) ; reliable(K) = .false. - k1_min = k1 ; k_found = .true. ; exit - endif - ! Check whether rt is at interface K. - if (k1 > 1) then ; if ((rt <= ppoly_i_E(k1,1)) .and. (rt >= ppoly_i_E(k1-1,2))) then - ! rt is at interface K1 - z_col_new(K) = z_col(K1) ; reliable(K) = .true. - k1_min = k1 ; k_found = .true. ; exit - elseif (unstable_int(K1) .and. (rt >= ru_min_int(k1)) .and. (rt <= ru_max_int(K1))) then - ! rt would be found at an unstable interface that it can not pass. - ! It is possible that this can never happen? - z_col_new(K) = z_col(K1) ; reliable(K) = .false. - k1_min = k1 ; k_found = .true. ; exit - endif ; endif - enddo - - if (.not.k_found) then - ! This should not happen unless k1_min = 1. - if (k1_min < 2) then - z_col_new(K) = z_col(k1_min) - else - z_col_new(K) = z_col(k1_min) - endif - endif - - else ! Search downward - do k1=K,nz - if ((rt < ppoly_i_E(k1,2)) .and. (rt > ppoly_i_E(k1,1))) then - ! rt is in layer k. - k_layer = k1 - k1_min = k1 ; k_found = .true. ; exit - elseif (unstable_lay(k1) .and. (rt >= ru_min_lay(k1)) .and. (rt <= ru_max_lay(K1))) then - ! rt would be found at unstable layer that it can not penetrate. - ! It is possible that this can never happen? - z_col_new(K) = z_col(K1) - reliable(K) = .false. - k1_min = k1 ; k_found = .true. ; exit - endif - if (k1 < nz) then ; if ((rt <= ppoly_i_E(k1+1,1)) .and. (rt >= ppoly_i_E(k1,2))) then - ! rt is at interface K1+1 - - z_col_new(K) = z_col(K1+1) ; reliable(K) = .true. - k1_min = k1+1 ; k_found = .true. ; exit - elseif (unstable_int(K1+1) .and. (rt >= ru_min_int(k1+1)) .and. (rt <= ru_max_int(K1+1))) then - ! rt would be found at an unstable interface that it can not pass. - ! It is possible that this can never happen? - z_col_new(K) = z_col(K1+1) - reliable(K) = .false. - k1_min = k1+1 ; k_found = .true. ; exit - endif ; endif - enddo - if (.not.k_found) then - z_col_new(K) = z_col(nz+1) - if (rt >= ppoly_i_E(nz,2)) then - reliable(K) = .true. - else - reliable(K) = .false. - endif - endif - endif - - if (k_layer > 0) then ! The new location is inside of layer k_layer. - ! Note that this is coded assuming that this layer is stably stratified. - if (.not.(ppoly_i_E(k1,2) > ppoly_i_E(k1,1))) call MOM_error(FATAL, & - "build_grid_SLight: Erroneously searching for an interface in an unstratified layer.") - - ! Use the false position method to find the location (degree <= 1) or the first guess. - zf = (rt - ppoly_i_E(k1,1)) / (ppoly_i_E(k1,2) - ppoly_i_E(k1,1)) - - if (ppoly_degree > 1) then ! Iterate to find the solution. - a(:) = 0.0 ; a(1) = ppoly_i_coefficients(k_layer,1) - rt - do m=2,ppoly_degree+1 ; a(m) = ppoly_i_coefficients(k_layer,m) ; enddo - ! Bracket the root. - zf1 = 0.0 ; rfn1 = a(1) - zf2 = 1.0 ; rfn2 = a(1) + (a(2) + (a(3) + (a(4) + a(5)))) - if (rfn1 * rfn2 > 0.0) call MOM_error(FATAL, "build_grid_SLight: Bad bracketing.") - - do itt=1,max_itt - rfn = a(1) + zf*(a(2) + zf*(a(3) + zf*(a(4) + zf*a(5)))) - ! Reset one of the ends of the bracket. - if (rfn * rfn1 > 0.0) then - zf1 = zf ; rfn1 = rfn - else - zf2 = zf ; rfn2 = rfn - endif - if (rfn1 == rfn2) exit - - drfn_dzf = (a(2) + zf*(2.0*a(3) + zf*(3.0*a(4) + zf*4.0*a(5)))) - sgn = 1.0 ; if (drfn_dzf < 0.0) sgn = -1.0 - - if ((sgn*(zf - rfn) >= zf1 * abs(drfn_dzf)) .and. & - (sgn*(zf - rfn) <= zf2 * abs(drfn_dzf))) then - delta_zf = -rfn / drfn_dzf - zf = zf + delta_zf - else ! Newton's method goes out of bounds, so use a false position method estimate - zf_prev = zf - zf = ( rfn2 * zf1 - rfn1 * zf2 ) / (rfn2 - rfn1) - delta_zf = zf - zf_prev - endif - - if (abs(delta_zf) < tol) exit - enddo - endif - z_col_new(K) = z_col(k_layer) + zf * z_sgn * h_col(k_layer) - reliable(K) = .true. - endif - - endif - - enddo - z_col_new(nz+1) = z_col(nz+1) ; reliable(nz+1) = .true. - -end subroutine rho_interfaces_col - -end module coord_slight diff --git a/src/ALE/regrid_consts.F90 b/src/ALE/regrid_consts.F90 index 9fe638dd5b..0c5ccf268f 100644 --- a/src/ALE/regrid_consts.F90 +++ b/src/ALE/regrid_consts.F90 @@ -16,8 +16,6 @@ module regrid_consts integer, parameter :: REGRIDDING_SIGMA = 4 !< Sigma coordinates identifier integer, parameter :: REGRIDDING_ARBITRARY = 5 !< Arbitrary coordinates identifier integer, parameter :: REGRIDDING_HYCOM1 = 6 !< Simple HyCOM coordinates without BBL -integer, parameter :: REGRIDDING_SLIGHT = 7 !< Identifier for stretched coordinates in the - !! lightest water, isopycnal below integer, parameter :: REGRIDDING_SIGMA_SHELF_ZSTAR = 8 !< Identifiered for z* coordinates at the bottom, !! sigma-near the top integer, parameter :: REGRIDDING_ADAPTIVE = 9 !< Adaptive coordinate mode identifier @@ -31,7 +29,6 @@ module regrid_consts character(len=*), parameter :: REGRIDDING_ARBITRARY_STRING = "ARB" !< Arbitrary coordinates character(len=*), parameter :: REGRIDDING_HYCOM1_STRING = "HYCOM1" !< Hycom string character(len=*), parameter :: REGRIDDING_HYBGEN_STRING = "HYBGEN" !< Hybgen string -character(len=*), parameter :: REGRIDDING_SLIGHT_STRING = "SLIGHT" !< Hybrid S-rho string character(len=*), parameter :: REGRIDDING_SIGMA_SHELF_ZSTAR_STRING = "SIGMA_SHELF_ZSTAR" !< Hybrid z*/sigma character(len=*), parameter :: REGRIDDING_ADAPTIVE_STRING = "ADAPTIVE" !< Adaptive coordinate string character(len=*), parameter :: DEFAULT_COORDINATE_MODE = REGRIDDING_LAYER_STRING !< Default coordinate mode @@ -63,7 +60,6 @@ function coordinateMode(string) case (trim(REGRIDDING_SIGMA_STRING)); coordinateMode = REGRIDDING_SIGMA case (trim(REGRIDDING_HYCOM1_STRING)); coordinateMode = REGRIDDING_HYCOM1 case (trim(REGRIDDING_HYBGEN_STRING)); coordinateMode = REGRIDDING_HYBGEN - case (trim(REGRIDDING_SLIGHT_STRING)); coordinateMode = REGRIDDING_SLIGHT case (trim(REGRIDDING_ARBITRARY_STRING)); coordinateMode = REGRIDDING_ARBITRARY case (trim(REGRIDDING_SIGMA_SHELF_ZSTAR_STRING)); coordinateMode = REGRIDDING_SIGMA_SHELF_ZSTAR case (trim(REGRIDDING_ADAPTIVE_STRING)); coordinateMode = REGRIDDING_ADAPTIVE @@ -85,7 +81,6 @@ function coordinateUnitsI(coordMode) case (REGRIDDING_SIGMA); coordinateUnitsI = "Non-dimensional" case (REGRIDDING_HYCOM1); coordinateUnitsI = "m" case (REGRIDDING_HYBGEN); coordinateUnitsI = "m" - case (REGRIDDING_SLIGHT); coordinateUnitsI = "m" case (REGRIDDING_ADAPTIVE); coordinateUnitsI = "m" case default ; call MOM_error(FATAL, "coordinateUnts: "//& "Unrecognized coordinate mode.") @@ -121,7 +116,6 @@ logical function state_dependent_int(mode) case (REGRIDDING_SIGMA); state_dependent_int = .false. case (REGRIDDING_HYCOM1); state_dependent_int = .true. case (REGRIDDING_HYBGEN); state_dependent_int = .true. - case (REGRIDDING_SLIGHT); state_dependent_int = .true. case (REGRIDDING_ADAPTIVE); state_dependent_int = .true. case default ; call MOM_error(FATAL, "state_dependent: "//& "Unrecognized choice of coordinate.") diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90 index 99e6c386c6..9152efa9ec 100644 --- a/src/diagnostics/MOM_obsolete_params.F90 +++ b/src/diagnostics/MOM_obsolete_params.F90 @@ -73,6 +73,14 @@ subroutine find_obsolete_params(param_file) call obsolete_real(param_file, "ZSTAR_RIGID_SURFACE_THRESHOLD") call obsolete_logical(param_file, "HENYEY_IGW_BACKGROUND_NEW") + call obsolete_real(param_file, "SLIGHT_DZ_SURFACE") + call obsolete_int(param_file, "SLIGHT_NZ_SURFACE_FIXED") + call obsolete_real(param_file, "SLIGHT_SURFACE_AVG_DEPTH") + call obsolete_real(param_file, "SLIGHT_NLAY_TO_INTERIOR") + call obsolete_logical(param_file, "SLIGHT_FIX_HALOCLINES") + call obsolete_real(param_file, "HALOCLINE_FILTER_LENGTH") + call obsolete_real(param_file, "HALOCLINE_STRAT_TOL") + ! Test for inconsistent parameter settings. split = .true. ; test_logic = .false. call read_param(param_file,"SPLIT",split) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index cd5682d2d9..ff0eda6325 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -327,10 +327,6 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe call build_rho_column(get_rho_CS(remap_cs%regrid_cs), GV%ke, & GV%Z_to_H*(G%bathyT(i,j)+G%Z_ref), h(i,j,:), T(i,j,:), S(i,j,:), & eqn_of_state, zInterfaces, h_neglect, h_neglect_edge) - elseif (remap_cs%vertical_coord == coordinateMode('SLIGHT')) then -! call build_slight_column(remap_cs%regrid_cs,remap_cs%remap_cs, nz, & -! GV%Z_to_H*(G%bathyT(i,j)+G%Z_ref), sum(h(i,j,:)), zInterfaces) - call MOM_error(FATAL,"diag_remap_update: SLIGHT coordinate not coded for diagnostics yet!") elseif (remap_cs%vertical_coord == coordinateMode('HYCOM1')) then ! call build_hycom1_column(remap_cs%regrid_cs, nz, & ! GV%Z_to_H*(G%bathyT(i,j)+G%Z_ref), sum(h(i,j,:)), zInterfaces)