Skip to content

Issue649 addd diagnostic fsd #661

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

Open
docguibou wants to merge 8 commits into
base: develop
Choose a base branch
from
Open

Issue649 addd diagnostic fsd #661

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
b940e76
Adding a diagnostic FSD, compiles with OASIS
Nov 28, 2024
37c9456
Fixing bugs in diagnostic FSD. Compiles and run.
Nov 28, 2024
7937ea4
removing files added by mistake
Nov 28, 2024
9015d5f
Added new floe size outputs + rewrote computation from prognostic FSD
Nov 28, 2024
3b31cad
various fixes recommended by Tim W.
Jan 6, 2025
7e76947
Update model/finiteelement.cpp
docguibou Jan 6, 2025
6716f34
Update model/finiteelement.cpp
docguibou Jan 6, 2025
75472e3
Update model/finiteelement.cpp
docguibou Jan 6, 2025
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
214 changes: 164 additions & 50 deletions model/finiteelement.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4266,7 +4266,122 @@ FiniteElement::updateSigmaDamage(double const dt)
}//loop over elements
}//updateSigmaDamage

//------------------------------------------------------------------------------------------------------
//! Helper function to compute the n-th moment of the FSD.
//! Called by `computeParametersDiagnosticFSD`
double FiniteElement::computeCharacteristicDiameter(double dmin, double dmax, const double gamma, int n)
{
// Helper function to compute the k-th moment of the diameter distribution
auto moment_diameter = [=](int k) -> double {
double exponent = 1 - gamma + k; // Compute the power-law exponent
if (std::abs(exponent) < 1e-6)
{
// Special case: floe_exponent = moment_exponent - 1 (logarithmic divergence)
return std::log(dmax / dmin);
}
return (std::pow(dmax, exponent) - std::pow(dmin, exponent)) / exponent;
};

// Compute the numerator (n-th moment) and denominator (m-th moment)
const int m = n-1 ;
double numerator = moment_diameter(n);
double denominator = moment_diameter(m);

// Return the ratio
return numerator / denominator;
}

//------------------------------------------------------------------------------------------------------
//! Computes floe size paremetres for a FSD largely inspired from Denton et Timmermans 2022 and brainstormed with C. Rousset
//! The idea is to assume a single exponent power law for the FSD and have a fixed upper limit.
//! The exponent varies with concentration and as suggested by D&T.
//! Called by `updateIceDiagnostics`
void
FiniteElement::computeParametersDiagnosticFSD(const int i, const double upper_lim_floe_size, const double poisson)
{
// Largely inspired from Denton et Timmermans 2022 and brainstormed with C. Rousset
// The idea is to assume a single exponent power law for the FSD and have a fixed upper limit.
// The exponent varies with concentration and as suggested by D&T
// The FSD characteristic floe sizes are sensitive to choice of lower and upper limits
// Denton et Timmermans 2022 use a number density FSD as a function of floe area with an exponent alpha (f(a)=Ca**-alpha)).
// One can how that replacing area by floe size gives an exponent gamma=2*alpha-1
// Careful that D&T22 work with m=-alpha, I prefer thinking with a positive exponent for the power-law
const double gamma = 2.*(2.-0.26*M_conc[i])-1 ; // 2*alpha(conc)-1
double sea_ice_thickness = M_thick[i] ;
if (M_ice_cat_type == setup::IceCategoryType::YOUNG_ICE)
sea_ice_thickness += M_h_young[i];///ctot ;
const double dmin = 0.5 * std::pow ( std::pow(PI,4)*5.5e9*std::pow(sea_ice_thickness,3) /
(48*physical::rhow*physical::g * (1-std::pow(poisson,2)) )
, 0.25 ) ;
// TODO : Introduce proper variables for characteristic floe sizes
// Compute moments
D_dchar[i] = this->computeCharacteristicDiameter(dmin, upper_lim_floe_size, gamma, 3.0);
D_dlength[i] = this->computeCharacteristicDiameter(dmin, upper_lim_floe_size, gamma, 2.0);
D_dnum[i] = this->computeCharacteristicDiameter(dmin, upper_lim_floe_size, gamma, 1.0);
D_dmean[i] = D_dnum[i] ;
D_dmax[i] = D_dchar[i] ;
}
#ifdef OASIS
//! calculate parameters summarising the prognostic FSD.
//! Called by `updateIceDiagnostics`
void
Comment thread
docguibou marked this conversation as resolved.
FiniteElement::computeParametersPrognosticFSD(const int i)
{
double M_1_mech= 0. ;
double M_0 = 0, M_1 = 0., M_neg1 = 0., M_neg2 = 0.;
if (M_conc_mech_fsd[M_num_fsd_bins-1][i]<=M_dmax_c_threshold*D_conc[i])
{
//If more than 90% of sea ice that is in the element is "broken"
Copy link
Copy Markdown
Contributor

@tdcwilliams tdcwilliams Jan 3, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can you add an explanation for the ordering? It looks like you start at the largest floe size bin and decrease - what is the reason for that? (I would expect the other way so if 90% of the ice is in the smallest bin you'd stop there with a small Dmax, otherwise keep increasing

Copy link
Copy Markdown
Contributor Author

@docguibou docguibou Jan 6, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So happy you're back! Hmmm I haven't changed the ordering actually. I suppose my reasoning at the time was that there is more elements with unbroken ice than with broken ice, so the loop would stop faster starting from the larger floes. If it really does not make sense to you, I can rewrite it.

Copy link
Copy Markdown
Contributor

@tdcwilliams tdcwilliams Jan 6, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unless I am misunderstanding, the current Dmax comes from \int_{Dmax}^\infty = 0.9, but reversing gives \int_0^{Dmax} = 0.9? These are quite different. But perhaps we could think about it and if it is important and change it in the other branch

Copy link
Copy Markdown
Contributor Author

@docguibou docguibou Jan 6, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I defined M_dmax_c_threshold as (1-0.9) to account for that, that's ok no? Super confusing with these comments though.

double ctot_fsd = M_conc_mech_fsd[M_num_fsd_bins-1][i];
for(int j=M_num_fsd_bins-2;j>-1;j--)
{
ctot_fsd += M_conc_mech_fsd[j][i] ;
if (ctot_fsd>M_dmax_c_threshold*D_conc[i])
{
D_dmax[i]=M_fsd_bin_up_limits[j] ;
break ;
}
}
for (int j = 0; j < M_num_fsd_bins; j++)
M_1_mech += M_conc_mech_fsd[j][i] * M_fsd_bin_centres[j];
}
else
{
D_dmax[i]=(M_conc_mech_fsd[M_num_fsd_bins-1][i]/D_conc[i] - M_dmax_c_threshold) / (1.-M_dmax_c_threshold)
* (M_fsd_unbroken_floe_size- M_fsd_bin_up_limits[M_num_fsd_bins-1]) + M_fsd_bin_up_limits[M_num_fsd_bins-1] ;
D_dmax[i]=std::min(D_dmax[i],M_fsd_unbroken_floe_size);

for (int j = 0; j < M_num_fsd_bins-1; j++)
M_1_mech += M_conc_mech_fsd[j][i] * M_fsd_bin_centres[j];
M_1_mech += M_conc_mech_fsd[M_num_fsd_bins-1][i] * D_dmax[i];
}
D_dmax[i] =std::min(D_dmax[i],M_fsd_unbroken_floe_size);
D_dmax[i] =std::max(D_dmax[i],M_fsd_min_floe_size);
// Loop over bins to compute moments for the thermo FSD
for (int j = 0; j < M_num_fsd_bins-1; j++)
{
double d_j = M_fsd_bin_centres[j];
M_1 += M_conc_fsd[j][i] * d_j;
M_neg1 += M_conc_fsd[j][i] / d_j; // M_{-1}
M_neg2 += M_conc_fsd[j][i] / (d_j * d_j); // M_{-2}
}
double d_n = (D_dmax[i]-M_fsd_bin_low_limits[M_num_fsd_bins-1])/2. + M_fsd_bin_low_limits[M_num_fsd_bins-1];
Comment thread
tdcwilliams marked this conversation as resolved.
M_1 += M_conc_fsd[M_num_fsd_bins-1][i] * d_n;
M_neg1 += M_conc_fsd[M_num_fsd_bins-1][i] / d_n; // M_{-1}
M_neg2 += M_conc_fsd[M_num_fsd_bins-1][i] / (d_n * d_n); // M_{-2}
M_0 = D_conc[i];
// No need to normalize the other moments as we do a ratio
// Now compute the characteristic floe sizes
D_dmean[i] = M_1_mech / M_0; // Characteristic floe size associated with the number of floes
Copy link
Copy Markdown
Contributor

@tdcwilliams tdcwilliams Jan 6, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

looks like D_dchar but for mech FSD?

Copy link
Copy Markdown
Contributor Author

@docguibou docguibou Jan 6, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is! I am reluctant to change that as it is what we published, but to be consistent with assumptions from scattering parameterizations, I should have not used the areal FSD. However, that would be assuming that absolutely all floes contribute to scattering, which is likely to be wrong. I can discuss that with you at some point.

D_dchar[i] = M_1 / M_0; // Mean floe size associated with area (D_dmean), same as representative floe size defined by Roach et al. 2018
D_dlength[i] = M_0 / M_neg1; // Characteristic floe size associated with floe edge length (perimetre)
D_dnum[i] = M_neg1 / M_neg2; // Characteristic floe size associated with the number of floes

D_dmean[i] = std::max(D_dmean[i],M_fsd_min_floe_size);
D_dlength[i] = std::max(D_dlength[i],M_fsd_min_floe_size);
D_dchar[i] = std::max(D_dchar[i],M_fsd_min_floe_size);
D_dnum[i] = std::max(D_dnum[i],M_fsd_min_floe_size);
}
//------------------------------------------------------------------------------------------------------
//! Redistribute the floes in the FSD after break_up if prob[i]>0
//! v 0.1 of the function, should be done properly with adjustable parameters in a namelist
Expand Down Expand Up @@ -4380,8 +4495,6 @@ FiniteElement::redistributeFSD()
M_conc_fsd[j][i] -= broken_area ;

//!Define a redistributor beta
//! with uniform redistribution of floes for sizes below the broken categories
//! (for any D such as Dmin < D < D_broken_cat, it generates an equal number of floes)
for (int k=0; k<=j ; k++)
{
double beta = (std::pow(M_fsd_bin_up_limits[k],3)- std::pow(M_fsd_bin_low_limits[k],3))
Expand Down Expand Up @@ -7328,6 +7441,12 @@ FiniteElement::initModelVariables()
M_variables_elt.push_back(&D_dmax);
D_dmean = ModelVariable(ModelVariable::variableID::D_dmean);
M_variables_elt.push_back(&D_dmean);
D_dchar = ModelVariable(ModelVariable::variableID::D_dchar);
M_variables_elt.push_back(&D_dchar);
D_dlength = ModelVariable(ModelVariable::variableID::D_dlength);
M_variables_elt.push_back(&D_dlength);
D_dnum = ModelVariable(ModelVariable::variableID::D_dnum);
M_variables_elt.push_back(&D_dnum);

//! - 2) loop over M_variables_elt in order to sort them
//! for restart/regrid/export
Expand Down Expand Up @@ -7865,10 +7984,23 @@ FiniteElement::initOASIS()
void
FiniteElement::updateIceDiagnostics()
{
// diagnostic FSD parametres
const double upper_lim_floe_size = 3000 ;
const double poisson = 0.3 ;
bool use_diagnostic_FSD = true ;
#ifdef OASIS
if (M_num_fsd_bins>0)
use_diagnostic_FSD = false ;
#endif
D_conc.resize(M_num_elements); //! \param D_conc (double) Total concentration (diagnostic)
D_thick.resize(M_num_elements); //! \param D_thick (double) Total thickness (diagnostic)
D_snow_thick.resize(M_num_elements); //! \param D_snow_thick (double) Total snow thickness (diagnostic)
D_tsurf.resize(M_num_elements); //! \param D_tsurf (double) Mean surface temperature (diagnostic)
D_dmax.resize(M_num_elements) ;
D_dmean.resize(M_num_elements) ;
D_dlength.resize(M_num_elements);
D_dchar.resize(M_num_elements) ;
D_dnum.resize(M_num_elements) ;
for(int k=0; k<2; k++)
D_sigma[k].resize(M_num_elements);

Expand Down Expand Up @@ -7904,57 +8036,24 @@ FiniteElement::updateIceDiagnostics()
double const dyN = shape_coeff[j+3];
D_divergence[i] += dxN*u + dyN*v;
}

// FSD relative parameters
D_dmean[i] = 0. ;
D_dmax[i] = 0. ;
#if defined (OASIS)
if (M_num_fsd_bins>0)
// TODO Here compute characteristic floe sizes from either diagnostic or prognostic FSD
if (D_conc[i]<1e-12)
{
if (D_conc[i]>0)
{
D_dmax[i]=1000. ;
D_dmean[i]=1000.;
if (M_conc_mech_fsd[M_num_fsd_bins-1][i]<=M_dmax_c_threshold*D_conc[i])
{
//If more than 90% of sea ice that is in the element is "broken"
double ctot_fsd = M_conc_mech_fsd[M_num_fsd_bins-1][i];
for(int j=M_num_fsd_bins-2;j>-1;j--)
{
ctot_fsd += M_conc_mech_fsd[j][i] ;
if (ctot_fsd>M_dmax_c_threshold*D_conc[i])
{
D_dmax[i]=M_fsd_bin_up_limits[j] ;
break ;
}
}
D_dmean[i]=0.;
for(int j=0;j<M_num_fsd_bins;j++)
D_dmean[i] += M_conc_mech_fsd[j][i] * M_fsd_bin_centres[j]/D_conc[i];
}
else
{
D_dmax[i]=(M_conc_mech_fsd[M_num_fsd_bins-1][i]/D_conc[i] - M_dmax_c_threshold) / (1.-M_dmax_c_threshold)
* (M_fsd_unbroken_floe_size- M_fsd_bin_up_limits[M_num_fsd_bins-1]) + M_fsd_bin_up_limits[M_num_fsd_bins-1] ;
D_dmax[i]=std::min(D_dmax[i],M_fsd_unbroken_floe_size);

D_dmean[i]=0.;
for(int j=0;j<M_num_fsd_bins-1;j++)
D_dmean[i] += M_conc_mech_fsd[j][i] * M_fsd_bin_centres[j]/D_conc[i] ;
D_dmean[i]+= M_conc_mech_fsd[M_num_fsd_bins-1][i] * D_dmax[i] /D_conc[i] ;
}
D_dmax[i] =std::min(D_dmax[i],M_fsd_unbroken_floe_size);
D_dmax[i] =std::max(D_dmax[i],M_fsd_min_floe_size);
D_dmean[i]=std::min(D_dmax[i],D_dmean[i]);
D_dmean[i]=std::max(D_dmean[i],M_fsd_min_floe_size);
}
else
{
D_dmax[i] = 0. ;
D_dmean[i] = 0. ;
}
D_dmax[i] = 0. ;
D_dmean[i] = 0. ;
D_dlength[i] = 0.;
D_dchar[i] = 0.;
D_dnum[i] = 0.;
}
else
{
#ifdef OASIS
if (M_num_fsd_bins>0)
this->computeParametersPrognosticFSD(i);
#endif
if (use_diagnostic_FSD)
this->computeParametersDiagnosticFSD(i, upper_lim_floe_size, poisson) ;
}
}
}

Expand Down Expand Up @@ -8889,6 +8988,18 @@ FiniteElement::updateMeans(GridOutput& means, double time_factor)
for (int i=0; i<M_local_nelements; i++)
it->data_mesh[i] += D_dmean[i]*time_factor;
break;
case (GridOutput::variableID::dchar):
for (int i=0; i<M_local_nelements; i++)
it->data_mesh[i] += D_dchar[i]*time_factor;
break;
case (GridOutput::variableID::dlength):
for (int i=0; i<M_local_nelements; i++)
it->data_mesh[i] += D_dlength[i]*time_factor;
break;
case (GridOutput::variableID::dnum):
for (int i=0; i<M_local_nelements; i++)
it->data_mesh[i] += D_dnum[i]*time_factor;
break;

// Coupling variables (not covered elsewhere)
// NB: reversed sign convention!
Expand Down Expand Up @@ -9108,6 +9219,9 @@ FiniteElement::initMoorings()
("saltflux", GridOutput::variableID::saltflux)
("dmax",GridOutput::variableID::dmax)
("dmean",GridOutput::variableID::dmean)
("dchar",GridOutput::variableID::dchar)
("dlength",GridOutput::variableID::dlength)
("dnum",GridOutput::variableID::dnum)
// Forcing
("tair", GridOutput::variableID::tair)
("sphuma", GridOutput::variableID::sphuma)
Expand Down
12 changes: 9 additions & 3 deletions model/finiteelement.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -301,13 +301,16 @@ class FiniteElement
void initUpdateGhosts();
int globalNumToprocId(int global_num);

// FSD related functions
void computeParametersDiagnosticFSD(const int i, const double upper_lim_floe_size, const double poisson);
double computeCharacteristicDiameter(const double dmin, const double dmax, const double gamma, const int n) ;
#ifdef OASIS
bool M_couple_waves;
bool M_recv_wave_stress;
// FSD related functions
void initFsd();
void redistributeFSD();
void updateFSD();
void computeParametersPrognosticFSD(const int i);
std::vector<double> computeWaveBreakingProb();
double computeLateralAreaFSD(const int cpt);
double computeLeadFractionFSD(const int cpt);
Expand Down Expand Up @@ -823,8 +826,11 @@ class FiniteElement
ModelVariable D_fwflux; // Fresh-water flux at ocean surface [kg/m2/s]
ModelVariable D_fwflux_ice; // Fresh-water flux at ocean surface due to ice processes [kg/m2/s]
ModelVariable D_brine; // Brine release into the ocean [kg/m2/s]
ModelVariable D_dmax; //max floe size [m]
ModelVariable D_dmean; //mean floe size [m]
ModelVariable D_dmax; //max floe size from mechanical FSD (for exchange with wave model, Boutin et al. 2021) [m]
ModelVariable D_dmean; //mean floe size from mechanical FSD (for exchange with wave model, Boutin et al. 2021) [m]
ModelVariable D_dchar; //characterisic floe size (Roach et al., 2018) associated with floe area from (observed) FSD [m]
ModelVariable D_dlength; // floe size associated with floe edge length from (observed) FSD [m]
ModelVariable D_dnum; // floe size considering floe number from (observed) FSD [m]
ModelVariable D_tau_ow; // Ocean atmosphere drag coefficient - still needs to be multiplied with the wind [Pa/s/m] (for the coupled ice-ocean system)
ModelVariable D_evap; // Evaporation out of the ocean [kg/m2/s]
ModelVariable D_rain; // Rain into the ocean [kg/m2/s]
Expand Down
28 changes: 28 additions & 0 deletions model/gridoutput.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -218,6 +218,9 @@ class GridOutput
// WIM variables
dmax = 300,
dmean = 301,
dchar = 302,
dlength = 303,
dnum = 304,

// Coupling variables not already covered elsewhere
taux = 901,
Expand Down Expand Up @@ -715,6 +718,7 @@ class GridOutput
longName = "Ice-atmosphere thermodynamic drag";
stdName = "ice-atmosphere_thermodynamic_drag";
Units = "1/s";
break;
case (variableID::meltpond_fraction):
name = "meltpond_fraction";
longName = "Meltpond fraction";
Expand Down Expand Up @@ -812,6 +816,30 @@ class GridOutput
cell_methods = "area: mean where sea_ice";
break;

case (variableID::dchar):
name = "dchar";
longName = "Mean floe size associated with floe area";
stdName = "char_floe_size";
Units = "m";
cell_methods = "area: mean where sea_ice";
break;

case (variableID::dnum):
name = "dnum";
longName = "Mean floe size associated with number of floes";
stdName = "mean_num_floe_size";
Units = "m";
cell_methods = "area: mean where sea_ice";
break;

case (variableID::dlength):
name = "dlength";
longName = "Mean floe size associated with floe perimeter";
stdName = "mean_length_floe_size";
Units = "m";
cell_methods = "area: mean where sea_ice";
break;

//forcing variables
case (variableID::tair):
name = "t2m";
Expand Down
Loading