diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 0d325292f5..813d7110d1 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -20,7 +20,7 @@ module MOM_regridding use regrid_consts, only : state_dependent, coordinateUnits use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE use regrid_consts, only : REGRIDDING_LAYER, REGRIDDING_ZSTAR -use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SIGMA +use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SCALAR, REGRIDDING_SIGMA use regrid_consts, only : REGRIDDING_ARBITRARY, REGRIDDING_SIGMA_SHELF_ZSTAR use regrid_consts, only : REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, REGRIDDING_ADAPTIVE use regrid_interp, only : interp_CS_type, set_interp_scheme, set_interp_extrap, set_interp_answer_date @@ -28,6 +28,7 @@ module MOM_regridding use coord_zlike, only : init_coord_zlike, zlike_CS, set_zlike_params, build_zstar_column, end_coord_zlike use coord_sigma, only : init_coord_sigma, sigma_CS, set_sigma_params, build_sigma_column, end_coord_sigma use coord_rho, only : init_coord_rho, rho_CS, set_rho_params, build_rho_column, end_coord_rho +use coord_scalar, only : init_coord_scalar, scalar_CS, set_scalar_params, build_scalar_column, end_coord_scalar use coord_rho, only : old_inflate_layers_1d use coord_hycom, only : init_coord_hycom, hycom_CS, set_hycom_params, build_hycom1_column, end_coord_hycom use coord_adapt, only : init_coord_adapt, adapt_CS, set_adapt_params, build_adapt_column, end_coord_adapt @@ -118,11 +119,15 @@ module MOM_regridding !! Higher values use more robust forms of the same remapping expressions. integer :: remap_answer_date = 99991231 + !> If true, use histogram procedure to map between source and target grids for diagnostics + logical :: histogram_extensive_diags = .false. + logical :: use_hybgen_unmix = .false. !< If true, use the hybgen unmixing code before remapping type(zlike_CS), pointer :: zlike_CS => null() !< Control structure for z-like coordinate generator type(sigma_CS), pointer :: sigma_CS => null() !< Control structure for sigma coordinate generator type(rho_CS), pointer :: rho_CS => null() !< Control structure for rho coordinate generator + type(scalar_CS), pointer :: scalar_CS => null() !< Control structure for scalar coordinate generator type(hycom_CS), pointer :: hycom_CS => null() !< Control structure for hybrid coordinate generator type(adapt_CS), pointer :: adapt_CS => null() !< Control structure for adaptive coordinate generator type(hybgen_regrid_CS), pointer :: hybgen_CS => NULL() !< Control structure for hybgen regridding @@ -141,7 +146,8 @@ module MOM_regridding 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 +public get_zlike_CS, get_sigma_CS, get_rho_CS, get_scalar_CS +public check_if_histogram_extensive_diags !> Documentation for coordinate options character(len=*), parameter, public :: regriddingCoordinateModeDoc = & @@ -150,6 +156,7 @@ module MOM_regridding " SIGMA_SHELF_ZSTAR - stretched geopotential z* ignoring shelf\n"//& " SIGMA - terrain following coordinates\n"//& " RHO - continuous isopycnal\n"//& + " SCALAR - any scalar variable ** for diagnostic grids only ** \n"//& " HYCOM1 - HyCOM-like hybrid coordinate\n"//& " HYBGEN - Hybrid coordinate from the Hycom hybgen code\n"//& " ADAPTIVE - optimize for smooth neutral density surfaces" @@ -210,6 +217,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m integer :: regrid_answer_date ! The vintage of the regridding expressions to use. real :: tmpReal ! A temporary variable used in setting other variables [various] real :: P_Ref ! The coordinate variable reference pression [R L2 T-2 ~> Pa] + logical :: histogram_extensive_diags ! For diagnostic grids, extensive diags are remapped using a histogramming procedure real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z). real :: dz_extra ! The thickness of an added layer to append to the woa09_dz profile when ! maximum_depth is large [m] (not in Z). @@ -226,6 +234,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m real, dimension(:), allocatable :: dz_max ! Thicknesses used to find maximum interface depths ! [H ~> m or kg m-2] or other units real, dimension(:), allocatable :: rho_target ! Target density used in HYBRID mode [kg m-3] + real, dimension(:), allocatable :: scalar_target ! Target scalar used in SCALAR mode [kg m-3] ! Thicknesses [m] that give level centers approximately corresponding to table 2 of WOA09 ! These are approximate because the WOA09 depths are not smoothly spaced. Levels ! 1, 4, 5, 9, 12, 24, and 36 are 2.5, 2.5, 1.25 12.5, 37.5 and 62.5 m deeper than WOA09 @@ -417,6 +426,8 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m expected_units = 'nondim' ; alt_units = expected_units elseif (CS%regridding_scheme == REGRIDDING_RHO) then expected_units = 'kg m-3' ; alt_units = expected_units + elseif (CS%regridding_scheme == REGRIDDING_SCALAR) then + expected_units = 'degC' ; alt_units = expected_units else expected_units = 'meters' ; alt_units = 'm' endif @@ -432,6 +443,9 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m if (CS%regridding_scheme == REGRIDDING_RHO) then allocate(rho_target(ke+1)) call MOM_read_data(trim(fileName), trim(varName), rho_target) + elseif (CS%regridding_scheme == REGRIDDING_SCALAR) then + allocate(scalar_target(ke+1)) + call MOM_read_data(trim(fileName), trim(varName), scalar_target) else allocate(dz(ke)) allocate(z_max(ke+1)) @@ -587,7 +601,12 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m allocate( CS%coordinateResolution(CS%nk), source=-1.E30 ) if (state_dependent(CS%regridding_scheme)) then ! Target values - allocate( CS%target_density(CS%nk+1), source=-1.E30*US%kg_m3_to_R ) + !! gmac : adding if statement to allow selection of different units, may not be necessary + if (coordinateMode(coord_mode) == REGRIDDING_RHO) then + allocate( CS%target_density(CS%nk+1), source=-1.E30*US%kg_m3_to_R ) + elseif (coordinateMode(coord_mode) == REGRIDDING_SCALAR) then + allocate( CS%target_density(CS%nk+1), source=-1.E30*US%degC_to_C ) + endif endif if (allocated(dz)) then @@ -595,6 +614,8 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m call setCoordinateResolution(dz, CS, scale=1.0) elseif (coordinateMode(coord_mode) == REGRIDDING_RHO) then call setCoordinateResolution(dz, CS, scale=US%kg_m3_to_R) + elseif (coordinateMode(coord_mode) == REGRIDDING_SCALAR) then + call setCoordinateResolution(dz, CS, scale=US%degC_to_C) elseif (coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then call setCoordinateResolution(dz, CS, scale=GV%m_to_H) CS%coord_scale = GV%H_to_m @@ -604,9 +625,11 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m endif endif - ! set coord_scale for RHO regridding independent of allocation status of dz + ! set coord_scale for RHO and SCALAR regridding independent of allocation status of dz if (coordinateMode(coord_mode) == REGRIDDING_RHO) then CS%coord_scale = US%R_to_kg_m3 + elseif (coordinateMode(coord_mode) == REGRIDDING_SCALAR) then + CS%coord_scale = US%C_to_degC endif ! ensure CS%ref_pressure is rescaled properly @@ -623,6 +646,17 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m 'RHO target densities for interfaces', units=coordinateUnits(coord_mode)) endif + if (allocated(scalar_target)) then + call set_target_densities(CS, US%degC_to_C*scalar_target) + deallocate(scalar_target) + + ! \todo This line looks like it would overwrite the target densities set just above? + elseif (coordinateMode(coord_mode) == REGRIDDING_SCALAR) then + call set_target_densities_from_GV(GV, US, CS) + call log_param(param_file, mdl, "!TARGET_DENSITIES", US%C_to_degC*CS%target_density(:), & + 'SCALAR target densities for interfaces', units=coordinateUnits(coord_mode)) + endif + ! initialise coordinate-specific control structure call initCoord(CS, GV, US, coord_mode, param_file) @@ -645,7 +679,13 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m "When interpolating potential density profiles we can add "//& "some artificial compressibility solely to make homogeneous "//& "regions appear stratified.", units="nondim", default=0.) - call set_regrid_params(CS, compress_fraction=tmpReal, ref_pressure=P_Ref) + call get_param(param_file, mdl, create_coord_param(param_prefix, "HISTOGRAM_EXTENSIVE_DIAGS", param_suffix), & + histogram_extensive_diags, & + "If true, extensive diagnostics are remapped using a histogram procedure. "//& + "This is therefore suitable for coordinates that are non-monotonic "//& + "in the vertical dimension. This should only be set True for **diagnostic**"//& + "coordinates.", units="nondim", default=.false.) + call set_regrid_params(CS, compress_fraction=tmpReal, ref_pressure=P_Ref, histogram_extensive_diags=histogram_extensive_diags) endif if (main_parameters) then @@ -830,6 +870,7 @@ subroutine end_regridding(CS) if (associated(CS%zlike_CS)) call end_coord_zlike(CS%zlike_CS) if (associated(CS%sigma_CS)) call end_coord_sigma(CS%sigma_CS) if (associated(CS%rho_CS)) call end_coord_rho(CS%rho_CS) + if (associated(CS%scalar_CS)) call end_coord_scalar(CS%scalar_CS) if (associated(CS%hycom_CS)) call end_coord_hycom(CS%hycom_CS) if (associated(CS%adapt_CS)) call end_coord_adapt(CS%adapt_CS) if (associated(CS%hybgen_CS)) call end_hybgen_regrid(CS%hybgen_CS) @@ -1998,7 +2039,7 @@ function uniformResolution(nk,coordMode,maxDepth,rhoLight,rhoHeavy) REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_ADAPTIVE ) uniformResolution(:) = maxDepth / real(nk) - case ( REGRIDDING_RHO ) + case ( REGRIDDING_RHO, REGRIDDING_SCALAR ) uniformResolution(:) = (rhoHeavy - rhoLight) / real(nk) case ( REGRIDDING_SIGMA ) @@ -2031,7 +2072,9 @@ subroutine initCoord(CS, GV, US, coord_mode, param_file) case (REGRIDDING_SIGMA) call init_coord_sigma(CS%sigma_CS, CS%nk, CS%coordinateResolution) case (REGRIDDING_RHO) - call init_coord_rho(CS%rho_CS, CS%nk, CS%ref_pressure, CS%target_density, CS%interp_CS) + call init_coord_rho(CS%rho_CS, CS%nk, CS%ref_pressure, CS%target_density, CS%histogram_extensive_diags, CS%interp_CS) + case (REGRIDDING_SCALAR) + call init_coord_scalar(CS%scalar_CS, CS%nk, CS%ref_pressure, CS%target_density, CS%histogram_extensive_diags, CS%interp_CS) case (REGRIDDING_HYCOM1) call init_coord_hycom(CS%hycom_CS, CS%nk, CS%coordinateResolution, CS%target_density, & CS%interp_CS) @@ -2288,6 +2331,16 @@ function getCoordinateInterfaces( CS, undo_scaling ) call MOM_error(FATAL, 'MOM_regridding, getCoordinateInterfaces: '//& 'target densities not set!') + if (unscale) then + getCoordinateInterfaces(:) = CS%coord_scale * CS%target_density(:) + else + getCoordinateInterfaces(:) = CS%target_density(:) + endif + elseif (CS%regridding_scheme == REGRIDDING_SCALAR) then + if (.not. CS%target_density_set) & + call MOM_error(FATAL, 'MOM_regridding, getCoordinateInterfaces: '//& + 'target densities not set!') + if (unscale) then getCoordinateInterfaces(:) = CS%coord_scale * CS%target_density(:) else @@ -2330,6 +2383,8 @@ function getCoordinateUnits( CS ) getCoordinateUnits = 'fraction' case ( REGRIDDING_RHO ) getCoordinateUnits = 'kg/m3' + case ( REGRIDDING_SCALAR ) + getCoordinateUnits = 'degC' case ( REGRIDDING_ARBITRARY ) getCoordinateUnits = 'unknown' case default @@ -2356,6 +2411,8 @@ function getCoordinateShortName( CS ) getCoordinateShortName = 'sigma' case ( REGRIDDING_RHO ) getCoordinateShortName = 'rho' + case ( REGRIDDING_SCALAR ) + getCoordinateShortName = 'scalar' case ( REGRIDDING_ARBITRARY ) getCoordinateShortName = 'coordinate' case ( REGRIDDING_HYCOM1 ) @@ -2374,7 +2431,7 @@ end function getCoordinateShortName !> Can be used to set any of the parameters for MOM_regridding. subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_grid_weight, & interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, & - compress_fraction, ref_pressure, & + compress_fraction, ref_pressure, histogram_extensive_diags, & integrate_downward_for_e, remap_answers_2018, remap_answer_date, regrid_answer_date, & adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin, adaptDrho0) type(regridding_CS), intent(inout) :: CS !< Regridding control structure @@ -2388,6 +2445,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri real, optional, intent(in) :: compress_fraction !< Fraction of compressibility to add to potential density [nondim] real, optional, intent(in) :: ref_pressure !< The reference pressure for density-dependent !! coordinates [R L2 T-2 ~> Pa] + logical, optional, intent(in) :: histogram_extensive_diags !< If true, use histogrammign procedure for remapping extensive diagnostics logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface positions downward !! from the top. logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions @@ -2426,6 +2484,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri if (present(compress_fraction)) CS%compressibility_fraction = compress_fraction if (present(ref_pressure)) CS%ref_pressure = ref_pressure if (present(integrate_downward_for_e)) CS%integrate_downward_for_e = integrate_downward_for_e + if (present(histogram_extensive_diags)) CS%histogram_extensive_diags = histogram_extensive_diags if (present(remap_answers_2018)) then if (remap_answers_2018) then CS%remap_answer_date = 20181231 @@ -2445,10 +2504,19 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri case (REGRIDDING_RHO) if (present(min_thickness)) call set_rho_params(CS%rho_CS, min_thickness=min_thickness) if (present(ref_pressure)) call set_rho_params(CS%rho_CS, ref_pressure=ref_pressure) + if (present(histogram_extensive_diags)) call set_rho_params(CS%rho_CS, histogram_extensive_diags=histogram_extensive_diags) if (present(integrate_downward_for_e)) & call set_rho_params(CS%rho_CS, integrate_downward_for_e=integrate_downward_for_e) if (associated(CS%rho_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) & call set_rho_params(CS%rho_CS, interp_CS=CS%interp_CS) + case (REGRIDDING_SCALAR) + if (present(min_thickness)) call set_scalar_params(CS%scalar_CS, min_thickness=min_thickness) + if (present(ref_pressure)) call set_scalar_params(CS%scalar_CS, ref_pressure=ref_pressure) + if (present(histogram_extensive_diags)) call set_scalar_params(CS%scalar_CS, histogram_extensive_diags=histogram_extensive_diags) + if (present(integrate_downward_for_e)) & + call set_scalar_params(CS%scalar_CS, integrate_downward_for_e=integrate_downward_for_e) + if (associated(CS%scalar_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) & + call set_scalar_params(CS%scalar_CS, interp_CS=CS%interp_CS) case (REGRIDDING_HYCOM1) if (associated(CS%hycom_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) & call set_hycom_params(CS%hycom_CS, interp_CS=CS%interp_CS) @@ -2498,6 +2566,14 @@ function get_rho_CS(CS) get_rho_CS = CS%rho_CS end function get_rho_CS +!> This returns a copy of the scalar_CS stored in the regridding control structure. +function get_scalar_CS(CS) + type(regridding_CS), intent(in) :: CS !< Regridding control structure + type(scalar_CS) :: get_scalar_CS + + get_scalar_CS = CS%scalar_CS +end function get_scalar_CS + !------------------------------------------------------------------------------ !> Return coordinate-derived thicknesses for fixed coordinate systems function getStaticThickness( CS, SSH, depth ) @@ -2626,6 +2702,14 @@ integer function rho_function1( string, rho_target ) end function rho_function1 +subroutine check_if_histogram_extensive_diags(CS,histogram_extensive_diags) + type(regridding_CS), intent(in) :: CS + logical, intent(inout) :: histogram_extensive_diags + + histogram_extensive_diags = CS%histogram_extensive_diags + +end subroutine check_if_histogram_extensive_diags + !> \namespace mom_regridding !! !! A vertical grid is defined solely by the cell thicknesses, \f$h\f$. diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 7257319edb..59898d98da 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -84,6 +84,7 @@ module MOM_remapping public initialize_remapping, end_remapping, remapping_set_param, extract_member_remapping_CS public remapping_unit_tests, build_reconstructions_1d, average_value_ppoly public interpolate_column, reintegrate_column, dzFromH1H2 +public histogram_column ! The following are private parameter constants integer, parameter :: REMAPPING_PCM = 0 !< O(h^1) remapping scheme @@ -1385,6 +1386,32 @@ subroutine reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, uh_dest) end subroutine reintegrate_column +!> Conservatively histogram data, uh_dest, weights mapping source grid to target grid +subroutine histogram_column(nsrc, uh_src, ndest, uh_dest, histogram_weights) + integer, intent(in) :: nsrc !< Number of source cells + real, dimension(nsrc), intent(in) :: uh_src !< Values at source cell interfaces [A H] + integer, intent(in) :: ndest !< Number of destination cells + real, dimension(ndest), intent(inout) :: uh_dest !< Interpolated value at destination cell interfaces [A H] + real, dimension(nsrc,ndest), intent(in) :: histogram_weights !< Weights mapping source to destination grid + + ! Local variables + integer :: k_dest, k_src + real, dimension(nsrc) :: weights + real :: uh_dest_sum + + uh_dest(:) = 0.0 + do k_dest = 1,ndest ! Loop through destination grid layers + weights = histogram_weights(:,k_dest) + uh_dest_sum = 0.0 + do k_src = 1,nsrc ! Loop through source grid + ! cumulatively sum field multiplied by the weight associated with each grid + uh_dest_sum = uh_dest_sum + uh_src(k_src)*weights(k_src) + enddo + uh_dest(k_dest) = uh_dest_sum + enddo + +end subroutine histogram_column + !> Returns the average value of a reconstruction within a single source cell, i0, !! between the non-dimensional positions xa and xb (xa<=xb) with dimensional !! separation dh. diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90 index c967687dc8..9184c3f477 100644 --- a/src/ALE/coord_rho.F90 +++ b/src/ALE/coord_rho.F90 @@ -7,6 +7,7 @@ module coord_rho use MOM_remapping, only : remapping_CS, remapping_core_h use MOM_EOS, only : EOS_type, calculate_density use regrid_interp, only : interp_CS_type, build_and_interpolate_grid, DEGREE_MAX +use regrid_interp, only : build_histogram_weights implicit none ; private @@ -29,6 +30,9 @@ module coord_rho !> Nominal density of interfaces [R ~> kg m-3] real, allocatable, dimension(:) :: target_density + !> If true, use a histogramming procedure for extensive diagnostics (rather than the reintegrate procedure) + logical :: histogram_extensive_diags = .false. + !> Interpolation control structure type(interp_CS_type) :: interp_CS end type rho_CS @@ -38,11 +42,12 @@ module coord_rho contains !> Initialise a rho_CS with pointers to parameters -subroutine init_coord_rho(CS, nk, ref_pressure, target_density, interp_CS) +subroutine init_coord_rho(CS, nk, ref_pressure, target_density, histogram_extensive_diags, interp_CS) type(rho_CS), pointer :: CS !< Unassociated pointer to hold the control structure integer, intent(in) :: nk !< Number of layers in the grid real, intent(in) :: ref_pressure !< Coordinate reference pressure [R L2 T-2 ~> Pa] real, dimension(:), intent(in) :: target_density !< Nominal density of interfaces [R ~> kg m-3] + logical, intent(in) :: histogram_extensive_diags !< Boolean to select how to deal with extensive diagnostics type(interp_CS_type), intent(in) :: interp_CS !< Controls for interpolation if (associated(CS)) call MOM_error(FATAL, "init_coord_rho: CS already associated!") @@ -52,6 +57,7 @@ subroutine init_coord_rho(CS, nk, ref_pressure, target_density, interp_CS) CS%nk = nk CS%ref_pressure = ref_pressure CS%target_density(:) = target_density(:) + CS%histogram_extensive_diags = histogram_extensive_diags CS%interp_CS = interp_CS end subroutine init_coord_rho @@ -67,7 +73,7 @@ subroutine end_coord_rho(CS) end subroutine end_coord_rho !> This subroutine can be used to set the parameters for the coord_rho module -subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, interp_CS, ref_pressure) +subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, histogram_extensive_diags, interp_CS, ref_pressure) type(rho_CS), pointer :: CS !< Coordinate control structure real, optional, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2] logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface @@ -75,6 +81,7 @@ subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, interp_CS !! from the bottom upward, as does the rest of the model. real, optional, intent(in) :: ref_pressure !< The reference pressure for density-dependent !! coordinates [R L2 T-2 ~> Pa] + logical, optional, intent(in) :: histogram_extensive_diags !< If true, use histogram approach to for outputing extensive diags type(interp_CS_type), optional, intent(in) :: interp_CS !< Controls for interpolation @@ -82,6 +89,7 @@ subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, interp_CS if (present(min_thickness)) CS%min_thickness = min_thickness if (present(integrate_downward_for_e)) CS%integrate_downward_for_e = integrate_downward_for_e + if (present(histogram_extensive_diags)) CS%histogram_extensive_diags = histogram_extensive_diags if (present(interp_CS)) CS%interp_CS = interp_CS if (present(ref_pressure)) CS%ref_pressure = ref_pressure end subroutine set_rho_params @@ -91,7 +99,7 @@ end subroutine set_rho_params !! 1. Density profiles are calculated on the source grid. !! 2. Positions of target densities (for interfaces) are found by interpolation. subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & - z_rigid_top, eta_orig, h_neglect, h_neglect_edge) + histogram_weights, z_rigid_top, eta_orig, h_neglect, h_neglect_edge) type(rho_CS), intent(in) :: CS !< coord_rho control structure integer, intent(in) :: nz !< Number of levels on source grid (i.e. length of h, T, S) real, intent(in) :: depth !< Depth of ocean bottom (positive downward) [H ~> m or kg m-2] @@ -109,6 +117,7 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & !! of cell reconstructions [H ~> m or kg m-2] real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose !! of edge value calculations [H ~> m or kg m-2] + real, optional, dimension(nz,CS%nk), intent (inout) :: histogram_weights ! Matrix of weights mapping source grid cells to target grid layers ! Local variables integer :: k, count_nonzero_layers @@ -141,6 +150,12 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & h_nv, xTmp, CS%target_density, CS%nk, h_new, & x1, h_neglect, h_neglect_edge) + if (CS%histogram_extensive_diags) then + ! Based on source column density profile, derive weights that map source to target grid + call build_histogram_weights(CS%interp_CS, densities, nz, h, CS%target_density, CS%nk, & + histogram_weights, h_neglect, h_neglect_edge) + endif + ! Inflate vanished layers call old_inflate_layers_1d(CS%min_thickness, CS%nk, h_new) diff --git a/src/ALE/coord_scalar.F90 b/src/ALE/coord_scalar.F90 new file mode 100644 index 0000000000..0576338613 --- /dev/null +++ b/src/ALE/coord_scalar.F90 @@ -0,0 +1,437 @@ +!> Regrid columns for a continuous scalar coordinate +module coord_scalar + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_error_handler, only : MOM_error, FATAL +use MOM_remapping, only : remapping_CS, remapping_core_h +use MOM_EOS, only : EOS_type, calculate_density +use regrid_interp, only : interp_CS_type, build_and_interpolate_grid, DEGREE_MAX +use regrid_interp, only : build_histogram_weights + +implicit none ; private + +!> Control structure containing required parameters for the scalar coordinate +type, public :: scalar_CS ; private + + !> Number of layers + integer :: nk + + !> Minimum thickness allowed for layers, often in [H ~> m or kg m-2] + real :: min_thickness = 0. + + !> Reference pressure for density calculations [R L2 T-2 ~> Pa] + real :: ref_pressure + + !> If true, integrate for interface positions from the top downward. + !! If false, integrate from the bottom upward, as does the rest of the model. + logical :: integrate_downward_for_e = .false. + + !> Nominal density of interfaces [R ~> kg m-3] + real, allocatable, dimension(:) :: target_density + + !> If true, use a histogramming procedure for extensive diagnostics (rather than the reintegrate procedure) + logical :: histogram_extensive_diags = .false. + + !> Interpolation control structure + type(interp_CS_type) :: interp_CS +end type scalar_CS + +public init_coord_scalar, set_scalar_params, build_scalar_column, old_inflate_layers_1d, end_coord_scalar + +contains + +!> Initialise a scalar_CS with pointers to parameters +subroutine init_coord_scalar(CS, nk, ref_pressure, target_density, histogram_extensive_diags, interp_CS) + type(scalar_CS), pointer :: CS !< Unassociated pointer to hold the control structure + integer, intent(in) :: nk !< Number of layers in the grid + real, intent(in) :: ref_pressure !< Coordinate reference pressure [R L2 T-2 ~> Pa] + real, dimension(:), intent(in) :: target_density !< Nominal density of interfaces [R ~> kg m-3] + logical, intent(in) :: histogram_extensive_diags !< Boolean to select how to deal with extensive diagnostics + type(interp_CS_type), intent(in) :: interp_CS !< Controls for interpolation + + if (associated(CS)) call MOM_error(FATAL, "init_coord_scalar: CS already associated!") + allocate(CS) + allocate(CS%target_density(nk+1)) + + CS%nk = nk + CS%ref_pressure = ref_pressure + CS%target_density(:) = target_density(:) + CS%histogram_extensive_diags = histogram_extensive_diags + CS%interp_CS = interp_CS + +end subroutine init_coord_scalar + +!> This subroutine deallocates memory in the control structure for the coord_scalar module +subroutine end_coord_scalar(CS) + type(scalar_CS), pointer :: CS !< Coordinate control structure + + ! nothing to do + if (.not. associated(CS)) return + deallocate(CS%target_density) + deallocate(CS) +end subroutine end_coord_scalar + +!> This subroutine can be used to set the parameters for the coord_scalar module +subroutine set_scalar_params(CS, min_thickness, integrate_downward_for_e, histogram_extensive_diags, interp_CS, ref_pressure) + type(scalar_CS), pointer :: CS !< Coordinate control structure + real, optional, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2] + logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface + !! positions from the top downward. If false, integrate + !! from the bottom upward, as does the rest of the model. + real, optional, intent(in) :: ref_pressure !< The reference pressure for density-dependent + !! coordinates [R L2 T-2 ~> Pa] + logical, optional, intent(in) :: histogram_extensive_diags !< If true, use histogram approach to for outputing extensive diags + + type(interp_CS_type), optional, intent(in) :: interp_CS !< Controls for interpolation + + if (.not. associated(CS)) call MOM_error(FATAL, "set_scalar_params: CS not associated") + + if (present(min_thickness)) CS%min_thickness = min_thickness + if (present(integrate_downward_for_e)) CS%integrate_downward_for_e = integrate_downward_for_e + if (present(histogram_extensive_diags)) CS%histogram_extensive_diags = histogram_extensive_diags + if (present(interp_CS)) CS%interp_CS = interp_CS + if (present(ref_pressure)) CS%ref_pressure = ref_pressure +end subroutine set_scalar_params + +!> Build a scalar coordinate column +!! +!! 1. Density profiles are calculated on the source grid. +!! 2. Positions of target densities (for interfaces) are found by interpolation. +subroutine build_scalar_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & + histogram_weights, z_rigid_top, eta_orig, h_neglect, h_neglect_edge) + type(scalar_CS), intent(in) :: CS !< coord_scalar control structure + integer, intent(in) :: nz !< Number of levels on source grid (i.e. length of h, T, S) + real, intent(in) :: depth !< Depth of ocean bottom (positive downward) [H ~> m or kg m-2] + real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(nz), intent(in) :: T !< Temperature for source column [C ~> degC] + real, dimension(nz), intent(in) :: S !< Salinity for source column [S ~> ppt] + type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure + real, dimension(CS%nk+1), & + intent(inout) :: z_interface !< Absolute positions of interfaces [H ~> m or kg m-2] + real, optional, intent(in) :: z_rigid_top !< The height of a rigid top (positive upward in the same + !! units as depth) [H ~> m or kg m-2] + real, optional, intent(in) :: eta_orig !< The actual original height of the top in the same + !! units as depth) [H ~> m or kg m-2] + real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose + !! of cell reconstructions [H ~> m or kg m-2] + real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose + !! of edge value calculations [H ~> m or kg m-2] + real, optional, dimension(nz,CS%nk), intent (inout) :: histogram_weights ! Matrix of weights mapping source grid cells to target grid layers + + ! Local variables + integer :: k, count_nonzero_layers + integer, dimension(nz) :: mapping + real, dimension(nz) :: pres ! Pressures used to calculate density [R L2 T-2 ~> Pa] + real, dimension(nz) :: h_nv ! Thicknesses of non-vanishing layers [H ~> m or kg m-2] + real, dimension(nz) :: densities ! Layer density [R ~> kg m-3] + real, dimension(nz+1) :: xTmp ! Temporary positions [H ~> m or kg m-2] + real, dimension(CS%nk) :: h_new ! New thicknesses [H ~> m or kg m-2] + real, dimension(CS%nk+1) :: x1 ! Interface heights [H ~> m or kg m-2] + + ! Construct source column with vanished layers removed (stored in h_nv) + call copy_finite_thicknesses(nz, h, CS%min_thickness, count_nonzero_layers, h_nv, mapping) + + if (count_nonzero_layers > 1) then + xTmp(1) = 0.0 + do k = 1,count_nonzero_layers + xTmp(k+1) = xTmp(k) + h_nv(k) + enddo + + ! Compute densities on source column + pres(:) = CS%ref_pressure + ! call calculate_density(T, S, pres, densities, eqn_of_state) + do k = 1,count_nonzero_layers + densities(k) = T(mapping(k)) + enddo + + ! Based on source column density profile, interpolate to generate a new grid + call build_and_interpolate_grid(CS%interp_CS, densities, count_nonzero_layers, & + h_nv, xTmp, CS%target_density, CS%nk, h_new, & + x1, h_neglect, h_neglect_edge) + + if (CS%histogram_extensive_diags) then + ! Based on source column density profile, derive weights that map source to target grid + call build_histogram_weights(CS%interp_CS, densities, nz, h, CS%target_density, CS%nk, & + histogram_weights, h_neglect, h_neglect_edge) + endif + + ! Inflate vanished layers + call old_inflate_layers_1d(CS%min_thickness, CS%nk, h_new) + + ! Comment: The following adjustment of h_new, and re-calculation of h_new via x1 needs to be removed + x1(1) = 0.0 ; do k = 1,CS%nk ; x1(k+1) = x1(k) + h_new(k) ; enddo + do k = 1,CS%nk + h_new(k) = x1(k+1) - x1(k) + enddo + + else ! count_nonzero_layers <= 1 + if (nz == CS%nk) then + h_new(:) = h(:) ! This keeps old behavior + else + h_new(:) = 0. + h_new(1) = h(1) + endif + endif + + ! Return interface positions + if (CS%integrate_downward_for_e) then + ! Remapping is defined integrating from zero + z_interface(1) = 0. + do k = 1,CS%nk + z_interface(k+1) = z_interface(k) - h_new(k) + enddo + else + ! The rest of the model defines grids integrating up from the bottom + z_interface(CS%nk+1) = -depth + do k = CS%nk,1,-1 + z_interface(k) = z_interface(k+1) + h_new(k) + enddo + endif + +end subroutine build_scalar_column + +!### build_scalar_column_iteratively is never used or called. + +!> Iteratively build a scalar coordinate column +!! +!! The algorithm operates as follows within each column: +!! +!! 1. Given T & S within each layer, the layer densities are computed. +!! 2. Based on these layer densities, a global density profile is reconstructed +!! (this profile is monotonically increasing and may be discontinuous) +!! 3. The new grid interfaces are determined based on the target interface +!! densities. +!! 4. T & S are remapped onto the new grid. +!! 5. Return to step 1 until convergence or until the maximum number of +!! iterations is reached, whichever comes first. +subroutine build_scalar_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_state, & + zInterface, h_neglect, h_neglect_edge, dev_tol) + type(scalar_CS), intent(in) :: CS !< Regridding control structure + type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options + integer, intent(in) :: nz !< Number of levels + real, intent(in) :: depth !< Depth of ocean bottom [Z ~> m] + real, dimension(nz), intent(in) :: h !< Layer thicknesses in Z coordinates [Z ~> m] + real, dimension(nz), intent(in) :: T !< T for column [C ~> degC] + real, dimension(nz), intent(in) :: S !< S for column [S ~> ppt] + type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure + real, dimension(nz+1), intent(inout) :: zInterface !< Absolute positions of interfaces [Z ~> m] + real, optional, intent(in) :: h_neglect !< A negligibly small width for the + !! purpose of cell reconstructions + !! in the same units as h [Z ~> m] + real, optional, intent(in) :: h_neglect_edge !< A negligibly small width + !! for the purpose of edge value calculations + !! in the same units as h [Z ~> m] + real, optional, intent(in) :: dev_tol !< The tolerance for the deviation between + !! successive grids for determining when the + !! iterative solver has converged [Z ~> m] + + ! Local variables + real, dimension(nz+1) :: x0, x1, xTmp ! Temporary interface heights [Z ~> m] + real, dimension(nz) :: pres ! The pressure used in the equation of state [R L2 T-2 ~> Pa]. + real, dimension(nz) :: densities ! Layer densities [R ~> kg m-3] + real, dimension(nz) :: T_tmp, S_tmp ! A temporary profile of temperature [C ~> degC] and salinity [S ~> ppt]. + real, dimension(nz) :: h0, h1, hTmp ! Temporary thicknesses [Z ~> m] + real :: deviation ! When iterating to determine the final grid, this is the + ! deviation between two successive grids [Z ~> m]. + real :: deviation_tol ! Deviation tolerance between succesive grids in + ! regridding iterations [Z ~> m] + real :: threshold ! The minimum thickness for a layer to be considered to exist [Z ~> m] + integer, dimension(nz) :: mapping ! The indices of the massive layers in the initial column. + integer :: k, m, count_nonzero_layers + + ! Maximum number of regridding iterations + integer, parameter :: NB_REGRIDDING_ITERATIONS = 1 + + threshold = CS%min_thickness + pres(:) = CS%ref_pressure + T_tmp(:) = T(:) + S_tmp(:) = S(:) + h0(:) = h(:) + + ! Start iterations to build grid + m = 1 + deviation_tol = 1.0e-15*depth ; if (present(dev_tol)) deviation_tol = dev_tol + + do m=1,NB_REGRIDDING_ITERATIONS + + ! Construct column with vanished layers removed + call copy_finite_thicknesses(nz, h0, threshold, count_nonzero_layers, hTmp, mapping) + if ( count_nonzero_layers <= 1 ) then + h1(:) = h0(:) + exit ! stop iterations here + endif + + xTmp(1) = 0.0 + do k = 1,count_nonzero_layers + xTmp(k+1) = xTmp(k) + hTmp(k) + enddo + + ! Compute densities within current water column + call calculate_density(T_tmp, S_tmp, pres, densities, eqn_of_state) + + do k = 1,count_nonzero_layers + densities(k) = densities(mapping(k)) + enddo + + ! One regridding iteration + ! Based on global density profile, interpolate to generate a new grid + call build_and_interpolate_grid(CS%interp_CS, densities, count_nonzero_layers, & + hTmp, xTmp, CS%target_density, nz, h1, x1, h_neglect, h_neglect_edge) + + call old_inflate_layers_1d( CS%min_thickness, nz, h1 ) + x1(1) = 0.0 ; do k = 1,nz ; x1(k+1) = x1(k) + h1(k) ; enddo + + ! Remap T and S from previous grid to new grid + do k = 1,nz + h1(k) = x1(k+1) - x1(k) + enddo + + call remapping_core_h(remapCS, nz, h0, S, nz, h1, S_tmp) + + call remapping_core_h(remapCS, nz, h0, T, nz, h1, T_tmp) + + ! Compute the deviation between two successive grids + deviation = 0.0 + x0(1) = 0.0 + x1(1) = 0.0 + do k = 2,nz + x0(k) = x0(k-1) + h0(k-1) + x1(k) = x1(k-1) + h1(k-1) + deviation = deviation + (x0(k)-x1(k))**2 + enddo + deviation = sqrt( deviation / (nz-1) ) + + if ( deviation <= deviation_tol ) exit + + ! Copy final grid onto start grid for next iteration + h0(:) = h1(:) + enddo ! end regridding iterations + + if (CS%integrate_downward_for_e) then + zInterface(1) = 0. + do k = 1,nz + zInterface(k+1) = zInterface(k) - h1(k) + ! Adjust interface position to accommodate inflating layers + ! without disturbing the interface above + enddo + else + ! The rest of the model defines grids integrating up from the bottom + zInterface(nz+1) = -depth + do k = nz,1,-1 + zInterface(k) = zInterface(k+1) + h1(k) + ! Adjust interface position to accommodate inflating layers + ! without disturbing the interface above + enddo + endif + +end subroutine build_scalar_column_iteratively + +!> Copy column thicknesses with vanished layers removed +subroutine copy_finite_thicknesses(nk, h_in, thresh, nout, h_out, mapping) + integer, intent(in) :: nk !< Number of layer for h_in, T_in, S_in + real, dimension(nk), intent(in) :: h_in !< Thickness of input column [H ~> m or kg m-2] or [Z ~> m] + real, intent(in) :: thresh !< Thickness threshold defining vanished + !! layers [H ~> m or kg m-2] or [Z ~> m] + integer, intent(out) :: nout !< Number of non-vanished layers + real, dimension(nk), intent(out) :: h_out !< Thickness of output column [H ~> m or kg m-2] or [Z ~> m] + integer, dimension(nk), intent(out) :: mapping !< Index of k-out corresponding to k-in + ! Local variables + integer :: k, k_thickest + real :: thickness_in_vanished ! Summed thicknesses in discarded layers [H ~> m or kg m-2] or [Z ~> m] + real :: thickest_h_out ! Thickness of the thickest layer [H ~> m or kg m-2] or [Z ~> m] + + ! Build up new grid + nout = 0 + thickness_in_vanished = 0.0 + thickest_h_out = h_in(1) + k_thickest = 1 + do k = 1, nk + mapping(k) = nout ! Note k>=nout always + h_out(k) = 0. ! Make sure h_out is set everywhere + if (h_in(k) > thresh) then + ! For non-vanished layers + nout = nout + 1 + mapping(nout) = k + h_out(nout) = h_in(k) + if (h_out(nout) > thickest_h_out) then + thickest_h_out = h_out(nout) + k_thickest = nout + endif + else + ! Add up mass in vanished layers + thickness_in_vanished = thickness_in_vanished + h_in(k) + endif + enddo + + ! No finite layers + if (nout <= 1) return + + ! Adjust for any lost volume in vanished layers + h_out(k_thickest) = h_out(k_thickest) + thickness_in_vanished + +end subroutine copy_finite_thicknesses + +!------------------------------------------------------------------------------ +!> Inflate vanished layers to finite (nonzero) width +subroutine old_inflate_layers_1d( min_thickness, nk, h ) + + ! Argument + real, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2] or other units + integer, intent(in) :: nk !< Number of layers in the grid + real, dimension(:), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2] or other units + + ! Local variable + integer :: k + integer :: k_found + integer :: count_nonzero_layers + real :: delta ! An increase to a layer to increase it to the minimum thickness in the + ! same units as h, often [H ~> m or kg m-2] + real :: correction ! The accumulated correction that will be applied to the thickest layer + ! to give mass conservation in the same units as h, often [H ~> m or kg m-2] + real :: maxThickness ! The thickness of the thickest layer in the same units as h, often [H ~> m or kg m-2] + + ! Count number of nonzero layers + count_nonzero_layers = 0 + do k = 1,nk + if ( h(k) > min_thickness ) then + count_nonzero_layers = count_nonzero_layers + 1 + endif + enddo + + ! If all layer thicknesses are greater than the threshold, exit routine + if ( count_nonzero_layers == nk ) return + + ! If all thicknesses are zero, inflate them all and exit + if ( count_nonzero_layers == 0 ) then + do k = 1,nk + h(k) = min_thickness + enddo + return + endif + + ! Inflate zero layers + correction = 0.0 + do k = 1,nk + if ( h(k) <= min_thickness ) then + delta = min_thickness - h(k) + correction = correction + delta + h(k) = h(k) + delta + endif + enddo + + ! Modify thicknesses of nonzero layers to ensure volume conservation + maxThickness = h(1) + k_found = 1 + do k = 1,nk + if ( h(k) > maxThickness ) then + maxThickness = h(k) + k_found = k + endif + enddo + + h(k_found) = h(k_found) - correction + +end subroutine old_inflate_layers_1d + +end module coord_scalar diff --git a/src/ALE/regrid_consts.F90 b/src/ALE/regrid_consts.F90 index 0c5ccf268f..7c42d1f8b9 100644 --- a/src/ALE/regrid_consts.F90 +++ b/src/ALE/regrid_consts.F90 @@ -20,11 +20,13 @@ module regrid_consts !! sigma-near the top integer, parameter :: REGRIDDING_ADAPTIVE = 9 !< Adaptive coordinate mode identifier integer, parameter :: REGRIDDING_HYBGEN = 10 !< Hybgen coordinates identifier +integer, parameter :: REGRIDDING_SCALAR = 11 !< Scalar coordinates identifier character(len=*), parameter :: REGRIDDING_LAYER_STRING = "LAYER" !< Layer string character(len=*), parameter :: REGRIDDING_ZSTAR_STRING_OLD = "Z*" !< z* string (legacy name) character(len=*), parameter :: REGRIDDING_ZSTAR_STRING = "ZSTAR" !< z* string character(len=*), parameter :: REGRIDDING_RHO_STRING = "RHO" !< Rho string +character(len=*), parameter :: REGRIDDING_SCALAR_STRING = "SCALAR" !< Scalar string character(len=*), parameter :: REGRIDDING_SIGMA_STRING = "SIGMA" !< Sigma string character(len=*), parameter :: REGRIDDING_ARBITRARY_STRING = "ARB" !< Arbitrary coordinates character(len=*), parameter :: REGRIDDING_HYCOM1_STRING = "HYCOM1" !< Hycom string @@ -57,6 +59,7 @@ function coordinateMode(string) case (trim(REGRIDDING_ZSTAR_STRING)); coordinateMode = REGRIDDING_ZSTAR case (trim(REGRIDDING_ZSTAR_STRING_OLD)); coordinateMode = REGRIDDING_ZSTAR case (trim(REGRIDDING_RHO_STRING)); coordinateMode = REGRIDDING_RHO + case (trim(REGRIDDING_SCALAR_STRING)); coordinateMode = REGRIDDING_SCALAR case (trim(REGRIDDING_SIGMA_STRING)); coordinateMode = REGRIDDING_SIGMA case (trim(REGRIDDING_HYCOM1_STRING)); coordinateMode = REGRIDDING_HYCOM1 case (trim(REGRIDDING_HYBGEN_STRING)); coordinateMode = REGRIDDING_HYBGEN @@ -78,6 +81,7 @@ function coordinateUnitsI(coordMode) case (REGRIDDING_ZSTAR); coordinateUnitsI = "m" case (REGRIDDING_SIGMA_SHELF_ZSTAR); coordinateUnitsI = "m" case (REGRIDDING_RHO); coordinateUnitsI = "kg m^-3" + case (REGRIDDING_SCALAR); coordinateUnitsI = "degC" case (REGRIDDING_SIGMA); coordinateUnitsI = "Non-dimensional" case (REGRIDDING_HYCOM1); coordinateUnitsI = "m" case (REGRIDDING_HYBGEN); coordinateUnitsI = "m" @@ -113,6 +117,7 @@ logical function state_dependent_int(mode) case (REGRIDDING_ZSTAR); state_dependent_int = .false. case (REGRIDDING_SIGMA_SHELF_ZSTAR); state_dependent_int = .false. case (REGRIDDING_RHO); state_dependent_int = .true. + case (REGRIDDING_SCALAR); state_dependent_int = .true. case (REGRIDDING_SIGMA); state_dependent_int = .false. case (REGRIDDING_HYCOM1); state_dependent_int = .true. case (REGRIDDING_HYBGEN); state_dependent_int = .true. diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index 6e0be9ebba..77c8090e99 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -39,6 +39,7 @@ module regrid_interp public regridding_set_ppolys, build_and_interpolate_grid public set_interp_scheme, set_interp_extrap, set_interp_answer_date +public get_histogram_weights, build_histogram_weights ! List of interpolation schemes integer, parameter :: INTERPOLATION_P1M_H2 = 0 !< O(h^2) @@ -357,6 +358,36 @@ subroutine build_and_interpolate_grid(CS, densities, n0, h0, x0, target_values, n1, h1, x1, answer_date=CS%answer_date) end subroutine build_and_interpolate_grid +!> Build array of weights mapping source grid to target grid +subroutine build_histogram_weights(CS, densities, n0, h0, target_values, & + n1, histogram_weights, h_neglect, h_neglect_edge) +type(interp_CS_type), intent(in) :: CS !< A control structure for regrid_interp +integer, intent(in) :: n0 !< The number of points on the input grid +integer, intent(in) :: n1 !< The number of points on the output grid +real, dimension(n0), intent(in) :: h0 !< Initial cell widths [H] +real, dimension(n0), intent(in) :: densities !< Input cell densities [R ~> kg m-3] +real, dimension(n1+1), intent(in) :: target_values !< Target values of interfaces [R ~> kg m-3] +real, dimension(n0,n1), intent(inout):: histogram_weights !< Matrix of weights mapping source grid cells + !! to target grid layers +real, optional, intent(in) :: h_neglect !< A negligibly small width for the + !! purpose of cell reconstructions [H] + !! in the same units as h0. +real, optional, intent(in) :: h_neglect_edge !< A negligibly small width + !! for the purpose of edge value calculations [H] + !! in the same units as h0. + +real, dimension(n0,2) :: ppoly0_E ! Polynomial edge values [R ~> kg m-3] +real, dimension(n0,2) :: ppoly0_S ! Polynomial edge slopes [R H-1] +real, dimension(n0,DEGREE_MAX+1) :: ppoly0_C ! Polynomial interpolant coeficients on the local 0-1 grid [R ~> kg m-3] +integer :: degree + +call regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, ppoly0_C, & +degree, h_neglect, h_neglect_edge) +call get_histogram_weights(n0, ppoly0_E, ppoly0_C, target_values, degree, & +n1, histogram_weights, answer_date=CS%answer_date) + +end subroutine build_histogram_weights + !> Given a target value, find corresponding coordinate for given polynomial !! !! Here, 'ppoly' is assumed to be a piecewise discontinuous polynomial of degree @@ -507,6 +538,290 @@ function get_polynomial_coordinate( N, h, x_g, edge_values, ppoly_coefs, & x_tgt = x_g(k_found) + xi0 * h(k_found) end function get_polynomial_coordinate +! For a given target value, find all locations of its interface in the column. +! Return an array containing the index of the interface in each cell. +! If the interface index is -1, the target value is less than the value in the cell +! If the interface index is 2, the target value is greater than the values in the cell +! If the interface index is 0, the target value falls between the edges of the cell obove and this one +! If the interface value is between 0 and 1, this is the normalized position of the interface in the cell, +! determined via NR solver. +function get_interface_indices( N, edge_values, ppoly_coefs, & + target_value, degree, answer_date ) result ( interfaces ) + + ! Arguments + integer, intent(in) :: N !< Number of grid cells + real, dimension(N,2), intent(in) :: edge_values !< Edge values of interpolating polynomials [A] + real, dimension(N,DEGREE_MAX+1), intent(in) :: ppoly_coefs !< Coefficients of interpolating polynomials [A] + real, intent(in) :: target_value !< Target value to find position for [A] + integer, intent(in) :: degree !< Degree of the interpolating polynomials + integer, intent(in) :: answer_date !< The vintage of the expressions to use + real, dimension(N) :: interfaces !< An array indicating location of target interfaces in relation to each cell + + ! Local variables + real :: xi0 ! normalized target coordinate [nondim] + real, dimension(DEGREE_MAX) :: a ! polynomial coefficients [A] + real :: numerator ! The numerator of an expression [A] + real :: denominator ! The denominator of an expression [A] + real :: delta ! Newton-Raphson increment [nondim] +! real :: x ! global target coordinate [nondim] + real :: eps ! offset used to get away from boundaries [nondim] + real :: grad ! gradient during N-R iterations [A] + integer :: i, k, iter ! loop indices + integer :: k_found ! index of target cell + character(len=320) :: mesg + logical :: use_2018_answers ! If true use older, less accurate expressions. + + eps = NR_OFFSET + k_found = -1 + use_2018_answers = (answer_date < 20190101) + + !!! THIS IS A HACK AND COULD BE DONE BETTER + !!! but because of discontinuous values at cell boundaries, I need to do the first cell here separately + k = 1 + if ( target_value <= edge_values(k,1) ) then + interfaces(k) = -1 + elseif (target_value >= edge_values(k,2) ) then + interfaces(k) = 2 + else + ! Interface is within cell, so use Newton-Raphson iterations to find it + + ! Reset all polynomial coefficients to 0 and copy those pertaining to + ! the found cell + a(:) = 0.0 + do i = 1,degree+1 + a(i) = ppoly_coefs(k,i) + enddo + + ! Guess the middle of the cell to start Newton-Raphson iterations + xi0 = 0.5 + + ! Newton-Raphson iterations + do iter = 1,NR_ITERATIONS + + if (use_2018_answers) then + numerator = a(1) + a(2)*xi0 + a(3)*xi0*xi0 + a(4)*xi0*xi0*xi0 + & + a(5)*xi0*xi0*xi0*xi0 - target_value + denominator = a(2) + 2*a(3)*xi0 + 3*a(4)*xi0*xi0 + 4*a(5)*xi0*xi0*xi0 + else ! These expressions are mathematicaly equivalent but more accurate. + numerator = (a(1) - target_value) + xi0*(a(2) + xi0*(a(3) + xi0*(a(4) + a(5)*xi0))) + denominator = a(2) + xi0*(2.*a(3) + xi0*(3.*a(4) + 4.*a(5)*xi0)) + endif + + delta = -numerator / denominator + + xi0 = xi0 + delta + + ! Check whether new estimate is out of bounds. If the new estimate is + ! indeed out of bounds, we manually set it to be equal to the overtaken + ! bound with a small offset towards the interior when the gradient of + ! the function at the boundary is zero (in which case, the Newton-Raphson + ! algorithm does not converge). + if ( xi0 < 0.0 ) then + xi0 = 0.0 + grad = a(2) + if ( grad == 0.0 ) xi0 = xi0 + eps + endif + + if ( xi0 > 1.0 ) then + xi0 = 1.0 + if (use_2018_answers) then + grad = a(2) + 2*a(3) + 3*a(4) + 4*a(5) + else ! These expressions are mathematicaly equivalent but more accurate. + grad = a(2) + (2.*a(3) + (3.*a(4) + 4.*a(5))) + endif + if ( grad == 0.0 ) xi0 = xi0 - eps + endif + + ! break if converged or too many iterations taken + if ( abs(delta) < NR_TOLERANCE ) exit + enddo ! end Newton-Raphson iterations + + interfaces(k) = xi0 + + endif + + !!! Now do the rest of the column + + do k = 2,N + if ( ( target_value < edge_values(k,1) ) .AND. ( target_value < edge_values(k,2) ) ) then + ! target less than whole cell + interfaces(k) = -1 + elseif ( (target_value > edge_values(k,1) ) .AND. (target_value > edge_values(k,1) ) ) then + ! target greater than whole cell + interfaces(k) = 2 + elseif ( ( target_value >= edge_values(k-1,2) ) .AND. ( target_value <= edge_values(k,1) ) ) then + ! target greater than cell above but less than current cell + interfaces(k) = 0 + elseif ( ( target_value <= edge_values(k-1,2) ) .AND. ( target_value >= edge_values(k,1) ) ) then + ! target less than cell above but greater than current cell + interfaces(k) = 0 + else + ! Interface is within cell, so use Newton-Raphson iterations to find it + + ! Reset all polynomial coefficients to 0 and copy those pertaining to + ! the found cell + a(:) = 0.0 + do i = 1,degree+1 + a(i) = ppoly_coefs(k,i) + enddo + + ! Guess the middle of the cell to start Newton-Raphson iterations + xi0 = 0.5 + + ! Newton-Raphson iterations + do iter = 1,NR_ITERATIONS + + if (use_2018_answers) then + numerator = a(1) + a(2)*xi0 + a(3)*xi0*xi0 + a(4)*xi0*xi0*xi0 + & + a(5)*xi0*xi0*xi0*xi0 - target_value + denominator = a(2) + 2*a(3)*xi0 + 3*a(4)*xi0*xi0 + 4*a(5)*xi0*xi0*xi0 + else ! These expressions are mathematicaly equivalent but more accurate. + numerator = (a(1) - target_value) + xi0*(a(2) + xi0*(a(3) + xi0*(a(4) + a(5)*xi0))) + denominator = a(2) + xi0*(2.*a(3) + xi0*(3.*a(4) + 4.*a(5)*xi0)) + endif + + delta = -numerator / denominator + + xi0 = xi0 + delta + + ! Check whether new estimate is out of bounds. If the new estimate is + ! indeed out of bounds, we manually set it to be equal to the overtaken + ! bound with a small offset towards the interior when the gradient of + ! the function at the boundary is zero (in which case, the Newton-Raphson + ! algorithm does not converge). + if ( xi0 < 0.0 ) then + xi0 = 0.0 + grad = a(2) + if ( grad == 0.0 ) xi0 = xi0 + eps + endif + + if ( xi0 > 1.0 ) then + xi0 = 1.0 + if (use_2018_answers) then + grad = a(2) + 2*a(3) + 3*a(4) + 4*a(5) + else ! These expressions are mathematicaly equivalent but more accurate. + grad = a(2) + (2.*a(3) + (3.*a(4) + 4.*a(5))) + endif + if ( grad == 0.0 ) xi0 = xi0 - eps + endif + + ! break if converged or too many iterations taken + if ( abs(delta) < NR_TOLERANCE ) exit + enddo ! end Newton-Raphson iterations + + interfaces(k) = xi0 + + endif + enddo + + ! print *, "target_value", target_value + ! print *, "edge_values(:,1)", edge_values(:,1) + ! print *, "edge_values(:,2)", edge_values(:,2) + ! print *, "interfaces", interfaces + +end function get_interface_indices + +!> Given target values (e.g., density), build new grid based on polynomial +!! +!! Given the grid 'grid0' and the piecewise polynomial interpolant +!! 'ppoly0' (possibly discontinuous), the coordinates of the new grid 'grid1' +!! are determined by finding the corresponding target interface densities. +subroutine get_histogram_weights( n0, ppoly0_E, ppoly0_coefs, & + target_values, degree, n1, histogram_weights, answer_date ) +integer, intent(in) :: n0 !< Number of points on source grid +integer, intent(in) :: n1 !< Number of points on target grid +real, dimension(n0,2), intent(in) :: ppoly0_E !< Edge values of interpolating polynomials [A] +real, dimension(n0,DEGREE_MAX+1), & +intent(in) :: ppoly0_coefs !< Coefficients of interpolating polynomials [A] +real, dimension(n1+1), intent(in) :: target_values !< Target values of interfaces [A] +integer, intent(in) :: degree !< Degree of interpolating polynomials +integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use +real, dimension(n0,n1), intent(inout) :: histogram_weights !< Matrix of weights mapping source grid cells + !< to target grid layers + +! Local variables +real, dimension(n0) :: iindices_l !< array to hold interface indices of target values +real, dimension(n0) :: iindices_u !< array to hold interface indices of target values +real :: index_u !< dummy array to hold upper interface index +real :: index_l !< dummy array to hold lower interface index +real :: edge_value_shallow !< dummy array to hold edge value at shallow edge of cell on source grid +real :: edge_value_deep !< dummy array to hold edge value at deep edge of cell on source grid +integer :: k0 ! loop index +integer :: k1 ! loop index +real :: tl ! current interface target density [A] +real :: tu ! dummy variable for interface target value above t +real, dimension(n0) :: weightsum ! dummy array for checking whether all of a bin is accounted for + +do k1=1,n1 + tl = target_values(k1) + tu = target_values(k1+1) + iindices_l(:) = get_interface_indices( n0, ppoly0_E, ppoly0_coefs, tl, degree, answer_date ) + iindices_u(:) = get_interface_indices( n0, ppoly0_E, ppoly0_coefs, tu, degree, answer_date ) + do k0=1,n0 + edge_value_shallow = ppoly0_E(k0,1) + edge_value_deep = ppoly0_E(k0,2) + index_l = iindices_l(k0) + index_u = iindices_u(k0) + if ( ( index_l == -1 ) .AND. ( index_u == -1 ) ) then + ! Both interfaces have values lower than this cell + histogram_weights(k0,k1) = 0 + elseif ( (index_l == 2 ) .AND. ( index_u == 2 ) ) then + ! Both interfaces have values greater than this cell + histogram_weights(k0,k1) = 0 + elseif ( (index_l == -1 ) .AND. ( index_u == 2 ) ) then + ! The upper interface is greater, the lower interface is less + ! Therefore all of the cell is within the bin + histogram_weights(k0,k1) = 1 + !! It should never occur the other way round + !! - that the upper interface target value is lower than the cell, + !! and the lower interface target value is greater than the cell, + !! since this would imply that the lower interface target is greater than the upper interface target + elseif ( ( (index_l >= 0 ) .AND. (index_l < 1) ) .AND. ( (index_u >=0 ) .AND. (index_u < 1) ) ) then + ! Both interfaces are within the cell + ! Their orientation doesn't matter, so the weight is just the magnitude of the difference + histogram_weights(k0,k1) = abs(index_l - index_u) + elseif ( ( ( (index_l >= 0 ) .AND. (index_l < 1) ) ) .AND. ( (index_u == -1 ) .OR. ( index_u == 2 ) ) ) then + ! Lower interface is within cell (and upper interface is not) + if ( edge_value_shallow > edge_value_deep ) then + ! Values decreasing in cell + ! Therefore include shallower portion of cell + histogram_weights(k0,k1) = index_l + else + ! Values increasing in cell + ! Therefore include deeper portion of cell + histogram_weights(k0,k1) = 1 - index_l + endif + elseif ( ( ( (index_u >= 0 ) .AND. (index_u < 1) ) ) .AND. ( (index_l == -1 ) .OR. ( index_l == 2 ) ) ) then + ! Upper interface is within cell (and lower interface is not) + if ( edge_value_shallow > edge_value_deep ) then + ! Values decreasing in cell + ! Therefore include deeper portion of cell + histogram_weights(k0,k1) = 1 - index_u + else + ! Values increasing in cell + ! Therefore include shallower portion of cell + histogram_weights(k0,k1) = index_u + endif + endif + enddo + + ! print *, "tl", tl + ! print *, "tu", tu + ! print *, "iindices_l", iindices_l + ! print *, "iindices_u", iindices_u + ! print *, "histogram_weights", histogram_weights(:,k1) + +enddo + +! weightsum=0.0 +! do k1 = 1,n1 +! weightsum(:) = weightsum(:) + histogram_weights(:,k1) +! enddo + +! print *, "weightsum", weightsum + +end subroutine get_histogram_weights + !> Numeric value of interpolation_scheme corresponding to scheme name integer function interpolation_scheme(interp_scheme) character(len=*), intent(in) :: interp_scheme !< Name of the interpolation scheme diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 3bb73e4c57..2c946e8434 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -17,6 +17,7 @@ module MOM_diag_mediator use MOM_diag_remap, only : diag_remap_ctrl, diag_remap_update, diag_remap_calc_hmask use MOM_diag_remap, only : diag_remap_init, diag_remap_end, diag_remap_do_remap use MOM_diag_remap, only : vertically_reintegrate_diag_field, vertically_interpolate_diag_field +use MOM_diag_remap, only : vertically_histogram_diag_field use MOM_diag_remap, only : horizontally_average_diag_field, diag_remap_get_axes_info use MOM_diag_remap, only : diag_remap_configure_axes, diag_remap_axes_configured use MOM_diag_remap, only : diag_remap_diag_registration_closed, diag_remap_set_active @@ -34,6 +35,7 @@ module MOM_diag_mediator use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type +use MOM_regridding, only : check_if_histogram_extensive_diags implicit none ; private @@ -510,6 +512,7 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical) ! Allocate these arrays since the size of the diagnostic array is now known allocate(diag_cs%diag_remap_cs(i)%h(G%isd:G%ied,G%jsd:G%jed, diag_cs%diag_remap_cs(i)%nz)) allocate(diag_cs%diag_remap_cs(i)%h_extensive(G%isd:G%ied,G%jsd:G%jed, diag_cs%diag_remap_cs(i)%nz)) + allocate(diag_cs%diag_remap_cs(i)%hweights3d(G%isd:G%ied,G%jsd:G%jed, GV%ke, diag_cs%diag_remap_cs(i)%nz)) ! This vertical coordinate has been configured so can be used. if (diag_remap_axes_configured(diag_cs%diag_remap_cs(i))) then @@ -1565,6 +1568,8 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) dz_diag, & ! Layer vertical extents for remapping [Z ~> m] dz_begin ! Layer vertical extents for remapping extensive quantities [Z ~> m] + logical :: histogram_extensive_diags + if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator) ! For intensive variables only, we can choose to use a different diagnostic grid to map to @@ -1612,36 +1617,63 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) staggered_in_y = diag%axes%is_v_point .or. diag%axes%is_q_point if (diag%v_extensive .and. .not.diag%axes%is_native) then - ! The field is vertically integrated and needs to be re-gridded + ! The field is vertically integrated and needs to either be re-gridded + ! or binned into a new coordinate using a histogram approach if (present(mask)) then call MOM_error(FATAL,"post_data_3d: no mask for regridded field.") endif - if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) - allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz)) - if (diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%Z_based_coord) then - call vertically_reintegrate_diag_field( & - diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), diag_cs%G, & - dz_begin, diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%h_extensive, & - staggered_in_x, staggered_in_y, diag%axes%mask3d, field, remapped_field) - else - call vertically_reintegrate_diag_field( & - diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), diag_cs%G, & - diag_cs%h_begin, diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%h_extensive, & - staggered_in_x, staggered_in_y, diag%axes%mask3d, field, remapped_field) - endif - if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) - if (associated(diag%axes%mask3d)) then - ! Since 3d masks do not vary in the vertical, just use as much as is - ! needed. - call post_data_3d_low(diag, remapped_field, diag_cs, is_static, & - mask=diag%axes%mask3d) + call check_if_histogram_extensive_diags(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%regrid_cs, histogram_extensive_diags) + if ( histogram_extensive_diags ) then + ! Take histogramming approach + if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) + allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz)) + ! Shouldn't ever be specified for a Z-based coordinate, so can call directly, unlike below + call vertically_histogram_diag_field( & + diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), diag_cs%G, & + diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%hweights3d, & + staggered_in_x, staggered_in_y, diag%axes%mask3d, field, remapped_field) + if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) + if (associated(diag%axes%mask3d)) then + ! Since 3d masks do not vary in the vertical, just use as much as is + ! needed. + call post_data_3d_low(diag, remapped_field, diag_cs, is_static, & + mask=diag%axes%mask3d) + else + call post_data_3d_low(diag, remapped_field, diag_cs, is_static) + endif + if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) + deallocate(remapped_field) + if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) else - call post_data_3d_low(diag, remapped_field, diag_cs, is_static) + ! The field is vertically integrated and needs to be re-gridded + ! Take regridding approach + if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) + allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz)) + if (diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%Z_based_coord) then + call vertically_reintegrate_diag_field( & + diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), diag_cs%G, & + dz_begin, diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%h_extensive, & + staggered_in_x, staggered_in_y, diag%axes%mask3d, field, remapped_field) + else + call vertically_reintegrate_diag_field( & + diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), diag_cs%G, & + diag_cs%h_begin, diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%h_extensive, & + staggered_in_x, staggered_in_y, diag%axes%mask3d, field, remapped_field) + endif + if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) + if (associated(diag%axes%mask3d)) then + ! Since 3d masks do not vary in the vertical, just use as much as is + ! needed. + call post_data_3d_low(diag, remapped_field, diag_cs, is_static, & + mask=diag%axes%mask3d) + else + call post_data_3d_low(diag, remapped_field, diag_cs, is_static) + endif + if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) + deallocate(remapped_field) + if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) endif - if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) - deallocate(remapped_field) - if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) elseif (diag%axes%needs_remapping) then ! Remap this field to another vertical coordinate. if (present(mask)) then @@ -3546,7 +3578,7 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S, update_intensiv diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h) else call diag_remap_update(diag_cs%diag_remap_cs(m), diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, & - diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h) + diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h, diag_cs%diag_remap_cs(m)%hweights3d) endif enddo endif @@ -3558,7 +3590,7 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S, update_intensiv diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h_extensive) else call diag_remap_update(diag_cs%diag_remap_cs(m), diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, & - diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h_extensive) + diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h_extensive, diag_cs%diag_remap_cs(m)%hweights3d) endif enddo endif diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 38553a4351..065b77dbd1 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -41,15 +41,17 @@ module MOM_diag_remap use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : EOS_type use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h -use MOM_remapping, only : interpolate_column, reintegrate_column +use MOM_remapping, only : interpolate_column, reintegrate_column, histogram_column use MOM_regridding, only : regridding_CS, initialize_regridding, end_regridding use MOM_regridding, only : set_regrid_params, get_regrid_size +use MOM_regridding, only : check_if_histogram_extensive_diags use MOM_regridding, only : getCoordinateInterfaces, set_h_neglect, set_dz_neglect -use MOM_regridding, only : get_zlike_CS, get_sigma_CS, get_rho_CS +use MOM_regridding, only : get_zlike_CS, get_sigma_CS, get_rho_CS, get_scalar_CS use regrid_consts, only : coordinateMode use coord_zlike, only : build_zstar_column use coord_sigma, only : build_sigma_column use coord_rho, only : build_rho_column +use coord_scalar, only : build_scalar_column implicit none ; private @@ -64,6 +66,7 @@ module MOM_diag_remap public diag_remap_diag_registration_closed public vertically_reintegrate_diag_field public vertically_interpolate_diag_field +public vertically_histogram_diag_field public horizontally_average_diag_field !> Represents remapping of diagnostics to a particular vertical coordinate. @@ -90,6 +93,7 @@ module MOM_diag_remap !! vertical extents in [Z ~> m], depending on the setting of Z_based_coord. real, dimension(:,:,:), allocatable :: h_extensive !< Remap grid thicknesses in [H ~> m or kg m-2] or !! vertical extents in [Z ~> m] for remapping extensive variables + real, dimension(:,:,:,:), allocatable :: hweights3d !< Mapping of histogram weights integer :: interface_axes_id = 0 !< Vertical axes id for remapping at interfaces integer :: layer_axes_id = 0 !< Vertical axes id for remapping on layers logical :: om4_remap_via_sub_cells !< Use the OM4-era ramap_via_sub_cells @@ -204,6 +208,9 @@ subroutine diag_remap_configure_axes(remap_cs, GV, US, param_file) elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then units = 'kg m-3' longname = 'Target Potential Density' + elseif (remap_cs%vertical_coord == coordinateMode('SCALAR')) then + units = 'degC' + longname = 'Target Scalar Values' else units = 'meters' longname = 'Depth' @@ -263,7 +270,7 @@ function diag_remap_axes_configured(remap_cs) !! height or layer thicknesses changes. In the case of density-based !! coordinates then technically we should also regenerate the !! target grid whenever T/S change. -subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_target) +subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_target, hweights3d) type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diagnostic coordinate control structure type(ocean_grid_type), pointer :: G !< The ocean's grid type type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -279,6 +286,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe real, dimension(SZI_(G),SZJ_(G),remap_cs%nz), & intent(inout) :: h_target !< The new diagnostic thicknesses in [H ~> m or kg m-2] !! or [Z ~> m], depending on the value of remap_cs%Z_based_coord + real, optional, dimension(SZI_(G),SZJ_(G),SZK_(GV),remap_cs%nz) :: hweights3d ! Local variables real, dimension(remap_cs%nz + 1) :: zInterfaces ! Interface positions [H ~> m or kg m-2] or [Z ~> m] @@ -287,7 +295,9 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe real :: h_tot(SZI_(G),SZJ_(G)) ! The total thickness of the water column [H ~> m or kg m-2] or [Z ~> m] real :: Z_unit_scale ! A conversion factor from Z-units the internal work units in this routine, ! in units of [H Z-1 ~> 1 or kg m-3] or [nondim], depending on remap_cs%Z_based_coord. - integer :: i, j, k, is, ie, js, je, nz + integer :: i, j, k, is, ie, js, je, nz, k0, k1 + real, dimension(SZK_(GV),remap_cs%nz) :: histogram_weights + logical :: histogram_extensive_diags ! Note that coordinateMode('LAYER') is never 'configured' so will always return here. if (.not. remap_cs%configured) return @@ -360,8 +370,38 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe ! This function call can work with 5 arguments in units of [Z ~> m] or [H ~> kg m-2]. call build_rho_column(get_rho_CS(remap_cs%regrid_cs), GV%ke, & bottom_depth(i,j), h(i,j,:), T(i,j,:), S(i,j,:), & - eqn_of_state, zInterfaces, h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) + eqn_of_state, zInterfaces, histogram_weights=histogram_weights, h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) + do k=1,nz ; h_target(i,j,k) = zInterfaces(K) - zInterfaces(K+1) ; enddo + + ! Fill out 4d array of histogram weights + call check_if_histogram_extensive_diags(remap_cs%regrid_cs,histogram_extensive_diags) + if ( histogram_extensive_diags ) then + do k0 = 1,GV%ke + do k1 = 1,nz + hweights3d(i,j,k0,k1) = histogram_weights(k0,k1) + enddo + enddo + endif + endif ; enddo ; enddo + elseif (remap_cs%vertical_coord == coordinateMode('SCALAR')) then + do j=js-1,je+1 ; do i=is-1,ie+1 ; if (G%mask2dT(i,j) > 0.0) then + ! This function call can work with 5 arguments in units of [Z ~> m] or [H ~> kg m-2]. + call build_scalar_column(get_scalar_CS(remap_cs%regrid_cs), GV%ke, & + bottom_depth(i,j), h(i,j,:), T(i,j,:), S(i,j,:), & + eqn_of_state, zInterfaces, histogram_weights=histogram_weights, h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) + + do k=1,nz ; h_target(i,j,k) = zInterfaces(K) - zInterfaces(K+1) ; enddo + + ! Fill out 4d array of histogram weights + call check_if_histogram_extensive_diags(remap_cs%regrid_cs,histogram_extensive_diags) + if ( histogram_extensive_diags ) then + do k0 = 1,GV%ke + do k1 = 1,nz + hweights3d(i,j,k0,k1) = histogram_weights(k0,k1) + enddo + enddo + endif endif ; enddo ; enddo elseif (remap_cs%vertical_coord == coordinateMode('HYCOM1')) then call MOM_error(FATAL,"diag_remap_update: HYCOM1 coordinate not coded for diagnostics yet!") @@ -657,6 +697,89 @@ subroutine vertically_reintegrate_field(remap_cs, G, isdf, jsdf, h, h_target, st end subroutine vertically_reintegrate_field +!> Vertically histogram a diagnostic field to alternative vertical grid. +subroutine vertically_histogram_diag_field(remap_cs, G, hweights3d, staggered_in_x, staggered_in_y, & + mask, field, histogrammed_field) + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + real, dimension(:,:,:,:), intent(in) :: hweights3d !< Weights used to map from grid to histogram + logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points + logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points + real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim]. Note that because this + !! is a pointer it retains its declared indexing conventions. + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(:,:,:), intent(out) :: histogrammed_field !< Field argument remapped to alternative coordinate [A] + + ! Local variables + integer :: isdf, jsdf !< The starting i- and j-indices in memory for field + + call assert(remap_cs%initialized, 'vertically_histogram_diag_field: remap_cs not initialized.') + call assert(size(field, 3) == size(hweights3d, 3), 'vertically_histogram_diag_field: Remap field and thickness z-axes do not match.') + + isdf = G%isd ; if (staggered_in_x) Isdf = G%IsdB + jsdf = G%jsd ; if (staggered_in_y) Jsdf = G%JsdB + + if (associated(mask)) then + call vertically_histogram_field(remap_cs, G, isdf, jsdf, hweights3d, staggered_in_x, staggered_in_y, & + field, histogrammed_field, mask(:,:,1)) + else + call vertically_histogram_field(remap_cs, G, isdf, jsdf, hweights3d, staggered_in_x, staggered_in_y, & + field, histogrammed_field) + endif + +end subroutine vertically_histogram_diag_field + +!> The internal routine to vertically histogram a diagnostic field to +!! an alternative vertical grid. +subroutine vertically_histogram_field(remap_cs, G, isdf, jsdf, hweights3d, staggered_in_x, staggered_in_y, & + field, histogrammed_field, mask) + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + integer, intent(in) :: isdf !< The starting i-index in memory for field + integer, intent(in) :: jsdf !< The starting j-index in memory for field + real, dimension(G%isd:,G%jsd:,:,:), intent(in) :: hweights3d !< The weights for histogramming from source to target + logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points + logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points + real, dimension(isdf:,jsdf:,:), & + intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(isdf:,jsdf:,:), & + intent(out) :: histogrammed_field !< Field argument remapped to alternative coordinate [A] + real, dimension(isdf:,jsdf:), & + optional, intent(in) :: mask !< A mask for the field [nondim] + + ! Local variables + integer :: i, j ! Grid index + integer :: nz_src, nz_dest + + nz_src = size(field,3) + nz_dest = remap_cs%nz + histogrammed_field(:,:,:) = 0. + + if (staggered_in_x .and. .not. staggered_in_y) then + ! U-points + call assert(.false., 'vertically_histogram_diag_field: U point histogramming is not coded yet.') + elseif (staggered_in_y .and. .not. staggered_in_x) then + ! V-points + call assert(.false., 'vertically_histogram_diag_field: V point histogramming is not coded yet.') + elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then + ! H-points + if (present(mask)) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (mask(i,J) > 0.0) then + call histogram_column(nz_src, field(i,j,:), & + nz_dest, histogrammed_field(i,j,:), hweights3d(i,j,:,:)) + endif ; enddo ; enddo + else + do j=G%jsc,G%jec ; do i=G%isc,G%iec + call histogram_column(nz_src, field(i,j,:), & + nz_dest, histogrammed_field(i,j,:), hweights3d(i,j,:,:)) + enddo ; enddo + endif + else + call assert(.false., 'vertically_histogram_diag_field: Q point remapping is not coded yet.') + endif + +end subroutine vertically_histogram_field + !> Vertically interpolate diagnostic field to alternative vertical grid. subroutine vertically_interpolate_diag_field(remap_cs, G, h, staggered_in_x, staggered_in_y, & mask, field, interpolated_field)