From 915cfe225e5e3c42103f20501c60e6aa90cd613c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 26 Jan 2024 18:10:41 -0500 Subject: [PATCH] +ALE_remap_scalar with arbitrary thickness units This commit add the new optional arguments h_neglect and h_neglect_edge to ALE_remap_scalar to allow for the thicknesses used in this routine to be provided in any self-consistent units, including [Z ~> m], instead of just [H ~> m or kg m-2]. To help make use of this new capability, this commit also adds the new functions set_h_neglect and set_dz_neglect to the MOM_regridding module. build_grid_rho and build_grid_HyCOM1 have been refactored to use set_h_neglect in place of the corresponding duplicated code blocks. This commit also adds the new optional argument h_in_Z_units to MOM_initialize_tracer_from_Z, which in turn uses this new capability for ALE_remap_scalar to use vertical layer extents (in Z units) rather than thicknesses (in H units). Although there are new optional arguments to public interfaces, they are not yet being exercised with this commit so no answers are changed. Moreover, even if they were being exercised, all Boussinesq solutions would give identical answers. --- src/ALE/MOM_ALE.F90 | 36 +++++++---- src/ALE/MOM_regridding.F90 | 60 ++++++++++++++----- .../MOM_tracer_initialization_from_Z.F90 | 56 +++++++++++------ 3 files changed, 108 insertions(+), 44 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 77ee1192a2..543d77a0f3 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -1260,16 +1260,17 @@ end subroutine mask_near_bottom_vel !! h_dst must be dimensioned as a model array with GV%ke layers while h_src can !! have an arbitrary number of layers specified by nk_src. subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap, & - answers_2018, answer_date ) + answers_2018, answer_date, h_neglect, h_neglect_edge) type(remapping_CS), intent(in) :: CS !< Remapping control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure integer, intent(in) :: nk_src !< Number of levels on source grid real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: h_src !< Level thickness of source grid - !! [H ~> m or kg m-2] + !! [H ~> m or kg m-2] or other units + !! if H_neglect is provided real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: s_src !< Scalar on source grid, in arbitrary units [A] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(in) :: h_dst !< Level thickness of destination grid - !! [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(in) :: h_dst !< Level thickness of destination grid in the + !! same units as h_src, often [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(inout) :: s_dst !< Scalar on destination grid, in the same !! arbitrary units as s_src [A] logical, optional, intent(in) :: all_cells !< If false, only reconstruct for @@ -1283,10 +1284,16 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c !! use more robust forms of the same expressions. integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use !! for remapping + real, optional, intent(in) :: h_neglect !< A negligibly small thickness used in + !! remapping cell reconstructions, in the same + !! units as h_src, often [H ~> m or kg m-2] + real, optional, intent(in) :: h_neglect_edge !< A negligibly small thickness used in + !! remapping edge value calculations, in the same + !! units as h_src, often [H ~> m or kg m-2] ! Local variables integer :: i, j, k, n_points real :: dx(GV%ke+1) ! Change in interface position [H ~> m or kg m-2] - real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2] + real :: h_neg, h_neg_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2] logical :: ignore_vanished_layers, use_remapping_core_w, use_2018_remap ignore_vanished_layers = .false. @@ -1297,12 +1304,17 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c use_2018_remap = .true. ; if (present(answers_2018)) use_2018_remap = answers_2018 if (present(answer_date)) use_2018_remap = (answer_date < 20190101) - if (.not.use_2018_remap) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then - h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 + if (present(h_neglect)) then + h_neg = h_neglect + h_neg_edge = h_neg ; if (present(h_neglect_edge)) h_neg_edge = h_neglect_edge else - h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 + if (.not.use_2018_remap) then + h_neg = GV%H_subroundoff ; h_neg_edge = GV%H_subroundoff + elseif (GV%Boussinesq) then + h_neg = GV%m_to_H*1.0e-30 ; h_neg_edge = GV%m_to_H*1.0e-10 + else + h_neg = GV%kg_m2_to_H*1.0e-30 ; h_neg_edge = GV%kg_m2_to_H*1.0e-10 + endif endif !$OMP parallel do default(shared) firstprivate(n_points,dx) @@ -1318,10 +1330,10 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c if (use_remapping_core_w) then call dzFromH1H2( n_points, h_src(i,j,1:n_points), GV%ke, h_dst(i,j,:), dx ) call remapping_core_w(CS, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), & - GV%ke, dx, s_dst(i,j,:), h_neglect, h_neglect_edge) + GV%ke, dx, s_dst(i,j,:), h_neg, h_neg_edge) else call remapping_core_h(CS, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), & - GV%ke, h_dst(i,j,:), s_dst(i,j,:), h_neglect, h_neglect_edge) + GV%ke, h_dst(i,j,:), s_dst(i,j,:), h_neg, h_neg_edge) endif else s_dst(i,j,:) = 0. diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 8ef0679358..904164c8e7 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -144,6 +144,7 @@ module MOM_regridding public getCoordinateResolution, getCoordinateInterfaces public getCoordinateUnits, getCoordinateShortName, getStaticThickness public DEFAULT_COORDINATE_MODE +public set_h_neglect, set_dz_neglect public get_zlike_CS, get_sigma_CS, get_rho_CS !> Documentation for coordinate options @@ -1416,13 +1417,7 @@ subroutine build_rho_grid( G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS, #endif logical :: ice_shelf - if (CS%remap_answer_date >= 20190101) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then - h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 - else - h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 - endif + h_neglect = set_h_neglect(GV, CS%remap_answer_date, h_neglect_edge) nz = GV%ke ice_shelf = present(frac_shelf_h) @@ -1575,13 +1570,7 @@ subroutine build_grid_HyCOM1( G, GV, US, h, nom_depth_H, tv, h_new, dzInterface, real :: z_top_col, totalThickness logical :: ice_shelf - if (CS%remap_answer_date >= 20190101) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then - h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 - else - h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 - endif + h_neglect = set_h_neglect(GV, CS%remap_answer_date, h_neglect_edge) if (.not.CS%target_density_set) call MOM_error(FATAL, "build_grid_HyCOM1 : "//& "Target densities must be set before build_grid_HyCOM1 is called.") @@ -2095,6 +2084,49 @@ subroutine write_regrid_file( CS, GV, filepath ) end subroutine write_regrid_file +!> Set appropriate values for the negligible thicknesses used for remapping based on an answer date. +function set_h_neglect(GV, remap_answer_date, h_neglect_edge) result(h_neglect) + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + integer, intent(in) :: remap_answer_date !< The vintage of the expressions to use + !! for remapping. Values below 20190101 recover the + !! remapping answers from 2018. Higher values use more + !! robust forms of the same remapping algorithms. + real, intent(out) :: h_neglect_edge !< A negligibly small thickness used in + !! remapping edge value calculations [H ~> m or kg m-2] + real :: h_neglect !< A negligibly small thickness used in + !! remapping cell reconstructions [H ~> m or kg m-2] + + if (remap_answer_date >= 20190101) then + h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff + elseif (GV%Boussinesq) then + h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 + else + h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 + endif +end function set_h_neglect + +!> Set appropriate values for the negligible vertical layer extents used for remapping based on an answer date. +function set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) result(dz_neglect) + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + integer, intent(in) :: remap_answer_date !< The vintage of the expressions to use + !! for remapping. Values below 20190101 recover the + !! remapping answers from 2018. Higher values use more + !! robust forms of the same remapping algorithms. + real, intent(out) :: dz_neglect_edge !< A negligibly small vertical layer extent + !! used in remapping edge value calculations [Z ~> m] + real :: dz_neglect !< A negligibly small vertical layer extent + !! used in remapping cell reconstructions [Z ~> m] + + if (remap_answer_date >= 20190101) then + dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff + elseif (GV%Boussinesq) then + dz_neglect = US%m_to_Z*1.0e-30 ; dz_neglect_edge = US%m_to_Z*1.0e-10 + else + dz_neglect = GV%kg_m2_to_H * (GV%H_to_m*US%m_to_Z) * 1.0e-30 + dz_neglect_edge = GV%kg_m2_to_H * (GV%H_to_m*US%m_to_Z) * 1.0e-10 + endif +end function set_dz_neglect !------------------------------------------------------------------------------ !> Query the fixed resolution data diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index 808430df2c..5a172b5d97 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -3,20 +3,21 @@ module MOM_tracer_initialization_from_Z ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_debugging, only : hchksum -use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end -use MOM_cpu_clock, only : CLOCK_ROUTINE, CLOCK_LOOP -use MOM_domains, only : pass_var +use MOM_debugging, only : hchksum +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock, only : CLOCK_ROUTINE, CLOCK_LOOP +use MOM_domains, only : pass_var use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint -use MOM_file_parser, only : get_param, param_file_type, log_version -use MOM_grid, only : ocean_grid_type +use MOM_file_parser, only : get_param, param_file_type, log_version +use MOM_grid, only : ocean_grid_type use MOM_horizontal_regridding, only : myStats, horiz_interp_and_extrap_tracer use MOM_interface_heights, only : dz_to_thickness_simple -use MOM_remapping, only : remapping_CS, initialize_remapping -use MOM_unit_scaling, only : unit_scale_type -use MOM_verticalGrid, only : verticalGrid_type -use MOM_ALE, only : ALE_remap_scalar +use MOM_regridding, only : set_dz_neglect +use MOM_remapping, only : remapping_CS, initialize_remapping +use MOM_unit_scaling, only : unit_scale_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_ALE, only : ALE_remap_scalar implicit none ; private @@ -36,12 +37,13 @@ module MOM_tracer_initialization_from_Z !> Initializes a tracer from a z-space data file, including any lateral regridding that is needed. subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_nam, & src_var_unit_conversion, src_var_record, homogenize, & - useALEremapping, remappingScheme, src_var_gridspec ) + useALEremapping, remappingScheme, src_var_gridspec, h_in_Z_units ) type(ocean_grid_type), intent(inout) :: 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 thickness [H ~> m or kg m-2]. + intent(in) :: h !< Layer thicknesses, in [H ~> m or kg m-2] or + !! [Z ~> m] depending on the value of h_in_Z_units. real, dimension(:,:,:), pointer :: tr !< Pointer to array to be initialized [CU ~> conc] type(param_file_type), intent(in) :: PF !< parameter file character(len=*), intent(in) :: src_file !< source filename @@ -54,12 +56,18 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ character(len=*), optional, intent(in) :: remappingScheme !< remapping scheme to use. character(len=*), optional, intent(in) :: src_var_gridspec !< Source variable name in a gridspec file. !! This is not implemented yet. + logical, optional, intent(in) :: h_in_Z_units !< If present and true, the input grid + !! thicknesses are in the units of height + !! ([Z ~> m]) instead of the usual units of + !! thicknesses ([H ~> m or kg m-2]) + ! Local variables real :: land_fill = 0.0 ! A value to use to replace missing values [CU ~> conc] real :: convert ! A conversion factor into the model's internal units [CU conc-1 ~> 1] integer :: recnum character(len=64) :: remapScheme logical :: homog, useALE + logical :: h_is_in_Z_units ! This include declares and sets the variable "version". # include "version_variable.h" @@ -84,6 +92,10 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ type(verticalGrid_type) :: GV_loc ! A temporary vertical grid structure real :: missing_value ! A value indicating that there is no valid input data at this point [CU ~> conc] + real :: dz_neglect ! A negligibly small vertical layer extent used in + ! remapping cell reconstructions [Z ~> m] + real :: dz_neglect_edge ! A negligibly small vertical layer extent used in + ! remapping edge value calculations [Z ~> m] integer :: nPoints ! The number of valid input data points in a column integer :: id_clock_routine, id_clock_ALE integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. @@ -143,6 +155,8 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ convert = 1.0 if (PRESENT(src_var_unit_conversion)) convert = src_var_unit_conversion + h_is_in_Z_units = .false. ; if (present(h_in_Z_units)) h_is_in_Z_units = h_in_Z_units + call horiz_interp_and_extrap_tracer(src_file, src_var_nam, recnum, & G, tr_z, mask_z, z_in, z_edges_in, missing_value, & scale=convert, homogenize=homog, m_to_Z=US%m_to_Z, answer_date=hor_regrid_answer_date) @@ -185,12 +199,18 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ dzSrc(i,j,:) = h1(:) enddo ; enddo - ! Equation of state data is not available, so a simpler rescaling will have to suffice, - ! but it might be problematic in non-Boussinesq mode. - GV_loc = GV ; GV_loc%ke = kd - call dz_to_thickness_simple(dzSrc, hSrc, G, GV_loc, US) - - call ALE_remap_scalar(remapCS, G, GV, kd, hSrc, tr_z, h, tr, all_cells=.false., answer_date=remap_answer_date ) + if (h_is_in_Z_units) then + dz_neglect = set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) + call ALE_remap_scalar(remapCS, G, GV, kd, hSrc, tr_z, h, tr, all_cells=.false., answer_date=remap_answer_date, & + H_neglect=dz_neglect, H_neglect_edge=dz_neglect_edge) + else + ! Equation of state data is not available, so a simpler rescaling will have to suffice, + ! but it might be problematic in non-Boussinesq mode. + GV_loc = GV ; GV_loc%ke = kd + call dz_to_thickness_simple(dzSrc, hSrc, G, GV_loc, US) + + call ALE_remap_scalar(remapCS, G, GV, kd, hSrc, tr_z, h, tr, all_cells=.false., answer_date=remap_answer_date ) + endif deallocate( hSrc ) deallocate( dzSrc )