From 218a965c5f8ed1d8545c894aa99bf321e6d8b337 Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Thu, 29 Aug 2024 21:30:45 +0100 Subject: [PATCH 01/19] adds procedures for finding all interfaces in a column based on target values, and deriving histogram weights from those interfaces --- src/ALE/regrid_interp.F90 | 193 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 193 insertions(+) diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index b3100fe8ae..d7a24558b8 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -503,6 +503,199 @@ 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 indices 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) + + do k = 1,N + if ( target_value <= edge_values(k,1) ) then + interfaces(k) = -1 + elseif (target_value >= edge_values(k,2) ) then + interfaces(k) = 2 + elseif ( ( target_value >= edge_values(k-1,2) ) .AND. ( target_value <= edge_values(k,1) ) ) then + 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 + +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, 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 + + + +! 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 upper interface index +real :: index_u_prior !< dummy array to hold upper interface index +real :: index_l_prior !< dummy array to hold upper interface index +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 + +do k1=1,n1 + tl = target_value(k) + tu = target_values(k+1) + iindices_l(:) = get_interface_indices( N, ppoly0_E, ppoly0_coefs, tl, degree, answer_date ) + iindices_u(:) = get_interface_indices( N, ppoly0_E, ppoly0_coefs, tu, degree, answer_date ) + !!! Note that I need a different approach for getting weights in very first cell + !!! and then will start this loop from k=2 + do k0=1,n0 + index_l = iindices_l(k0) + index_l_prior = iindices_l(k0-1) + index_u = iindices_u(k0) + index_u_prior = iindices_u(k0-1) + 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 rounnd + !! - that the upper interface is lower than the cell, + !! and the lower interface is greater than the cell, + !! since this would imply that the lower interface is greater than the upper interface + 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 + if ( index_l_prior < 2 ) then + ! and is less than or within the previous cell + ! Therefore include shallower portion of cell + histogram_weights(k0,k1) = index_l + elseif ( index_l_prior == 2 ) then + ! and is greater than the previous 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 + ! Lower interface is within cell + if ( index_u_prior < 2 ) then + ! and is greater than or within the previous cell + ! Therefore include deeper portion of cell + histogram_weights(k0,k1) = 1 - index_u + elseif ( index_u_prior == 2 ) then + ! and is less than the previous cell + ! Therefore include shallower portion of cell + histogram_weights(k0,k1) = index_u + endif + endif + enddo +enddo + +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 From fe91d89e8af9ec5b57eb247daee468f70ac49fe1 Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Fri, 30 Aug 2024 10:02:04 +0100 Subject: [PATCH 02/19] allows handling of first grid cell of source grid --- src/ALE/regrid_interp.F90 | 109 ++++++++++++++++++++++++++++++-------- 1 file changed, 88 insertions(+), 21 deletions(-) diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index d7a24558b8..bffd5f6864 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -540,7 +540,74 @@ function get_interface_indices( N, edge_values, ppoly_coefs, & k_found = -1 use_2018_answers = (answer_date < 20190101) - do k = 1,N + !!! 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) ) then interfaces(k) = -1 elseif (target_value >= edge_values(k,2) ) then @@ -630,26 +697,26 @@ subroutine get_histogram_weights( n0, ppoly0_E, ppoly0_coefs, & 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 upper interface index -real :: index_u_prior !< dummy array to hold upper interface index -real :: index_l_prior !< 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 do k1=1,n1 - tl = target_value(k) + tl = target_values(k) tu = target_values(k+1) iindices_l(:) = get_interface_indices( N, ppoly0_E, ppoly0_coefs, tl, degree, answer_date ) iindices_u(:) = get_interface_indices( N, ppoly0_E, ppoly0_coefs, tu, degree, answer_date ) !!! Note that I need a different approach for getting weights in very first cell !!! and then will start this loop from k=2 do k0=1,n0 + edge_value_shallow = ppoly0_E(k0,1) + edge_value_deep = ppoly0_E(k0,2) index_l = iindices_l(k0) - index_l_prior = iindices_l(k0-1) index_u = iindices_u(k0) - index_u_prior = iindices_u(k0-1) if ( ( index_l == -1 ) .AND. ( index_u == -1 ) ) then ! Both interfaces have values lower than this cell histogram_weights(k0,k1) = 0 @@ -660,33 +727,33 @@ subroutine get_histogram_weights( n0, ppoly0_E, ppoly0_coefs, & ! 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 rounnd - !! - that the upper interface is lower than the cell, - !! and the lower interface is greater than the cell, - !! since this would imply that the lower interface is greater than the upper interface + !! 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 - if ( index_l_prior < 2 ) then - ! and is less than or within the previous cell + ! 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 - elseif ( index_l_prior == 2 ) then - ! and is greater than the previous cell + 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 - ! Lower interface is within cell - if ( index_u_prior < 2 ) then - ! and is greater than or within the previous cell + ! 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 - elseif ( index_u_prior == 2 ) then - ! and is less than the previous cell + else + ! Values increasing in cell ! Therefore include shallower portion of cell histogram_weights(k0,k1) = index_u endif From 226460eda57f76cf9577a7fdd727def66df16824 Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Fri, 30 Aug 2024 12:12:12 +0100 Subject: [PATCH 03/19] implements evaluation of weights in diag_remap_update and creates logical for histogramming of extensive variables --- src/ALE/coord_rho.F90 | 20 ++++++++++++++--- src/ALE/regrid_interp.F90 | 38 +++++++++++++++++++++++++++----- src/framework/MOM_diag_remap.F90 | 14 +++++++++++- 3 files changed, 62 insertions(+), 10 deletions(-) diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90 index 3ed769f4e4..b8ab054d0b 100644 --- a/src/ALE/coord_rho.F90 +++ b/src/ALE/coord_rho.F90 @@ -29,6 +29,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 +41,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 +56,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 +72,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 +80,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] + real, 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 +88,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 +98,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 +116,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 +149,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, 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/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index bffd5f6864..238444739c 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -353,6 +353,34 @@ 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, 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) :: 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 @@ -520,7 +548,7 @@ function get_interface_indices( N, edge_values, ppoly_coefs, & 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 indices in relation to each cell + real, dimension(N) :: interfaces !< An array indicating location of target interfaces in relation to each cell ! Local variables real :: xi0 ! normalized target coordinate [nondim] @@ -681,7 +709,7 @@ end function get_interface_indices !! '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, answer_date ) + 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] @@ -690,8 +718,8 @@ subroutine get_histogram_weights( n0, ppoly0_E, ppoly0_coefs, & 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 @@ -710,8 +738,6 @@ subroutine get_histogram_weights( n0, ppoly0_E, ppoly0_coefs, & tu = target_values(k+1) iindices_l(:) = get_interface_indices( N, ppoly0_E, ppoly0_coefs, tl, degree, answer_date ) iindices_u(:) = get_interface_indices( N, ppoly0_E, ppoly0_coefs, tu, degree, answer_date ) - !!! Note that I need a different approach for getting weights in very first cell - !!! and then will start this loop from k=2 do k0=1,n0 edge_value_shallow = ppoly0_E(k0,1) edge_value_deep = ppoly0_E(k0,2) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index e8e6a756e9..6c26bf5375 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -263,7 +263,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 +279,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] @@ -288,6 +289,8 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe 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 + 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,7 +363,16 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe 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) + do k=1,nz ; h_target(i,j,k) = zInterfaces(K) - zInterfaces(K+1) ; enddo + call check_if_histogram_extensive_diags(remap_cs%regrid_cs,histogram_extensive_diags) + if ( histogram_extensive_diags ) then + 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, & + histogram_weights=histogram_weights, h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) + hweights3d(i,j,:,:) = histogram_weights + 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!") From e69be9eee85cdab01d6286ea0008a8b2da2c594c Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Fri, 30 Aug 2024 12:37:12 +0100 Subject: [PATCH 04/19] adds histogram weights to diagnostic control structure and allows updating in diag mediator --- src/framework/MOM_diag_mediator.F90 | 4 ++-- src/framework/MOM_diag_remap.F90 | 20 +++++++++++--------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index c28e2e5896..d6c695598f 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -3544,7 +3544,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 @@ -3556,7 +3556,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 6c26bf5375..43f521ed2d 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -90,6 +90,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 @@ -288,7 +289,7 @@ 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 @@ -362,16 +363,17 @@ 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 - call check_if_histogram_extensive_diags(remap_cs%regrid_cs,histogram_extensive_diags) - if ( histogram_extensive_diags ) then - 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, & - histogram_weights=histogram_weights, h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) - hweights3d(i,j,:,:) = histogram_weights + + ! Fill out 4d array of histogram weights + if ( allocated(hweights3d) ) then + do k0 = 1,GV%ke + do k1 = 1,nk + hweights3d(i,j,k0,k1) = histogram_weights(k0,k1) + enddo + enddo endif endif ; enddo ; enddo elseif (remap_cs%vertical_coord == coordinateMode('HYCOM1')) then From 4fcb60e08755a36b18e4e1661de262c14891f38d Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Fri, 30 Aug 2024 13:31:22 +0100 Subject: [PATCH 05/19] changes approach to deciding whether to allocate to 4d array --- src/framework/MOM_diag_mediator.F90 | 1 + src/framework/MOM_diag_remap.F90 | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index d6c695598f..bb4e5fbf81 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -510,6 +510,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, G%ksd:G%ked, 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 diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 43f521ed2d..0499930b2c 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -368,7 +368,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe do k=1,nz ; h_target(i,j,k) = zInterfaces(K) - zInterfaces(K+1) ; enddo ! Fill out 4d array of histogram weights - if ( allocated(hweights3d) ) then + if ( get_rho_CS(remap_cs%regrid_cs)%histogram_extensive_diags ) then do k0 = 1,GV%ke do k1 = 1,nk hweights3d(i,j,k0,k1) = histogram_weights(k0,k1) From eaa8e3d47b60d634443843b86d597ced9f4bfbb9 Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Fri, 30 Aug 2024 13:33:50 +0100 Subject: [PATCH 06/19] makes functions available elsewhere --- src/ALE/regrid_interp.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index 238444739c..e2ad848c94 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) From 6916c81ef57090c126888ea71e3b50e070d4c574 Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Fri, 30 Aug 2024 14:34:12 +0100 Subject: [PATCH 07/19] fixes bugs to get compiled; adds capacity to specify histogram logical in parameter file --- src/ALE/MOM_regridding.F90 | 23 +++++++++++++++++++++-- src/ALE/coord_rho.F90 | 5 +++-- src/ALE/regrid_interp.F90 | 17 +++++++++-------- src/framework/MOM_diag_mediator.F90 | 2 +- src/framework/MOM_diag_remap.F90 | 8 +++++--- 5 files changed, 39 insertions(+), 16 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 0d325292f5..35efbd4643 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -118,6 +118,9 @@ 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 @@ -142,6 +145,7 @@ module MOM_regridding public DEFAULT_COORDINATE_MODE public set_h_neglect, set_dz_neglect public get_zlike_CS, get_sigma_CS, get_rho_CS +public check_if_histogram_extensive_diags !> Documentation for coordinate options character(len=*), parameter, public :: regriddingCoordinateModeDoc = & @@ -210,6 +214,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). @@ -645,7 +650,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 @@ -2031,7 +2042,7 @@ 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_HYCOM1) call init_coord_hycom(CS%hycom_CS, CS%nk, CS%coordinateResolution, CS%target_density, & CS%interp_CS) @@ -2626,6 +2637,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/coord_rho.F90 b/src/ALE/coord_rho.F90 index b8ab054d0b..22f1946b49 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 @@ -80,7 +81,7 @@ subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, histogram !! 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] - real, optional, intent(in) :: histogram_extensive_diags !< If true, use histogram approach to for outputing extensive diags + 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 @@ -151,7 +152,7 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & 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, CS%target_density, CS%nk, & + call build_histogram_weights(CS%interp_CS, densities, nz, h, CS%target_density, CS%nk, & histogram_weights, h_neglect, h_neglect_edge) endif diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index e2ad848c94..3e99168a84 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -355,11 +355,12 @@ subroutine build_and_interpolate_grid(CS, densities, n0, h0, x0, target_values, end subroutine build_and_interpolate_grid !> Build array of weights mapping source grid to target grid -subroutine build_histogram_weights(CS, densities, n0, target_values, & +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 @@ -735,10 +736,10 @@ subroutine get_histogram_weights( n0, ppoly0_E, ppoly0_coefs, & real :: tu ! dummy variable for interface target value above t do k1=1,n1 - tl = target_values(k) - tu = target_values(k+1) - iindices_l(:) = get_interface_indices( N, ppoly0_E, ppoly0_coefs, tl, degree, answer_date ) - iindices_u(:) = get_interface_indices( N, ppoly0_E, ppoly0_coefs, tu, degree, answer_date ) + 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) @@ -758,11 +759,11 @@ subroutine get_histogram_weights( n0, ppoly0_E, ppoly0_coefs, & !! - 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 + 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 + 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 @@ -773,7 +774,7 @@ subroutine get_histogram_weights( n0, ppoly0_E, ppoly0_coefs, & ! 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 + 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 diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index bb4e5fbf81..33e9a49860 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -510,7 +510,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, G%ksd:G%ked, 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 diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 0499930b2c..908a03ca8c 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -44,6 +44,7 @@ module MOM_diag_remap use MOM_remapping, only : interpolate_column, reintegrate_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 regrid_consts, only : coordinateMode @@ -90,7 +91,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 + 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 @@ -368,9 +369,10 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe do k=1,nz ; h_target(i,j,k) = zInterfaces(K) - zInterfaces(K+1) ; enddo ! Fill out 4d array of histogram weights - if ( get_rho_CS(remap_cs%regrid_cs)%histogram_extensive_diags ) then + 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,nk + do k1 = 1,nz hweights3d(i,j,k0,k1) = histogram_weights(k0,k1) enddo enddo From 0cedaa584305f877a04f26213d0ca3c69f65551f Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Fri, 30 Aug 2024 16:42:04 +0100 Subject: [PATCH 08/19] fixes some bugs, implements capacity to do histogramming (though doesn't yet actually do it) at the point of posting data --- src/ALE/MOM_regridding.F90 | 5 ++- src/framework/MOM_diag_mediator.F90 | 65 ++++++++++++++++------------- 2 files changed, 40 insertions(+), 30 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 35efbd4643..3339314710 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -654,7 +654,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m 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**" + "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 @@ -2385,7 +2385,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 @@ -2399,6 +2399,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 diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 33e9a49860..56a6cfcf3d 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -34,6 +34,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 @@ -1566,6 +1567,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 @@ -1613,36 +1616,42 @@ 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 - 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) + 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 + print *, "Eeep!" + ! Here is where I do the histogram approach! 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) + ! The field is vertically integrated and needs to be re-gridded + 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) + 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 From f7d51fb4909aaa784da2699d7812e9c77478e0a4 Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Sat, 31 Aug 2024 11:43:09 +0100 Subject: [PATCH 09/19] includes function call to histogram diags in post_data_3d --- src/framework/MOM_diag_mediator.F90 | 29 +++++++++++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 56a6cfcf3d..a468b48152 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -1616,12 +1616,37 @@ 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 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 + 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 - print *, "Eeep!" - ! Here is where I do the histogram approach! + ! 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 ! The field is vertically integrated and needs to be re-gridded + ! Take regridding approach if (present(mask)) then call MOM_error(FATAL,"post_data_3d: no mask for regridded field.") endif From 429edc6714a8bf43d6e33f78c92beabd27e8eb05 Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Sat, 31 Aug 2024 12:05:44 +0100 Subject: [PATCH 10/19] adds functions for vertical histogramming of diagnostic fields --- src/framework/MOM_diag_remap.F90 | 81 ++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 908a03ca8c..202b6134d9 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -681,6 +681,87 @@ 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(h, 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 + + 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) From 6ce14184f931659f3e9005091294d0f81435bbf6 Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Sat, 31 Aug 2024 12:20:55 +0100 Subject: [PATCH 11/19] adds functions to do the vertical histogram based on weights --- src/ALE/MOM_remapping.F90 | 26 ++++++++++++++++++++++++++ src/framework/MOM_diag_remap.F90 | 2 +- 2 files changed, 27 insertions(+), 1 deletion(-) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 4495e4a170..7b363ec30c 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -78,6 +78,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 @@ -1215,6 +1216,31 @@ 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 + real, dimension(k_src) :: 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/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 202b6134d9..0171e8c68b 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -41,7 +41,7 @@ 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 From b3739f4540bba860e54684e0538a38ccd61689fe Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Sat, 31 Aug 2024 12:50:17 +0100 Subject: [PATCH 12/19] miscellaneous bug fixes --- src/ALE/MOM_remapping.F90 | 3 ++- src/framework/MOM_diag_mediator.F90 | 5 +---- src/framework/MOM_diag_remap.F90 | 13 ++++++++----- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 7b363ec30c..3f7ebbe015 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -1225,7 +1225,8 @@ subroutine histogram_column(nsrc, uh_src, ndest, uh_dest, histogram_weights) real, dimension(nsrc,ndest), intent(in) :: histogram_weights !< Weights mapping source to destination grid ! Local variables - real, dimension(k_src) :: weights + integer :: k_dest, k_src + real, dimension(nsrc) :: weights real :: uh_dest_sum uh_dest(:) = 0.0 diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index a468b48152..23c78a7f50 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 @@ -1647,10 +1648,6 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) else ! The field is vertically integrated and needs to be re-gridded ! Take regridding 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 diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 0171e8c68b..99604fa478 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -65,6 +65,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. @@ -690,7 +691,7 @@ subroutine vertically_histogram_diag_field(remap_cs, G, hweights3d, staggered_in 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. + !! 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] @@ -698,8 +699,7 @@ subroutine vertically_histogram_diag_field(remap_cs, G, hweights3d, staggered_in 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(h, 3), & - 'vertically_histogram_diag_field: Remap field and thickness z-axes do not match.') + 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 @@ -708,10 +708,10 @@ subroutine vertically_histogram_diag_field(remap_cs, G, hweights3d, staggered_in 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, & + 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 @@ -734,7 +734,10 @@ subroutine vertically_histogram_field(remap_cs, G, isdf, jsdf, hweights3d, stagg ! 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 From 9c580f5e1685b9cac5e01747c4057b7916d021b0 Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Sat, 31 Aug 2024 14:42:33 +0100 Subject: [PATCH 13/19] fixes issue with regrid parameters not being correctly set --- src/ALE/MOM_regridding.F90 | 2 ++ src/framework/MOM_diag_mediator.F90 | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 3339314710..bc84aea62c 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -2438,6 +2438,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 @@ -2457,6 +2458,7 @@ 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))) & diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 23c78a7f50..36cb2f8f73 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -1624,7 +1624,9 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) endif call check_if_histogram_extensive_diags(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%regrid_cs, histogram_extensive_diags) + print *, "histogram_extensive_diags", histogram_extensive_diags if ( histogram_extensive_diags ) then + print *, "Eeep! inside logical for histogram check" ! 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)) From f75d2382a47626c3f51a8e0e5a153ecd7f64a232 Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Sat, 31 Aug 2024 14:44:00 +0100 Subject: [PATCH 14/19] removes print statements --- src/framework/MOM_diag_mediator.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 36cb2f8f73..23c78a7f50 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -1624,9 +1624,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) endif call check_if_histogram_extensive_diags(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%regrid_cs, histogram_extensive_diags) - print *, "histogram_extensive_diags", histogram_extensive_diags if ( histogram_extensive_diags ) then - print *, "Eeep! inside logical for histogram check" ! 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)) From c500d6c9ea24a1bdacdfcda953fcc3e0e643028a Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Sat, 31 Aug 2024 15:33:01 +0100 Subject: [PATCH 15/19] implements scalar coordinate (currently temperature) machinery into regridding and diagnostic remapping --- src/ALE/MOM_regridding.F90 | 72 +++++++++++++++++++++++++++++--- src/ALE/regrid_consts.F90 | 5 +++ src/framework/MOM_diag_remap.F90 | 25 ++++++++++- 3 files changed, 96 insertions(+), 6 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index bc84aea62c..5fc6a59c8d 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 @@ -126,6 +127,7 @@ module MOM_regridding 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 rho 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 @@ -144,7 +146,7 @@ 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 @@ -154,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" @@ -231,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 @@ -422,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 @@ -437,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)) @@ -592,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 @@ -600,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 @@ -609,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 @@ -628,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) @@ -841,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) @@ -2009,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 ) @@ -2043,6 +2073,8 @@ subroutine initCoord(CS, GV, US, coord_mode, param_file) 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%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) @@ -2299,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 @@ -2341,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 @@ -2367,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 ) @@ -2463,6 +2509,14 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri 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_rho_params(CS%scalar_CS, min_thickness=min_thickness) + if (present(ref_pressure)) call set_rho_params(CS%scalar_CS, ref_pressure=ref_pressure) + if (present(histogram_extensive_diags)) call set_rho_params(CS%scalar_CS, histogram_extensive_diags=histogram_extensive_diags) + if (present(integrate_downward_for_e)) & + call set_rho_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_rho_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) @@ -2512,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 ) 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/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 99604fa478..4d934b3ada 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -46,11 +46,12 @@ module MOM_diag_remap 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 @@ -207,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' @@ -369,6 +373,25 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe 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 From ad018340b420074095a0f6bae69a91862b356892 Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Sat, 31 Aug 2024 15:37:00 +0100 Subject: [PATCH 16/19] very simply replicates rho coordinate but replaces with temperature rather than density --- src/ALE/coord_scalar.F90 | 437 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 437 insertions(+) create mode 100644 src/ALE/coord_scalar.F90 diff --git a/src/ALE/coord_scalar.F90 b/src/ALE/coord_scalar.F90 new file mode 100644 index 0000000000..037b7dcd4c --- /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, h_neglect, h_neglect_edge) + + call remapping_core_h(remapCS, nz, h0, T, nz, h1, T_tmp, h_neglect, h_neglect_edge) + + ! 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 From 2525f79f95122b4173f80774c4f023bf296da66d Mon Sep 17 00:00:00 2001 From: Graeme MacGilchrist Date: Mon, 7 Oct 2024 13:18:08 +0100 Subject: [PATCH 17/19] historical changes --- src/ALE/MOM_regridding.F90 | 10 +++++----- src/ALE/regrid_interp.F90 | 31 +++++++++++++++++++++++++++++-- 2 files changed, 34 insertions(+), 7 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 5fc6a59c8d..b9042e6c24 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -2510,13 +2510,13 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri 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_rho_params(CS%scalar_CS, min_thickness=min_thickness) - if (present(ref_pressure)) call set_rho_params(CS%scalar_CS, ref_pressure=ref_pressure) - if (present(histogram_extensive_diags)) call set_rho_params(CS%scalar_CS, histogram_extensive_diags=histogram_extensive_diags) + 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_rho_params(CS%scalar_CS, integrate_downward_for_e=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_rho_params(CS%scalar_CS, interp_CS=CS%interp_CS) + 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) diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index 3e99168a84..5614653b48 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -381,6 +381,7 @@ subroutine build_histogram_weights(CS, densities, n0, h0, target_values, & 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 @@ -638,11 +639,17 @@ function get_interface_indices( N, edge_values, ppoly_coefs, & !!! Now do the rest of the column do k = 2,N - if ( target_value <= edge_values(k,1) ) then + 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,2) ) then + 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 @@ -703,6 +710,11 @@ function get_interface_indices( N, edge_values, ppoly_coefs, & 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 @@ -734,6 +746,7 @@ subroutine get_histogram_weights( n0, ppoly0_E, ppoly0_coefs, & 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) @@ -787,8 +800,22 @@ subroutine get_histogram_weights( n0, ppoly0_E, ppoly0_coefs, & 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 From 8f2f030e08a7fb59b32e11685504998243a7d14a Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Tue, 14 Jan 2025 11:36:33 -0800 Subject: [PATCH 18/19] Eliminated h_neglect argument to remapping_core_h Follows from https://github.com/NOAA-GFDL/MOM6/pull/705 --- src/ALE/coord_scalar.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ALE/coord_scalar.F90 b/src/ALE/coord_scalar.F90 index 037b7dcd4c..0576338613 100644 --- a/src/ALE/coord_scalar.F90 +++ b/src/ALE/coord_scalar.F90 @@ -287,9 +287,9 @@ subroutine build_scalar_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_ h1(k) = x1(k+1) - x1(k) enddo - call remapping_core_h(remapCS, nz, h0, S, nz, h1, S_tmp, h_neglect, h_neglect_edge) + call remapping_core_h(remapCS, nz, h0, S, nz, h1, S_tmp) - call remapping_core_h(remapCS, nz, h0, T, nz, h1, T_tmp, h_neglect, h_neglect_edge) + call remapping_core_h(remapCS, nz, h0, T, nz, h1, T_tmp) ! Compute the deviation between two successive grids deviation = 0.0 From 749f112b2bad031805c0e88d7a6fad24f333d250 Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Tue, 4 Feb 2025 11:50:29 -0800 Subject: [PATCH 19/19] Improve comments --- src/ALE/MOM_regridding.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index b9042e6c24..813d7110d1 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -127,7 +127,7 @@ module MOM_regridding 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 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 @@ -681,7 +681,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m "regions appear stratified.", units="nondim", default=0.) 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"//& + "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.)