Skip to content

Clean up iceshelf framework #816

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

Merged
xylar merged 2 commits into from
May 8, 2024
Merged

Clean up iceshelf framework #816

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
177c265
Clean up adjust_ssh
xylar Apr 7, 2024
f161a60
Save the total ssh and landIcePressure change
xylar Apr 7, 2024
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
77 changes: 50 additions & 27 deletions compass/ocean/iceshelf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,11 +79,21 @@ def adjust_ssh(variable, iteration_count, step, update_pio=True,
step.update_namelist_pio('namelist.ocean')
partition(ntasks, config, logger)

for iterIndex in range(iteration_count):
logger.info(f" * Iteration {iterIndex + 1}/{iteration_count}")
with xarray.open_dataset('adjusting_init0.nc') as ds:
ds = ds.isel(Time=0)
orig_ssh = ds.ssh
orig_land_ice_pressure = ds.landIcePressure

in_filename = f'adjusting_init{iterIndex}.nc'
out_filename = f'adjusting_init{iterIndex + 1}.nc'
on_a_sphere = ds.attrs['on_a_sphere'].lower() == 'yes'

modify_mask = numpy.logical_and(ds.maxLevelCell > 0,
ds.modifyLandIcePressureMask == 1)

for iter_index in range(iteration_count):
logger.info(f" * Iteration {iter_index + 1}/{iteration_count}")

in_filename = f'adjusting_init{iter_index}.nc'
out_filename = f'adjusting_init{iter_index + 1}.nc'
symlink(in_filename, 'adjusting_init.nc')

logger.info(" * Running forward model")
Expand All @@ -99,37 +109,32 @@ def adjust_ssh(variable, iteration_count, step, update_pio=True,

ds = ds.isel(Time=0)

on_a_sphere = ds.attrs['on_a_sphere'].lower() == 'yes'

initSSH = ds.ssh
init_ssh = ds.ssh
if 'minLevelCell' in ds:
minLevelCell = ds.minLevelCell - 1
else:
minLevelCell = xarray.zeros_like(ds.maxLevelCell)

with xarray.open_dataset('output_ssh.nc') as ds_ssh:
# get the last time entry
ds_ssh = ds_ssh.isel(Time=ds_ssh.sizes['Time'] - 1)
finalSSH = ds_ssh.ssh
ds_ssh = ds_ssh.isel(Time=-1)
final_ssh = ds_ssh.ssh
topDensity = ds_ssh.density.isel(nVertLevels=minLevelCell)

mask = numpy.logical_and(ds.maxLevelCell > 0,
ds.modifyLandIcePressureMask == 1)

deltaSSH = mask * (finalSSH - initSSH)
delta_ssh = modify_mask * (final_ssh - init_ssh)

# then, modify the SSH or land-ice pressure
if variable == 'ssh':
ssh = finalSSH.expand_dims(dim='Time', axis=0)
ssh = final_ssh.expand_dims(dim='Time', axis=0)
ds_out['ssh'] = ssh
# also update the landIceDraft variable, which will be used to
# compensate for the SSH due to land-ice pressure when
# computing sea-surface tilt
ds_out['landIceDraft'] = ssh
# we also need to stretch layerThickness to be compatible with
# the new SSH
stretch = ((finalSSH + ds.bottomDepth) /
(initSSH + ds.bottomDepth))
stretch = ((final_ssh + ds.bottomDepth) /
(init_ssh + ds.bottomDepth))
ds_out['layerThickness'] = ds_out.layerThickness * stretch
landIcePressure = ds.landIcePressure.values
else:
Expand All @@ -140,15 +145,15 @@ def adjust_ssh(variable, iteration_count, step, update_pio=True,
# land-ice pressure is too large, the sign of the second term
# makes sense.
gravity = constants['SHR_CONST_G']
deltaLandIcePressure = topDensity * gravity * deltaSSH
deltaLandIcePressure = topDensity * gravity * delta_ssh

landIcePressure = numpy.maximum(
0.0, ds.landIcePressure + deltaLandIcePressure)

ds_out['landIcePressure'] = \
landIcePressure.expand_dims(dim='Time', axis=0)

finalSSH = initSSH
final_ssh = init_ssh

if convert_to_cdf5:
name, ext = os.path.splitext(out_filename)
Expand All @@ -161,14 +166,14 @@ def adjust_ssh(variable, iteration_count, step, update_pio=True,
subprocess.check_call(args)

# Write the largest change in SSH and its lon/lat to a file
with open(f'maxDeltaSSH_{iterIndex:03d}.log', 'w') as \
with open(f'maxDeltaSSH_{iter_index:03d}.log', 'w') as \
log_file:

mask = landIcePressure > 0.
iCell = numpy.abs(deltaSSH.where(mask)).argmax().values
icell = numpy.abs(delta_ssh.where(mask)).argmax().values

ds_cell = ds.isel(nCells=iCell)
deltaSSHMax = deltaSSH.isel(nCells=iCell).values
ds_cell = ds.isel(nCells=icell)
delta_ssh_max = delta_ssh.isel(nCells=icell).values

if on_a_sphere:
coords = (f'lon/lat: '
Expand All @@ -177,21 +182,39 @@ def adjust_ssh(variable, iteration_count, step, update_pio=True,
else:
coords = (f'x/y: {1e-3 * ds_cell.xCell.values:f} '
f'{1e-3 * ds_cell.yCell.values:f}')
string = (f'deltaSSHMax: '
f'{deltaSSHMax:g}, {coords}')
string = (f'delta_ssh_max: '
f'{delta_ssh_max:g}, {coords}')
logger.info(f' {string}')
log_file.write(f'{string}\n')
string = (f'ssh: {finalSSH.isel(nCells=iCell).values:g}, '
string = (f'ssh: {final_ssh.isel(nCells=icell).values:g}, '
f'landIcePressure: '
f'{landIcePressure.isel(nCells=iCell).values:g}')
f'{landIcePressure.isel(nCells=icell).values:g}')
logger.info(f' {string}')
log_file.write(f'{string}\n')

if delta_ssh_threshold is not None:
if abs(deltaSSHMax) < delta_ssh_threshold:
if abs(delta_ssh_max) < delta_ssh_threshold:
break

logger.info(" - Complete\n")

if out_filename is not None:
shutil.copy(out_filename, 'adjusted_init.nc')

with xarray.open_dataset('adjusted_init.nc') as ds:
ds = ds.isel(Time=0)
final_ssh = ds.ssh
final_land_ice_pressure = ds.landIcePressure
delta_ssh = final_ssh - orig_ssh
masked_delta_ssh = modify_mask * delta_ssh
delta_land_ice_pressure = \
final_land_ice_pressure - orig_land_ice_pressure
masked_delta_land_ice_pressure = \
modify_mask * delta_land_ice_pressure
ds_out = xarray.Dataset()
ds_out['delta_ssh'] = delta_ssh
ds_out['masked_delta_ssh'] = masked_delta_ssh
ds_out['delta_land_ice_pressure'] = delta_land_ice_pressure
ds_out['masked_delta_land_ice_pressure'] = \
masked_delta_land_ice_pressure
write_netcdf(ds_out, 'total_delta.nc')