From eeb27ec6c4cea403925e74aa16422d755762a853 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 26 Jan 2024 18:10:41 -0500 Subject: [PATCH 1/3] +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 ) From e56a5a9c914e30fbdaceee6541b5134f9cfa5953 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 27 Jan 2024 13:10:18 -0500 Subject: [PATCH 2/3] (*)MOM_temp_salt_init_from_Z Z-unit tracer remap Revise MOM_temp_salt_initialize_from_Z in cases when Z_INIT_REMAP_GENERAL is False to call ALE_remap_scalar with vertical layer extents (in Z units) rather than layer thicknesses (in H units). When in fully non-Boussinesq mode, this same routine uses dz_to_thickness (using the full equation of state) rather than dz_to_thickness_simple to initialize the layer thicknesses. Boussinesq answers are bitwise identical, but answers can change in some fully non-Boussinesq cases. --- .../MOM_state_initialization.F90 | 36 ++++++++++++++----- 1 file changed, 27 insertions(+), 9 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 7dfced262b..855a6f2aa0 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -92,6 +92,7 @@ module MOM_state_initialization use MOM_ALE, only : ALE_remap_scalar, ALE_regrid_accelerated, TS_PLM_edge_values use MOM_regridding, only : regridding_CS, set_regrid_params, getCoordinateResolution use MOM_regridding, only : regridding_main, regridding_preadjust_reqs, convective_adjustment +use MOM_regridding, only : set_dz_neglect use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer, homogenize_field use MOM_oda_incupd, only: oda_incupd_CS, initialize_oda_incupd_fixed, initialize_oda_incupd @@ -2483,6 +2484,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just real, dimension(:,:,:), allocatable :: h1 ! Thicknesses on the input grid [H ~> m or kg m-2]. real, dimension(:,:,:), allocatable :: dz_interface ! Change in position of interface due to ! regridding [H ~> m or kg m-2] + 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] real :: zTopOfCell, zBottomOfCell ! Heights in Z units [Z ~> m]. type(regridding_CS) :: regridCS ! Regridding parameters and work arrays type(remapping_CS) :: remapCS ! Remapping parameters and work arrays @@ -2768,6 +2773,11 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just frac_shelf_h=frac_shelf_h ) deallocate( dz_interface ) + + call ALE_remap_scalar(remapCS, G, GV, nkd, h1, tmpT1dIn, h, tv%T, all_cells=remap_full_column, & + old_remap=remap_old_alg, answer_date=remap_answer_date ) + call ALE_remap_scalar(remapCS, G, GV, nkd, h1, tmpS1dIn, h, tv%S, all_cells=remap_full_column, & + old_remap=remap_old_alg, answer_date=remap_answer_date ) else ! This is the old way of initializing to z* coordinates only allocate( hTarget(nz) ) @@ -2788,16 +2798,24 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just enddo ; enddo deallocate( hTarget ) - ! This is a simple conversion of the target grid to thickness units that may not be - ! appropriate in non-Boussinesq mode. - call dz_to_thickness_simple(dz, h, G, GV, US) + dz_neglect = set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) + call ALE_remap_scalar(remapCS, G, GV, nkd, dz1, tmpT1dIn, dz, tv%T, all_cells=remap_full_column, & + old_remap=remap_old_alg, answer_date=remap_answer_date, & + H_neglect=dz_neglect, H_neglect_edge=dz_neglect_edge) + call ALE_remap_scalar(remapCS, G, GV, nkd, dz1, tmpS1dIn, dz, tv%S, all_cells=remap_full_column, & + old_remap=remap_old_alg, answer_date=remap_answer_date, & + H_neglect=dz_neglect, H_neglect_edge=dz_neglect_edge) + + if (GV%Boussinesq .or. GV%semi_Boussinesq) then + ! This is a simple conversion of the target grid to thickness units that is not + ! appropriate in non-Boussinesq mode. + call dz_to_thickness_simple(dz, h, G, GV, US) + else + ! Convert dz into thicknesses in units of H using the equation of state as appropriate. + call dz_to_thickness(dz, tv, h, G, GV, US) + endif endif - call ALE_remap_scalar(remapCS, G, GV, nkd, h1, tmpT1dIn, h, tv%T, all_cells=remap_full_column, & - old_remap=remap_old_alg, answer_date=remap_answer_date ) - call ALE_remap_scalar(remapCS, G, GV, nkd, h1, tmpS1dIn, h, tv%S, all_cells=remap_full_column, & - old_remap=remap_old_alg, answer_date=remap_answer_date ) - deallocate( dz1 ) deallocate( h1 ) deallocate( tmpT1dIn ) @@ -2879,7 +2897,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just ks, G, GV, US, PF, just_read) endif - ! Now convert thicknesses to units of H. + ! Now convert dz into thicknesses in units of H. call dz_to_thickness(dz, tv, h, G, GV, US) endif ! useALEremapping From 46ef67c443c072aa239d9a74d46939e0929b6525 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 28 Jan 2024 10:16:54 -0500 Subject: [PATCH 3/3] *+non-Boussinesq revisions to MOM_generic_tracer Revised initialize_MOM_generic_tracer to use thickness_to_dz to get the layer vertical extents and then provide these to MOM_initialize_tracer_from_Z to read in initial generic tracer concentrations from Z-space files. The previous approach inappropriately used an simple multiplicative rescaling (via a call to dz_to_thickness_simple in MOM_initialize_tracer_from_Z) by a factor that includes the Boussinesq reference density when in non-Boussinesq mode. A new thermo_vars_type arguments was added to initialize_MOM_generic_tracer to allow for this change. Also revised MOM_generic_tracer_column_physics to use thickness_to_dz instead of a simple multiplicative rescaling to get the layer vertical extents (in m) that are used in calls to generic_tracer_source. The multiplicative factor that was used previously (GV%H_to_m) includes the Boussinesq reference density and hence is inappropriate in non-Boussinesq mode; using thickness_to_dz avoids this. Also added comments documenting the meaning and units of about 30 real variables in the MOM_generic_tracer routines. There is a new mandatory argument to initialize_MOM_generic_tracer. All Boussinseq mode answers are bitwise identical, but in non-Boussinesq mode mode generic tracer answers are changed by avoiding the use of the Boussinesq reference density in several places. --- src/tracer/MOM_generic_tracer.F90 | 144 ++++++++++++++----------- src/tracer/MOM_tracer_flow_control.F90 | 2 +- 2 files changed, 84 insertions(+), 62 deletions(-) diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 131110e6b2..7f550d8de5 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -38,6 +38,7 @@ module MOM_generic_tracer use MOM_forcing_type, only : forcing, optics_type use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type + use MOM_interface_heights, only : thickness_to_dz use MOM_io, only : file_exists, MOM_read_data, slasher use MOM_open_boundary, only : ocean_OBC_type use MOM_open_boundary, only : register_obgc_segments, fill_obgc_segments @@ -75,8 +76,10 @@ module MOM_generic_tracer character(len = 200) :: IC_file !< The file in which the generic tracer initial values can !! be found, or an empty string for internal initialization. logical :: Z_IC_file !< If true, the generic_tracer IC_file is in Z-space. The default is false. - real :: tracer_IC_val = 0.0 !< The initial value assigned to tracers. - real :: tracer_land_val = -1.0 !< The values of tracers used where land is masked out. + real :: tracer_IC_val = 0.0 !< The initial value assigned to tracers, in + !! concentration units [conc] + real :: tracer_land_val = -1.0 !< The values of tracers used where land is masked out, in + !! concentration units [conc] logical :: tracers_may_reinit !< If true, tracers may go through the !! initialization code if they are not found in the restart files. @@ -102,6 +105,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) type(tracer_registry_type), pointer :: tr_Reg !< Pointer to the control structure for the tracer !! advection and diffusion module. type(MOM_restart_CS), target, intent(inout) :: restart_CS !< MOM restart control struct + ! Local variables logical :: register_MOM_generic_tracer logical :: obc_has @@ -113,14 +117,17 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) ! These can be overridden later in via the field manager? integer :: ntau, axes(3) - type(g_tracer_type), pointer :: g_tracer,g_tracer_next - character(len=fm_string_len) :: g_tracer_name,longname,units - character(len=fm_string_len) :: obc_src_file_name,obc_src_field_name - real :: lfac_in,lfac_out - real, dimension(:,:,:,:), pointer :: tr_field - real, dimension(:,:,:), pointer :: tr_ptr - real, dimension(HI%isd:HI%ied, HI%jsd:HI%jed,GV%ke) :: grid_tmask - integer, dimension(HI%isd:HI%ied, HI%jsd:HI%jed) :: grid_kmt + type(g_tracer_type), pointer :: g_tracer, g_tracer_next + character(len=fm_string_len) :: g_tracer_name, longname,units + character(len=fm_string_len) :: obc_src_file_name, obc_src_field_name + real :: lfac_in ! Multiplicative factor used in setting the tracer-specific inverse length + ! scales associated with inflowing tracer reservoirs at OBCs [nondim] + real :: lfac_out ! Multiplicative factor used in setting the tracer-specific inverse length + ! scales associated with outflowing tracer reservoirs at OBCs [nondim] + real, dimension(:,:,:,:), pointer :: tr_field ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(:,:,:), pointer :: tr_ptr ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)) :: grid_tmask ! A 3-d copy of G%mask2dT [nondim] + integer, dimension(SZI_(HI),SZJ_(HI)) :: grid_kmt ! A 2-d array of nk register_MOM_generic_tracer = .false. if (associated(CS)) then @@ -141,7 +148,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) ! Read all relevant parameters and write them to the model log. call log_version(param_file, sub_name, version, "") call get_param(param_file, sub_name, "GENERIC_TRACER_IC_FILE", CS%IC_file, & - "The file in which the generic trcer initial values can "//& + "The file in which the generic tracer initial values can "//& "be found, or an empty string for internal initialization.", & default=" ") if ((len_trim(CS%IC_file) > 0) .and. (scan(CS%IC_file,'/') == 0)) then @@ -169,7 +176,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) !Fields cannot be diag registered as they are allocated and have to registered later. grid_tmask(:,:,:) = 0.0 - grid_kmt(:,:) = 0.0 + grid_kmt(:,:) = 0 axes(:) = -1 ! @@ -222,23 +229,26 @@ end function register_MOM_generic_tracer !> Register OBC segments for generic tracers subroutine register_MOM_generic_tracer_segments(CS, GV, OBC, tr_Reg, param_file) - type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, - !! where, and what open boundary conditions are used. - type(tracer_registry_type), pointer :: tr_Reg !< Pointer to the control structure for the tracer - !! advection and diffusion module. - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, + !! where, and what open boundary conditions are used. + type(tracer_registry_type), pointer :: tr_Reg !< Pointer to the control structure for the tracer + !! advection and diffusion module. + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + ! Local variables logical :: obc_has ! This include declares and sets the variable "version". # include "version_variable.h" - character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer_segments' type(g_tracer_type), pointer :: g_tracer,g_tracer_next character(len=fm_string_len) :: g_tracer_name - character(len=fm_string_len) :: obc_src_file_name,obc_src_field_name - real :: lfac_in,lfac_out + character(len=fm_string_len) :: obc_src_file_name, obc_src_field_name + real :: lfac_in ! Multiplicative factor used in setting the tracer-specific inverse length + ! scales associated with inflowing tracer reservoirs at OBCs [nondim] + real :: lfac_out ! Multiplicative factor used in setting the tracer-specific inverse length + ! scales associated with outflowing tracer reservoirs at OBCs [nondim] if (.NOT. associated(OBC)) return !Get the tracer list @@ -266,6 +276,7 @@ subroutine register_MOM_generic_tracer_segments(CS, GV, OBC, tr_Reg, param_file) enddo end subroutine register_MOM_generic_tracer_segments + !> Initialize phase II: Initialize required variables for generic tracers !! There are some steps of initialization that cannot be done in register_MOM_generic_tracer !! This is the place and time to do them: @@ -275,15 +286,17 @@ end subroutine register_MOM_generic_tracer_segments !! !! This subroutine initializes the NTR tracer fields in tr(:,:,:,:) !! and it sets up the tracer output. - subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, diag, OBC, CS, & - sponge_CSp, ALE_sponge_CSp) + subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, tv, param_file, diag, OBC, & + CS, sponge_CSp, ALE_sponge_CSp) logical, intent(in) :: restart !< .true. if the fields have already been !! read from a restart file. type(time_type), target, intent(in) :: day !< Time of the start of the run. type(ocean_grid_type), intent(inout) :: 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)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + 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 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters type(diag_ctrl), target, intent(in) :: diag !< Regulates diagnostic output. type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, @@ -298,10 +311,11 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, integer :: i, j, k, isc, iec, jsc, jec, nk type(g_tracer_type), pointer :: g_tracer,g_tracer_next character(len=fm_string_len) :: g_tracer_name - real, dimension(:,:,:,:), pointer :: tr_field - real, dimension(:,:,:), pointer :: tr_ptr - real, dimension(G%isd:G%ied, G%jsd:G%jed, 1:GV%ke) :: grid_tmask - integer, dimension(G%isd:G%ied, G%jsd:G%jed) :: grid_kmt + real, dimension(:,:,:,:), pointer :: tr_field ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(:,:,:), pointer :: tr_ptr ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Layer vertical extent [Z ~> m] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: grid_tmask ! A 3-d copy of G%mask2dT [nondim] + integer, dimension(SZI_(G),SZJ_(G)) :: grid_kmt ! A 2-d array of nk !! 2010/02/04 Add code to re-initialize Generic Tracers if needed during a model simulation !! By default, restart cpio should not contain a Generic Tracer IC file and step below will be skipped. @@ -316,6 +330,8 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, !For each tracer name get its fields g_tracer=>CS%g_tracer_list + call thickness_to_dz(h, tv, dz, G, GV, US) + do if (INDEX(CS%IC_file, '_NULL_') /= 0) then call MOM_error(WARNING, "The name of the IC_file "//trim(CS%IC_file)//& @@ -335,12 +351,11 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, "initializing generic tracer "//trim(g_tracer_name)//& " using MOM_initialize_tracer_from_Z ") - call MOM_initialize_tracer_from_Z(h, tr_ptr, G, GV, US, param_file, & - src_file = g_tracer%src_file, & - src_var_nam = g_tracer%src_var_name, & - src_var_unit_conversion = g_tracer%src_var_unit_conversion,& - src_var_record = g_tracer%src_var_record, & - src_var_gridspec = g_tracer%src_var_gridspec ) + call MOM_initialize_tracer_from_Z(dz, tr_ptr, G, GV, US, param_file, & + src_file=g_tracer%src_file, src_var_nam=g_tracer%src_var_name, & + src_var_unit_conversion=g_tracer%src_var_unit_conversion, & + src_var_record=g_tracer%src_var_record, src_var_gridspec=g_tracer%src_var_gridspec, & + h_in_Z_units=.true.) !Check/apply the bounds for each g_tracer do k=1,nk ; do j=jsc,jec ; do i=isc,iec @@ -466,8 +481,9 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables type(optics_type), intent(in) :: optics !< The structure containing optical properties. - real, optional, intent(in) :: evap_CFL_limit !< Limits how much water can be fluxed out of - !! the top layer Stored previously in diabatic CS. + real, optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can + !! be fluxed out of the top layer in a timestep [nondim] + ! Stored previously in diabatic CS. real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which fluxes !! can be applied [H ~> m or kg m-2] ! Stored previously in diabatic CS. @@ -479,14 +495,17 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, type(g_tracer_type), pointer :: g_tracer, g_tracer_next character(len=fm_string_len) :: g_tracer_name - real, dimension(:,:), pointer :: stf_array,trunoff_array,runoff_tracer_flux_array + real, dimension(:,:), pointer :: stf_array ! The surface flux of the tracer [conc kg m-2 s-1] + real, dimension(:,:), pointer :: trunoff_array ! The tracer concentration in the river runoff [conc] + real, dimension(:,:), pointer :: runoff_tracer_flux_array ! The runoff tracer flux [conc kg m-2 s-1] - real :: surface_field(SZI_(G),SZJ_(G)) + real :: surface_field(SZI_(G),SZJ_(G)) ! The surface value of some field, here only used for salinity [S ~> ppt] real :: dz_ml(SZI_(G),SZJ_(G)) ! The mixed layer depth in the MKS units used for generic tracers [m] - real :: sosga + real :: sosga ! The global mean surface salinity [ppt] - real, dimension(G%isd:G%ied,G%jsd:G%jed,GV%ke) :: rho_dzt, dzt - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: rho_dzt ! Layer mass per unit area [kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dzt ! Layer vertical extents [m] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! A work array of thicknesses [H ~> m or kg m-2] integer :: i, j, k, isc, iec, jsc, jec, nk isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke @@ -536,14 +555,15 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, ! rho_dzt(:,:,:) = GV%H_to_kg_m2 * GV%Angstrom_H - do k = 1, nk ; do j = jsc, jec ; do i = isc, iec !{ + do k=1,nk ; do j=jsc,jec ; do i=isc,iec rho_dzt(i,j,k) = GV%H_to_kg_m2 * h_old(i,j,k) - enddo ; enddo ; enddo !} + enddo ; enddo ; enddo dzt(:,:,:) = 1.0 - do k = 1, nk ; do j = jsc, jec ; do i = isc, iec !{ - dzt(i,j,k) = GV%H_to_m * h_old(i,j,k) - enddo ; enddo ; enddo !} + call thickness_to_dz(h_old, tv, dzt, G, GV, US) + do k=1,nk ; do j=jsc,jec ; do i=isc,iec + dzt(i,j,k) = US%Z_to_m * dzt(i,j,k) + enddo ; enddo ; enddo dz_ml(:,:) = 0.0 do j=jsc,jec ; do i=isc,iec surface_field(i,j) = tv%S(i,j,1) @@ -639,8 +659,8 @@ function MOM_generic_tracer_stock(h, stocks, G, GV, CS, names, units, stock_inde ! Local variables type(g_tracer_type), pointer :: g_tracer, g_tracer_next - real, dimension(:,:,:,:), pointer :: tr_field - real, dimension(:,:,:), pointer :: tr_ptr + real, dimension(:,:,:,:), pointer :: tr_field ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(:,:,:), pointer :: tr_ptr ! A pointer to a generic tracer field, in concentration units [conc] character(len=128), parameter :: sub_name = 'MOM_generic_tracer_stock' integer :: m @@ -802,7 +822,7 @@ subroutine array_global_min_max(tr_array, tmask, isd, jsd, isc, iec, jsc, jec, n real :: tmax, tmin ! Maximum and minimum tracer values, in the same units as tr_array real :: tmax0, tmin0 ! First-guest values of tmax and tmin. integer :: itmax, jtmax, ktmax, itmin, jtmin, ktmin - real :: fudge ! A factor that is close to 1 that is used to find the location of the extrema. + real :: fudge ! A factor that is close to 1 that is used to find the location of the extrema [nondim]. ! arrays to enable vectorization integer :: iminarr(3), imaxarr(3) @@ -853,7 +873,7 @@ subroutine array_global_min_max(tr_array, tmask, isd, jsd, isc, iec, jsc, jec, n ! Now find the location of the global extrema. ! - ! Note that the fudge factor above guarantees that the location of max (min) is uinque, + ! Note that the fudge factor above guarantees that the location of max (min) is unique, ! since tmax0 (tmin0) has slightly different values on each processor. ! Otherwise, the function tr_array(i,j,k) could be equal to global max (min) at more ! than one point in space and this would be a much more difficult problem to solve. @@ -899,16 +919,16 @@ subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, GV, CS) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. -! Local variables - real :: sosga + ! Local variables + real :: sosga ! The global mean surface salinity [ppt] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV),1) :: rho0 ! An unused array of densities [kg m-3] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dzt ! Layer vertical extents [m] character(len=128), parameter :: sub_name = 'MOM_generic_tracer_surface_state' - real, dimension(G%isd:G%ied,G%jsd:G%jed,1:GV%ke,1) :: rho0 - real, dimension(G%isd:G%ied,G%jsd:G%jed,1:GV%ke) :: dzt !Set coupler values !nnz: fake rho0 - rho0=1.0 + rho0(:,:,:,:) = 1.0 dzt(:,:,:) = GV%H_to_m * h(:,:,:) @@ -937,7 +957,7 @@ subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, GV, CS) !Niki: The problem with calling diagnostic outputs here is that this subroutine is called every dt_cpld ! hence if dt_therm > dt_cpld we get output (and contribution to the mean) at times that tracers ! had not been updated. - ! Moving this to the end of column physics subrotuine fixes this issue. + ! Moving this to the end of column physics subroutine fixes this issue. end subroutine MOM_generic_tracer_surface_state @@ -976,7 +996,7 @@ end subroutine MOM_generic_flux_init subroutine MOM_generic_tracer_fluxes_accumulate(flux_tmp, weight) type(forcing), intent(in) :: flux_tmp !< A structure containing pointers to !! thermodynamic and tracer forcing fields. - real, intent(in) :: weight !< A weight for accumulating this flux + real, intent(in) :: weight !< A weight for accumulating this flux [nondim] call generic_tracer_coupler_accumulate(flux_tmp%tr_fluxes, weight) @@ -986,10 +1006,12 @@ end subroutine MOM_generic_tracer_fluxes_accumulate subroutine MOM_generic_tracer_get(name,member,array, CS) character(len=*), intent(in) :: name !< Name of requested tracer. character(len=*), intent(in) :: member !< The tracer element to return. - real, dimension(:,:,:), intent(out) :: array !< Array filled by this routine. - type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + real, dimension(:,:,:), intent(out) :: array !< Array filled by this routine, in arbitrary units [A] + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. - real, dimension(:,:,:), pointer :: array_ptr + ! Local variables + real, dimension(:,:,:), pointer :: array_ptr ! The tracer in the generic tracer structures, in + ! arbitrary units [A] character(len=128), parameter :: sub_name = 'MOM_generic_tracer_get' call g_tracer_get_pointer(CS%g_tracer_list,name,member,array_ptr) diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index c8ce2f5f75..6d035e1d27 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -340,7 +340,7 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag call initialize_CFC_cap(restart, day, G, GV, US, h, diag, OBC, CS%CFC_cap_CSp) if (CS%use_MOM_generic_tracer) & - call initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, diag, OBC, & + call initialize_MOM_generic_tracer(restart, day, G, GV, US, h, tv, param_file, diag, OBC, & CS%MOM_generic_tracer_CSp, sponge_CSp, ALE_sponge_CSp) if (CS%use_pseudo_salt_tracer) & call initialize_pseudo_salt_tracer(restart, day, G, GV, US, h, diag, OBC, CS%pseudo_salt_tracer_CSp, &