Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25149) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25150) @@ -5662,4 +5662,8 @@ this->SealevelriseEustaticIce(Sgi,peustatic,masks, oceanarea); + /*Compute hydrological contribution to sea level rise: */ + this->SealevelriseEustaticHydro(Sgi,peustatic,masks, oceanarea); + + } /*}}}*/ @@ -5777,4 +5781,70 @@ /*convolve:*/ for(int i=0;i(eustatic)); + *peustatic=eustatic; + return; +} +/*}}}*/ +void Tria::SealevelriseEustaticHydro(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble oceanarea){ /*{{{*/ + + /*diverse:*/ + int gsize; + IssmDouble area; + IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes eustatic + IssmDouble W; //change in water height thickness (Farrel and Clarke, Equ. 4) + bool notfullygrounded=false; + bool scaleoceanarea= false; + + /*elastic green function:*/ + IssmDouble* G=NULL; + + /*ice properties: */ + IssmDouble rho_water; + IssmDouble rho_freshwater; + + /*constants:*/ + IssmDouble constant=0; + + /*Initialize eustatic component: do not skip this step :):*/ + IssmDouble eustatic = 0.; + + /*early return if we are on an ice cap:*/ + if(masks->isiceonly[this->lid]){ *peustatic=0; return; } + + /*early return if we are fully floating:*/ + if (masks->isfullyfloating[this->lid]){ constant=0; *peustatic=0; return; } + + /*If we are here, we are on earth, not on ice: */ + + /*recover material parameters: */ + rho_water=FindParam(MaterialsRhoSeawaterEnum); + rho_freshwater=FindParam(MaterialsRhoFreshwaterEnum); + + /*recover ocean area scaling: */ + this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); + + /*retrieve precomputed G:*/ + this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); + + /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ + this->GetInput2Value(&area,AreaEnum); + + /*Retrieve water height at vertices: */ + Input2* deltathickness_input=this->GetInput2(SurfaceloadWaterHeightChangeEnum); + if (!deltathickness_input)_error_("SurfaceloadWaterHeightChangeEnum input needed to compute sea level rise!"); + deltathickness_input->GetInputAverage(&W); + + /*Compute eustatic component:*/ + _assert_(oceanarea>0.); + if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 + eustatic += rho_freshwater*area*phi*W/(oceanarea*rho_water); + + /*convert from m to kg/m^2:*/ + W=W*rho_freshwater*phi; + + /*convolve:*/ + for(int i=0;i *steric_rate_g = NULL; Vector *dynamic_rate_g = NULL; - Vector *hydro_rate_g = NULL; Vector *U_esa_rate= NULL; Vector *N_esa_rate= NULL; @@ -296,5 +295,4 @@ GetStericRate(&steric_rate_g,femmodel); GetDynamicRate(&dynamic_rate_g,femmodel); - GetVectorFromInputsx(&hydro_rate_g,femmodel,SurfaceloadWaterHeightChangeEnum,VertexSIdEnum); if(geodetic){ GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum); @@ -304,5 +302,5 @@ } - /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate) * dt + steric_rate + dynamic_rate + hydro_rate* dt*/ + /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate) * dt + steric_rate + dynamic_rate dt*/ if(geodetic){ SL->AXPY(N_gia_rate,dt); @@ -311,5 +309,4 @@ SL->AXPY(steric_rate_g,dt); SL->AXPY(dynamic_rate_g,dt); - SL->AXPY(hydro_rate_g,dt); /*compute new bedrock position: */ @@ -328,5 +325,4 @@ delete steric_rate_g; delete dynamic_rate_g; - delete hydro_rate_g; if(geodetic){ delete U_esa_rate; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 25149) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 25150) @@ -365,4 +365,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_water",MaterialsRhoSeawaterEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.materials.earth_density",MaterialsEarthDensityEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_freshwater",MaterialsRhoFreshwaterEnum)); break; } Index: /issm/trunk-jpl/src/m/classes/materials.m =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.m (revision 25149) +++ /issm/trunk-jpl/src/m/classes/materials.m (revision 25150) @@ -58,4 +58,5 @@ self.addprop('rho_ice'); self.addprop('rho_water'); + self.addprop('rho_freshwater'); self.addprop('earth_density'); otherwise @@ -138,4 +139,7 @@ % average density of the Earth (kg/m^3) self.earth_density=5512; + + %fresh water density (kg/m^3) + self.rho_freshwater=1000.; otherwise @@ -184,4 +188,5 @@ fielddisplay(self,'rho_water','ocean water density [kg/m^3]'); fielddisplay(self,'earth_density','mantle density [kg/m^3]'); + fielddisplay(self,'rho_freshwater','fresh water density [kg/m^3]'); otherwise @@ -233,4 +238,5 @@ md = checkfield(md,'fieldname','materials.rho_water','>',0); md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1); + md = checkfield(md,'fieldname','materials.rho_freshwater','>',0); otherwise @@ -280,4 +286,5 @@ WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double'); WriteData(fid,prefix,'object',self,'class','materials','fieldname','earth_density','format','Double'); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double'); otherwise @@ -334,4 +341,5 @@ writejsdouble(fid,[modelname '.materials.rho_water'],self.rho_water); writejsdouble(fid,[modelname '.materials.earth_density'],self.earth_density); + writejsdouble(fid,[modelname '.materials.rho_freshwater'],self.rho_freshwater); otherwise