diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 2a9ebc09c8..137f6cee9b 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -206,12 +206,12 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) call get_param(param_file, mdl, "REMAPPING_SCHEME", string, & "This sets the reconstruction scheme used "//& "for vertical remapping for all variables. "//& - "It can be one of the following schemes: "//& + "It can be one of the following schemes: \n"//& trim(remappingSchemesDoc), default=remappingDefaultScheme) call get_param(param_file, mdl, "VELOCITY_REMAPPING_SCHEME", vel_string, & "This sets the reconstruction scheme used for vertical remapping "//& "of velocities. By default it is the same as REMAPPING_SCHEME. "//& - "It can be one of the following schemes: "//& + "It can be one of the following schemes: \n"//& trim(remappingSchemesDoc), default=trim(string)) call get_param(param_file, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, & "If true, cell-by-cell reconstructions are checked for "//& diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 53072909a5..5e46b8d1f6 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -166,6 +166,7 @@ module MOM_regridding " P1M_H4 (2nd-order accurate)\n"//& " P1M_IH4 (2nd-order accurate)\n"//& " PLM (2nd-order accurate)\n"//& + " PPM_CW (3rd-order accurate)\n"//& " PPM_H4 (3rd-order accurate)\n"//& " PPM_IH4 (3rd-order accurate)\n"//& " P3M_IH4IH3 (4th-order accurate)\n"//& @@ -269,7 +270,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m "determine the new grid. These parameters are "//& "only relevant when REGRIDDING_COORDINATE_MODE is "//& "set to a function of state. Otherwise, it is not "//& - "used. It can be one of the following schemes: "//& + "used. It can be one of the following schemes: \n"//& trim(regriddingInterpSchemeDoc), default=trim(string2)) call set_regrid_params(CS, interp_scheme=string) @@ -582,6 +583,13 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m call set_regrid_params(CS, min_thickness=0.) endif + if (main_parameters .and. coordinateMode(coord_mode) == REGRIDDING_HYCOM1) then + call get_param(param_file, mdl, "HYCOM1_ONLY_IMPROVES", tmpLogical, & + "When regridding, an interface is only moved if this improves the fit to the target density.", & + default=.false.) + call set_hycom_params(CS%hycom_CS, only_improves=tmpLogical) + endif + CS%use_hybgen_unmix = .false. if (coordinateMode(coord_mode) == REGRIDDING_HYBGEN) then call get_param(param_file, mdl, "USE_HYBGEN_UNMIX", CS%use_hybgen_unmix, & @@ -865,7 +873,7 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, & call build_grid_arbitrary( G, GV, h, dzInterface, trickGnuCompiler, CS ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_HYCOM1 ) - call build_grid_HyCOM1( G, GV, G%US, h, tv, h_new, dzInterface, CS, frac_shelf_h ) + call build_grid_HyCOM1( G, GV, G%US, h, tv, h_new, dzInterface, remapCS, CS, frac_shelf_h ) case ( REGRIDDING_HYBGEN ) call hybgen_regrid(G, GV, G%US, h, tv, CS%hybgen_CS, dzInterface, PCM_cell) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) @@ -1515,12 +1523,13 @@ end subroutine build_rho_grid !! \remark { Based on Bleck, 2002: An ocean-ice general circulation model framed in !! hybrid isopycnic-Cartesian coordinates, Ocean Modelling 37, 55-88. !! http://dx.doi.org/10.1016/S1463-5003(01)00012-9 } -subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS, frac_shelf_h ) +subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, remapCS, CS, frac_shelf_h ) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Existing model thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure + type(remapping_CS), intent(in) :: remapCS !< The remapping control structure type(regridding_CS), intent(in) :: CS !< Regridding control structure real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< Changes in interface position @@ -1575,7 +1584,7 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS, frac_she ( 0.5 * ( z_col(K) + z_col(K+1) ) * (GV%H_to_RZ*GV%g_Earth) - tv%P_Ref ) enddo - call build_hycom1_column(CS%hycom_CS, tv%eqn_of_state, GV%ke, nominalDepth, & + call build_hycom1_column(CS%hycom_CS, remapCS, tv%eqn_of_state, GV%ke, nominalDepth, & h(i,j,:), tv%T(i,j,:), tv%S(i,j,:), p_col, & z_col, z_col_new, zScale=GV%Z_to_H, & h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 9ebc0601d2..eeb4590a08 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -8,12 +8,14 @@ module MOM_remapping use MOM_io, only : stdout, stderr use MOM_string_functions, only : uppercase use regrid_edge_values, only : edge_values_explicit_h4, edge_values_implicit_h4 +use regrid_edge_values, only : edge_values_explicit_h4cw use regrid_edge_values, only : edge_values_implicit_h4, edge_values_implicit_h6 use regrid_edge_values, only : edge_slopes_implicit_h3, edge_slopes_implicit_h5 use remapping_attic, only : remapping_attic_unit_tests use PCM_functions, only : PCM_reconstruction use PLM_functions, only : PLM_reconstruction, PLM_boundary_extrapolation use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation +use PPM_functions, only : PPM_monotonicity use PQM_functions, only : PQM_reconstruction, PQM_boundary_extrapolation_v1 use MOM_hybgen_remap, only : hybgen_plm_coefs, hybgen_ppm_coefs, hybgen_weno_coefs @@ -48,6 +50,7 @@ module MOM_remapping integer, parameter :: REMAPPING_PCM = 0 !< O(h^1) remapping scheme integer, parameter :: REMAPPING_PLM = 2 !< O(h^2) remapping scheme integer, parameter :: REMAPPING_PLM_HYBGEN = 3 !< O(h^2) remapping scheme +integer, parameter :: REMAPPING_PPM_CW =10 !< O(h^3) remapping scheme integer, parameter :: REMAPPING_PPM_H4 = 4 !< O(h^3) remapping scheme integer, parameter :: REMAPPING_PPM_IH4 = 5 !< O(h^3) remapping scheme integer, parameter :: REMAPPING_PPM_HYBGEN = 6 !< O(h^3) remapping scheme @@ -287,7 +290,7 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & local_remapping_scheme = REMAPPING_PCM elseif (n0<=3) then local_remapping_scheme = min( local_remapping_scheme, REMAPPING_PLM ) - elseif (n0<=4) then + elseif (n0<=4 .and. local_remapping_scheme /= REMAPPING_PPM_CW ) then local_remapping_scheme = min( local_remapping_scheme, REMAPPING_PPM_H4 ) endif select case ( local_remapping_scheme ) @@ -310,6 +313,15 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & if ( CS%boundary_extrapolation ) & call PLM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) iMethod = INTEGRATION_PLM + case ( REMAPPING_PPM_CW ) + ! identical to REMAPPING_PPM_HYBGEN + call edge_values_explicit_h4cw( n0, h0, u0, ppoly_r_E, h_neglect_edge ) + call PPM_monotonicity( n0, u0, ppoly_r_E ) + call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answer_date=CS%answer_date ) + if ( CS%boundary_extrapolation ) then + call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) + endif + iMethod = INTEGRATION_PPM case ( REMAPPING_PPM_H4 ) call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answer_date=CS%answer_date ) call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answer_date=CS%answer_date ) @@ -1283,6 +1295,9 @@ subroutine setReconstructionType(string,CS) case ("PLM_HYBGEN") CS%remapping_scheme = REMAPPING_PLM_HYBGEN degree = 1 + case ("PPM_CW") + CS%remapping_scheme = REMAPPING_PPM_CW + degree = 2 case ("PPM_H4") CS%remapping_scheme = REMAPPING_PPM_H4 degree = 2 diff --git a/src/ALE/PPM_functions.F90 b/src/ALE/PPM_functions.F90 index aa24806d68..805a70d502 100644 --- a/src/ALE/PPM_functions.F90 +++ b/src/ALE/PPM_functions.F90 @@ -13,7 +13,7 @@ module PPM_functions implicit none ; private -public PPM_reconstruction, PPM_boundary_extrapolation +public PPM_reconstruction, PPM_boundary_extrapolation, PPM_monotonicity !> A tiny width that is so small that adding it to cell widths does not !! change the value due to a computational representation. It is used @@ -127,6 +127,35 @@ subroutine PPM_limiter_standard( N, h, u, edge_values, h_neglect, answer_date ) end subroutine PPM_limiter_standard +!> Adjusts edge values using the original monotonicity constraint (Colella & Woodward, JCP 1984) +!! Based on hybgen_ppm_coefs +subroutine PPM_monotonicity( N, u, edge_values ) + integer, intent(in) :: N !< Number of cells + real, dimension(:), intent(in) :: u !< cell average properties (size N) [A] + real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A] + + ! Local variables + integer :: k ! Loop index + real :: a6,da ! scalar temporaries + + ! Loop on interior cells to impose monotonicity + ! Eq. 1.10 of (Colella & Woodward, JCP 84) + do k = 2,N-1 + if (((u(k+1)-u(k))*(u(k)-u(k-1)) <= 0.)) then !local extremum + edge_values(k,1) = u(k) + edge_values(k,2) = u(k) + else + da = edge_values(k,2)-edge_values(k,1) + a6 = 6.0*u(k) - 3.0*(edge_values(k,1)+edge_values(k,2)) + if (da*a6 > da*da) then !peak in right half of zone + edge_values(k,1) = 3.0*u(k) - 2.0*edge_values(k,2) + elseif (da*a6 < -da*da) then !peak in left half of zone + edge_values(k,2) = 3.0*u(k) - 2.0*edge_values(k,1) + endif + endif + enddo ! end loop on interior cells + +end subroutine PPM_monotonicity !------------------------------------------------------------------------------ !> Reconstruction by parabolas within boundary cells diff --git a/src/ALE/coord_hycom.F90 b/src/ALE/coord_hycom.F90 index 5a3ffaff52..aa2715eb42 100644 --- a/src/ALE/coord_hycom.F90 +++ b/src/ALE/coord_hycom.F90 @@ -4,8 +4,10 @@ module coord_hycom ! 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 +use regrid_interp, only : interp_CS_type, build_and_interpolate_grid, regridding_set_ppolys +use regrid_interp, only : DEGREE_MAX implicit none ; private @@ -27,6 +29,9 @@ module coord_hycom !> Maximum thicknesses of layers [H ~> m or kg m-2] real, allocatable, dimension(:) :: max_layer_thickness + !> If true, an interface only moves if it improves the density fit + logical :: only_improves = .false. + !> Interpolation control structure type(interp_CS_type) :: interp_CS end type hycom_CS @@ -69,10 +74,11 @@ subroutine end_coord_hycom(CS) end subroutine end_coord_hycom !> This subroutine can be used to set the parameters for the coord_hycom module -subroutine set_hycom_params(CS, max_interface_depths, max_layer_thickness, interp_CS) +subroutine set_hycom_params(CS, max_interface_depths, max_layer_thickness, only_improves, interp_CS) type(hycom_CS), pointer :: CS !< Coordinate control structure real, dimension(:), optional, intent(in) :: max_interface_depths !< Maximum depths of interfaces [H ~> m or kg m-2] real, dimension(:), optional, intent(in) :: max_layer_thickness !< Maximum thicknesses of layers [H ~> m or kg m-2] + logical, optional, intent(in) :: only_improves !< If true, an interface only moves if it improves the density fit type(interp_CS_type), optional, intent(in) :: interp_CS !< Controls for interpolation if (.not. associated(CS)) call MOM_error(FATAL, "set_hycom_params: CS not associated") @@ -91,13 +97,16 @@ subroutine set_hycom_params(CS, max_interface_depths, max_layer_thickness, inter CS%max_layer_thickness(:) = max_layer_thickness(:) endif + if (present(only_improves)) CS%only_improves = only_improves + if (present(interp_CS)) CS%interp_CS = interp_CS end subroutine set_hycom_params !> Build a HyCOM coordinate column -subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, & +subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, p_col, & z_col, z_col_new, zScale, h_neglect, h_neglect_edge) type(hycom_CS), intent(in) :: CS !< Coordinate control structure + type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure integer, intent(in) :: nz !< Number of levels real, intent(in) :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2]) @@ -116,8 +125,17 @@ subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, & ! Local variables integer :: k - real, dimension(nz) :: rho_col ! Layer densities in a column [R ~> kg m-3] - real, dimension(CS%nk) :: h_col_new ! New layer thicknesses + real, dimension(nz) :: rho_col ! Layer densities in a column [R ~> kg m-3] + real, dimension(CS%nk) :: h_col_new ! New layer thicknesses [H ~> m or kg m-2] + real, dimension(CS%nk) :: r_col_new ! New layer densities [R ~> kg m-3] + real, dimension(CS%nk) :: T_col_new ! New layer temperatures [C ~> degC] + real, dimension(CS%nk) :: S_col_new ! New layer salinities [S ~> ppt] + real, dimension(CS%nk) :: p_col_new ! New layer pressure [R L2 T-2 ~> Pa] + real, dimension(CS%nk+1) :: RiA_ini ! Initial nk+1 interface density anomaly w.r.t. the + ! interface target densities [R ~> kg m-3] + real, dimension(CS%nk+1) :: RiA_new ! New interface density anomaly w.r.t. the + ! interface target densities [R ~> kg m-3] + real :: z_1, z_nz ! mid point of 1st and last layers [H ~> m or kg m-2] real :: z_scale ! A scaling factor from the input thicknesses to the target thicknesses, ! perhaps 1 or a factor in [H Z-1 ~> 1 or kg m-3] real :: stretching ! z* stretching, converts z* to z [nondim]. @@ -130,18 +148,43 @@ subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, & z_scale = 1.0 ; if (present(zScale)) z_scale = zScale - ! Work bottom recording potential density - call calculate_density(T, S, p_col, rho_col, eqn_of_state) - ! This ensures the potential density profile is monotonic - ! although not necessarily single valued. - do k = nz-1, 1, -1 - rho_col(k) = min( rho_col(k), rho_col(k+1) ) - enddo + if (CS%only_improves .and. nz == CS%nk) then + call build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, CS%nk, depth, & + h, T, S, p_col, rho_col, RiA_ini, h_neglect, h_neglect_edge) + else + ! Work bottom recording potential density + call calculate_density(T, S, p_col, rho_col, eqn_of_state) + ! This ensures the potential density profile is monotonic + ! although not necessarily single valued. + do k = nz-1, 1, -1 + rho_col(k) = min( rho_col(k), rho_col(k+1) ) + enddo + endif ! Interpolates for the target interface position with the rho_col profile ! Based on global density profile, interpolate to generate a new grid call build_and_interpolate_grid(CS%interp_CS, rho_col, nz, h(:), z_col, & CS%target_density, CS%nk, h_col_new, z_col_new, h_neglect, h_neglect_edge) + if (CS%only_improves .and. nz == CS%nk) then + ! Only move an interface if it improves the density fit + z_1 = 0.5 * ( z_col(1) + z_col(2) ) + z_nz = 0.5 * ( z_col(nz) + z_col(nz+1) ) + do k = 1,CS%nk + p_col_new(k) = p_col(1) + ( 0.5 * ( z_col_new(K) + z_col_new(K+1) ) - z_1 ) / ( z_nz - z_1 ) * & + ( p_col(nz) - p_col(1) ) + enddo + ! Remap from original h and T,S to get T,S_col_new + call remapping_core_h(remapCS, nz, h(:), T, CS%nk, h_col_new, T_col_new, h_neglect, h_neglect_edge) + call remapping_core_h(remapCS, nz, h(:), S, CS%nk, h_col_new, S_col_new, h_neglect, h_neglect_edge) + call build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, CS%nk, depth, & + h_col_new, T_col_new, S_col_new, p_col_new, r_col_new, RiA_new, h_neglect, h_neglect_edge) + do k= 2,CS%nk + if ( abs(RiA_ini(K)) <= abs(RiA_new(K)) .and. z_col(K) > z_col_new(K-1) .and. & + z_col(K) < z_col_new(K+1)) then + z_col_new(K) = z_col(K) + endif + enddo + endif !only_improves ! Sweep down the interfaces and make sure that the interface is at least ! as deep as a nominal target z* grid @@ -165,4 +208,59 @@ subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, & enddo ; endif end subroutine build_hycom1_column +!> Calculate interface density anomaly w.r.t. the target. +subroutine build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, nz, depth, h, T, S, p_col, & + R, RiAnom, h_neglect, h_neglect_edge) + type(hycom_CS), intent(in) :: CS !< Coordinate control structure + type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options + type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure + integer, intent(in) :: nz !< Number of levels + real, intent(in) :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2]) + real, dimension(nz), intent(in) :: T !< Temperature of column [C ~> degC] + real, dimension(nz), intent(in) :: S !< Salinity of column [S ~> ppt] + real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(nz), intent(in) :: p_col !< Layer pressure [R L2 T-2 ~> Pa] + !! to desired units for zInterface, perhaps GV%Z_to_H. + real, dimension(nz), intent(out) :: R !< Layer density [R ~> kg m-3] + real, dimension(nz+1), intent(out) :: RiAnom !< The interface density anomaly + !! w.r.t. the interface target + !! densities [R ~> kg m-3] + real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose of + !! cell reconstruction [H ~> m or kg m-2] + real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose of + !! edge value calculation [H ~> m or kg m-2] + ! Local variables + integer :: degree,k + real, dimension(nz) :: rho_col ! Layer densities in a column [R ~> kg m-3] + real, dimension(nz,2) :: ppoly_E ! Polynomial edge values [R ~> kg m-3] + real, dimension(nz,2) :: ppoly_S ! Polynomial edge slopes [R H-1] + real, dimension(nz,DEGREE_MAX+1) :: ppoly_C ! Polynomial interpolant coeficients on the local 0-1 grid [R ~> kg m-3] + + ! Work bottom recording potential density + call calculate_density(T, S, p_col, rho_col, eqn_of_state) + ! This ensures the potential density profile is monotonic + ! although not necessarily single valued. + do k = nz-1, 1, -1 + rho_col(k) = min( rho_col(k), rho_col(k+1) ) + enddo + + call regridding_set_ppolys(CS%interp_CS, rho_col, nz, h, ppoly_E, ppoly_S, ppoly_C, & + degree, h_neglect, h_neglect_edge) + + R(1) = rho_col(1) + RiAnom(1) = ppoly_E(1,1) - CS%target_density(1) + do k= 2,nz + R(k) = rho_col(k) + if (ppoly_E(k-1,2) > CS%target_density(k)) then + RiAnom(k) = ppoly_E(k-1,2) - CS%target_density(k) !interface is heavier than target + elseif (ppoly_E(k,1) < CS%target_density(k)) then + RiAnom(k) = ppoly_E(k,1) - CS%target_density(k) !interface is lighter than target + else + RiAnom(k) = 0.0 !interface spans the target + endif + enddo + RiAnom(nz+1) = ppoly_E(nz,2) - CS%target_density(nz+1) + +end subroutine build_hycom1_target_anomaly + end module coord_hycom diff --git a/src/ALE/regrid_edge_values.F90 b/src/ALE/regrid_edge_values.F90 index 3f59fac60f..9b574348af 100644 --- a/src/ALE/regrid_edge_values.F90 +++ b/src/ALE/regrid_edge_values.F90 @@ -14,7 +14,7 @@ module regrid_edge_values ! The following routines are visible to the outside world ! ----------------------------------------------------------------------------- public bound_edge_values, average_discontinuous_edge_values, check_discontinuous_edge_values -public edge_values_explicit_h2, edge_values_explicit_h4 +public edge_values_explicit_h2, edge_values_explicit_h4, edge_values_explicit_h4cw public edge_values_implicit_h4, edge_values_implicit_h6 public edge_slopes_implicit_h3, edge_slopes_implicit_h5 @@ -357,6 +357,106 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answer_date ) end subroutine edge_values_explicit_h4 +!> Compute h4 edge values (explicit fourth order accurate) +!! in the same units as u. +!! +!! From (Colella & Woodward, JCP, 1984) and based on hybgen_ppm_coefs. +!! +!! Compute edge values based on fourth-order explicit estimates. +!! These estimates are based on a cubic interpolant spanning four cells +!! and evaluated at the location of the middle edge. An interpolant spanning +!! cells i-2, i-1, i and i+1 is evaluated at edge i-1/2. The estimate for +!! each edge is unique. +!! +!! i-2 i-1 i i+1 +!! ..--o------o------o------o------o--.. +!! i-1/2 +!! +!! For this fourth-order scheme, at least four cells must exist. +subroutine edge_values_explicit_h4cw( N, h, u, edge_val, h_neglect ) + integer, intent(in) :: N !< Number of cells + real, dimension(N), intent(in) :: h !< cell widths [H] + real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] + real, dimension(N,2), intent(inout) :: edge_val !< Returned edge values [A]; the second index + !! is for the two edges of each cell. + real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + + ! Local variables + real :: dp(N) ! Input grid layer thicknesses, but with a minimum thickness [H ~> m or kg m-2] + real :: hNeglect ! A negligible thickness in the same units as h + real :: da ! Difference between the unlimited scalar edge value estimates [A] + real :: a6 ! Scalar field differences that are proportional to the curvature [A] + real :: slk, srk ! Differences between adjacent cell averages of scalars [A] + real :: sck ! Scalar differences across a cell. + real :: au(N) ! Scalar field difference across each cell [A] + real :: al(N), ar(N) ! Scalar field at the left and right edges of a cell [A] + real :: h112(N+1), h122(N+1) ! Combinations of thicknesses [H ~> m or kg m-2] + real :: I_h12(N+1) ! Inverses of combinations of thickesses [H-1 ~> m-1 or m2 kg-1] + real :: h2_h123(N) ! A ratio of a layer thickness of the sum of 3 adjacent thicknesses [nondim] + real :: I_h0123(N) ! Inverse of the sum of 4 adjacent thicknesses [H-1 ~> m-1 or m2 kg-1] + real :: h01_h112(N+1) ! A ratio of sums of adjacent thicknesses [nondim], 2/3 in the limit of uniform thicknesses. + real :: h23_h122(N+1) ! A ratio of sums of adjacent thicknesses [nondim], 2/3 in the limit of uniform thicknesses. + integer :: k + + hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect + + ! Set the thicknesses for very thin layers to some minimum value. + do k=1,N ; dp(k) = max(h(k), hNeglect) ; enddo + + !compute grid metrics + do k=2,N + h112(K) = 2.*dp(k-1) + dp(k) + h122(K) = dp(k-1) + 2.*dp(k) + I_h12(K) = 1.0 / (dp(k-1) + dp(k)) + enddo !k + do k=2,N-1 + h2_h123(k) = dp(k) / (dp(k) + (dp(k-1)+dp(k+1))) + enddo + do K=3,N-1 + I_h0123(K) = 1.0 / ((dp(k-2) + dp(k-1)) + (dp(k) + dp(k+1))) + + h01_h112(K) = (dp(k-2) + dp(k-1)) / (2.0*dp(k-1) + dp(k)) + h23_h122(K) = (dp(k) + dp(k+1)) / (dp(k-1) + 2.0*dp(k)) + enddo + + !Compute average slopes: Colella, Eq. (1.8) + au(1) = 0. + do k=2,N-1 + slk = u(k )-u(k-1) + srk = u(k+1)-u(k) + if (slk*srk > 0.) then + sck = h2_h123(k)*( h112(K)*srk*I_h12(K+1) + h122(K+1)*slk*I_h12(K) ) + au(k) = sign(min(abs(2.0*slk), abs(sck), abs(2.0*srk)), sck) + else + au(k) = 0. + endif + enddo !k + au(N) = 0. + + !Compute "first guess" edge values: Colella, Eq. (1.6) + al(1) = u(1) ! 1st layer PCM + ar(1) = u(1) ! 1st layer PCM + al(2) = u(1) ! 1st layer PCM + do K=3,N-1 + ! This is a 4th order explicit edge value estimate. + al(k) = (dp(k)*u(k-1) + dp(k-1)*u(k)) * I_h12(K) & + + I_h0123(K)*( 2.*dp(k)*dp(k-1)*I_h12(K)*(u(k)-u(k-1)) * & + ( h01_h112(K) - h23_h122(K) ) & + + (dp(k)*au(k-1)*h23_h122(K) - dp(k-1)*au(k)*h01_h112(K)) ) + ar(k-1) = al(k) + enddo !k + ar(N-1) = u(N) ! last layer PCM + al(N) = u(N) ! last layer PCM + ar(N) = u(N) ! last layer PCM + + !Set coefficients + do k=1,N + edge_val(k,1) = al(k) + edge_val(k,2) = ar(k) + enddo !k + +end subroutine edge_values_explicit_h4cw + !> Compute ih4 edge values (implicit fourth order accurate) !! in the same units as u. !! diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index 4d09daf6f3..e119ce9d53 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -7,11 +7,13 @@ module regrid_interp use MOM_string_functions, only : uppercase use regrid_edge_values, only : edge_values_explicit_h2, edge_values_explicit_h4 +use regrid_edge_values, only : edge_values_explicit_h4cw use regrid_edge_values, only : edge_values_implicit_h4, edge_values_implicit_h6 use regrid_edge_values, only : edge_slopes_implicit_h3, edge_slopes_implicit_h5 use PLM_functions, only : PLM_reconstruction, PLM_boundary_extrapolation use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation +use PPM_functions, only : PPM_monotonicity use PQM_functions, only : PQM_reconstruction, PQM_boundary_extrapolation_v1 use P1M_functions, only : P1M_interpolation, P1M_boundary_extrapolation @@ -45,6 +47,7 @@ module regrid_interp integer, parameter :: INTERPOLATION_P1M_H4 = 1 !< O(h^2) integer, parameter :: INTERPOLATION_P1M_IH4 = 2 !< O(h^2) integer, parameter :: INTERPOLATION_PLM = 3 !< O(h^2) +integer, parameter :: INTERPOLATION_PPM_CW =10 !< O(h^3) integer, parameter :: INTERPOLATION_PPM_H4 = 4 !< O(h^3) integer, parameter :: INTERPOLATION_PPM_IH4 = 5 !< O(h^3) integer, parameter :: INTERPOLATION_P3M_IH4IH3 = 6 !< O(h^4) @@ -144,6 +147,25 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & call PLM_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) endif + case ( INTERPOLATION_PPM_CW ) + if ( n0 >= 4 ) then + degree = DEGREE_2 + call edge_values_explicit_h4cw( n0, h0, densities, ppoly0_E, h_neglect_edge ) + call PPM_monotonicity( n0, densities, ppoly0_E ) + call PPM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date ) + if (extrapolate) then + call PPM_boundary_extrapolation( n0, h0, densities, ppoly0_E, & + ppoly0_coefs, h_neglect ) + endif + else + degree = DEGREE_1 + call edge_values_explicit_h2( n0, h0, densities, ppoly0_E ) + call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date ) + if (extrapolate) then + call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs ) + endif + endif + case ( INTERPOLATION_PPM_H4 ) if ( n0 >= 4 ) then degree = DEGREE_2 @@ -486,7 +508,7 @@ end function get_polynomial_coordinate !> 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 - !! Valid values include "P1M_H2", "P1M_H4", "P1M_IH2", "PLM", "PPM_H4", + !! Valid values include "P1M_H2", "P1M_H4", "P1M_IH2", "PLM", "PPM_CW", "PPM_H4", !! "PPM_IH4", "P3M_IH4IH3", "P3M_IH6IH5", "PQM_IH4IH3", and "PQM_IH6IH5" select case ( uppercase(trim(interp_scheme)) ) @@ -494,6 +516,7 @@ integer function interpolation_scheme(interp_scheme) case ("P1M_H4"); interpolation_scheme = INTERPOLATION_P1M_H4 case ("P1M_IH2"); interpolation_scheme = INTERPOLATION_P1M_IH4 case ("PLM"); interpolation_scheme = INTERPOLATION_PLM + case ("PPM_CW"); interpolation_scheme = INTERPOLATION_PPM_CW case ("PPM_H4"); interpolation_scheme = INTERPOLATION_PPM_H4 case ("PPM_IH4"); interpolation_scheme = INTERPOLATION_PPM_IH4 case ("P3M_IH4IH3"); interpolation_scheme = INTERPOLATION_P3M_IH4IH3 @@ -509,7 +532,7 @@ end function interpolation_scheme subroutine set_interp_scheme(CS, interp_scheme) type(interp_CS_type), intent(inout) :: CS !< A control structure for regrid_interp character(len=*), intent(in) :: interp_scheme !< Name of the interpolation scheme - !! Valid values include "P1M_H2", "P1M_H4", "P1M_IH2", "PLM", "PPM_H4", + !! Valid values include "P1M_H2", "P1M_H4", "P1M_IH2", "PLM", "PPM_CW", "PPM_H4", !! "PPM_IH4", "P3M_IH4IH3", "P3M_IH6IH5", "PQM_IH4IH3", and "PQM_IH6IH5" CS%interpolation_scheme = interpolation_scheme(interp_scheme)