Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add PPM_CM and HYCOM1_ONLY_IMPROVES #341

Merged
merged 2 commits into from
Mar 30, 2023
Merged

Add PPM_CM and HYCOM1_ONLY_IMPROVES #341

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a17fb86
Add PPM_CM and HYCOM1_ONLY_IMPROVES
awallcraft Mar 19, 2023
3fdc31e
Merge branch 'dev/gfdl' into PPMCW
marshallward Mar 30, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/ALE/MOM_ALE.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 "//&
Expand Down
17 changes: 13 additions & 4 deletions src/ALE/MOM_regridding.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"//&
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
17 changes: 16 additions & 1 deletion src/ALE/MOM_remapping.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 )
Expand All @@ -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 )
Expand Down Expand Up @@ -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
Expand Down
31 changes: 30 additions & 1 deletion src/ALE/PPM_functions.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading