WSL/SLF GitLab Repository

Commit a4b0a56d authored by Nander Wever's avatar Nander Wever
Browse files

Implementing the hydrostatic balance calculation for sea ice simulations with Alpine3D.

parent a49b182f
......@@ -97,7 +97,7 @@ SnowpackInterface::SnowpackInterface(const mio::Config& io_cfg, const size_t& nb
const std::string& grids_requirements,
const bool is_restart_in)
: run_info(), io(io_cfg), pts(prepare_pts(vec_pts)),dem(dem_in),
is_restart(is_restart_in), useCanopy(false), enable_lateral_flow(false), a3d_view(false),
is_restart(is_restart_in), useCanopy(false), enable_lateral_flow(false), enforce_hydrostatic_balance(false), a3d_view(false),
do_io_locally(true), station_name(),glacier_katabatic_flow(false), snow_production(false), snow_grooming(false),
Tsoil_idx(), grids_start(0), grids_days_between(0), ts_start(0.), ts_days_between(0.), prof_start(0.), prof_days_between(0.),
grids_write(true), ts_write(false), prof_write(false), snow_write(false), snow_poi_written(false), glacier_from_grid(false),
......@@ -175,6 +175,9 @@ SnowpackInterface::SnowpackInterface(const mio::Config& io_cfg, const size_t& nb
//check if lateral flow is enabled
sn_cfg.getValue("LATERAL_FLOW", "Alpine3D", enable_lateral_flow, IOUtils::nothrow);
//check if lateral flow is enabled
sn_cfg.getValue("ENFORCE_HYDROSTATIC_BALANCE", "Alpine3D", enforce_hydrostatic_balance, IOUtils::nothrow);
//check if A3D view should be used for grids
sn_cfg.getValue("A3D_VIEW", "Output", a3d_view, IOUtils::nothrow);
......@@ -314,6 +317,7 @@ SnowpackInterface& SnowpackInterface::operator=(const SnowpackInterface& source)
glaciers = source.glaciers;
techSnow = source.techSnow;
enable_lateral_flow = source.enable_lateral_flow;
enforce_hydrostatic_balance = source.enforce_hydrostatic_balance;
a3d_view = source.a3d_view;
sn_cfg = source.sn_cfg;
......@@ -344,10 +348,11 @@ SnowpackInterface& SnowpackInterface::operator=(const SnowpackInterface& source)
std::string SnowpackInterface::getGridsRequirements() const
{
std::string ret = "ELEV";
if (glacier_katabatic_flow) {
return "GLACIER TSS HS";
ret += " GLACIER TSS HS";
}
return "";
return ret;
}
mio::Config SnowpackInterface::readAndTweakConfig(const mio::Config& io_cfg, const bool have_pts)
......@@ -827,6 +832,13 @@ void SnowpackInterface::calcNextStep()
// timing
timer.restart();
const std::string variant = sn_cfg.get("VARIANT", "SnowpackAdvanced", "");
double sealevel = IOUtils::nodata;
double totSeaIceMass = IOUtils::nodata;
if (variant == "SEAICE") {
sealevel = calcHydrostaticBalance(0., totSeaIceMass, true);
}
size_t errCount = 0;
const mio::Grid2DObject tmp_psum(psum, mpi_offset, 0, mpi_nx, dimy);
const mio::Grid2DObject tmp_psum_ph(psum_ph, mpi_offset, 0, mpi_nx, dimy);
......@@ -859,6 +871,10 @@ void SnowpackInterface::calcNextStep()
calcLateralFlow();
}
if (variant == "SEAICE") {
sealevel = calcHydrostaticBalance(sealevel, totSeaIceMass, false);
}
//Retrieve special points data and write files
if (!pts.empty()) write_special_points();
......@@ -1154,6 +1170,7 @@ SN_SNOWSOIL_DATA SnowpackInterface::getIcePixel(const double glacier_height, con
if (MPIControl::instance().master() || do_io_locally) {
const bool useSoil = sn_cfg.get("SNP_SOIL", "Snowpack");
const std::string variant = sn_cfg.get("VARIANT", "SnowpackAdvanced", "");
const std::string coordsys = sn_cfg.get("COORDSYS", "Input");
const std::string coordparam = sn_cfg.get("COORDPARAM", "Input", "");
Coords llcorner_out( dem.llcorner );
......@@ -1246,6 +1263,11 @@ SN_SNOWSOIL_DATA SnowpackInterface::getIcePixel(const double glacier_height, con
if (is_special_point) { //create SMET files for special points
write_SMET_header(snowPixel.meta, landuse(ix, iy));
}
if (variant == "SEAICE") {
snowPixel.Seaice->ForcedSeaLevel = snowPixel.Cdata.height;
snowPixel.Cdata.height = 0.;
}
}
}
if (ii == MPIControl::instance().master_rank() || do_io_locally) {
......@@ -1407,3 +1429,97 @@ void SnowpackInterface::calcLateralFlow()
}
return;
}
double SnowpackInterface::calcHydrostaticBalance(const double& sealevel_in, double& tot_mass_in, const bool& begin)
{
std::vector<SnowStation*> snow_pixel;
// Retrieve snow stations
for (size_t ii = 0; ii < workers.size(); ii++) {
workers[ii]->getLateralFlow(snow_pixel);
}
// Calculate new hydrostatic balance
size_t ix=0; // The source cell x coordinate
size_t errCount = 0;
double tot_mass = 0.;
double tot_OceanSalinity = 0.;
double tot_SeaLevel = 0.;
size_t tot_mass_n = 0;
size_t tot_skipped = 0;
//#pragma omp parallel for schedule(dynamic, 1) reduction(+: errCount)
for (size_t ii = 0; ii < workers.size(); ii++) { // Cycle over all workers
try {
for (size_t jj = 0; jj < worker_deltax[ii]; jj++) { // Cycle over x range per worker
for (size_t iy = 0; iy < dimy; iy++) { // Cycle over y
ix = worker_startx[ii] + jj;
const size_t index_SnowStation_src = ix * dimy + iy; // Index of source cell of water
if (snow_pixel[index_SnowStation_src] != NULL) { // Make sure it is not a NULL pointer (in case of skipped cells)
tot_mass += snow_pixel[index_SnowStation_src]->swe;
tot_OceanSalinity += snow_pixel[index_SnowStation_src]->Seaice->OceanSalinity;
tot_SeaLevel += snow_pixel[index_SnowStation_src]->Seaice->ForcedSeaLevel;
tot_mass_n++;
} else {
tot_skipped++;
}
}
}
} catch(const std::exception& e) {
++errCount;
cout << e.what() << std::endl;
}
}
if (errCount>0) {
//something wrong took place, quitting. At least we tried writing the special points out
std::abort(); //force core dump
}
if (tot_mass_n == 0) {
return IOUtils::nodata;
}
tot_mass /= tot_mass_n;
tot_OceanSalinity /= tot_mass_n;
tot_SeaLevel /= tot_mass_n;
const double tmp_sealevel = (tot_mass) / (Constants::density_water + SeaIce::betaS * tot_OceanSalinity); // Enforcing hydrostatic balance
const double delta_sl1 = (sealevel_in - tot_SeaLevel);
const double delta_sl2 = (tot_mass - tot_mass_in) / (Constants::density_water + SeaIce::betaS * tot_OceanSalinity); // Delta sea level, keeping offset of hydrostatic balance
if (begin) {
tot_mass_in = tot_mass;
return ((enforce_hydrostatic_balance) ? (tmp_sealevel) : (tot_SeaLevel));
}
printf("[i] HydroStatic Balance from %d pixels (%d skipped): sealevel=%f --> sealevel=%f --> sealevel=%f\n", int(tot_mass_n), int(tot_skipped), sealevel_in, tot_SeaLevel, tmp_sealevel);
// Apply hydrostatic adjustment
#pragma omp parallel for schedule(dynamic, 1) reduction(+: errCount)
for (size_t ii = 0; ii < workers.size(); ii++) { // Cycle over all workers
try {
for (size_t jj = 0; jj < worker_deltax[ii]; jj++) { // Cycle over x range per worker
for (size_t iy = 0; iy < dimy; iy++) { // Cycle over y
ix = worker_startx[ii] + jj;
const size_t index_SnowStation_src = ix * dimy + iy; // Index of source cell of water
if (snow_pixel[index_SnowStation_src] != NULL) { // Make sure it is not a NULL pointer (in case of skipped cells)
if (snow_pixel[index_SnowStation_src]->Seaice->ForcedSeaLevel != IOUtils::nodata) {
if (enforce_hydrostatic_balance) {
snow_pixel[index_SnowStation_src]->Seaice->ForcedSeaLevel += (tmp_sealevel - tot_SeaLevel);
} else {
snow_pixel[index_SnowStation_src]->Seaice->ForcedSeaLevel += delta_sl1 + delta_sl2;
}
}
}
}
}
} catch(const std::exception& e) {
++errCount;
cout << e.what() << std::endl;
}
}
if (errCount>0) {
//something wrong took place, quitting. At least we tried writing the special points out
std::abort(); //force core dump
}
return tmp_sealevel;
}
......@@ -187,6 +187,7 @@ class Runoff; // forward declaration, cyclic header include
const std::vector<SurfaceFluxes*>& surface_flux);
void write_special_points();
void calcLateralFlow();
double calcHydrostaticBalance(const double& sealevel_in, double& tot_mass_in, const bool& begin);
RunInfo run_info;
......@@ -198,7 +199,7 @@ class Runoff; // forward declaration, cyclic header include
const mio::DEMObject dem;
// Config dependent information
bool is_restart, useCanopy, enable_lateral_flow, a3d_view;
bool is_restart, useCanopy, enable_lateral_flow, enforce_hydrostatic_balance, a3d_view;
bool do_io_locally; // if false all I/O will only be done on the master process
std::string station_name; // value for the key OUTPUT::EXPERIMENT
bool glacier_katabatic_flow, snow_production, snow_grooming;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment