diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 9bc71dd15f..77959d2990 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -1923,7 +1923,6 @@ subroutine convective_adjustment(G, GV, h, tv) end subroutine convective_adjustment - !------------------------------------------------------------------------------ !> Return a uniform resolution vector in the units of the coordinate function uniformResolution(nk,coordMode,maxDepth,rhoLight,rhoHeavy) diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90 index 5299c74b1b..b367157ecc 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 coord_utils, only : build_scalar_column implicit none ; private @@ -87,7 +88,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, & +subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, ksort, & 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) @@ -98,6 +99,7 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & type(EOS_type), pointer :: eqn_of_state !< Equation of state structure real, dimension(CS%nk+1), & intent(inout) :: z_interface !< Absolute positions of interfaces + real, dimension(nz), intent(out) :: ksort !< k-indices for monotonically increasing column [ppt] 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 @@ -116,56 +118,13 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & ! 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) = densities(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) - - ! 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 + ! Compute densities on source column + pres(:) = CS%ref_pressure + call calculate_density(T, S, pres, densities, eqn_of_state) + + ! Build column + call build_scalar_column(CS, nz, depth, h, densities, z_interface, ksort, & + h_neglect, h_neglect_edge) end subroutine build_rho_column diff --git a/src/ALE/coord_utils.F90 b/src/ALE/coord_utils.F90 new file mode 100644 index 0000000000..01911383a5 --- /dev/null +++ b/src/ALE/coord_utils.F90 @@ -0,0 +1,240 @@ +!> Regrid columns for the continuous isopycnal (rho) coordinate +module coord_utils + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_error_handler, only : MOM_error, FATAL +use regrid_interp, only : interp_CS_type, build_and_interpolate_grid, DEGREE_MAX + +implicit none ; private + +public build_scalar_column, sort_scalar_k_1d, copy_finite_thicknesses, old_inflate_layers_1d + +!> Build a scalar coordinate column +!! +!! 1. Positions of target scalar values (for interfaces) are found by interpolation. +subroutine build_scalar_column(CS, nz, depth, h, phi, z_interface, ksort, & + h_neglect, h_neglect_edge) + type(rho_CS), intent(in) :: CS !< coord_* control structure + integer, intent(in) :: nz !< Number of levels on source grid (i.e. length of h, phi) + real, intent(in) :: depth !< Depth of ocean bottom (positive in m) + real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(nz), intent(in) :: phi !< Scalar for source column + real, dimension(CS%nk+1), & + intent(inout) :: z_interface !< Absolute positions of interfaces + real, dimension(nz), intent(out) :: ksort !< k-indices for monotonically increasing column + 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] + + ! Local variables + integer :: k, count_nonzero_layers + integer, dimension(nz) :: mapping + real, dimension(nz) :: h_nv ! Thicknesses of non-vanishing layers [H ~> m or kg m-2] + 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 + + ! 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 + + do k = 1,count_nonzero_layers + phi(k) = phi(mapping(k)) + enddo + + ! Determine k-indices for a monotonically increasing column + ! Do not change the scalar distribution in the column + call sort_scalar_k_1d(nz, phi, ksort) + + ! Based on source column density profile, interpolate to generate a new grid + call build_and_interpolate_grid(CS%interp_CS, phi, count_nonzero_layers, & + h_nv, xTmp, CS%target_density, CS%nk, h_new, & + x1, h_neglect, h_neglect_edge) + + ! 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 + +!------------------------------------------------------------------------------ +!> Return the index of a sorted array of scalar values +subroutine sort_scalar_k_1d(nz, phi, ksort) + integer, intent(in) :: nz !< number of levels on source grid + real, dimension(nz), & + intent(in) :: phi !< Array of scalar quantity to be sorted + integer, dimension(nz), & + intent(out) :: ksort !< An array of indicies for a + !! monotonically increasing scalar +!------------------------------------------------------------------------------ +! Check each water column to see if a given scalar is monotonically increasing. +! If not, return an array of the sorted indices (bubble sort algorithm). +! No need to return the sorted scalar array itself. +!------------------------------------------------------------------------------ + + ! Local variables + integer :: k + real :: P0, P1 ! temperatures + logical :: monotonic + + ! Repeat swapping of indices until complete + do + monotonic = .true. + do k = 1,nz-1 + ! Gather information of scalar value in current and next cells + P0 = phi(k) ; P1 = phi(k+1) + ! If the scalar value of the current cell is larger than the scalar + ! below it, we swap the cell indices + if ( P0 > P1 ) then + ksort(k) = k+1 ; ksort(k+1) = k + monotonic = .false. + endif + enddo ! k + + if ( monotonic ) exit + enddo + +end subroutine sort_scalar_k_1d + +!> 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] + integer, intent(in) :: nk !< Number of layers in the grid + real, dimension(:), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2] + + ! Local variable + integer :: k + integer :: k_found + integer :: count_nonzero_layers + real :: delta + real :: correction + real :: maxThickness + + ! 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_utils