Skip to content

Commit

Permalink
Add Jan's calving to sandbox
Browse files Browse the repository at this point in the history
Co-authored-by: Jan Malles <[email protected]>
  • Loading branch information
fmaussion and jmalles committed Feb 4, 2024
1 parent dd3c1b3 commit eae1b85
Show file tree
Hide file tree
Showing 2 changed files with 22 additions and 9 deletions.
21 changes: 16 additions & 5 deletions oggm/sandbox/calving.py
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@


class CalvingFluxBasedModel_v1(FlowlineModel):
"""
"""Jan Malles' implementation of calving as in his paper.
"""

def __init__(self, flowlines, mb_model=None, y0=0., glen_a=None,
Expand Down Expand Up @@ -754,6 +754,10 @@ def get_diagnostics(self, fl_id=-1):
df['ice_velocity'] = (var[1:nx+1] + var[:nx])/2 + _u_slide
df['surface_ice_velocity'] = (((var[1:nx+1] + var[:nx])/2 *
self._surf_vel_fac) + _u_slide)
var = self.u_drag[fl_id]
df['deformation_velocity'] = (var[1:nx+1] + var[:nx])/2
var = self.u_slide[fl_id]
df['sliding_velocity'] = (var[1:nx+1] + var[:nx])/2
var = self.shapefac_stag[fl_id]
df['shape_fac'] = (var[1:nx+1] + var[:nx])/2

Expand All @@ -765,7 +769,8 @@ def get_diagnostics(self, fl_id=-1):

class CalvingFluxBasedModel_v2(FlowlineModel):
"""This is the version currently in progress - written by Fabien based
on Jan's code.
on Jan's code. Code is prettier but does not exactly do the same at the
moment.
"""

def __init__(self, flowlines, mb_model=None, y0=0., glen_a=None,
Expand Down Expand Up @@ -1038,7 +1043,7 @@ def step(self, dt):
# Check that the "last_above_wl" is not just the last in a
# "lake" (over-deepening) which is followed by land again.
# If after last_above_wl we still have ice, we don't calve.
calving_is_happening = fl.thick[last_above_wl + 1] > 0
calving_is_happening = not fl.thick[last_above_wl + 1] > 0

# Determine water depth at the front
h = fl.thick[last_above_wl]
Expand Down Expand Up @@ -1086,7 +1091,7 @@ def step(self, dt):

# Add "stretching stress" to basal shear/driving stress
if calving_is_happening:
stress[stretch_first:stretch_last] += stretch_factor * (pull_last / stretch_dist)
stress[stretch_first:stretch_last+1] += stretch_factor * (pull_last / stretch_dist)

# Compute velocities
# Deformation is the usual formulation with changed stress
Expand All @@ -1095,7 +1100,9 @@ def step(self, dt):
# Sliding is increased where there is water
# Determine height above buoyancy
eff_water_depth_stag = (rho_ocean / self.rho) * water_depth_stag
# Avoid dividing by zero where thick equals water depth
z_a_b = utils.clip_min(thick_stag - eff_water_depth_stag, 0.01)
z_a_b[thick_stag == 0] = 1 # Stress is zero there
us_stag[:] = (stress ** N / z_a_b) * self.fs

# Force velocity beyond grounding line to be zero in order to
Expand Down Expand Up @@ -1138,9 +1145,9 @@ def step(self, dt):

# Temporary check for above
z_a_b = thick_stag
z_a_b[thick_stag == 0] = 1
vel = (stress ** N / z_a_b) * self.fs
assert np.allclose(us_stag, vel)

u_stag[:] = ud_stag + us_stag

# Staggered flux rate
Expand Down Expand Up @@ -1406,6 +1413,10 @@ def get_diagnostics(self, fl_id=-1):
df['ice_flux'] = (var[1:nx+1] + var[:nx])/2
var = self.u_stag[fl_id]
df['ice_velocity'] = (var[1:nx+1] + var[:nx])/2
var = self.ud_stag[fl_id]
df['deformation_velocity'] = (var[1:nx+1] + var[:nx])/2
var = self.us_stag[fl_id]
df['sliding_velocity'] = (var[1:nx+1] + var[:nx])/2

# Surface vel is deformation corrected by factor + sliding
var = self.ud_stag[fl_id]
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -156,11 +156,13 @@ def test_flux_gate_with_calving(self):
if do_plot: # pragma: no cover
plt.plot(model.fls[-1].surface_h, label=model_class.__name__)

np.testing.assert_allclose(vols[1], vols[2])
np.testing.assert_allclose(vols[1], vols[0], rtol=0.01) # Normal!
# Jan's is different than the original one
np.testing.assert_allclose(vols[1], vols[0], rtol=0.01)
np.testing.assert_allclose(calvs[1], calvs[0], rtol=0.12)

np.testing.assert_allclose(calvs[1], calvs[2])
np.testing.assert_allclose(calvs[1], calvs[0], rtol=0.12) # Normal!
# Now between my implementation and Jans
# np.testing.assert_allclose(vols[1], vols[2])
# np.testing.assert_allclose(calvs[1], calvs[2])

if do_plot: # pragma: no cover
plt.legend()
Expand Down

0 comments on commit eae1b85

Please sign in to comment.