Changeset 25958 for issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

Timestamp:

01/27/21 17:49:18 (4 years ago)

Author:

Eric.Larour

Message:

CHG: switched to rates on solidearth surafce loads.

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (12 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -5523 +5523 @@
 /*}}}*/
 void    Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
+
+
         /*early return if we are not on an ice cap OR ocean:*/
         if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){
@@ -5599 +5601 @@
         }
         else if(masks->isiceonly[this->lid]){
-                IssmDouble rho_ice, I;
-
-                /*recover material parameters: */
+                IssmDouble rho_ice, dIdt, dt;
+                
+
+                /*recover parameters: */
                 rho_ice=FindParam(MaterialsRhoIceEnum);
+                dt=FindParam(TimesteppingTimeStepEnum);
 
                 /*Compute ice thickness change: */
-                Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessChangeEnum);
+                Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessRateEnum);
                 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!");
-                deltathickness_input->GetInputAverage(&I);
-
-                dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
-                dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
-                dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
+                deltathickness_input->GetInputAverage(&dIdt);
+
+
+                dI_list[0] = -4*PI*(rho_ice*dIdt*dt*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
+                dI_list[1] = -4*PI*(rho_ice*dIdt*dt*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
+                dI_list[2] = +4*PI*(rho_ice*dIdt*dt*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
         }
 
@@ -5819 +5824 @@
         IssmDouble area;
         IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic
-        IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
+        IssmDouble Idot;  //ice thickness rate (Farrel and Clarke, Equ. 4)
         bool notfullygrounded=false;
         bool scaleoceanarea= false;
         bool computerigid= false;
         int  glfraction=1;
+        IssmDouble dt=0;
 
         /*output: */
@@ -5868 +5874 @@
         rho_ice=FindParam(MaterialsRhoIceEnum);
         rho_water=FindParam(MaterialsRhoSeawaterEnum);
+        dt=FindParam(TimesteppingTimeStepEnum);
         this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
         this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
@@ -5892 +5899 @@
 
         /*Retrieve ice thickness at vertices: */
-        Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessChangeEnum);
+        Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessRateEnum);
         if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!");
 
         /*/Average ice thickness over grounded area of the element only: {{{*/
-        if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
+        if(!notfullygrounded)deltathickness_input->GetInputAverage(&Idot);
         else{
                 IssmDouble total_weight=0;
@@ -5909 +5916 @@
                 /* Start  looping on the number of gaussian points and average over these gaussian points: */
                 total_weight=0;
-                I=0;
+                Idot=0;
                 while(gauss->next()){
-                        IssmDouble Ig=0;
-                        deltathickness_input->GetInputValue(&Ig,gauss);
-                        I+=Ig*gauss->weight;
+                        IssmDouble Idotg=0;
+                        deltathickness_input->GetInputValue(&Idotg,gauss);
+                        Idot+=Idotg*gauss->weight;
                         total_weight+=gauss->weight;
                 }
-                I=I/total_weight;
+                Idot=Idot/total_weight;
                 delete gauss;
         }
@@ -5924 +5931 @@
         _assert_(oceanarea>0.);
         if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
-        bslcice = rho_ice*area*phi*I/(oceanarea*rho_water);
+        bslcice = rho_ice*area*phi*Idot*dt/(oceanarea*rho_water);
         _assert_(!xIsNan<IssmDouble>(bslcice));
 
         if(computerigid){
                 /*convert from m to kg/m^2:*/
-                I=I*rho_ice*phi;
+                Idot=Idot*rho_ice*phi;
 
                 /*convolve:*/
-                for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*I;
+                for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*Idot*dt;
         }
 
@@ -5950 +5957 @@
         IssmDouble area;
         IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic
-        IssmDouble W;  //change in water height thickness (Farrel and Clarke, Equ. 4)
+        IssmDouble Wdot;  //change in water height thickness (Farrel and Clarke, Equ. 4)
+        IssmDouble dt=0;
         bool notfullygrounded=false;
         bool scaleoceanarea= false;
@@ -5983 +5991 @@
         rho_water=FindParam(MaterialsRhoSeawaterEnum);
         rho_freshwater=FindParam(MaterialsRhoFreshwaterEnum);
+        dt=FindParam(TimesteppingTimeStepEnum);
         this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
         this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
@@ -5993 +6002 @@
 
         /*Retrieve water height at vertices: */
-        Input* deltathickness_input=this->GetInput(SurfaceloadWaterHeightChangeEnum);
-        if (!deltathickness_input)_error_("SurfaceloadWaterHeightChangeEnum input needed to compute sea level change!");
-        deltathickness_input->GetInputAverage(&W);
+        Input* deltathickness_input=this->GetInput(SurfaceloadWaterHeightRateEnum);
+        if (!deltathickness_input)_error_("SurfaceloadWaterHeightRateEnum input needed to compute sea level change!");
+        deltathickness_input->GetInputAverage(&Wdot);
 
         /*Compute barystatic component:*/
         _assert_(oceanarea>0.);
         if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
-        bslchydro = rho_freshwater*area*phi*W/(oceanarea*rho_water);
+        bslchydro = rho_freshwater*area*phi*Wdot*dt/(oceanarea*rho_water);
         _assert_(!xIsNan<IssmDouble>(bslchydro));
 
         if(computeelastic){
                 /*convert from m to kg/m^2:*/
-                W=W*rho_freshwater*phi;
+                Wdot=Wdot*rho_freshwater*phi;
 
                 /*convolve:*/
-                for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*W;
+                for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*Wdot*dt;
         }
 
@@ -6118 +6127 @@
         /*diverse:*/
         int gsize;
-        IssmDouble I, S, BP;            //change in relative ice thickness and sea level
+        IssmDouble Idot, S, BP;         //change in relative ice thickness and sea level
         IssmDouble rho_ice,rho_water;
         int horiz;
@@ -6182 +6191 @@
         }
         else if (masks->isiceonly[this->lid]){
+                IssmDouble dt=0;
+                dt=FindParam(TimesteppingTimeStepEnum);
 
                 /*Compute ice thickness change: */
-                Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessChangeEnum);
+                Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessRateEnum);
                 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!");
-                deltathickness_input->GetInputAverage(&I);
+                deltathickness_input->GetInputAverage(&Idot);
 
                 /*convert to kg/m^2*/
-                I=I*rho_ice;
+                Idot=Idot*rho_ice;
 
                 for(int i=0;i<gsize;i++){
-                        Up[i]+=I*GU[i];
+                        Up[i]+=Idot*dt*GU[i];
                         if(horiz){
-                                North[i]+=I*GN[i];
-                                East[i]+=I*GE[i];
+                                North[i]+=Idot*dt*GN[i];
+                                East[i]+=Idot*dt*GE[i];
                         }
                 }