Skip to content

Commit

Permalink
Merge branch 'Hallberg-NOAA-ALE_remap_fix' into dev/gfdl
Browse files Browse the repository at this point in the history
  • Loading branch information
@adcroft
adcroft committed Dec 3, 2019
2 parents bfa303d + 4158a7e commit c9a7f98
Show file tree
Hide file tree
Showing 6 changed files with 417 additions and 355 deletions.
43 changes: 29 additions & 14 deletions src/ALE/MOM_ALE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,19 +63,23 @@ module MOM_ALE

!> ALE control structure
type, public :: ALE_CS ; private
logical :: remap_uv_using_old_alg !< If true, uses the old "remapping via a delta z"
!! method. If False, uses the new method that
!! remaps between grids described by h.
logical :: remap_uv_using_old_alg !< If true, uses the old "remapping via a delta z"
!! method. If False, uses the new method that
!! remaps between grids described by h.

real :: regrid_time_scale !< The time-scale used in blending between the current (old) grid
!! and the target (new) grid [T ~> s]

type(regridding_CS) :: regridCS !< Regridding parameters and work arrays
type(remapping_CS) :: remapCS !< Remapping parameters and work arrays

integer :: nk !< Used only for queries, not directly by this module
integer :: nk !< Used only for queries, not directly by this module

logical :: remap_after_initialization !< Indicates whether to regrid/remap after initializing the state.
logical :: remap_after_initialization !< Indicates whether to regrid/remap after initializing the state.

logical :: answers_2018 !< If true, use the order of arithmetic and expressions for remapping
!! that recover the answers from the end of 2018. Otherwise, use more
!! robust and accurate forms of mathematically equivalent expressions.

logical :: show_call_tree !< For debugging

Expand Down Expand Up @@ -145,6 +149,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
character(len=40) :: mdl = "MOM_ALE" ! This module's name.
character(len=80) :: string ! Temporary strings
real :: filter_shallow_depth, filter_deep_depth
logical :: default_2018_answers
logical :: check_reconstruction
logical :: check_remapping
logical :: force_bounds_in_subcell
Expand Down Expand Up @@ -192,11 +197,19 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
call get_param(param_file, mdl, "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &
"If true, values at the interfaces of boundary cells are "//&
"extrapolated instead of piecewise constant", default=.false.)
call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
"This sets the default value for the various _2018_ANSWERS parameters.", &
default=.true.)
call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", CS%answers_2018, &
"If true, use the order of arithmetic and expressions that recover the "//&
"answers from the end of 2018. Otherwise, use updated and more robust "//&
"forms of the same expressions.", default=default_2018_answers)
call initialize_remapping( CS%remapCS, string, &
boundary_extrapolation=remap_boundary_extrap, &
check_reconstruction=check_reconstruction, &
check_remapping=check_remapping, &
force_bounds_in_subcell=force_bounds_in_subcell)
force_bounds_in_subcell=force_bounds_in_subcell, &
answers_2018=CS%answers_2018)

call get_param(param_file, mdl, "REMAP_AFTER_INITIALIZATION", CS%remap_after_initialization, &
"If true, applies regridding and remapping immediately after "//&
Expand All @@ -220,7 +233,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
"REGRID_FILTER_SHALLOW_DEPTH the filter weights adopt a cubic profile.", &
units="m", default=0., scale=GV%m_to_H)
call set_regrid_params(CS%regridCS, depth_of_time_filter_shallow=filter_shallow_depth, &
depth_of_time_filter_deep=filter_deep_depth)
depth_of_time_filter_deep=filter_deep_depth)
call get_param(param_file, mdl, "REGRID_USE_OLD_DIRECTION", local_logical, &
"If true, the regridding ntegrates upwards from the bottom for "//&
"interface positions, much as the main model does. If false "//&
Expand Down Expand Up @@ -1089,13 +1102,13 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_ext

! Local variables
integer :: i, j, k
real :: hTmp(GV%ke)
real :: tmp(GV%ke)
real :: hTmp(GV%ke) ! A 1-d copy of h [H ~> m or kg m-2]
real :: tmp(GV%ke) ! A 1-d copy of a column of temperature [degC] or salinity [ppt]
real, dimension(CS%nk,2) :: &
ppol_E !Edge value of polynomial
ppol_E ! Edge value of polynomial in [degC] or [ppt]
real, dimension(CS%nk,3) :: &
ppol_coefs !Coefficients of polynomial
real :: h_neglect, h_neglect_edge
ppol_coefs ! Coefficients of polynomial, all in [degC] or [ppt]
real :: h_neglect, h_neglect_edge ! Tiny thicknesses [H ~> m or kg m-2]

!### Try replacing both of these with GV%H_subroundoff
if (GV%Boussinesq) then
Expand All @@ -1116,7 +1129,8 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_ext
ppol_E(:,:) = 0.0
ppol_coefs(:,:) = 0.0
!### Try to replace the following value of h_neglect with GV%H_subroundoff.
call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=h_neglect_edge )
call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=h_neglect_edge, &
answers_2018=CS%answers_2018 )
call PPM_reconstruction( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
if (bdry_extrap) &
call PPM_boundary_extrapolation( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
Expand All @@ -1131,7 +1145,8 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_ext
ppol_coefs(:,:) = 0.0
tmp(:) = tv%T(i,j,:)
!### Try to replace the following value of h_neglect with GV%H_subroundoff.
call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=1.0e-10*GV%m_to_H )
call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=1.0e-10*GV%m_to_H, &
answers_2018=CS%answers_2018 )
call PPM_reconstruction( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
if (bdry_extrap) &
call PPM_boundary_extrapolation(GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
Expand Down
37 changes: 23 additions & 14 deletions src/ALE/MOM_remapping.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,8 @@ module MOM_remapping
logical :: check_remapping = .false.
!> If true, the intermediate values used in remapping are forced to be bounded.
logical :: force_bounds_in_subcell = .false.
!> If true use older, less acccurate expressions.
logical :: answers_2018 = .true.
end type

! The following routines are visible to the outside world
Expand Down Expand Up @@ -84,13 +86,14 @@ module MOM_remapping

!> Set parameters within remapping object
subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation, &
check_reconstruction, check_remapping, force_bounds_in_subcell)
check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
type(remapping_CS), intent(inout) :: CS !< Remapping control structure
character(len=*), optional, intent(in) :: remapping_scheme !< Remapping scheme to use
logical, optional, intent(in) :: boundary_extrapolation !< Indicate to extrapolate in boundary cells
logical, optional, intent(in) :: check_reconstruction !< Indicate to check reconstructions
logical, optional, intent(in) :: check_remapping !< Indicate to check results of remapping
logical, optional, intent(in) :: force_bounds_in_subcell !< Force subcells values to be bounded
logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.

if (present(remapping_scheme)) then
call setReconstructionType( remapping_scheme, CS )
Expand All @@ -107,6 +110,9 @@ subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation, &
if (present(force_bounds_in_subcell)) then
CS%force_bounds_in_subcell = force_bounds_in_subcell
endif
if (present(answers_2018)) then
CS%answers_2018 = answers_2018
endif
end subroutine remapping_set_param

subroutine extract_member_remapping_CS(CS, remapping_scheme, degree, boundary_extrapolation, check_reconstruction, &
Expand Down Expand Up @@ -392,31 +398,31 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, &
endif
iMethod = INTEGRATION_PLM
case ( REMAPPING_PPM_H4 )
call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge )
call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 )
call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect )
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_IH4 )
call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge )
call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 )
call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect )
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_PQM_IH4IH3 )
call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge )
call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_S, h_neglect )
call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 )
call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_S, h_neglect, answers_2018=CS%answers_2018 )
call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefs, h_neglect )
if ( CS%boundary_extrapolation ) then
call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, &
ppoly_r_coefs, h_neglect )
endif
iMethod = INTEGRATION_PQM
case ( REMAPPING_PQM_IH6IH5 )
call edge_values_implicit_h6( n0, h0, u0, ppoly_r_E, h_neglect_edge )
call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_S, h_neglect )
call edge_values_implicit_h6( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 )
call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_S, h_neglect, answers_2018=CS%answers_2018 )
call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefs, h_neglect )
if ( CS%boundary_extrapolation ) then
call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, &
Expand Down Expand Up @@ -1537,19 +1543,20 @@ end subroutine dzFromH1H2

!> Constructor for remapping control structure
subroutine initialize_remapping( CS, remapping_scheme, boundary_extrapolation, &
check_reconstruction, check_remapping, force_bounds_in_subcell)
check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
! Arguments
type(remapping_CS), intent(inout) :: CS !< Remapping control structure
character(len=*), intent(in) :: remapping_scheme !< Remapping scheme to use
logical, optional, intent(in) :: boundary_extrapolation !< Indicate to extrapolate in boundary cells
logical, optional, intent(in) :: check_reconstruction !< Indicate to check reconstructions
logical, optional, intent(in) :: check_remapping !< Indicate to check results of remapping
logical, optional, intent(in) :: force_bounds_in_subcell !< Force subcells values to be bounded
logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.

! Note that remapping_scheme is mandatory fir initialize_remapping()
! Note that remapping_scheme is mandatory for initialize_remapping()
call remapping_set_param(CS, remapping_scheme=remapping_scheme, boundary_extrapolation=boundary_extrapolation, &
check_reconstruction=check_reconstruction, check_remapping=check_remapping, &
force_bounds_in_subcell=force_bounds_in_subcell)
force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=answers_2018)

end subroutine initialize_remapping

Expand Down Expand Up @@ -1615,13 +1622,15 @@ logical function remapping_unit_tests(verbose)
data h2 /6*0.5/ ! 6 uniform layers with total depth of 3
type(remapping_CS) :: CS !< Remapping control structure
real, allocatable, dimension(:,:) :: ppoly0_E, ppoly0_S, ppoly0_coefs
logical :: answers_2018 ! If true use older, less acccurate expressions.
integer :: i
real :: err, h_neglect, h_neglect_edge
logical :: thisTest, v

v = verbose
h_neglect = hNeglect_dflt
h_neglect_edge = 1.0e-10
answers_2018 = .true.

write(*,*) '==== MOM_remapping: remapping_unit_tests ================='
remapping_unit_tests = .false. ! Normally return false
Expand All @@ -1643,7 +1652,7 @@ logical function remapping_unit_tests(verbose)
remapping_unit_tests = remapping_unit_tests .or. thisTest

thisTest = .false.
call initialize_remapping(CS, 'PPM_H4')
call initialize_remapping(CS, 'PPM_H4', answers_2018=answers_2018)
if (verbose) write(*,*) 'h0 (test data)'
if (verbose) call dumpGrid(n0,h0,x0,u0)

Expand All @@ -1667,7 +1676,7 @@ logical function remapping_unit_tests(verbose)
ppoly0_S(:,:) = 0.0
ppoly0_coefs(:,:) = 0.0

call edge_values_explicit_h4( n0, h0, u0, ppoly0_E, h_neglect=1e-10 )
call edge_values_explicit_h4( n0, h0, u0, ppoly0_E, h_neglect=1e-10, answers_2018=answers_2018 )
call PPM_reconstruction( n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect )
call PPM_boundary_extrapolation( n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect )
u1(:) = 0.
Expand Down Expand Up @@ -1798,7 +1807,7 @@ logical function remapping_unit_tests(verbose)
test_answer(v, 3, ppoly0_coefs(:,2), (/0.,4.,0./), 'Non-uniform line PLM: P1')

call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), ppoly0_E, &
h_neglect=1e-10 )
h_neglect=1e-10, answers_2018=answers_2018 )
! The next two tests currently fail due to roundoff.
thisTest = test_answer(v, 5, ppoly0_E(:,1), (/0.,2.,4.,6.,8./), 'Line H4: left edges')
thisTest = test_answer(v, 5, ppoly0_E(:,2), (/2.,4.,6.,8.,10./), 'Line H4: right edges')
Expand All @@ -1814,7 +1823,7 @@ logical function remapping_unit_tests(verbose)
test_answer(v, 5, ppoly0_coefs(:,3), (/0.,0.,0.,0.,0./), 'Line PPM: P2')

call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,1.,7.,19.,37./), ppoly0_E, &
h_neglect=1e-10 )
h_neglect=1e-10, answers_2018=answers_2018 )
! The next two tests currently fail due to roundoff.
thisTest = test_answer(v, 5, ppoly0_E(:,1), (/3.,0.,3.,12.,27./), 'Parabola H4: left edges')
thisTest = test_answer(v, 5, ppoly0_E(:,2), (/0.,3.,12.,27.,48./), 'Parabola H4: right edges')
Expand Down
Loading

0 comments on commit c9a7f98

Please sign in to comment.