Skip to content

Commit

Permalink
PICO: Added some in-line description of the functions and a test case
Browse files Browse the repository at this point in the history
  • Loading branch information
@talbrecht
talbrecht committed Apr 1, 2021
1 parent dc5f8f9 commit 9c32530
Show file tree
Hide file tree
Showing 2 changed files with 68 additions and 23 deletions.
35 changes: 12 additions & 23 deletions src/coupler/ocean/PicoGeometry.cc
View file
Expand Up @@ -132,18 +132,15 @@ void PicoGeometry::update(const IceModelVec2S &bed_elevation, const IceModelVec2
compute_distances_cf(m_ocean_mask, m_ice_rises, exclude_ice_rises, m_distance_cf);
}

// these two could be done at the same time
// computing ice_shelf_mask and box_mask could be done at the same time
{
compute_ice_shelf_mask(m_ice_rises, m_lake_mask, m_ice_shelves);
m_n_shelves = m_ice_shelves.max() + 1;

get_basin_neighbors(cell_type, m_basin_mask, m_n_basin_neighbors);

std::vector<int> cfs_in_basins_per_shelf(m_n_shelves*m_n_basins,0);
std::vector<int> most_shelf_cells_in_basin(m_n_shelves, 0);
identify_calving_front_connection(cell_type, m_basin_mask, m_ice_shelves, most_shelf_cells_in_basin, cfs_in_basins_per_shelf);

//split_ice_shelves(cell_type, m_basin_mask, m_n_basin_neighbors, m_ice_shelves);
split_ice_shelves(cell_type, m_basin_mask, m_n_basin_neighbors, most_shelf_cells_in_basin, cfs_in_basins_per_shelf, m_ice_shelves);

compute_continental_shelf_mask(bed_elevation, m_ice_rises, continental_shelf_depth, m_continental_shelf);
Expand Down Expand Up @@ -505,7 +502,8 @@ void PicoGeometry::compute_ocean_mask(const IceModelVec2CellType &cell_type, Ice
}

/*!
* FIXME
* Find the two neighboring basins by checking for the basin boundaries on the ice free ocean.
* Could there be more than two neighbors? Should we better identify the intersection at the coastline?
*/
void PicoGeometry::get_basin_neighbors(const IceModelVec2CellType &cell_type,
const IceModelVec2Int &basin_mask,
Expand All @@ -524,7 +522,6 @@ void PicoGeometry::get_basin_neighbors(const IceModelVec2CellType &cell_type,
auto M = cell_type.int_star(i, j);
int bn = 0; //neighbor basin id

// finding the basin boundary on the ice free ocean
if (cell_type.as_int(i, j) == MASK_ICE_FREE_OCEAN and
((M.n == MASK_ICE_FREE_OCEAN and B.n != b) or
(M.s == MASK_ICE_FREE_OCEAN and B.s != b) or
Expand All @@ -547,16 +544,17 @@ void PicoGeometry::get_basin_neighbors(const IceModelVec2CellType &cell_type,
}
}
}
// compute global sums
//GlobalSum(m_grid->com, result.data(), 2*m_n_basins);
for (int b = 1; b < 2*m_n_basins; ++b) {
result[b] = GlobalSum(m_grid->com, result[b]);
if (b%2==0)
m_log->message(2, "PICO, get basin neighbors with b=%d, b1=%d and b2=%d \n",b,result[b],result[b+1]);
m_log->message(2, "PICO, get basin neighbors of b=%d: b1=%d and b2=%d \n",b/2,result[b],result[b+1]);
}
}

/*!
* FIXME
* Find all basins b, in which the ice shelf s has a calving front with potential ocean water intrusion.
* Find the basin bmax, in which the ice shelf s has the most cells
*/
void PicoGeometry::identify_calving_front_connection(const IceModelVec2CellType &cell_type,
const IceModelVec2Int &basin_mask,
Expand All @@ -575,7 +573,6 @@ void PicoGeometry::identify_calving_front_connection(const IceModelVec2CellType
int b = basin_mask.as_int(i, j);
n_shelf_cells_per_basin[s*m_n_basins+b]++;

// find all basins b, in which the ice shelf s has a calving front with potential ocean water intrusion
if (cell_type.as_int(i, j) == MASK_FLOATING) {
auto M = cell_type.int_star(i, j);
if (M.n == MASK_ICE_FREE_OCEAN or M.e == MASK_ICE_FREE_OCEAN or M.s == MASK_ICE_FREE_OCEAN or M.w == MASK_ICE_FREE_OCEAN) {
Expand All @@ -587,6 +584,7 @@ void PicoGeometry::identify_calving_front_connection(const IceModelVec2CellType
}

//GlobalSum(m_grid->com, cfs_in_basins_per_shelf.data(), m_n_shelves*m_n_basins);
//GlobalSum(m_grid->com, n_shelf_cells_per_basin.data(), m_n_shelves*m_n_basins);
for (int s = 0; s < m_n_shelves; s++) {
int n_shelf_cells_per_basin_max = 0;
for (int b = 0; b < m_n_basins; b++) {
Expand All @@ -603,7 +601,7 @@ void PicoGeometry::identify_calving_front_connection(const IceModelVec2CellType


/*!
* FIXME
* Find all ice shelves s that spread across non-neighboring basins with calving fronts in those basins and add an ice shelf mask number.
*/
void PicoGeometry::split_ice_shelves(const IceModelVec2CellType &cell_type,
const IceModelVec2Int &basin_mask,
Expand All @@ -616,14 +614,6 @@ void PicoGeometry::split_ice_shelves(const IceModelVec2CellType &cell_type,

IceModelVec::AccessList list{ &cell_type, &basin_mask, &shelf_mask, &m_tmp};

// test if vector exists
//for (int b = 1; b < m_n_basins; ++b) {
// m_log->message(2, "PICO, use basin neighbors with b=%d, b1=%d and b2=%d \n",b,n_basin_neighbors[2*b],n_basin_neighbors[2*b+1]);
//}

//m_n_shelves = shelf_mask.max() + 1;

// now find all ice shelves spread across non-neighboring basins with calving fronts in those basins and add an ice shelf mask number

std::vector<int> n_shelf_cells_to_split(m_n_shelves*m_n_basins,0);
for (Points p(*m_grid); p; p.next()) {
Expand All @@ -634,11 +624,11 @@ void PicoGeometry::split_ice_shelves(const IceModelVec2CellType &cell_type,
int b0 = most_shelf_cells_in_basin[s];
if (b != b0 and b != n_basin_neighbors[2*b0] and b != n_basin_neighbors[2*b0+1] and cfs_in_basins_per_shelf[s*m_n_basins+b] > 0) {
n_shelf_cells_to_split[s*m_n_basins+b]++;
//m_log->message(2, "PICO, test split s=%d with b0=%d and b=%d at %d,%d \n",s,b0,b,i,j);
}
}
}

//GlobalSum(m_grid->com, n_shelf_cells_to_split.data(), m_n_shelves*m_n_basins);
std::vector<int> add_shelf_instance(m_n_shelves*m_n_basins,0);
int m_shelf_numbers_to_add = 0;
for (int s = 0; s < m_n_shelves; s++) {
Expand All @@ -648,7 +638,8 @@ void PicoGeometry::split_ice_shelves(const IceModelVec2CellType &cell_type,
if (n_shelf_cells_to_split[s*m_n_basins+b] > 0) {
m_shelf_numbers_to_add += 1;
add_shelf_instance[s*m_n_basins+b] = m_n_shelves + m_shelf_numbers_to_add;
m_log->message(2, "\nPICO, split ice shelf s=%d with b0=%d and b=%d and n=%d and si=%d\n",s,b0,b,n_shelf_cells_to_split[s*m_n_basins+b],add_shelf_instance[s*m_n_basins+b]);
m_log->message(3, "\nPICO, split ice shelf s=%d with bmax=%d and b=%d and n=%d and si=%d\n",
s,b0,b,n_shelf_cells_to_split[s*m_n_basins+b],add_shelf_instance[s*m_n_basins+b]);
}
}
}
Expand All @@ -659,9 +650,7 @@ void PicoGeometry::split_ice_shelves(const IceModelVec2CellType &cell_type,
int b = basin_mask.as_int(i, j);
int s = shelf_mask.as_int(i, j);
if (add_shelf_instance[s*m_n_basins+b] > 0) {
//shelf_mask.as_int(i, j) = add_shelf_instance[s*m_n_basins+b];
m_tmp(i, j) = add_shelf_instance[s*m_n_basins+b];
//m_log->message(2, "PICO, split ice shelf s=%d and b=%d and si=%d at %d,%d\n",s,b,add_shelf_instance[s*m_n_basins+b],i,j);
}
}
}
Expand Down
56 changes: 56 additions & 0 deletions test/test_pico_split/run_test.sh
View file
@@ -0,0 +1,56 @@
#!/bin/bash

# Copyright (C) 2021 PISM authors
# created by torsten.albrecht@pik-potsdam.de

###############################################################################
# Test simulation of Antarctic ice sheet model to check on PICO adjustments
# concerning large ice shelves crossing basin boundaries.
###############################################################################


INNAME=bedmap2_schmidtko14_50km.nc

# adding SMB
cp $INNAME input_test01.nc
ncap2 -O -s "climatic_mass_balance=0*salinity_ocean" input_test01.nc input_test01.nc
ncatted -O -a units,climatic_mass_balance,o,c,"kg m-2 year-1" \
-a long_name,climatic_mass_balance,o,c,"surface mass balance (accumulation/ablation) rate" \
-a standard_name,climatic_mass_balance,o,c,"land_ice_surface_specific_mass_balance_flux" \
input_test01.nc


#SIXTEENKMGRID="-Mx 381 -My 381 -Lz 6000 -Lbz 2000 -Mz 81 -Mbz 21"
FIFTYKMGRID="-Mx 120 -My 120 -Lz 6000 -Lbz 2000 -Mz 81 -Mbz 21"

stressbalance="-pik -stress_balance ssa+sia -ssa_method fd"
pico="-ocean pico -ocean_pico_file $INNAME -gamma_T 1.0e-5 -overturning_coeff 0.8e6 -exclude_icerises -continental_shelf_depth -2500"
surface="-surface pik"
regrid_opts="-bootstrap $FIFTYKMGRID" # -regrid_file $INNAME"


# merge Amundsen Sea and Ross ice shelves
ncap2 -O -s "where((basins==12 || basins==14) && mask==2 && topg<-600 && lon<-90) thk=topg/(-910.0/1028.0)-10.0" input_test01.nc input_test02.nc

# advance Ross Ice Shelf into Amundsen Sea basin, but without connection
ncap2 -O -s "where((basins==12 || basins==14) && mask==2 && topg<-600 && lon<-90 && lat<-77.4) thk=topg/(-910.0/1028.0)-10.0" input_test01.nc input_test03.nc

# merge Ross Ie Shelf with neighboring basin to check on weighted average
ncap2 -O -s "where((basins==12 || basins==13) && mask==2 && lon<-90) topg=-600.0" input_test01.nc input_test04.nc
ncap2 -O -s "where((basins==12 || basins==13) && mask==2 && topg<=-600 && lon<-90) thk=topg/(-910.0/1028.0)-10.0" input_test04.nc input_test04.nc


for testcase in "test01" "test02" "test03" "test04"; do
echo $testcase
../../bin/pismr -i input_${testcase}.nc $regrid_opts $stressbalance $surface $pico -y 0.001 -o ${testcase}.nc

# assert that pico inout temperatures in basins 12, 13 and 14 are 271.553311083714, 272.781500154734 and 273.596660429239 K
ncks -H -d x,59 -d y,37 -v pico_temperature_box0 ${testcase}.nc | grep -C 1 'pico_temperature_box0 =' | grep -v '^$' | awk 'END{print}'
ncks -H -d x,41 -d y,34 -v pico_temperature_box0 ${testcase}.nc | grep -C 1 'pico_temperature_box0 =' | grep -v '^$' | awk 'END{print}'
ncks -H -d x,28 -d y,49 -v pico_temperature_box0 ${testcase}.nc | grep -C 1 'pico_temperature_box0 =' | grep -v '^$' | awk 'END{print}'
# in test04 the values in basin 12 and 13 should be averaged to 271.673605540195 K

done



0 comments on commit 9c32530

Please sign in to comment.