Skip to content

Commit

Permalink
+ALE_remap_scalar with arbitrary thickness units
Browse files Browse the repository at this point in the history 
  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.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Feb 2, 2024
1 parent 76f0668 commit eeb27ec
Show file tree
Hide file tree
Showing 3 changed files with 108 additions and 44 deletions.
36 changes: 24 additions & 12 deletions src/ALE/MOM_ALE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.
Expand All @@ -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)
Expand All @@ -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.
Expand Down
60 changes: 46 additions & 14 deletions src/ALE/MOM_regridding.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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.")
Expand Down Expand Up @@ -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
Expand Down
56 changes: 38 additions & 18 deletions src/initialization/MOM_tracer_initialization_from_Z.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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"
Expand All @@ -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.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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 )
Expand Down

0 comments on commit eeb27ec

Please sign in to comment.