Skip to content

Do not zero N under floating ice in Interface #794

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.

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

Open
wants to merge 4 commits into
base: landice/develop
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
52 changes: 52 additions & 0 deletions Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -454,6 +454,58 @@ intel-nersc:
"OPENMP = $(OPENMP)" \
"CPPFLAGS = $(MODEL_FORMULATION) -D_MPI" )

intel-cori-knl:
( $(MAKE) all \
"FC_PARALLEL = ftn" \
"CC_PARALLEL = cc" \
"CXX_PARALLEL = CC" \
"FC_SERIAL = ftn" \
"CC_SERIAL = cc" \
"CXX_SERIAL = CC" \
"FFLAGS_PROMOTION = -real-size 64" \
"FFLAGS_OPT = -O3 -convert big_endian -free -align array64byte" \
"CFLAGS_OPT = -O3" \
"CXXFLAGS_OPT = -O3, -xMIC-AVX512" \
"LDFLAGS_OPT = -O3" \
"FFLAGS_OMP = -qopenmp" \
"CFLAGS_OMP = -qopenmp" \
"FFLAGS_DEBUG = -real-size 64 -g -convert big_endian -free -CU -CB -check all -gen-interfaces -warn interfaces -traceback" \
"CFLAGS_DEBUG = -g -traceback" \
"CXXFLAGS_DEBUG = -g -traceback, -xMIC-AVX512" \
"LDFLAGS_DEBUG = -g -traceback" \
"BUILD_TARGET = $(@)" \
"CORE = $(CORE)" \
"DEBUG = $(DEBUG)" \
"USE_PAPI = $(USE_PAPI)" \
"OPENMP = $(OPENMP)" \
"CPPFLAGS = $(MODEL_FORMULATION) -D_MPI" )

intel-cori-haswell:
( $(MAKE) all \
"FC_PARALLEL = ftn" \
"CC_PARALLEL = cc" \
"CXX_PARALLEL = CC" \
"FC_SERIAL = ftn" \
"CC_SERIAL = cc" \
"CXX_SERIAL = CC" \
"FFLAGS_PROMOTION = -real-size 64" \
"FFLAGS_OPT = -O3 -convert big_endian -free -align array64byte" \
"CFLAGS_OPT = -O3" \
"CXXFLAGS_OPT = -O3, -xCORE-AVX2" \
"LDFLAGS_OPT = -O3" \
"FFLAGS_OMP = -qopenmp" \
"CFLAGS_OMP = -qopenmp" \
"FFLAGS_DEBUG = -real-size 64 -g -convert big_endian -free -CU -CB -check all -gen-interfaces -warn interfaces -traceback" \
"CFLAGS_DEBUG = -g -traceback" \
"CXXFLAGS_DEBUG = -g -traceback, -xCORE-AVX2" \
"LDFLAGS_DEBUG = -g -traceback" \
"BUILD_TARGET = $(@)" \
"CORE = $(CORE)" \
"DEBUG = $(DEBUG)" \
"USE_PAPI = $(USE_PAPI)" \
"OPENMP = $(OPENMP)" \
"CPPFLAGS = $(MODEL_FORMULATION) -D_MPI" )

bluegene:
( $(MAKE) all \
"FC_PARALLEL = mpixlf95_r" \
Expand Down
4 changes: 4 additions & 0 deletions src/core_landice/Registry.xml
Original file line number Diff line number Diff line change
Expand Up @@ -100,6 +100,10 @@
description="If true, nonconvergence of the velocity solver results in a fatal error. If false, it is a warning."
possible_values=".true. or .false."
/>
<nml_option name="config_extend_effective_pressure" type="logical" default_value=".true." units="unitless"
description="If true, the field effectivePressure is extended by one cell before being passed to the velocity solver."
possible_values=".true. or .false."
/>
<nml_option name="config_effective_pressure_max" type="real" default_value="1.0e36" units="Pa"
description="An optional maximum limit on effective pressure passed to the velocity solver."
possible_values="Any positive, real number."
Expand Down
2 changes: 0 additions & 2 deletions src/core_landice/mode_forward/Interface_velocity_solver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1252,8 +1252,6 @@ void importFields(std::vector<std::pair<int, int> >& marineBdyExtensionMap, dou
}

bool floating = rho_ice * thicknessData[iv] + rho_ocean * bedTopographyData[iv] < 0;
if (floating && (effecPress_F != 0))
effecPressData[iv] = 0.0;
}
}

Expand Down
51 changes: 51 additions & 0 deletions src/core_landice/mode_forward/mpas_li_velocity.F
Original file line number Diff line number Diff line change
Expand Up @@ -808,13 +808,20 @@ subroutine calculate_beta(block, err)
type (mpas_pool_type), pointer :: thermalPool
type (mpas_pool_type), pointer :: meshPool
real (kind=RKIND), dimension(:), pointer :: betaSolve, beta
real (kind=RKIND), dimension(:), pointer :: effectivePressure, effectivePressureLimited
real (kind=RKIND), dimension(:), pointer :: basalTemperature, basalPmpTemperature
real (kind=RKIND), dimension(:), pointer :: muFriction
integer, dimension(:), pointer :: cellMask
logical, pointer :: config_use_glp
logical, pointer :: config_beta_thawed_only
logical, pointer :: config_extend_effective_pressure
real (kind=RKIND), pointer :: config_effective_pressure_max
logical, pointer :: hydroActive
real (kind=RKIND) :: betaAccum
integer :: nBetaValues
real (kind=RKIND) :: NAccum
real (kind=RKIND) :: Nmin
integer :: nNValues
integer :: iCell, iCell2
integer :: neighbor
integer, pointer :: nCells
Expand All @@ -825,15 +832,25 @@ subroutine calculate_beta(block, err)

call mpas_pool_get_config(liConfigs, 'config_use_glp', config_use_glp)
call mpas_pool_get_config(liConfigs, 'config_beta_thawed_only', config_beta_thawed_only)
call mpas_pool_get_config(liConfigs, 'config_extend_effective_pressure', config_extend_effective_pressure)
call mpas_pool_get_config(liConfigs, 'config_effective_pressure_max', config_effective_pressure_max)

call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)
call mpas_pool_get_subpool(block % structs, 'thermal', thermalPool)
call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)

call mpas_pool_get_array(velocityPool, 'betaSolve', betaSolve)
call mpas_pool_get_array(velocityPool, 'beta', beta)
call mpas_pool_get_array(velocityPool, 'muFriction', muFriction)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
call mpas_pool_get_array(velocityPool, 'effectivePressureLimited', effectivePressureLimited)

call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)

betaSolve = beta

Expand All @@ -846,6 +863,40 @@ subroutine calculate_beta(block, err)
end where
endif

! Effective pressure adjustments for velocity solver
effectivePressureLimited(:) = effectivePressure(:)

if (config_extend_effective_pressure) then
do iCell = 1, nCells
if ( .not. li_mask_is_grounded_ice(cellMask(iCell)) ) then
nNValues = 0
NAccum = 0.0_RKIND
Nmin = 1.0e20
do iCell2 = 1, nEdgesOnCell(iCell)
neighbor = cellsOnCell(iCell2, iCell)
if (li_mask_is_grounded_ice(cellMask(neighbor))) then
NAccum = NAccum + effectivePressureLimited(neighbor)
nNValues = nNValues + 1
Nmin = min(Nmin, effectivePressureLimited(neighbor))
endif
enddo
if (nNValues > 0) then
! effectivePressureLimited(iCell) = NAccum / real(nNValues, kind=RKIND) ! take average of neighboring effectivePressure values
! < TODO: could do this in log space
effectivePressureLimited(iCell) = Nmin
endif
endif ! if iCell is floating
enddo ! cell loop
endif

! Apply limits to N
do iCell = 1, nCells
effectivePressureLimited(iCell) = max(0.0_RKIND, effectivePressureLimited(iCell)) ! Never pass a negative N
effectivePressureLimited(iCell) = min(config_effective_pressure_max, effectivePressureLimited(iCell))
enddo
effectivePressureLimited(:) = effectivePressureLimited(:) * muFriction(:)


!--------------------------------------------------------------------
end subroutine calculate_beta

Expand Down
9 changes: 0 additions & 9 deletions src/core_landice/mode_forward/mpas_li_velocity_external.F
Original file line number Diff line number Diff line change
Expand Up @@ -449,7 +449,6 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
normalVelocity, uReconstructX, uReconstructY, uReconstructZ
real (kind=RKIND), dimension(:,:), pointer :: temperature
real (kind=RKIND), dimension(:), pointer :: stiffnessFactor
real (kind=RKIND), dimension(:), pointer :: effectivePressure
real (kind=RKIND), dimension(:), pointer :: effectivePressureLimited
real (kind=RKIND), dimension(:), pointer :: muFriction
real (kind=RKIND), pointer :: deltat
Expand All @@ -460,7 +459,6 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
logical, pointer :: config_output_external_velocity_solver_data
logical, pointer :: config_nonconvergence_error
real (kind=RKIND), pointer :: config_ice_density
real (kind=RKIND), pointer :: config_effective_pressure_max
integer, pointer :: anyDynamicVertexMaskChanged
integer, pointer :: dirichletMaskChanged
integer, pointer :: nEdges
Expand All @@ -486,7 +484,6 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
config_output_external_velocity_solver_data)
call mpas_pool_get_config(liConfigs, 'config_nonconvergence_error', config_nonconvergence_error)
call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density)
call mpas_pool_get_config(liConfigs, 'config_effective_pressure_max', config_effective_pressure_max)

! Mesh variables
call mpas_pool_get_array(meshPool, 'layerThicknessFractions', layerThicknessFractions)
Expand Down Expand Up @@ -524,19 +521,13 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)

! Hydro variables
call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
call mpas_pool_get_array(velocityPool, 'effectivePressureLimited', effectivePressureLimited)

#if defined(USE_EXTERNAL_L1L2) || defined(USE_EXTERNAL_FIRSTORDER) || defined(USE_EXTERNAL_STOKES)
! Capture Albany output
call interface_redirect_stdout(timestepNumber)
#endif

! Create version of effectivePressure field with limits applied
do iCell = 1, nCells
effectivePressureLimited(iCell) = max(0.0_RKIND, effectivePressure(iCell)) ! Never pass a negative N
effectivePressureLimited(iCell) = min(config_effective_pressure_max, effectivePressureLimited(iCell))
enddo

! ==================================================================
! External dycore calls to be made only when vertex mask changes
Expand Down