From 3f58f8a3abe355fa818d4f5d888f3c13c95a6070 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 3 Feb 2022 00:24:17 -0500 Subject: [PATCH 1/9] read_variable_2d modified to accept 3 or 4 dims The read_variable_2d function was previously configured to only run if the start and nread arrays matched the size of the field they were accessing. This was incompatible with the history of the function, which had previously required a fourth time axis (of one record), then was later modified to not require this axis. As a result, there are now files in use both with and without a time axis. This patch relaxes this check to ensure that the read is quasi-2d, i.e. the first two axes can read a segment of a 2d field, but will now reshape the start and nread arrays to match the field being read. Some additional checks are also added to ensure that it only reads one 2d slice. --- src/framework/MOM_io.F90 | 107 ++++++++++++++++++++++++++++++++------- 1 file changed, 89 insertions(+), 18 deletions(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 2ea19df183..8928d2e56b 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -892,9 +892,17 @@ end subroutine read_variable_1d_int !> Read a 2d array from a netCDF input file and save to a variable. !! -!! Start and nread ranks may exceed var, but must match the rank of the -!! variable in the netCDF file. This allows for reading slices of larger -!! arrays. +!! Start and nread lenths may exceed var rank. This allows for reading slices +!! of larger arrays. +!! +!! Previous versions of the model required a time axis on IO fields. This +!! constraint was dropped in later versions. As a result, versions both with +!! and without a time axis now exist. In order to support all such versions, +!! we use a reshaped version of start and nread in order to read the variable +!! as it exists in the file. +!! +!! Certain constraints are still applied to start and nread in order to ensure +!! that varname is a valid 2d array, or contains valid 2d slices. !! !! I/O occurs only on the root PE, and data is broadcast to other ranks. !! Due to potentially large memory communication and storage, this subroutine @@ -908,11 +916,40 @@ subroutine read_variable_2d(filename, varname, var, start, nread, ncid_in) integer, optional, intent(in) :: ncid_in !< netCDF ID of an opened file. !! If absent, the file is opened and closed within this routine. - integer :: ncid, varid, ndims, rc - character(len=*), parameter :: hdr = "read_variable_2d" + integer :: ncid, varid + integer :: field_ndims, dim_len + integer, allocatable :: field_dimids(:), field_shape(:) + integer, allocatable :: field_start(:), field_nread(:) + integer :: i, rc + character(len=*), parameter :: hdr = "read_variable_2d: " character(len=128) :: msg - logical :: size_mismatch + ! Validate shape of start and nread + if (present(start)) then + if (size(start) < 2) & + call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) & + // " start must have at least two dimensions.") + endif + + if (present(nread)) then + if (size(nread) < 2) & + call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) & + // " nread must have at least two dimensions.") + + if (any(nread(3:) > 1)) & + call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) & + // " nread may only read a single level in higher dimensions.") + endif + + ! Since start and nread may be reshaped, we cannot rely on netCDF to ensure + ! that their lengths are equivalent, and must do it here. + if (present(start) .and. present(nread)) then + if (size(start) /= size(nread)) & + call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) & + // " start and nread must have the same length.") + endif + + ! Open and read `varname` from `filename` if (is_root_pe()) then if (present(ncid_in)) then ncid = ncid_in @@ -923,23 +960,57 @@ subroutine read_variable_2d(filename, varname, var, start, nread, ncid_in) call get_varid(varname, ncid, filename, varid, match_case=.false.) if (varid < 0) call MOM_error(FATAL, "Unable to get netCDF varid for "//trim(varname)//& " in "//trim(filename)) - ! Verify that start(:) and nread(:) ranks match variable's dimension count - rc = nf90_inquire_variable(ncid, varid, ndims=ndims) + + ! Query for the dimensionality of the input field + rc = nf90_inquire_variable(ncid, varid, ndims=field_ndims) if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //& - " Difficulties reading "//trim(varname)//" from "//trim(filename)) + ": Difficulties reading "//trim(varname)//" from "//trim(filename)) + + ! Confirm that field is at least 2d + if (field_ndims < 2) & + call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) // " " // & + trim(varname) // " from " // trim(filename) // " is not a 2D field.") + + ! If start and nread are present, then reshape them to match field dims + if (present(start) .or. present(nread)) then + allocate(field_shape(field_ndims)) + allocate(field_dimids(field_ndims)) + + rc = nf90_inquire_variable(ncid, varid, dimids=field_dimids) + if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //& + ": Difficulties reading "//trim(varname)//" from "//trim(filename)) + + do i = 1, field_ndims + rc = nf90_inquire_dimension(ncid, field_dimids(i), len=dim_len) + if (rc /= NF90_NOERR) & + call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) & + // ": Difficulties reading dimensions from " // trim(filename)) + field_shape(i) = dim_len + enddo - size_mismatch = .false. - if (present(start)) size_mismatch = size_mismatch .or. size(start) /= ndims - if (present(nread)) size_mismatch = size_mismatch .or. size(nread) /= ndims + ! Reshape start(:) and nreads(:) in case ranks differ + allocate(field_start(field_ndims)) + field_start(:) = 1 + if (present(start)) then + dim_len = min(size(start), size(field_start)) + field_start(:dim_len) = start(:dim_len) + endif + + allocate(field_nread(field_ndims)) + field_nread(:2) = field_shape(:2) + field_nread(3:) = 1 + if (present(nread)) field_shape(:2) = nread(:2) - if (size_mismatch) then - write (msg, '("'// hdr //': size(start) ", i0, " and/or size(nread) ", & - i0, " do not match ndims ", i0)') size(start), size(nread), ndims - call MOM_error(FATAL, trim(msg)) + rc = nf90_get_var(ncid, varid, var, field_start, field_nread) + + deallocate(field_start) + deallocate(field_nread) + deallocate(field_shape) + deallocate(field_dimids) + else + rc = nf90_get_var(ncid, varid, var) endif - ! NOTE: We could check additional information here (type, size, ...) - rc = nf90_get_var(ncid, varid, var, start, nread) if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //& " Difficulties reading "//trim(varname)//" from "//trim(filename)) From 56401b637942c171968cba3af309fcafa7019c63 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 22 Dec 2021 09:34:18 -0500 Subject: [PATCH 2/9] +Add MOM_check_scaling.F90 and MOM_scaling_check.F90 Added two new modules, the MOM6-specific MOM_check_scaling.F90 and the generic framework module MOM_scaling_check.F90, to assess the uniqueness of the unit scaling factors for all of the variables used by MOM6. If there are overlaps in scaling factors for different units, this also identifies and suggests alternate scaling factors with less overlaps. This commit includes the introduction of the new publicly visible routines check_scaling_factors(), scales_to_powers() and check_MOM6_scaling_factors. This new capability does not do anything for sufficiently low levels of model verbosity, and it is silent if the scaling factors are unique, or if less than 2 dimensions are being rescaled. All answers and output files are bitwise identical, but there can be additional messages to stdout. --- src/core/MOM.F90 | 3 + src/core/MOM_check_scaling.F90 | 221 +++++++++++++++++ src/framework/MOM_scaling_check.F90 | 353 ++++++++++++++++++++++++++++ 3 files changed, 577 insertions(+) create mode 100644 src/core/MOM_check_scaling.F90 create mode 100644 src/framework/MOM_scaling_check.F90 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c36c0545e1..4a82eac903 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -56,6 +56,7 @@ module MOM use MOM_ALE_sponge, only : rotate_ALE_sponge, update_ALE_sponge_field use MOM_barotropic, only : Barotropic_CS use MOM_boundary_update, only : call_OBC_register, OBC_register_end, update_OBC_CS +use MOM_check_scaling, only : check_MOM6_scaling_factors use MOM_coord_initialization, only : MOM_initialize_coord use MOM_diabatic_driver, only : diabatic, diabatic_driver_init, diabatic_CS, extract_diabatic_member use MOM_diabatic_driver, only : adiabatic, adiabatic_driver_init, diabatic_driver_end @@ -2285,6 +2286,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif ! dG_in is retained for now so that it can be used with write_ocean_geometry_file() below. + if (is_root_PE()) call check_MOM6_scaling_factors(CS%GV, US) + call callTree_waypoint("grids initialized (initialize_MOM)") call MOM_timing_init(CS) diff --git a/src/core/MOM_check_scaling.F90 b/src/core/MOM_check_scaling.F90 new file mode 100644 index 0000000000..b707b65bc9 --- /dev/null +++ b/src/core/MOM_check_scaling.F90 @@ -0,0 +1,221 @@ +!> This module is used to check the scaling factors used by the MOM6 ocean model +module MOM_check_scaling + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, assert, MOM_get_verbosity +use MOM_scaling_check, only : check_scaling_factors, scales_to_powers +use MOM_unit_scaling, only : unit_scale_type +use MOM_verticalGrid, only : verticalGrid_type + +implicit none ; private + +! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional +! consistency testing. These are noted in comments with units like Z, H, L, and T, along with +! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units +! vary with the Boussinesq approximation, the Boussinesq variant is given first. + +public check_MOM6_scaling_factors + +contains + +!> Evaluate whether the dimensional scaling factors provide unique tests for all of the combinations +!! of dimensions that are used in MOM6 (or perhaps widely used), and if they are not unique, explore +!! whether another combination of scaling factors can be found that is unique or has less common +!! cases with coinciding scaling. +subroutine check_MOM6_scaling_factors(GV, US) + type(verticalGrid_type), pointer :: GV !< The container for vertical grid data + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + + ! Local variables + integer, parameter :: ndims = 6 ! The number of rescalable dimensional factors. + real, dimension(ndims) :: scales ! An array of scaling factors for each of the basic units. + integer, dimension(ndims) :: scale_pow2 ! The powers of 2 that give each element of scales. + character(len=2), dimension(ndims) :: key + ! character(len=128) :: mesg, msg_frag + integer, allocatable :: weights(:) + character(len=80), allocatable :: descriptions(:) + ! logical :: verbose, very_verbose + integer :: n, ns, max_pow + + ! Set the names and scaling factors of the dimensions being rescaled. + key(:) = ["Z", "H", "L", "T", "R", "Q"] + scales(:) = (/ US%Z_to_m, GV%H_to_MKS, US%L_to_m, US%T_to_s, US%R_to_kg_m3, US%Q_to_J_kg /) + call scales_to_powers(scales, scale_pow2) + max_pow = 40 ! 60 + + ! The first call is just to find out how many elements are in the list of scaling combinations. + call compose_dimension_list(ns, descriptions, weights) + + allocate(descriptions(ns)) + do n=1,ns ; descriptions(n) = "" ; enddo + allocate(weights(ns), source=0) + ! This call records all the list of powers, the descriptions, and their weights. + call compose_dimension_list(ns, descriptions, weights) + + call check_scaling_factors("MOM6", descriptions, weights, key, scale_pow2, max_pow) + + deallocate(weights) + deallocate(descriptions) + +end subroutine check_MOM6_scaling_factors + + +!> This routine composes a list of the commonly used dimensional scaling factors in the MOM6 +!! code, along with weights reflecting the frequency of their occurrence in the MOM6 code or +!! other considerations of how likely the variables are be used. +subroutine compose_dimension_list(ns, des, wts) + integer, intent(out) :: ns !< The running sum of valid descriptions + character(len=*), allocatable, intent(inout) :: des(:) !< The unit descriptions that have been converted + integer, allocatable, intent(inout) :: wts(:) !< A list of the integer weights for each scaling factor, + !! perhaps the number of times it occurs in the MOM6 code. + + ns = 0 + ! Accumulate a list of units used in MOM6, in approximate descending order of frequency of occurrence. + call add_scaling(ns, des, wts, "[H ~> m or kg m-2]", 1239) ! Layer thicknesses + call add_scaling(ns, des, wts, "[Z ~> m]", 660) ! Depths and vertical distance + call add_scaling(ns, des, wts, "[L T-1 ~> m s-1]", 506) ! Horizontal velocities + call add_scaling(ns, des, wts, "[R ~> kg m-3]", 356) ! Densities + call add_scaling(ns, des, wts, "[T-1 ~> s-1]", 247) ! Rates + call add_scaling(ns, des, wts, "[T ~> s]", 237) ! Time intervals + call add_scaling(ns, des, wts, "[R L2 T-2 ~> Pa]", 231) ! Dynamic pressure + ! call add_scaling(ns, des, wts, "[R L2 T-2 ~> J m-3]") ! Energy density + call add_scaling(ns, des, wts, "[Z2 T-1 ~> m2 s-1]", 181) ! Vertical viscosities and diffusivities + call add_scaling(ns, des, wts, "[H L2 ~> m3 or kg]", 174) ! Cell volumes or masses + call add_scaling(ns, des, wts, "[H L2 T-1 ~> m3 s-1 or kg s-1]", 163) ! Volume or mass transports + call add_scaling(ns, des, wts, "[L T-2 ~> m s-2]", 136) ! Horizontal accelerations + call add_scaling(ns, des, wts, "[L ~> m]", 107) ! Horizontal distances + call add_scaling(ns, des, wts, "[Z T-1 ~> m s-1]", 104) ! Friction velocities and viscous coupling + call add_scaling(ns, des, wts, "[H-1 ~> m-1 or m2 kg-1]", 89) ! Inverse cell thicknesses + call add_scaling(ns, des, wts, "[L2 T-2 ~> m2 s-2]", 88) ! Resolved kinetic energy per unit mass + call add_scaling(ns, des, wts, "[R Z3 T-2 ~> J m-2]", 85) ! Integrated turbulent kinetic energy density + call add_scaling(ns, des, wts, "[L2 T-1 ~> m2 s-1]", 78) ! Horizontal viscosity or diffusivity + call add_scaling(ns, des, wts, "[T-2 ~> s-2]", 69) ! Squared shears and buoyancy frequency + call add_scaling(ns, des, wts, "[H L ~> m2 or kg m-1]", 68) ! Lateral cell face areas + call add_scaling(ns, des, wts, "[L2 ~> m2]", 67) ! Horizontal areas + + call add_scaling(ns, des, wts, "[R-1 ~> m3 kg-1]", 61) ! Specific volumes + call add_scaling(ns, des, wts, "[Q R Z T-1 ~> W m-2]", 62) ! Vertical heat fluxes + call add_scaling(ns, des, wts, "[Z-1 ~> m-1]", 60) ! Inverse vertical distances + call add_scaling(ns, des, wts, "[L2 Z-1 T-2 ~> m s-2]", 57) ! Gravitational acceleration + call add_scaling(ns, des, wts, "[R Z T-1 ~> kg m-2 s-1]", 52) ! Vertical mass fluxes + call add_scaling(ns, des, wts, "[H T-1 ~> m s-1 or kg m-2 s-1]", 51) ! Vertical thickness fluxes + call add_scaling(ns, des, wts, "[R Z3 T-3 ~> W m-2]", 45) ! Integrated turbulent kinetic energy sources + call add_scaling(ns, des, wts, "[R Z ~> kg m-2]", 42) ! Layer or column mass loads + call add_scaling(ns, des, wts, "[Z3 T-3 ~> m3 s-3]", 33) ! Integrated turbulent kinetic energy sources + call add_scaling(ns, des, wts, "[H2 ~> m2 or kg2 m-4]", 35) ! Squared layer thicknesses + call add_scaling(ns, des, wts, "[Z2 T-2 ~> m2 s-2]", 33) ! Turbulent kinetic energy + call add_scaling(ns, des, wts, "[L-1 ~> m-1]", 32) ! Inverse horizontal distances + call add_scaling(ns, des, wts, "[R L Z T-2 ~> Pa]", 27) ! Wind stresses + call add_scaling(ns, des, wts, "[T2 L-2 ~> s2 m-2]", 33) ! Inverse velocities squared + call add_scaling(ns, des, wts, "[R Z L2 T-2 ~> J m-2]", 25) ! Integrated energy + ! call add_scaling(ns, des, wts, "[R L2 Z T-2 ~> Pa m]") ! Depth integral of pressures (25) + call add_scaling(ns, des, wts, "[Z L2 T-2 ~> m3 s-2]", 25) ! Integrated energy + call add_scaling(ns, des, wts, "[H R ~> kg m-2 or kg2 m-5]", 24) ! Layer-integrated density + call add_scaling(ns, des, wts, "[L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1]", 20) ! pbce or gtot + call add_scaling(ns, des, wts, "[L-1 T-1 ~> m-1 s-1]", 19) ! Laplacian of velocity + + call add_scaling(ns, des, wts, "[L4 T-1 ~> m4 s-1]", 18) ! Biharmonic viscosity + call add_scaling(ns, des, wts, "[Z L T-1 ~> m2 s-1]", 17) ! Layer integrated velocities + call add_scaling(ns, des, wts, "[Z L-1 ~> nondim]", 15) ! Slopes + call add_scaling(ns, des, wts, "[Z L2 ~> m3]", 14) ! Diagnostic volumes + call add_scaling(ns, des, wts, "[H L T-1 ~> m2 s-1 or kg m-1 s-1]", 12) ! Layer integrated velocities + call add_scaling(ns, des, wts, "[L2 T-3 ~> m2 s-3]", 14) ! Buoyancy flux or MEKE sources [L2 T-3 ~> W kg-1] + call add_scaling(ns, des, wts, "[Z2 ~> m2]", 12) ! Squared vertical distances + call add_scaling(ns, des, wts, "[R Z L2 T-1 ~> kg s-1]", 12) ! Mass fluxes + call add_scaling(ns, des, wts, "[L-2 ~> m-2]", 12) ! Inverse areas + call add_scaling(ns, des, wts, "[L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]", 11) ! Gravitational acceleration over density + call add_scaling(ns, des, wts, "[Z T-2 ~> m s-2]", 10) ! Buoyancy differences or their derivatives + ! Could also add [Z T-2 degC-1 ~> m s-2 degC-1] or [Z T-2 ppt-1 ~> m s-2 ppt-1] + call add_scaling(ns, des, wts, "[R Z L2 T-3 ~> W m-2]", 10) ! Energy sources, including for MEKE + call add_scaling(ns, des, wts, "[L3 ~> m3]", 10) ! Metric dependent constants for viscosity + call add_scaling(ns, des, wts, "[Z-2 ~> m-2]", 9) ! Inverse of denominator in some weighted averages + call add_scaling(ns, des, wts, "[H-2 ~> m-2 or m4 kg-2]", 9) ! Mixed layer local work variables + call add_scaling(ns, des, wts, "[Z L2 T-1 ~> m3 s-1]", 9) ! Overturning (GM) streamfunction + call add_scaling(ns, des, wts, "[L2 Z-2 T-2 ~> s-2]", 9) ! Buoyancy frequency in some params. + call add_scaling(ns, des, wts, "[Q R Z ~> J m-2]", 8) ! time-integrated frazil heat flux + call add_scaling(ns, des, wts, "[Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]", 7) ! (TKE_to_Kd) + call add_scaling(ns, des, wts, "[Q degC-1 ~> J kg-1 degC-1]", 7) ! Heat capacity + + call add_scaling(ns, des, wts, "[R Z2 T-2 ~> J m-3]", 6) ! Potential energy height derivatives + call add_scaling(ns, des, wts, "[R Z3 T-2 H-1 ~> J m-3 or J kg-1]", 7) ! Partial derivatives of energy + call add_scaling(ns, des, wts, "[R L2 T-2 Z-1 ~> Pa m-1]", 7) ! Converts depth to pressure + call add_scaling(ns, des, wts, "[L4 Z-1 T-1 ~> m3 s-1]", 7) ! Rigidity of ice + call add_scaling(ns, des, wts, "[H L2 T-3 ~> m3 s-3]", 9) ! Kinetic energy diagnostics + call add_scaling(ns, des, wts, "[H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1]", 6) ! Layer potential vorticity + call add_scaling(ns, des, wts, "[R Z2 T-3 ~> W m-3]", 3) ! Kinetic energy dissipation rates + call add_scaling(ns, des, wts, "[Z2 L-2 ~> 1]", 1) ! Slopes squared + call add_scaling(ns, des, wts, "[Z H-1 ~> nondim or m3 kg-1]", 6) ! Thickness to height conversion + call add_scaling(ns, des, wts, "[Pa T2 R-1 L-2 ~> 1]", 6) ! Pressure conversion factor + ! Could also add [m T2 R-1 L-2 ~> m Pa-1] + ! Could also add [degC T2 R-1 L-2 ~> degC Pa-1] + call add_scaling(ns, des, wts, "[R H-1 ~> kg m-4 or m-1]", 5) ! Vertical density gradients + call add_scaling(ns, des, wts, "[R L4 T-4 ~> Pa m2 s-2]", 4) ! Integral in geopotential of pressure + call add_scaling(ns, des, wts, "[L Z-1 ~> nondim]", 4) ! Inverse slopes + call add_scaling(ns, des, wts, "[L-3 ~> m-3]", 4) ! Metric dependent constants for viscosity + call add_scaling(ns, des, wts, "[H T2 L-1 ~> s2 or kg s2 m-3]", 4) ! BT_cont_type face curvature fit + call add_scaling(ns, des, wts, "[H L-1 ~> nondim or kg m-3]", 4) ! BT_cont_type face curvature fit + call add_scaling(ns, des, wts, "[kg H-1 L-2 ~> kg m-3 or 1]", 20) ! Diagnostic conversions to mass + ! Could also add [m3 H-1 L-2 ~> 1 or m3 kg-1] + call add_scaling(ns, des, wts, "[Z T-2 R-1 ~> m4 s-2 kg-1]", 9) ! Gravitational acceleration over density + call add_scaling(ns, des, wts, "[R Z L4 T-3 ~> kg m2 s-3]", 9) ! MEKE fluxes + call add_scaling(ns, des, wts, "[R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]", 3) ! Thickness to pressure conversion + + call add_scaling(ns, des, wts, "[R-1 Z-1 ~> m2 kg-1]", 3) ! Inverse of column mass + call add_scaling(ns, des, wts, "[L4 ~> m4]", 3) ! Metric dependent constants for viscosity + call add_scaling(ns, des, wts, "[T-1 Z-1 ~> s-1 m-1]", 2) ! Barotropic PV, for some options + call add_scaling(ns, des, wts, "[R Z2 T-1 ~> J s m-3]", 2) ! River mixing term [R Z2 T-1 ~> Pa s] + call add_scaling(ns, des, wts, "[degC Q-1 ~> kg degC J-1]", 2) ! Inverse heat capacity + ! Could add call add_scaling(ns, des, wts, "[Q-1 ~> kg J-1]", 1) ! Inverse heat content + call add_scaling(ns, des, wts, "[L4 Z-2 T-1 ~> m2 s-1]", 2) ! Ice rigidity term + call add_scaling(ns, des, wts, "[R Z-1 ~> kg m-4]", 3) ! Vertical density gradient + call add_scaling(ns, des, wts, "[R Z L2 ~> kg]", 3) ! Depth and time integrated mass fluxes + call add_scaling(ns, des, wts, "[R L2 T-3 ~> W m-2]", 3) ! Depth integrated friction work + call add_scaling(ns, des, wts, "[ppt2 R-2 ~> ppt2 m6 kg-2]", 3) ! T / S gauge transformation + call add_scaling(ns, des, wts, "[R L-1 ~> kg m-4]", 2) ! Horizontal density gradient + ! Could add call add_scaling(ns, des, wts, "[H Z ~> m2 or kg m-1]", 2) ! Temporary variables + call add_scaling(ns, des, wts, "[Z3 R2 T-2 H-2 ~> kg2 m-5 s-2 or m s-2]", 2) ! Heating to PE change + call add_scaling(ns, des, wts, "[R2 L2 Z2 T-4 ~> Pa2]", 2) ! Squared wind stresses + call add_scaling(ns, des, wts, "[L-2 T-2 ~> m-2 s-2]", 2) ! Squared Laplacian of velocity + call add_scaling(ns, des, wts, "[T H Z-1 ~> s or s kg m-3]", 2) ! Time step times thickness conversion + call add_scaling(ns, des, wts, "[T H Z-1 R-1 ~> s m3 kg-1 or s]", 2) ! Time step over density with conversion + call add_scaling(ns, des, wts, "[H-3 ~> m-3 or m6 kg-3]", 1) ! A local term in ePBL + call add_scaling(ns, des, wts, "[H-4 ~> m-4 or m8 kg-4]", 1) ! A local term in ePBL + call add_scaling(ns, des, wts, "[H T Z-2 ~> s m-1 or kg s m-4]", 1) ! A local term in ePBL + + call add_scaling(ns, des, wts, "[H3 ~> m3 or kg3 m-6]", 1) ! Thickness cubed in a denominator + call add_scaling(ns, des, wts, "[H2 T-2 ~> m2 s-2 or kg2 m-4 s-2]", 1) ! Thickness times f squared + call add_scaling(ns, des, wts, "[H T2 R-1 Z-2 ~> m Pa-1 or s2 m-1]", 1) ! Pressure to thickness conversion + call add_scaling(ns, des, wts, "[L2 Z-2 ~> nondim]", 1) ! Inverse slope squared + call add_scaling(ns, des, wts, "[H R L2 T-2 ~> m Pa]", 1) ! Integral in thickness of pressure + call add_scaling(ns, des, wts, "[R T2 Z-1 ~> kg s2 m-4]", 1) ! Density divided by gravitational acceleration + +end subroutine compose_dimension_list + +!> Augment the count the valid unit descriptions, and add the provided description and its weight +!! to the end of the list if that list is allocated. +subroutine add_scaling(ns, descs, wts, scaling, weight) + integer, intent(inout) :: ns !< The running sum of valid descriptions. + character(len=*), allocatable, intent(inout) :: descs(:) !< The unit descriptions that have been converted + integer, allocatable, intent(inout) :: wts(:) !< A list of the integers for each scaling + character(len=*), intent(in) :: scaling !< The unit description that will be converted + integer, optional, intent(in) :: weight !< An optional weight or occurrence count + !! for this unit description, 1 by default. + + integer :: iend + + iend = index(scaling, "~>") + if (iend <= 1) then + call MOM_mesg("No scaling indicator ~> found for "//trim(scaling)) + else + ! Count and perhaps store this description and its weight. + ns = ns + 1 + if (allocated(descs)) descs(ns) = scaling + if (allocated(wts)) then + wts(ns) = 1 ; if (present(weight)) wts(ns) = weight + endif + endif + +end subroutine add_scaling + +end module MOM_check_scaling diff --git a/src/framework/MOM_scaling_check.F90 b/src/framework/MOM_scaling_check.F90 new file mode 100644 index 0000000000..e32e668b17 --- /dev/null +++ b/src/framework/MOM_scaling_check.F90 @@ -0,0 +1,353 @@ +!> This module is used to check the scaling factors used by the MOM6 ocean model +module MOM_scaling_check + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, assert, MOM_get_verbosity + +implicit none ; private + +! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional +! consistency testing. These are noted in comments with units like Z, H, L, and T, along with +! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units +! vary with the Boussinesq approximation, the Boussinesq variant is given first. + +public check_scaling_factors, scales_to_powers + +contains + +!> This subroutine does a checks whether the provided dimensional scaling factors give a unique +!! overall scaling for each of the combinations of units in description, and suggests a better +!! combination if it is not unique. However, this subroutine does nothing if the verbosity level +!! for this run is below 3. +subroutine check_scaling_factors(component, descs, weights, key, scales, max_powers) + character(len=*), intent(in) :: component !< The name of the component (e.g., MOM6) to use in messages + character(len=*), intent(in) :: descs(:) !< The descriptions for each combination of units + integer, intent(in) :: weights(:) !< A list of the weights for each described combination + character(len=*), intent(in) :: key(:) !< The key for the unit scaling + integer, intent(in) :: scales(:) !< The powers of 2 that give the scaling for each unit in key + integer, optional, intent(in) :: max_powers !< The maximum range of powers of 2 to search for + !! suggestions of better scaling factors, or 0 to avoid + !! suggesting improved factors. + + ! Local variables + integer, dimension(size(key)) :: next_scales, prev_scales, better_scales + character(len=512) :: mesg + character(len=64) :: msg_frag + integer, dimension(size(key), size(weights)) :: list + integer :: verbosity + logical :: same_key + integer :: orig_cost, test_cost, better_cost, prev_cost ! Various squared-weight mismatch costs. + integer :: better_dp ! The absolute change in powers with the better estimate. + integer :: ndims, ns, m, n, i, p, itt, max_itt, max_pow + + call assert((size(scales) == size(key)), "check_scaling_factors: Mismatched scales and key sizes.") + call assert((size(descs) == size(weights)), "check_scaling_factors: Mismatched descs and weights.") + + verbosity = MOM_get_verbosity() + ! Skip the rest of this routine if it would not write anything out. + if (verbosity < 3) return + + ndims = size(key) + ns = size(weights) + max_pow = 0 ; if (present(max_powers)) max_pow = max_powers + + list(:,:) = 0 + do n=1,ns + call encode_dim_powers(descs(n), key, list(:,n)) + enddo + + if (verbosity >= 7) then + write(mesg, '(I8)') ns + call MOM_mesg(trim(component)//": Extracted "//trim(adjustl(mesg))//" unit combinations from the list.") + mesg = "Dim Key: [" + do i=1,ndims ; mesg = trim(mesg)//" "//trim(key(i)) ; enddo + mesg = trim(mesg)//"]:" + call MOM_mesg(mesg) + do n=1,ns + call MOM_mesg(trim(component)//": Extracted ["//trim(int_array_msg(list(:,n)))//"] from "//trim(descs(n))) + enddo + + do n=1,ns ; do m=1,n-1 + same_key = .true. + do i=1,ndims ; if (list(i,n) /= list(i,m)) same_key = .false. ; enddo + if (same_key) then + call MOM_mesg(trim(component)//": The same powers occur for "//& + trim(descs(n))//" and "//trim(descs(m))//"." ) + endif + enddo ; enddo + endif + + orig_cost = non_unique_scales(scales, list, descs, weights, silent=(verbosity<4)) + + max_itt = 3*ndims ! Do up to 3 iterations for each rescalable dimension. + if (orig_cost /= 0) then + call MOM_mesg(trim(component)//": The dimensional scaling factors are not unique.") + prev_cost = orig_cost + prev_scales(:) = scales(:) + do itt=1,max_itt + ! Iterate to find a better solution. + better_scales(:) = prev_scales(:) + better_cost = prev_cost + better_dp = 0 + do i=1,ndims + if (scales(i) == 0) cycle ! DO not optimize unscaled dimensions. + next_scales(:) = prev_scales(:) + do p=-max_pow,max_pow + if ((p==0) .or. (p==prev_scales(i))) cycle + next_scales(i) = p + test_cost = non_unique_scales(next_scales, list, descs, weights, silent=.true.) + if ((test_cost < better_cost) .or. & + ((test_cost == better_cost) .and. (abs(p-prev_scales(i)) < better_dp))) then + ! This is a better scaling or has the same weighted mismatches but smaller + ! changes in rescaling factors, so it could be the next guess. + better_scales(:) = next_scales(:) + better_cost = test_cost + better_dp = abs(p - prev_scales(i)) + endif + enddo + enddo + if (better_cost < prev_cost) then + ! Store the new best guess and try again. + prev_scales(:) = better_scales(:) + prev_cost = better_cost + else ! No further optimization is possible. + exit + endif + if (better_cost == 0) exit + if (verbosity >= 7) then + write(mesg, '("Iteration ",I2," scaling cost reduced from ",I8," with original scales to ", I8)') & + itt, orig_cost, better_cost + call MOM_mesg(trim(component)//": "//trim(mesg)//" with revised scaling factors.") + endif + enddo + if (prev_cost < orig_cost) then + test_cost = non_unique_scales(prev_scales, list, descs, weights, silent=(verbosity<4)) + mesg = trim(component)//": Suggested improved scales: " + do i=1,ndims ; if ((prev_scales(i) /= scales(i)) .and. (scales(i) /= 0)) then + write(msg_frag, '(I3)') prev_scales(i) + mesg = trim(mesg)//" "//trim(key(i))//"_RESCALE_POWER = "//trim(adjustl(msg_frag)) + endif ; enddo + call MOM_mesg(mesg) + + write(mesg, '(I8)') orig_cost + write(msg_frag, '(I8)') test_cost + mesg = trim(component)//": Scaling overlaps reduced from "//trim(adjustl(mesg))//& + " with original scales to "//trim(adjustl(msg_frag))//" with suggested scales." + call MOM_mesg(mesg) + endif + + endif + +end subroutine check_scaling_factors + +!> Convert a unit scaling descriptor into an array of the dimensions of powers given in the key +subroutine encode_dim_powers(scaling, key, dim_powers) + + character(len=*), intent(in) :: scaling !< The unit description that will be converted + character(len=*), dimension(:), intent(in) :: key(:) !< The key for the unit scaling + integer, dimension(size(key)), intent(out) :: dim_powers !< The dimensions in scaling of each + !! element of they key. + + ! Local variables + character(len=:), allocatable :: actstr ! The full active remaining string to be parsed. + character(len=:), allocatable :: fragment ! The space-delimited fragment being parsed. + character(len=:), allocatable :: dimnm ! The probable dimension name + character(len=11) :: numbers ! The list of characters that could make up the exponent. + ! character(len=128) :: mesg + integer :: istart, iend, ieq, ief, ipow ! Positions in strings. + integer :: dp ! The power for this dimension. + integer :: ndim ! The number of dimensional scaling factors to consider. + integer :: n + + dim_powers(:) = 0 + + iend = index(scaling, "~>") - 1 + if (iend < 1) return + + ! Parse the key. + ndim = size(key) + numbers = "-0123456789" + + ! Strip away any leading square brace. + istart = index(scaling(:iend), "[") + 1 + ! If there is an "=" in the string, start after this. + ieq = index(scaling(istart:iend), "=", back=.true.) + if (ieq > 0) istart = istart + ieq + + ! Set up the active string to work on. + actstr = trim(adjustl(scaling(istart:iend))) + do ! Loop over each of the elements in the unit scaling descriptor. + if (len_trim(actstr) == 0) exit + ief = index(actstr, " ") - 1 + if (ief <= 0) ief = len_trim(actstr) + fragment = actstr(:ief) + ipow = scan(fragment, "-") + if (ipow == 0) ipow = scan(fragment, numbers) + + if (ipow == 0) then ! There is no exponent + dimnm = fragment + dp = 1 + ! call MOM_mesg("Parsing powerless fragment "//trim(fragment)//" from "//trim(scaling)) + else + if (verify(fragment(ipow:), numbers) == 0) then + read(fragment(ipow:),*) dp + dimnm = fragment(:ipow-1) + ! write(mesg, '(I3)') dp + ! call MOM_mesg("Parsed fragment "//trim(fragment)//" from "//trim(scaling)//& + ! " as "//trim(dimnm)//trim(adjustl(mesg))) + else + dimnm = fragment + dp = 1 + ! call MOM_mesg("Unparsed fragment "//trim(fragment)//" from "//trim(scaling)) + endif + endif + + do n=1,ndim + if (trim(dimnm) == trim(key(n))) then + dim_powers(n) = dim_powers(n) + dp + exit + endif + enddo + + ! Remove the leading fragment that has been parsed from actstr + actstr = trim(adjustl(actstr(ief+1:))) + enddo + +end subroutine encode_dim_powers + +!> Find the integer power of two that describe each of the scaling factors, or return 0 for +!! scaling factors that are not exceptionally close to an integer power of 2. +subroutine scales_to_powers(scale, pow2) + real, intent(in) :: scale(:) !< The scaling factor for each dimension + integer, intent(out) :: pow2(:) !< The exact powers of 2 for each scale, or 0 for non-exact powers of 2. + + real :: log2_sc ! The log base 2 of an element of scale + integer :: n, ndim + + ndim = size(scale) + + ! Find the integer power of two for the scaling factors, but skip the analysis of any factors + ! that are not close enough to being integer powers of 2. + do n=1,ndim + if (abs(scale(n)) > 0.0) then + log2_sc = log(abs(scale(n))) / log(2.0) + else + log2_sc = 0.0 + endif + if (abs(log2_sc - nint(log2_sc)) < 1.0e-6) then + ! This is close to an integer power of two. + pow2(n) = nint(log2_sc) + else + ! This is not being scaled by an integer power of 2, so return 0. + pow2(n) = 0 + endif + enddo + +end subroutine scales_to_powers + +!> Determine from the list of scaling factors and the unit combinations that are in use whether +!! all these combinations scale uniquely. +integer function non_unique_scales(scales, list, descs, weights, silent) + integer, intent(in) :: scales(:) !< The power of 2 that gives the scaling factor for each dimension + integer, intent(in) :: list(:,:) !< A list of the integers for each scaling + character(len=*), intent(in) :: descs(:) !< The unit descriptions that have been converted + integer, intent(in) :: weights(:) !< A list of the weights for each scaling + logical, optional, intent(in) :: silent !< If present and true, do not write any output. + + ! Local variables + integer, dimension(size(weights)) :: res_pow ! The net rescaling power for each combination. + integer, dimension(size(weights)) :: wt_merge ! The merged weights of scaling factors with common powers + ! for the dimensions being tested. + logical :: same_key, same_scales, verbose + character(len=256) :: mesg + integer :: nonzero_count ! The number of non-zero scaling factors + integer :: ndim ! The number of dimensional scaling factors to work with + integer :: i, n, m, ns + + verbose = .true. ; if (present(silent)) verbose = .not.silent + + ndim = size(scales) + ns = size(descs) + call assert((size(scales) == size(list, 1)), "non_unique_scales: Mismatched scales and list sizes.") + call assert((size(descs) == size(list, 2)), "non_unique_scales: Mismatched descs and list sizes.") + call assert((size(descs) == size(weights)), "non_unique_scales: Mismatched descs and weights.") + + ! Return .true. if all scaling powers are 0, or there is only one scaling factor in use. + nonzero_count = 0 ; do n=1,ndim ; if (scales(n) /= 0) nonzero_count = nonzero_count + 1 ; enddo + if (nonzero_count <= 1) return + + ! Figure out which unit combinations are unique for the set of dimensions and scaling factors + ! that are being tested, and combine the weights for scaling factors. + wt_merge(:) = weights(:) + do n=1,ns ; do m=1,n-1 + same_key = .true. + same_scales = .true. + do i=1,ndim + if (list(i,n) /= list(i,m)) same_key = .false. + if ((scales(i) /= 0) .and. (list(i,n) /= list(i,m))) same_scales = .false. + enddo + if (same_key .or. same_scales) then + if (wt_merge(n) > wt_merge(m)) then + wt_merge(n) = wt_merge(n) + wt_merge(m) + wt_merge(m) = 0 + else + wt_merge(m) = wt_merge(m) + wt_merge(n) + wt_merge(n) = 0 + endif + endif + if (wt_merge(n) == 0) exit ! Go to the next value of n. + enddo ; enddo + + do n=1,ns + res_pow(n) = 0 + do i=1,ndim + res_pow(n) = res_pow(n) + scales(i) * list(i,n) + enddo + enddo + + ! Determine the weighted cost of non-unique scaling factors. + non_unique_scales = 0 + do n=1,ns ; if (wt_merge(n) > 0) then ; do m=1,n-1 ; if (wt_merge(m) > 0) then + if (res_pow(n) == res_pow(m)) then + ! Use the product of the weights as the cost, as this should be vaguely proportional to + ! the likelihood that these factors would be combined in an expression. + non_unique_scales = min(non_unique_scales + wt_merge(n) * wt_merge(m), 99999999) + if (verbose) then + write(mesg, '(I8)') res_pow(n) + call MOM_mesg("The factors "//trim(descs(n))//" and "//trim(descs(m))//" both scale to "//& + trim(adjustl(mesg))//" for the given powers.") + + ! call MOM_mesg("Powers ["//trim(int_array_msg(list(:,n)))//"] and ["//& + ! trim(int_array_msg(list(:,m)))//"] with rescaling by ["//& + ! trim(int_array_msg(scales))//"]") + endif + endif + endif ; enddo ; endif ; enddo + +end function non_unique_scales + +!> Return a string the elements of an array of integers +function int_array_msg(array) + integer, intent(in) :: array(:) !< The array whose values are to be written. + character(len=16*size(array)) :: int_array_msg + + character(len=12) :: msg_frag + integer :: i, ni + ni = size(array) + + int_array_msg = "" + if (ni < 1) return + + do i=1,ni + write(msg_frag, '(I8)') array(i) + msg_frag = adjustl(msg_frag) + if (i == 1) then + int_array_msg = trim(msg_frag) + else + int_array_msg = trim(int_array_msg)//" "//trim(msg_frag) + endif + enddo +end function int_array_msg + +end module MOM_scaling_check From 75bf521e29e3434965b44cea288404bb0be341af Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 24 Jan 2022 17:04:25 -0500 Subject: [PATCH 3/9] +Move MOM_scaling_check.F90 to MOM_unique_scales.F90 Renamed the framework module MOM_scaling_check.F90 to MOM_unique_scales.F90 to help differentiate it from MOM_check_scaling.F90, and renamed the subroutine check_scaling_factors() as check_scaling_uniqueness(). Also added _Dimensional_consistency.dox to describe the dimensional consistency testing. This commit should address the issues raised in the review of MOM6 PR #49. All answers and output are bitwise identical. --- src/core/MOM_check_scaling.F90 | 6 +- ...caling_check.F90 => MOM_unique_scales.F90} | 13 +-- src/framework/_Dimensional_consistency.dox | 85 +++++++++++++++++++ 3 files changed, 95 insertions(+), 9 deletions(-) rename src/framework/{MOM_scaling_check.F90 => MOM_unique_scales.F90} (97%) create mode 100644 src/framework/_Dimensional_consistency.dox diff --git a/src/core/MOM_check_scaling.F90 b/src/core/MOM_check_scaling.F90 index b707b65bc9..55bd471fee 100644 --- a/src/core/MOM_check_scaling.F90 +++ b/src/core/MOM_check_scaling.F90 @@ -1,10 +1,10 @@ -!> This module is used to check the scaling factors used by the MOM6 ocean model +!> This module is used to check the dimensional scaling factors used by the MOM6 ocean model module MOM_check_scaling ! This file is part of MOM6. See LICENSE.md for the license. use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, assert, MOM_get_verbosity -use MOM_scaling_check, only : check_scaling_factors, scales_to_powers +use MOM_unique_scales, only : check_scaling_uniqueness, scales_to_powers use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type @@ -53,7 +53,7 @@ subroutine check_MOM6_scaling_factors(GV, US) ! This call records all the list of powers, the descriptions, and their weights. call compose_dimension_list(ns, descriptions, weights) - call check_scaling_factors("MOM6", descriptions, weights, key, scale_pow2, max_pow) + call check_scaling_uniqueness("MOM6", descriptions, weights, key, scale_pow2, max_pow) deallocate(weights) deallocate(descriptions) diff --git a/src/framework/MOM_scaling_check.F90 b/src/framework/MOM_unique_scales.F90 similarity index 97% rename from src/framework/MOM_scaling_check.F90 rename to src/framework/MOM_unique_scales.F90 index e32e668b17..730d11adb0 100644 --- a/src/framework/MOM_scaling_check.F90 +++ b/src/framework/MOM_unique_scales.F90 @@ -1,5 +1,6 @@ -!> This module is used to check the scaling factors used by the MOM6 ocean model -module MOM_scaling_check +!> This module provides tools that can be used to check the uniqueness of the dimensional +!! scaling factors used by the MOM6 ocean model or other models +module MOM_unique_scales ! This file is part of MOM6. See LICENSE.md for the license. @@ -12,7 +13,7 @@ module MOM_scaling_check ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units ! vary with the Boussinesq approximation, the Boussinesq variant is given first. -public check_scaling_factors, scales_to_powers +public check_scaling_uniqueness, scales_to_powers contains @@ -20,7 +21,7 @@ module MOM_scaling_check !! overall scaling for each of the combinations of units in description, and suggests a better !! combination if it is not unique. However, this subroutine does nothing if the verbosity level !! for this run is below 3. -subroutine check_scaling_factors(component, descs, weights, key, scales, max_powers) +subroutine check_scaling_uniqueness(component, descs, weights, key, scales, max_powers) character(len=*), intent(in) :: component !< The name of the component (e.g., MOM6) to use in messages character(len=*), intent(in) :: descs(:) !< The descriptions for each combination of units integer, intent(in) :: weights(:) !< A list of the weights for each described combination @@ -139,7 +140,7 @@ subroutine check_scaling_factors(component, descs, weights, key, scales, max_pow endif -end subroutine check_scaling_factors +end subroutine check_scaling_uniqueness !> Convert a unit scaling descriptor into an array of the dimensions of powers given in the key subroutine encode_dim_powers(scaling, key, dim_powers) @@ -350,4 +351,4 @@ function int_array_msg(array) enddo end function int_array_msg -end module MOM_scaling_check +end module MOM_unique_scales diff --git a/src/framework/_Dimensional_consistency.dox b/src/framework/_Dimensional_consistency.dox new file mode 100644 index 0000000000..0657724381 --- /dev/null +++ b/src/framework/_Dimensional_consistency.dox @@ -0,0 +1,85 @@ +/*! \page Dimensional_consistency Dimensional Consistency Testing + +\section section_Dimensional_consistency Dimensional Consistency Testing + + MOM6 uses a unique system for testing the dimensional consistency of all of +its expressions. The internal representations of dimensional variables are +rescaled by integer powers of 2 that depend on their units, with all input and +output being rescaled back to their original MKS units. By choosing different +powers of 2 for different units, the internal representations with different +units scale differently, so dimensionally inconsistent expressions will not +reproduce, but dimensionally inconsistent expressions give bitwise identical +results. So, for example, if horizontal lengths scale by a factor of 2^6=64, +and time is scaled by a factor of 2^4=16, horizontal velocities will scale by a +factor of 2^(6-4)=4. In this case, expressions that combine velocities, all +terms would scale by the same factor of 4. By contrast, if there were an +expression where a variable with units of length were added to one with the +units of a velocity, the results would scale inconsistently, and answers would +change with different scaling factors. + + What makes these integer powers of 2 special is the way that floating point +numbers are written as an O(1) mantissa times 2 raised to an integer exponent +between +/-1024. Multiplication by an integer power of 2 is just an integer +shift in the exponent, so as long as the model is not rescaled by an overly +large factor to encounter overflows and the model is not relying on automatic +underflows being converted to 0, all floating point operations can be carried +with one scale, and then rescaled to obtain identical answers. MOM6 has the +option to explicitly handle all relevant cases of underflows, and it can be +demonstrated to give identical answers when each of its units are scaled by +factors ranging from 2^-140 ~= 7.2e-43 to 2^140 ~= 1.4e42. + + When running with rescaling factors other than 2^0 = 1, there are some extra +array copies and multiplies of input fields or diagnostic output, so it is +slightly more efficient not to actively use the dimensional rescaling. For +production runs, we typically set all of the rescaling powers to 0, but for +debugging code problems, this rescaling can be an invaluable tool, especially +when combined with the very verbose runtime setting DEBUG=True in a MOM_input or +MOM_override file. Diffs of the output from runs with different scaling factors +readily highlights the earliest instances of differences, which can be used to +track down any dimensionally inconsistent expressions. Similarly, dimensional +inconsistencies in diagnostics is easily tracked down by comparing the output +from a pair of runs. + + All real variables in MOM6 should have comments describing their purpose, +along with their rescaled units and their mks counterparts with notation like +"! A velocity [L T-1 ~> m s-1]". If the units vary with the Boussinesq +approximation, the Boussinesq variant is given first. When variables are read +in, their dimensions are usually specified with a 'scale=' optional argument on +the MOM_get_param or MOM_read_data call, while the unscaling of diagnostics is +specified with a 'conversion=' factor. In both cases, these arguments it next +to a text string specifying the variable's units, which can then be check easily +for self-consistency. + + Currently in MOM6, the following dimensions have unique scaling, along with +the notation used to describe these variables in comments: + +\li Time, scaled by 2^T_RESCALE_POWER, denoted as [T ~> s] +\li Horizontal length, scaled by 2^L_RESCALE_POWER, denoted as [L ~> m] +\li Vertical height, scaled by 2^Z_RESCALE_POWER, denoted as [Z ~> m] +\li Vertical thickness, scaled by 2^H_RESCALE_POWER, denoted as [H ~> m or kg m-2] +\li Density, scaled by 2^R_RESCALE_POWER, denoted as [R ~> kg m-3] +\li Enthalpy (or heat content), scaled by 2^Q_RESCALE_POWER, denoted as [Q ~> J kg-1] + + These rescaling capabilities are also used by the SIS2 sea ice model, but it +does uses a non-Boussinesq mass scale of [R Z ~> kg m-2] for ice thicknesses, +rather than having a separate scaling factor (of [H ~> m or kg m-2]) that varies +between the Boussinesq and non-Boussinesq modes like MOM6 does. The actual +powers used in the scaling are specified separately for MOM6 and SIS2 and +need not be the same. + + Each of these units can be scaled in separate test runs, or all of them can be +rescaled simultaneously. In the latter case, MOM_unique_scales.F90 provides +tools to evaluate whether the specific combinations of units used by a model +scale by unique powers, and it can suggest scaling factors that provides unique +combinations of rescaling factors for the dimensions being tested, using a +cost-function based on the frequency with which units are used in the model (and +specified inside of MOM_check_scaling.F90), with a cost going as the product of +the frequency of units that resolve to the same scaling factor. + + A separate set of scaling factors could also be used for different chemical +tracer concentrations, for example. In this case, the tools in +MOM_unique_scales.F90 could still be used, but there would need to be a separate +equivalent of the unit_scaling_type with variables that are appropriate to the +units of the tracers. + +*/ From 64f432fcb5b962a890330788fc11c80572296fd1 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 15 Feb 2022 22:10:51 -0500 Subject: [PATCH 4/9] Diabatic driver: energetic_PBL -> ePBL, flag check Pointers to the diabatic driver's energetic PBL field are now only associated when `use_energetic_PBL` is true. The `energetic_PBL` field was also renamed to `ePBL` to avoid potential conflict with the `energetic_PBL` subroutine. Thanks to Alper Altuntas for detecting this issue and the proposed fix. Co-Authored-By: Alper Altuntas --- .../vertical/MOM_diabatic_driver.F90 | 21 +++++++++---------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 5eaca3c275..7b180f1d65 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -234,7 +234,7 @@ module MOM_diabatic_driver type(oda_incupd_CS), pointer :: oda_incupd_CSp => NULL() !< Control structure for a child module type(bulkmixedlayer_CS) :: bulkmixedlayer !< Bulk mixed layer control struct type(CVMix_conv_CS) :: CVMix_conv !< CVMix convection control struct - type(energetic_PBL_CS) :: energetic_PBL !< Energetic PBL control struct + type(energetic_PBL_CS) :: ePBL !< Energetic PBL control struct type(entrain_diffusive_CS) :: entrain_diffusive !< Diffusive entrainment control struct type(geothermal_CS) :: geothermal !< Geothermal control struct type(int_tide_CS) :: int_tide !< Internal tide control struct @@ -838,15 +838,15 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, & - CS%energetic_PBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) + CS%ePBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then - call energetic_PBL_get_MLD(CS%energetic_PBL, Hml(:,:), G, US) + call energetic_PBL_get_MLD(CS%ePBL, Hml(:,:), G, US) call pass_var(Hml, G%domain, halo=1) ! If visc%MLD exists, copy ePBL's MLD into it if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) elseif (associated(visc%MLD)) then - call energetic_PBL_get_MLD(CS%energetic_PBL, visc%MLD, G, US) + call energetic_PBL_get_MLD(CS%ePBL, visc%MLD, G, US) call pass_var(visc%MLD, G%domain, halo=1) endif @@ -1375,15 +1375,15 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, & - CS%energetic_PBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) + CS%ePBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then - call energetic_PBL_get_MLD(CS%energetic_PBL, Hml(:,:), G, US) + call energetic_PBL_get_MLD(CS%ePBL, Hml(:,:), G, US) call pass_var(Hml, G%domain, halo=1) ! If visc%MLD exists, copy ePBL's MLD into it if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) elseif (associated(visc%MLD)) then - call energetic_PBL_get_MLD(CS%energetic_PBL, visc%MLD, G, US) + call energetic_PBL_get_MLD(CS%ePBL, visc%MLD, G, US) call pass_var(visc%MLD, G%domain, halo=1) endif @@ -2589,14 +2589,13 @@ subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, evap_CFL_limit, if (present(opacity_CSp)) opacity_CSp => CS%opacity if (present(optics_CSp)) optics_CSp => CS%optics if (present(KPP_CSp)) KPP_CSp => CS%KPP_CSp - if (present(energetic_PBL_CSp)) energetic_PBL_CSp => CS%energetic_PBL + if (present(energetic_PBL_CSp) .and. CS%use_energetic_PBL) energetic_PBL_CSp => CS%ePBL if (present(diabatic_aux_CSp)) diabatic_aux_CSp => CS%diabatic_aux_CSp ! Constants within diabatic_CS if (present(evap_CFL_limit)) evap_CFL_limit = CS%evap_CFL_limit if (present(minimum_forcing_depth)) minimum_forcing_depth = CS%minimum_forcing_depth if (present(diabatic_halo)) diabatic_halo = CS%halo_TS_diff - end subroutine extract_diabatic_member !> Routine called for adiabatic physics @@ -3487,7 +3486,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di if (CS%use_bulkmixedlayer) & call bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS%bulkmixedlayer) if (CS%use_energetic_PBL) & - call energetic_PBL_init(Time, G, GV, US, param_file, diag, CS%energetic_PBL) + call energetic_PBL_init(Time, G, GV, US, param_file, diag, CS%ePBL) call regularize_layers_init(Time, G, GV, param_file, diag, CS%regularize_layers) @@ -3522,7 +3521,7 @@ subroutine diabatic_driver_end(CS) call diapyc_energy_req_end(CS%diapyc_en_rec_CSp) if (CS%use_energetic_PBL) & - call energetic_PBL_end(CS%energetic_PBL) + call energetic_PBL_end(CS%ePBL) call diabatic_aux_end(CS%diabatic_aux_CSp) From fc5253f73c34b2fafff0fe446d6e4cd75e317752 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 24 Jan 2022 14:07:37 -0500 Subject: [PATCH 5/9] (*)Correct memory declarations in MOM_regridding Corrected memory declarations in MOM_regridding.F90 in cases where the vertical size of the input and output grids are not the same. Although this is not know to have caused any particular problems, these inconsistencies could lead to segmentation faults in cases where the target grid (e.g., diagnostic output) is larger than the input grid (e.g., the model's native grid). In some cases, certain grid generation options were only written to work with the same size of input and output grids, and error handling has been added to these cases to gracefully bring down the model if they are used with different grid sizes. All answers are bitwise identical in the MOM6-examples test suite, but it is conceivable that this could correct subtle (memory-related) issues in some configurations. --- src/ALE/MOM_regridding.F90 | 108 +++++++++++++++++++++++-------------- 1 file changed, 67 insertions(+), 41 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 94d7852851..dafd165245 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -3,7 +3,7 @@ module MOM_regridding ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_error_handler, only : MOM_error, FATAL, WARNING +use MOM_error_handler, only : MOM_error, FATAL, WARNING, assert use MOM_file_parser, only : param_file_type, get_param, log_param use MOM_io, only : file_exists, field_exists, field_size, MOM_read_data use MOM_io, only : verify_variable_units, slasher @@ -776,11 +776,11 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_ type(regridding_CS), intent(in) :: CS !< Regridding control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after !! the last time step type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamical variables (T, S, ...) - real, dimension(SZI_(G),SZJ_(G), CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target coordinate - real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in position of each interface + real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target coordinate + real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in position of each interface real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: frac_shelf_h !< Fractional ice shelf coverage logical, optional, intent(in ) :: conv_adjust !< If true, do convective adjustment ! Local variables @@ -827,7 +827,7 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_ end select ! type of grid #ifdef __DO_SAFETY_CHECKS__ - call check_remapping_grid(G, GV, h, dzInterface,'in regridding_main') + if (CS%nk == GV%ke) call check_remapping_grid(G, GV, h, dzInterface,'in regridding_main') #endif end subroutine regridding_main @@ -1086,13 +1086,14 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth + real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth !! [H ~> m or kg m-2]. real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: frac_shelf_h !< Fractional !! ice shelf coverage [nondim]. ! Local variables real :: nominalDepth, minThickness, totalThickness, dh ! Depths and thicknesses [H ~> m or kg m-2] - real, dimension(SZK_(GV)+1) :: zOld, zNew ! Coordinate interface heights [H ~> m or kg m-2] + real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2] + real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2] integer :: i, j, k, nz logical :: ice_shelf @@ -1144,17 +1145,23 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) call filtered_grid_motion( CS, nz, zOld, zNew, dzInterface(i,j,:) ) #ifdef __DO_SAFETY_CHECKS__ - dh=max(nominalDepth,totalThickness) - if (abs(zNew(1)-zOld(1))>(nz-1)*0.5*epsilon(dh)*dh) then + dh = max(nominalDepth,totalThickness) + if (abs(zNew(1)-zOld(1)) > (nz-1)*0.5*epsilon(dh)*dh) then write(0,*) 'min_thickness=',CS%min_thickness write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness - write(0,*) 'dzInterface(1) = ',dzInterface(i,j,1),epsilon(dh),nz - do k=1,nz+1 + write(0,*) 'dzInterface(1) = ', dzInterface(i,j,1), epsilon(dh), nz, CS%nk + do k=1,min(nz,CS%nk)+1 write(0,*) k,zOld(k),zNew(k) enddo - do k=1,nz + do k=min(nz,CS%nk)+2,CS%nk+1 + write(0,*) k,zOld(nz+1),zNew(k) + enddo + do k=1,min(nz,CS%nk) write(0,*) k,h(i,j,k),zNew(k)-zNew(k+1),CS%coordinateResolution(k) enddo + do k=min(nz,CS%nk)+1,CS%nk + write(0,*) k, 0.0, zNew(k)-zNew(k+1), CS%coordinateResolution(k) + enddo call MOM_error( FATAL, & 'MOM_regridding, build_zstar_grid(): top surface has moved!!!' ) endif @@ -1183,14 +1190,15 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface ) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth + real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth !! [H ~> m or kg m-2] ! Local variables integer :: i, j, k integer :: nz real :: nominalDepth, totalThickness, dh - real, dimension(SZK_(GV)+1) :: zOld, zNew + real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2] + real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2] nz = GV%ke @@ -1227,12 +1235,18 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface ) write(0,*) 'min_thickness=',CS%min_thickness write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness write(0,*) 'dzInterface(1) = ',dzInterface(i,j,1),epsilon(dh),nz,CS%nk - do k=1,nz+1 + do k=1,min(nz,CS%nk)+1 write(0,*) k,zOld(k),zNew(k) enddo - do k=1,CS%nk + do k=min(nz,CS%nk)+2,CS%nk+1 + write(0,*) k,zOld(nz+1),zNew(k) + enddo + do k=1,min(nz,CS%nk) write(0,*) k,h(i,j,k),zNew(k)-zNew(k+1),totalThickness*CS%coordinateResolution(k),CS%coordinateResolution(k) enddo + do k=min(nz,CS%nk)+1,CS%nk + write(0,*) k,0.0,zNew(k)-zNew(k+1),totalThickness*CS%coordinateResolution(k),CS%coordinateResolution(k) + enddo call MOM_error( FATAL, & 'MOM_regridding, build_sigma_grid: top surface has moved!!!' ) endif @@ -1266,22 +1280,23 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel !------------------------------------------------------------------------------ ! Arguments + type(regridding_CS), intent(in) :: CS !< Regridding control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure - real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth + real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth !! [H ~> m or kg m-2] type(remapping_CS), intent(in) :: remapCS !< The remapping control structure - type(regridding_CS), intent(in) :: CS !< Regridding control structure real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice !! shelf coverage [nondim] ! Local variables - integer :: nz + integer :: nz ! The number of layers in the input grid integer :: i, j, k real :: nominalDepth ! Depth of the bottom of the ocean, positive downward [H ~> m or kg m-2] - real, dimension(SZK_(GV)+1) :: zOld, zNew ! Old and new interface heights [H ~> m or kg m-2] + real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2] + real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2] real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] real :: totalThickness ! Total thicknesses [H ~> m or kg m-2] #ifdef __DO_SAFETY_CHECKS__ @@ -1355,7 +1370,7 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel call filtered_grid_motion( CS, nz, zOld, zNew, dzInterface(i,j,:) ) #ifdef __DO_SAFETY_CHECKS__ - do k = 2,nz + do k=2,CS%nk if (zNew(k) > zOld(1)) then write(0,*) 'zOld=',zOld write(0,*) 'zNew=',zNew @@ -1380,12 +1395,18 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel write(0,*) 'min_thickness=',CS%min_thickness write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness write(0,*) 'zNew(1)-zOld(1) = ',zNew(1)-zOld(1),epsilon(dh),nz - do k=1,nz+1 + do k=1,min(nz,CS%nk)+1 write(0,*) k,zOld(k),zNew(k) enddo + do k=min(nz,CS%nk)+2,CS%nk+1 + write(0,*) k,zOld(nz+1),zNew(k) + enddo do k=1,nz write(0,*) k,h(i,j,k),zNew(k)-zNew(k+1) enddo + do k=min(nz,CS%nk)+1,CS%nk + write(0,*) k, 0.0, zNew(k)-zNew(k+1), CS%coordinateResolution(k) + enddo call MOM_error( FATAL, & 'MOM_regridding, build_rho_grid: top surface has moved!!!' ) endif @@ -1416,10 +1437,10 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS, frac_she !! coverage [nondim] ! Local variables - real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2] + 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(SZK_(GV)+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2] - real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [R L2 T-2 ~> Pa] + real, dimension(CS%nk+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2] integer :: i, j, k, nki real :: depth, nominalDepth real :: h_neglect, h_neglect_edge @@ -1496,10 +1517,10 @@ subroutine build_grid_adaptive(G, GV, US, h, tv, dzInterface, remapCS, CS) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth + type(regridding_CS), intent(in) :: CS !< Regridding control structure + real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth !! [H ~> m or kg m-2] type(remapping_CS), intent(in) :: remapCS !< The remapping control structure - type(regridding_CS), intent(in) :: CS !< Regridding control structure ! local variables integer :: i, j, k, nz ! indices and dimension lengths @@ -1512,6 +1533,9 @@ subroutine build_grid_adaptive(G, GV, US, h, tv, dzInterface, remapCS, CS) nz = GV%ke + call assert((GV%ke == CS%nk), "build_grid_adaptive is only written to work "//& + "with the same number of input and target layers.") + ! position surface at z = 0. zInt(:,:,1) = 0. @@ -1562,13 +1586,13 @@ subroutine build_grid_SLight(G, GV, US, h, tv, dzInterface, CS) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Existing model thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< Changes in interface position type(regridding_CS), intent(in) :: CS !< Regridding control structure + real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< Changes in interface position - real, dimension(SZK_(GV)+1) :: z_col ! Interface positions relative to the surface [H ~> m or kg m-2] - real, dimension(SZK_(GV)+1) :: z_col_new ! Interface positions relative to the surface [H ~> m or kg m-2] - real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2] - real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [R L2 T-2 ~> Pa] + real, 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] @@ -1585,8 +1609,10 @@ subroutine build_grid_SLight(G, GV, US, h, tv, dzInterface, CS) nz = GV%ke - if (.not.CS%target_density_set) call MOM_error(FATAL, "build_grid_SLight : "//& - "Target densities must be set before build_grid_SLight is called.") + 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 @@ -1691,13 +1717,13 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) !------------------------------------------------------------------------------ ! Arguments - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(in) :: h !< Original layer thicknesses [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface - !! depth [H ~> m or kg m-2] - real, intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2] - type(regridding_CS), intent(in) :: CS !< Regridding control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Original layer thicknesses [H ~> m or kg m-2] + type(regridding_CS), intent(in) :: CS !< Regridding control structure + real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface + !! depth [H ~> m or kg m-2] + real, intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2] ! Local variables integer :: i, j, k From c166358b21a5edb970571b31202c7185990913ea Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Tue, 1 Feb 2022 09:58:15 -0500 Subject: [PATCH 6/9] Add optional argument to FMS2 version of get_field_size - default behavior returns field size using fms1 format. --- config_src/infra/FMS2/MOM_io_infra.F90 | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index 62a43ab99b..f8999fa7b8 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/config_src/infra/FMS2/MOM_io_infra.F90 @@ -598,7 +598,7 @@ function field_exists(filename, field_name, domain, no_domain, MOM_domain) end function field_exists !> Given filename and fieldname, this subroutine returns the size of the field in the file -subroutine get_field_size(filename, fieldname, sizes, field_found, no_domain) +subroutine get_field_size(filename, fieldname, sizes, field_found, no_domain, fms1_format) character(len=*), intent(in) :: filename !< The name of the file to read character(len=*), intent(in) :: fieldname !< The name of the variable whose sizes are returned integer, dimension(:), intent(inout) :: sizes !< The sizes of the variable in each dimension @@ -607,16 +607,21 @@ subroutine get_field_size(filename, fieldname, sizes, field_found, no_domain) !! is a fatal error if the field is not found. logical, optional, intent(in) :: no_domain !< If present and true, do not check for file !! names with an appended tile number - + logical, optional, intent(in) :: fms1_format !< If true (default) , then for (Nx,Ny,Nt) data + !! return sizes=(Nx,Ny,1,Nt) if sizes if 4-dimensional ! Local variables type(FmsNetcdfFile_t) :: fileobj_read ! A handle to a non-domain-decomposed file for obtaining information ! about the exiting time axis entries in append mode. logical :: success ! If true, the file was opened successfully logical :: field_exists ! True if filename exists and field_name is in filename + logical :: fms1_fmt integer :: i, ndims + if (FMS2_reads) then field_exists = .false. + fms1_fmt=.true. + if (present(fms1_format)) fms1_fmt=fms1_format if (file_exists(filename)) then success = fms2_open_file(fileObj_read, trim(filename), "read") if (success) then @@ -627,6 +632,12 @@ subroutine get_field_size(filename, fieldname, sizes, field_found, no_domain) "get_field_size called with too few sizes for "//trim(fieldname)//" in "//trim(filename)) call get_variable_size(fileobj_read, fieldname, sizes(1:ndims)) do i=ndims+1,size(sizes) ; sizes(i) = 0 ; enddo + ! This preserves previous behavior when reading time-varying data without + ! a vertical extent. + if (fms1_fmt .and. size(sizes)==ndims+1) then + sizes(ndims+1)=sizes(ndims) + sizes(ndims)=1 + endif endif endif endif From e8416091589b98ae790c9c3a63e5c7980cc56c73 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Wed, 2 Feb 2022 17:10:12 -0500 Subject: [PATCH 7/9] remove unnecessary optional flag --- config_src/infra/FMS2/MOM_io_infra.F90 | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index f8999fa7b8..16b7e45bc6 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/config_src/infra/FMS2/MOM_io_infra.F90 @@ -598,7 +598,7 @@ function field_exists(filename, field_name, domain, no_domain, MOM_domain) end function field_exists !> Given filename and fieldname, this subroutine returns the size of the field in the file -subroutine get_field_size(filename, fieldname, sizes, field_found, no_domain, fms1_format) +subroutine get_field_size(filename, fieldname, sizes, field_found, no_domain) character(len=*), intent(in) :: filename !< The name of the file to read character(len=*), intent(in) :: fieldname !< The name of the variable whose sizes are returned integer, dimension(:), intent(inout) :: sizes !< The sizes of the variable in each dimension @@ -607,21 +607,16 @@ subroutine get_field_size(filename, fieldname, sizes, field_found, no_domain, fm !! is a fatal error if the field is not found. logical, optional, intent(in) :: no_domain !< If present and true, do not check for file !! names with an appended tile number - logical, optional, intent(in) :: fms1_format !< If true (default) , then for (Nx,Ny,Nt) data - !! return sizes=(Nx,Ny,1,Nt) if sizes if 4-dimensional ! Local variables type(FmsNetcdfFile_t) :: fileobj_read ! A handle to a non-domain-decomposed file for obtaining information ! about the exiting time axis entries in append mode. logical :: success ! If true, the file was opened successfully logical :: field_exists ! True if filename exists and field_name is in filename - logical :: fms1_fmt integer :: i, ndims if (FMS2_reads) then field_exists = .false. - fms1_fmt=.true. - if (present(fms1_format)) fms1_fmt=fms1_format if (file_exists(filename)) then success = fms2_open_file(fileObj_read, trim(filename), "read") if (success) then @@ -634,7 +629,7 @@ subroutine get_field_size(filename, fieldname, sizes, field_found, no_domain, fm do i=ndims+1,size(sizes) ; sizes(i) = 0 ; enddo ! This preserves previous behavior when reading time-varying data without ! a vertical extent. - if (fms1_fmt .and. size(sizes)==ndims+1) then + if (size(sizes)==ndims+1) then sizes(ndims+1)=sizes(ndims) sizes(ndims)=1 endif From 32e1ecf45f8d5064b3f7e986b1e9af226958a919 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Sun, 20 Feb 2022 08:21:17 -0700 Subject: [PATCH 8/9] Fixes issues with the GME code and get_param calls for Leith options (#65) * Eliminate GET_ALL_PARAMS in hor_visc_init Added do_not_log arguments to get_param calls in MOM_hor_visc.F90 that are only used conditionally, and eliminated the unlogged GET_ALL_PARAMS runtime parameter and get_all variable in hor_visc_init(). By design, all logging of parameters after this commit is identical to before, even for variables that are inactive and therefore should not be logged. In several places, there were some problems, mostly with the GME code, that have been noted in comments marked with '###'. Also cleaned up the code alignment and eliminated unneeded temporary variables in a few places in hor_visc(). All solutions are bitwise identical, and no output is changed. * Move call to get get_KH outside k-loop The call to get the 3-d GME diffusivity arrays and the subsequent blocking halo update was moved outside of the k-loop. * Increase loop range for calculation of GME fluxes * Makes GME filter rotationally invariant * Makes the GME filter rotationally invariant * Adds a runtime param to allow the user to control how many smoothing passes should be done. * Rearranges the get_param calls related to Leith The get_param calls for Leith were not in the correct location. This commit fixes that. * Adding halo updates * Fixes do loops indices and adds diagnostics * Adds option to save barotropic tension and strain; * Fixes many i and j loops indices associated with GME to avoid halo problems and unnecessary halo updates. With these changes, the model is now conserving mass and tracers when USE_GME = True. * Fixes issues related to the merging with dev/gfdl * Fixes calculation of FrictWork_GME The calculation now mimics the calculation of FrictWork and it includes the energy diffusion term. * Removes dependency of FrictWork_GME calculation on MEKE * Adding sh_xy_sq and sh_xx_sq in the OMP directives Co-authored-by: Robert Hallberg --- .../lateral/MOM_hor_visc.F90 | 254 ++++++++++-------- 1 file changed, 145 insertions(+), 109 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 0249f79c2d..78f48975e5 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -2,14 +2,14 @@ module MOM_hor_visc ! This file is part of MOM6. See LICENSE.md for the license. - -use MOM_checksums, only : hchksum, Bchksum +use MOM_checksums, only : hchksum, Bchksum, uvchksum use MOM_coms, only : min_across_PEs use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : post_product_u, post_product_sum_u use MOM_diag_mediator, only : post_product_v, post_product_sum_v use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_domains, only : pass_var, CORNER, pass_vector, AGRID, BGRID_NE +use MOM_domains, only : To_All, Scalar_Pair use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type @@ -185,6 +185,7 @@ module MOM_hor_visc ! 3D diagnostics hf_diffu(diffv) are commented because there is no clarity on proper remapping grid option. ! The code is retained for degugging purposes in the future. + integer :: num_smooth_gme !< number of smoothing passes for the GME fluxes. !>@{ !! Diagnostic id integer :: id_grid_Re_Ah = -1, id_grid_Re_Kh = -1 @@ -197,6 +198,8 @@ module MOM_hor_visc integer :: id_Ah_h = -1, id_Ah_q = -1 integer :: id_Kh_h = -1, id_Kh_q = -1 integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1 + integer :: id_dudx_bt = -1, id_dvdy_bt = -1 + integer :: id_dudy_bt = -1, id_dvdx_bt = -1 integer :: id_vort_xy_q = -1, id_div_xx_h = -1 integer :: id_sh_xy_q = -1, id_sh_xx_h = -1 integer :: id_FrictWork = -1, id_FrictWorkIntz = -1 @@ -452,12 +455,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, (CS%bound_Kh .and. .not.CS%better_bound_Kh) if (CS%use_GME) then + do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 boundary_mask_h(i,j) = (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1)) enddo ; enddo - do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 - boundary_mask_q(I,J) = (G%mask2dCv(i,J) * G%mask2dCv(i+1,J) * G%mask2dCu(I,j) * G%mask2dCu(I,j-1)) + do J=Jsq-2,Jeq+1 ; do I=Isq-2,Ieq+1 + boundary_mask_q(I,J) = (G%mask2dCv(i,J) * G%mask2dCv(i+1,J) * G%mask2dCu(I,j) * G%mask2dCu(I,j+1)) enddo ; enddo ! initialize diag. array with zeros @@ -468,81 +472,92 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Get barotropic velocities and their gradients call barotropic_get_tav(BT, ubtav, vbtav, G, US) + call pass_vector(ubtav, vbtav, G%Domain) - do j=js-1,je+2 ; do i=is-1,ie+2 + ! Calculate the barotropic horizontal tension + do j=Jsq-2,Jeq+2 ; do i=Isq-2,Ieq+2 dudx_bt(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * ubtav(I,j) - & G%IdyCu(I-1,j) * ubtav(I-1,j)) dvdy_bt(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * vbtav(i,J) - & G%IdxCv(i,J-1) * vbtav(i,J-1)) - enddo ; enddo - - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j) enddo ; enddo ! Components for the barotropic shearing strain - do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 dvdx_bt(I,J) = CS%DY_dxBu(I,J)*(vbtav(i+1,J)*G%IdyCv(i+1,J) & - vbtav(i,J)*G%IdyCv(i,J)) dudy_bt(I,J) = CS%DX_dyBu(I,J)*(ubtav(I,j+1)*G%IdxCu(I,j+1) & - ubtav(I,j)*G%IdxCu(I,j)) enddo ; enddo - call pass_vector(dudx_bt, dvdy_bt, G%Domain, stagger=BGRID_NE) - call pass_vector(dvdx_bt, dudy_bt, G%Domain, stagger=AGRID) + ! post barotropic tension and strain + if (CS%id_dudx_bt > 0) call post_data(CS%id_dudx_bt, dudx_bt, CS%diag) + if (CS%id_dvdy_bt > 0) call post_data(CS%id_dvdy_bt, dvdy_bt, CS%diag) + if (CS%id_dudy_bt > 0) call post_data(CS%id_dudy_bt, dudy_bt, CS%diag) + if (CS%id_dvdx_bt > 0) call post_data(CS%id_dvdx_bt, dvdx_bt, CS%diag) if (CS%no_slip) then - do J=js-1,Jeq ; do I=is-1,Ieq + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 sh_xy_bt(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_bt(I,J) + dudy_bt(I,J) ) enddo ; enddo else - do J=js-1,Jeq ; do I=is-1,Ieq + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 sh_xy_bt(I,J) = G%mask2dBu(I,J) * ( dvdx_bt(I,J) + dudy_bt(I,J) ) enddo ; enddo endif do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 - grad_vel_mag_bt_h(i,j) = boundary_mask_h(i,j) * (dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + & + grad_vel_mag_bt_h(i,j) = G%mask2dT(I,J) * boundary_mask_h(i,j) * (dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + & (0.25*((dvdx_bt(I,J)+dvdx_bt(I-1,J-1))+(dvdx_bt(I,J-1)+dvdx_bt(I-1,J))))**2 + & (0.25*((dudy_bt(I,J)+dudy_bt(I-1,J-1))+(dudy_bt(I,J-1)+dudy_bt(I-1,J))))**2) enddo ; enddo - do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 - grad_vel_mag_bt_q(I,J) = boundary_mask_q(I,J) * (dvdx_bt(I,J)**2 + dudy_bt(I,J)**2 + & + do J=Jsq-2,Jeq+1 ; do I=Isq-2,Ieq+1 + grad_vel_mag_bt_q(I,J) = G%mask2dBu(I,J) * boundary_mask_q(I,J) * (dvdx_bt(I,J)**2 + dudy_bt(I,J)**2 + & (0.25*((dudx_bt(i,j)+dudx_bt(i+1,j+1))+(dudx_bt(i,j+1)+dudx_bt(i+1,j))))**2 + & (0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1))+(dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2) enddo ; enddo - do j=js-1,je+1 ; do i=is-1,ie+1 + call pass_var(h, G%domain, halo=2) + + do j=js-2,je+2 ; do i=is-2,ie+2 htot(i,j) = 0.0 enddo ; enddo - do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 + do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2 htot(i,j) = htot(i,j) + GV%H_to_Z*h(i,j,k) enddo ; enddo ; enddo I_GME_h0 = 1.0 / CS%GME_h0 - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 if (grad_vel_mag_bt_h(i,j)>0) then - GME_effic_h(i,j) = CS%GME_efficiency * boundary_mask_h(i,j) * & + GME_effic_h(i,j) = CS%GME_efficiency * G%mask2dT(I,J) * & (MIN(htot(i,j) * I_GME_h0, 1.0)**2) else GME_effic_h(i,j) = 0.0 endif enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq + do J=Jsq-2,Jeq+1 ; do I=Isq-2,Ieq+1 if (grad_vel_mag_bt_q(I,J)>0) then h_arith_q = 0.25 * ((htot(i,j) + htot(i+1,j+1)) + (htot(i+1,j) + htot(i,j+1))) I_hq = 1.0 / h_arith_q h_harm_q = 0.25 * h_arith_q * ((htot(i,j)*I_hq + htot(i+1,j+1)*I_hq) + & (htot(i+1,j)*I_hq + htot(i,j+1)*I_hq)) - GME_effic_q(I,J) = CS%GME_efficiency * boundary_mask_q(I,J) * (MIN(h_harm_q * I_GME_h0, 1.0)**2) + GME_effic_q(I,J) = CS%GME_efficiency * G%mask2dBu(I,J) * (MIN(h_harm_q * I_GME_h0, 1.0)**2) else GME_effic_q(I,J) = 0.0 endif enddo ; enddo + call thickness_diffuse_get_KH(TD, KH_u_GME, KH_v_GME, G, GV) + + call pass_vector(KH_u_GME, KH_v_GME, G%domain, To_All+Scalar_Pair) + + if (CS%debug) & + call uvchksum("GME KH[u,v]_GME", KH_u_GME, KH_v_GME, G%HI, haloshift=2, scale=US%L_to_m**2*US%s_to_T) + endif ! use_GME !$OMP parallel do default(none) & @@ -555,7 +570,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP h_neglect, h_neglect3, FWfrac, inv_PI3, inv_PI6, H0_GME, & !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, & !$OMP div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, & - !$OMP TD, KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt & + !$OMP KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt & !$OMP ) & !$OMP private( & !$OMP i, j, k, n, & @@ -565,10 +580,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP vort_xy, vort_xy_dx, vort_xy_dy, div_xx, div_xx_dx, div_xx_dy, & !$OMP grad_div_mag_h, grad_div_mag_q, grad_vort_mag_h, grad_vort_mag_q, & !$OMP grad_vort, grad_vort_qg, grad_vort_mag_h_2d, grad_vort_mag_q_2d, & - !$OMP grad_vel_mag_h, grad_vel_mag_q, & + !$OMP grad_vel_mag_h, grad_vel_mag_q, sh_xx_sq, sh_xy_sq, & !$OMP grad_vel_mag_bt_h, grad_vel_mag_bt_q, grad_d2vel_mag_h, & !$OMP meke_res_fn, Shear_mag, Shear_mag_bc, vert_vort_mag, h_min, hrat_min, visc_bound_rem, & - !$OMP sh_xx_sq, sh_xy_sq, grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, & + !$OMP grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, & !$OMP Kh, Ah, AhSm, AhLth, local_strain, Sh_F_pow, & !$OMP dDel2vdx, dDel2udy, DY_dxCv, DX_dyCu, Del2vort_q, Del2vort_h, KE, & !$OMP h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff & @@ -1426,52 +1441,40 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif if (CS%use_GME) then - !### This call to get the 3-d GME diffusivity arrays and the subsequent blocking halo update - ! should occur outside of the k-loop, and perhaps the halo update should occur outside of - ! this routine altogether! - call thickness_diffuse_get_KH(TD, KH_u_GME, KH_v_GME, G, GV) - call pass_vector(KH_u_GME, KH_v_GME, G%Domain) - - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 GME_coeff = GME_effic_h(i,j) * 0.25 * & ((KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)) + (KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k))) GME_coeff = MIN(GME_coeff, CS%GME_limiter) if ((CS%id_GME_coeff_h>0) .or. find_FrictWork) GME_coeff_h(i,j,k) = GME_coeff str_xx_GME(i,j) = GME_coeff * sh_xx_bt(i,j) - enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq + do J=Jsq-2,Jeq+1 ; do I=Isq-2,Ieq+1 GME_coeff = GME_effic_q(I,J) * 0.25 * & ((KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)) + (KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k))) GME_coeff = MIN(GME_coeff, CS%GME_limiter) if (CS%id_GME_coeff_q>0) GME_coeff_q(I,J,k) = GME_coeff str_xy_GME(I,J) = GME_coeff * sh_xy_bt(I,J) - enddo ; enddo ! Applying GME diagonal term. This is linear and the arguments can be rescaled. - call smooth_GME(G, GME_flux_h=str_xx_GME) - call smooth_GME(G, GME_flux_q=str_xy_GME) + call smooth_GME(CS, G, GME_flux_h=str_xx_GME) + call smooth_GME(CS, G, GME_flux_q=str_xy_GME) do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 str_xx(i,j) = (str_xx(i,j) + str_xx_GME(i,j)) * (h(i,j,k) * CS%reduction_xx(i,j)) enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq - ! GME is applied below - if (CS%no_slip) then + ! GME is applied below + if (CS%no_slip) then + do J=js-1,Jeq ; do I=is-1,Ieq str_xy(I,J) = (str_xy(I,J) + str_xy_GME(I,J)) * (hq(I,J) * CS%reduction_xy(I,J)) - else + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq str_xy(I,J) = (str_xy(I,J) + str_xy_GME(I,J)) * (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J)) - endif - enddo ; enddo - - if (allocated(MEKE%GME_snk)) then - do j=js,je ; do i=is,ie - FrictWork_GME(i,j,k) = GME_coeff_h(i,j,k) * h(i,j,k) * GV%H_to_kg_m2 * grad_vel_mag_bt_h(i,j) enddo ; enddo endif @@ -1555,6 +1558,27 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, +(v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1) )) ) ) enddo ; enddo ; endif + if (CS%use_GME) then + do j=js,je ; do i=is,ie + ! Diagnose str_xx_GME*d_x u - str_yy_GME*d_y v + str_xy_GME*(d_y u + d_x v) + ! This is the old formulation that includes energy diffusion + FrictWork_GME(i,j,k) = GV%H_to_RZ * ( & + (str_xx_GME(i,j)*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & + -str_xx_GME(i,j)*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & + +0.25*((str_xy_GME(I,J)*( & + (u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J) & + +(v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J) ) & + +str_xy_GME(I-1,J-1)*( & + (u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1) & + +(v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1) )) & + +(str_xy_GME(I-1,J)*( & + (u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J) & + +(v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J) ) & + +str_xy_GME(I,J-1)*( & + (u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & + +(v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1) )) ) ) + enddo ; enddo ; endif + ! Make a similar calculation as for FrictWork above but accumulating into ! the vertically integrated MEKE source term, and adjusting for any ! energy loss seen as a reduction in the (biharmonic) frictional source term. @@ -1818,11 +1842,10 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "If true, use a Leith nonlinear eddy viscosity.", & default=.false., do_not_log=.not.CS%Laplacian) if (.not.CS%Laplacian) CS%Leith_Kh = .false. - ! This call duplicates one that occurs 26 lines later, and is probably unneccessary. - call get_param(param_file, mdl, "MODIFIED_LEITH", CS%Modified_Leith, & - "If true, add a term to Leith viscosity which is "//& - "proportional to the gradient of divergence.", & - default=.false., do_not_log=.not.CS%Laplacian) !### (.not.CS%Leith_Kh)? + call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, & + "The nondimensional Laplacian Leith constant, "//& + "often set to 1.0", units="nondim", default=0.0, & + fail_if_missing=CS%Leith_Kh, do_not_log=.not.CS%Leith_Kh) call get_param(param_file, mdl, "USE_MEKE", use_MEKE, & default=.false., do_not_log=.true.) call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", CS%res_scale_MEKE, & @@ -1831,27 +1854,6 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) do_not_log=.not.(CS%Laplacian.and.use_MEKE)) if (.not.(CS%Laplacian.and.use_MEKE)) CS%res_scale_MEKE = .false. - call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, & - "The nondimensional Laplacian Leith constant, "//& - "often set to 1.0", units="nondim", default=0.0, & - fail_if_missing=CS%Leith_Kh, do_not_log=.not.CS%Leith_Kh) - call get_param(param_file, mdl, "USE_QG_LEITH_VISC", CS%use_QG_Leith_visc, & - "If true, use QG Leith nonlinear eddy viscosity.", & - default=.false., do_not_log=.not.CS%Leith_Kh) - if (CS%use_QG_Leith_visc .and. .not. CS%Leith_Kh) call MOM_error(FATAL, & - "MOM_hor_visc.F90, hor_visc_init:"//& - "LEITH_KH must be True when USE_QG_LEITH_VISC=True.") - - !### The following two get_param_calls need to occur after Leith_Ah is read, but for now it replicates prior code. - CS%Leith_Ah = .false. - call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_Leith, & - "If true, include the beta term in the Leith nonlinear eddy viscosity.", & - default=CS%Leith_Kh, do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) ) - call get_param(param_file, mdl, "MODIFIED_LEITH", CS%modified_Leith, & - "If true, add a term to Leith viscosity which is "//& - "proportional to the gradient of divergence.", & - default=.false., do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) ) - call get_param(param_file, mdl, "BOUND_KH", CS%bound_Kh, & "If true, the Laplacian coefficient is locally limited "//& "to be stable.", default=.true., do_not_log=.not.CS%Laplacian) @@ -1944,6 +1946,21 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "typically 0.015 - 0.06.", units="nondim", default=0.0, & fail_if_missing=CS%Smagorinsky_Ah, do_not_log=.not.CS%Smagorinsky_Ah) + call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_Leith, & + "If true, include the beta term in the Leith nonlinear eddy viscosity.", & + default=CS%Leith_Kh, do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) ) + call get_param(param_file, mdl, "MODIFIED_LEITH", CS%modified_Leith, & + "If true, add a term to Leith viscosity which is "//& + "proportional to the gradient of divergence.", & + default=.false., do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) ) + call get_param(param_file, mdl, "USE_QG_LEITH_VISC", CS%use_QG_Leith_visc, & + "If true, use QG Leith nonlinear eddy viscosity.", & + default=.false., do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) ) + if (CS%use_QG_Leith_visc .and. .not. (CS%Leith_Kh .or. CS%Leith_Ah) ) then + call MOM_error(FATAL, "MOM_hor_visc.F90, hor_visc_init:"//& + "LEITH_KH or LEITH_AH must be True when USE_QG_LEITH_VISC=True.") + endif + call get_param(param_file, mdl, "BOUND_CORIOLIS", bound_Cor_def, default=.false.) call get_param(param_file, mdl, "BOUND_CORIOLIS_BIHARM", CS%bound_Coriolis, & "If true use a viscosity that increases with the square "//& @@ -1999,17 +2016,22 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "Use the split time stepping if true.", default=.true., do_not_log=.true.) if (CS%use_GME .and. .not.split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// & "cannot be used with SPLIT=False.") - call get_param(param_file, mdl, "GME_H0", CS%GME_h0, & - "The strength of GME tapers quadratically to zero when the bathymetric "//& - "depth is shallower than GME_H0.", & - units="m", scale=US%m_to_Z, default=1000.0, do_not_log=.not.CS%use_GME) - call get_param(param_file, mdl, "GME_EFFICIENCY", CS%GME_efficiency, & - "The nondimensional prefactor multiplying the GME coefficient.", & - units="nondim", default=1.0, do_not_log=.not.CS%use_GME) - call get_param(param_file, mdl, "GME_LIMITER", CS%GME_limiter, & - "The absolute maximum value the GME coefficient is allowed to take.", & - units="m2 s-1", scale=US%m_to_L**2*US%T_to_s, default=1.0e7, & - do_not_log=.not.CS%use_GME) + + if (CS%use_GME) then + call get_param(param_file, mdl, "GME_NUM_SMOOTHINGS", CS%num_smooth_gme, & + "Number of smoothing passes for the GME fluxes.", & + units="nondim", default=1) + call get_param(param_file, mdl, "GME_H0", CS%GME_h0, & + "The strength of GME tapers quadratically to zero when the bathymetric "//& + "depth is shallower than GME_H0.", & + units="m", scale=US%m_to_Z, default=1000.0) + call get_param(param_file, mdl, "GME_EFFICIENCY", CS%GME_efficiency, & + "The nondimensional prefactor multiplying the GME coefficient.", & + units="nondim", default=1.0) + call get_param(param_file, mdl, "GME_LIMITER", CS%GME_limiter, & + "The absolute maximum value the GME coefficient is allowed to take.", & + units="m2 s-1", scale=US%m_to_L**2*US%T_to_s, default=1.0e7) + endif if (CS%Laplacian .or. CS%biharmonic) then call get_param(param_file, mdl, "DT", dt, & @@ -2492,6 +2514,18 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) CS%min_grid_Kh = spacing(1.) * min_grid_sp_h2 * Idt endif if (CS%use_GME) then + CS%id_dudx_bt = register_diag_field('ocean_model', 'dudx_bt', diag%axesT1, Time, & + 'Zonal component of the barotropic shearing strain at h points', 's-1', & + conversion=US%s_to_T) + CS%id_dudy_bt = register_diag_field('ocean_model', 'dudy_bt', diag%axesB1, Time, & + 'Zonal component of the barotropic shearing strain at q points', 's-1', & + conversion=US%s_to_T) + CS%id_dvdy_bt = register_diag_field('ocean_model', 'dvdy_bt', diag%axesT1, Time, & + 'Meridional component of the barotropic shearing strain at h points', 's-1', & + conversion=US%s_to_T) + CS%id_dvdx_bt = register_diag_field('ocean_model', 'dvdx_bt', diag%axesB1, Time, & + 'Meridional component of the barotropic shearing strain at q points', 's-1', & + conversion=US%s_to_T) CS%id_GME_coeff_h = register_diag_field('ocean_model', 'GME_coeff_h', diag%axesTL, Time, & 'GME coefficient at h Points', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T) CS%id_GME_coeff_q = register_diag_field('ocean_model', 'GME_coeff_q', diag%axesBL, Time, & @@ -2501,7 +2535,8 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) endif CS%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,Time,& - 'Integral work done by lateral friction terms', & + 'Integral work done by lateral friction terms. If GME is turned on, this '//& + 'includes the GME contribution.', & 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) CS%id_FrictWorkIntz = register_diag_field('ocean_model','FrictWorkIntz',diag%axesT1,Time, & 'Depth integrated work done by lateral friction', & @@ -2531,7 +2566,8 @@ end subroutine align_aniso_tensor_to_grid !> Apply a 1-1-4-1-1 Laplacian filter one time on GME diffusive flux to reduce any !! horizontal two-grid-point noise -subroutine smooth_GME(G, GME_flux_h, GME_flux_q) +subroutine smooth_GME(CS, G, GME_flux_h, GME_flux_q) + type(hor_visc_CS), intent(in) :: CS !< Control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid real, dimension(SZI_(G),SZJ_(G)), optional, intent(inout) :: GME_flux_h!< GME diffusive flux !! at h points @@ -2542,15 +2578,18 @@ subroutine smooth_GME(G, GME_flux_h, GME_flux_q) real, dimension(SZIB_(G),SZJB_(G)) :: GME_flux_q_original real :: wc, ww, we, wn, ws ! averaging weights for smoothing integer :: i, j, k, s - do s=1,1 - ! Update halos + integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + do s=1,CS%num_smooth_gme if (present(GME_flux_h)) then - !### Work on a wider halo to eliminate this blocking send! - call pass_var(GME_flux_h, G%Domain) + ! Update halos if needed + if (s >= 2) call pass_var(GME_flux_h, G%Domain) GME_flux_h_original(:,:) = GME_flux_h(:,:) ! apply smoothing on GME - do j = G%jsc, G%jec - do i = G%isc, G%iec + do j = Jsq, Jeq+1 + do i = Isq, Ieq+1 ! skip land points if (G%mask2dT(i,j)==0.) cycle ! compute weights @@ -2558,24 +2597,22 @@ subroutine smooth_GME(G, GME_flux_h, GME_flux_q) we = 0.125 * G%mask2dT(i+1,j) ws = 0.125 * G%mask2dT(i,j-1) wn = 0.125 * G%mask2dT(i,j+1) - wc = 1.0 - (ww+we+wn+ws) - !### Add parentheses to make this rotationally invariant. + wc = 1.0 - ((ww+we)+(wn+ws)) GME_flux_h(i,j) = wc * GME_flux_h_original(i,j) & - + ww * GME_flux_h_original(i-1,j) & - + we * GME_flux_h_original(i+1,j) & - + ws * GME_flux_h_original(i,j-1) & - + wn * GME_flux_h_original(i,j+1) + + ((ww * GME_flux_h_original(i-1,j) & + + we * GME_flux_h_original(i+1,j)) & + + (ws * GME_flux_h_original(i,j-1) & + + wn * GME_flux_h_original(i,j+1))) enddo enddo endif - ! Update halos if (present(GME_flux_q)) then - !### Work on a wider halo to eliminate this blocking send! - call pass_var(GME_flux_q, G%Domain, position=CORNER, complete=.true.) + ! Update halos if needed + if (s >= 2) call pass_var(GME_flux_q, G%Domain, position=CORNER, complete=.true.) GME_flux_q_original(:,:) = GME_flux_q(:,:) ! apply smoothing on GME - do J = G%JscB, G%JecB - do I = G%IscB, G%IecB + do J = Jsq-1, Jeq + do I = Isq-1, Ieq ! skip land points if (G%mask2dBu(I,J)==0.) cycle ! compute weights @@ -2583,13 +2620,12 @@ subroutine smooth_GME(G, GME_flux_h, GME_flux_q) we = 0.125 * G%mask2dBu(I+1,J) ws = 0.125 * G%mask2dBu(I,J-1) wn = 0.125 * G%mask2dBu(I,J+1) - wc = 1.0 - (ww+we+wn+ws) - !### Add parentheses to make this rotationally invariant. + wc = 1.0 - ((ww+we)+(wn+ws)) GME_flux_q(I,J) = wc * GME_flux_q_original(I,J) & - + ww * GME_flux_q_original(I-1,J) & - + we * GME_flux_q_original(I+1,J) & - + ws * GME_flux_q_original(I,J-1) & - + wn * GME_flux_q_original(I,J+1) + + ((ww * GME_flux_q_original(I-1,J) & + + we * GME_flux_q_original(I+1,J)) & + + (ws * GME_flux_q_original(I,J-1) & + + wn * GME_flux_q_original(I,J+1))) enddo enddo endif From 149073fea0b9fa3b82880b59818710d9abbebbd5 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Sun, 20 Feb 2022 11:50:18 -0700 Subject: [PATCH 9/9] Remove hard-wired parameter in adjustEtaToFitBathymetry (#69) Subroutine adjustEtaToFitBathymetry had a hard-wired parameter (hTolerance = 0.1) controlling the tolerance when adjusting the thickness to fit the bathymetry. This patch adds an user-controlled parameter (THICKNESS_TOLERANCE), which replaces hTolerance. THICKNESS_TOLERANCE is only activated when ADJUST_THICKNESS=True. --- .../MOM_state_initialization.F90 | 32 +++++++++++++------ 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 22892817e6..8378644ea9 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -688,6 +688,8 @@ subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, f real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Interface heights, in depth units [Z ~> m]. integer :: inconsistent = 0 logical :: correct_thickness + real :: h_tolerance ! A parameter that controls the tolerance when adjusting the + ! thickness to fit the bathymetry [Z ~> m]. character(len=40) :: mdl = "initialize_thickness_from_file" ! This subroutine's name. character(len=200) :: filename, thickness_file, inputdir, mesg ! Strings for file/path integer :: i, j, k, is, ie, js, je, nz @@ -718,12 +720,18 @@ subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, f "If true, all mass below the bottom removed if the "//& "topography is shallower than the thickness input file "//& "would indicate.", default=.false., do_not_log=just_read) + if (correct_thickness) then + call get_param(param_file, mdl, "THICKNESS_TOLERANCE", h_tolerance, & + "A parameter that controls the tolerance when adjusting the "//& + "thickness to fit the bathymetry. Used when ADJUST_THICKNESS=True.", & + units="m", default=0.1, scale=US%m_to_Z, do_not_log=just_read) + endif if (just_read) return ! All run-time parameters have been read, so return. call MOM_read_data(filename, "eta", eta(:,:,:), G%Domain, scale=US%m_to_Z) if (correct_thickness) then - call adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta=G%Z_ref) + call adjustEtaToFitBathymetry(G, GV, US, eta, h, h_tolerance, dZ_ref_eta=G%Z_ref) else do k=nz,1,-1 ; do j=js,je ; do i=is,ie if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) then @@ -757,31 +765,29 @@ end subroutine initialize_thickness_from_file !! layers are contracted to ANGSTROM thickness (which may be 0). !! If the bottom most interface is above the topography then the entire column !! is dilated (expanded) to fill the void. -!! @remark{There is a (hard-wired) "tolerance" parameter such that the -!! criteria for adjustment must equal or exceed 10cm.} -subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta) +subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h, ht, dZ_ref_eta) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: eta !< Interface heights [Z ~> m]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, intent(in) :: ht !< Tolerance to exceed adjustment + !! criteria [Z ~> m] real, optional, intent(in) :: dZ_ref_eta !< The difference between the !! reference heights for bathyT and !! eta [Z ~> m], 0 by default. ! Local variables integer :: i, j, k, is, ie, js, je, nz, contractions, dilations - real :: hTolerance = 0.1 !< Tolerance to exceed adjustment criteria [Z ~> m] real :: dilate ! A factor by which the column is dilated [nondim] real :: dZ_ref ! The difference in the reference heights for G%bathyT and eta [Z ~> m] character(len=100) :: mesg is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - hTolerance = 0.1*US%m_to_Z dZ_ref = 0.0 ; if (present(dZ_ref_eta)) dZ_ref = dZ_ref_eta contractions = 0 do j=js,je ; do i=is,ie - if (-eta(i,j,nz+1) > (G%bathyT(i,j) + dZ_ref) + hTolerance) then + if (-eta(i,j,nz+1) > (G%bathyT(i,j) + dZ_ref) + ht) then eta(i,j,nz+1) = -(G%bathyT(i,j) + dZ_ref) contractions = contractions + 1 endif @@ -811,7 +817,7 @@ subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta) ! The whole column is dilated to accommodate deeper topography than ! the bathymetry would indicate. ! This should be... if ((G%mask2dt(i,j)*(eta(i,j,1)-eta(i,j,nz+1)) > 0.0) .and. & - if (-eta(i,j,nz+1) < (G%bathyT(i,j) + dZ_ref) - hTolerance) then + if (-eta(i,j,nz+1) < (G%bathyT(i,j) + dZ_ref) - ht) then dilations = dilations + 1 if (eta(i,j,1) <= eta(i,j,nz+1)) then do k=1,nz ; h(i,j,k) = (eta(i,j,1) + (G%bathyT(i,j) + dZ_ref)) / real(nz) ; enddo @@ -2402,6 +2408,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just real :: dilate ! A dilation factor to match topography [nondim] real :: missing_value_temp, missing_value_salt logical :: correct_thickness + real :: h_tolerance ! A parameter that controls the tolerance when adjusting the + ! thickness to fit the bathymetry [Z ~> m]. character(len=40) :: potemp_var, salin_var character(len=8) :: laynum @@ -2535,6 +2543,12 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just "If true, all mass below the bottom removed if the "//& "topography is shallower than the thickness input file "//& "would indicate.", default=.false., do_not_log=just_read) + if (correct_thickness) then + call get_param(PF, mdl, "THICKNESS_TOLERANCE", h_tolerance, & + "A parameter that controls the tolerance when adjusting the "//& + "thickness to fit the bathymetry. Used when ADJUST_THICKNESS=True.", & + units="m", default=0.1, scale=US%m_to_Z, do_not_log=just_read) + endif call get_param(PF, mdl, "FIT_TO_TARGET_DENSITY_IC", adjust_temperature, & "If true, all the interior layers are adjusted to "//& @@ -2732,7 +2746,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just Hmix_depth, eps_z, eps_rho, density_extrap_bug) if (correct_thickness) then - call adjustEtaToFitBathymetry(G, GV, US, zi, h, dZ_ref_eta=G%Z_ref) + call adjustEtaToFitBathymetry(G, GV, US, zi, h, h_tolerance, dZ_ref_eta=G%Z_ref) else do k=nz,1,-1 ; do j=js,je ; do i=is,ie if (zi(i,j,K) < (zi(i,j,K+1) + GV%Angstrom_Z)) then