Changeset 21752

Timestamp:

06/01/17 15:40:48 (8 years ago)

Author:

Eric.Larour

Message:

CHG: introduce steric rate in the SLR solution.

Location:

issm/trunk-jpl/src

Files:

• c/analyses/SealevelriseAnalysis.cpp (modified) (1 diff)
• c/cores/sealevelrise_core.cpp (modified) (7 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (3 diffs)
• m/classes/slr.m (modified) (6 diffs)

issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp ¶

@@ -36 +36 @@
         iomodel->FetchDataToInput(elements,"md.slr.deltathickness",SealevelriseDeltathicknessEnum);
         iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0);
+        iomodel->FetchDataToInput(elements,"md.slr.steric_rate",SealevelriseStericRateEnum);
 
 }/*}}}*/

issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp ¶

@@ -15 +15 @@
         Vector<IssmDouble> *Sg_absolute  = NULL; 
         Vector<IssmDouble> *Sg_eustatic  = NULL; 
+        Vector<IssmDouble> *steric_rate_g  = NULL; 
         Vector<IssmDouble> *U_radial  = NULL; 
         Vector<IssmDouble> *U_north   = NULL; 
@@ -33 +34 @@
         IssmDouble          *yy     = NULL;
         IssmDouble          *zz     = NULL;
+        IssmDouble          dt;
 
         /*Recover some parameters: */
@@ -40 +42 @@
         femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
         femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
-
+        
         /*first, recover lat,long and radius vectors from vertices: */
         VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); 
@@ -80 +82 @@
 
                 Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //ocean loading tems  (2nd and 5th terms on the RHS of Farrel and Clark)
-                
-                /*get results into elements:*/
-                //InputUpdateFromSolutionx(femmodel,Sg);                // from Eric 
-                InputUpdateFromVectorx(femmodel,Sg,SealevelEnum,VertexSIdEnum);
-
+
+                                
+        
                 /*compute other geodetic signatures, such as absolute sea level chagne, components of 3-D crustal motion: */
                 /*Initialize:*/
@@ -95 +95 @@
                 femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz); 
 
+                /*Now deal with steric ocean expansion by just shifting Sg by a spatial rate pattern : */
+                femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+                GetVectorFromInputsx(&steric_rate_g,femmodel,SealevelriseStericRateEnum,VertexPIdEnum);
+                Sg->AXPY(steric_rate_g,dt);
+                
                 /*compute: absolute sea level change = relative sea level change + vertical motion*/
                 Sg->Copy(Sg_absolute); Sg_absolute->AXPY(U_radial,1); 
@@ -102 +107 @@
                 InputUpdateFromVectorx(femmodel,U_north,SealevelNmotionEnum,VertexSIdEnum);     // north motion 
                 InputUpdateFromVectorx(femmodel,U_east,SealevelEmotionEnum,VertexSIdEnum);              // east motion 
-                InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelAbsoluteEnum,VertexSIdEnum);
+                InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelAbsoluteEnum,VertexSIdEnum); //absolute sea level
+                InputUpdateFromVectorx(femmodel,Sg,SealevelEnum,VertexSIdEnum); //relative sea level
                 
                 if(save_results){
@@ -119 +125 @@
                 delete U_east;
                 delete Sg_absolute;
+                delete steric_rate_g;
                 if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
         }

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h ¶

@@ -803 +803 @@
         SealevelriseAngularVelocityEnum,
         SealevelriseOceanAreaScalingEnum,
+        SealevelriseStericRateEnum,
         SealevelriseGElasticEnum,
         SealevelriseUElasticEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp ¶

@@ -783 +783 @@
                 case SealevelriseAngularVelocityEnum : return "SealevelriseAngularVelocity";
                 case SealevelriseOceanAreaScalingEnum : return "SealevelriseOceanAreaScaling";
+                case SealevelriseStericRateEnum : return "SealevelriseStericRate";
                 case SealevelriseGElasticEnum : return "SealevelriseGElastic";
                 case SealevelriseUElasticEnum : return "SealevelriseUElastic";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp ¶

@@ -801 +801 @@
               else if (strcmp(name,"SealevelriseAngularVelocity")==0) return SealevelriseAngularVelocityEnum;
               else if (strcmp(name,"SealevelriseOceanAreaScaling")==0) return SealevelriseOceanAreaScalingEnum;
+              else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum;
               else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum;
               else if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum;
@@ -874 +875 @@
               else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
               else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
-              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
+              if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
+              else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
               else if (strcmp(name,"StringParam")==0) return StringParamEnum;
               else if (strcmp(name,"Seg")==0) return SegEnum;
@@ -997 +998 @@
               else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
               else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
-              else if (strcmp(name,"P1P1")==0) return P1P1Enum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
+              if (strcmp(name,"P1P1")==0) return P1P1Enum;
+              else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
               else if (strcmp(name,"MINI")==0) return MINIEnum;
               else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;

issm/trunk-jpl/src/m/classes/slr.m ¶

@@ -24 +24 @@
                 rotation       = 0;
                 ocean_area_scaling = 0;
+                steric_rate    = 0; %rate of ocean expansion from steric effects. 
                 degacc         = 0;
                 requested_outputs      = {};
@@ -68 +69 @@
                 %numerical discretization accuracy
                 self.degacc=.01;
-                
+
+                %steric: 
+                self.steric_rate=0;
+        
+        
                 %output default:
                 self.requested_outputs={'default'};
@@ -93 +98 @@
                         md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
                         md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1);
+                        md = checkfield(md,'fieldname','slr.steric_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
                         md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10);
                         md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1);
@@ -134 +140 @@
                         fielddisplay(self,'rotation','earth rotational potential perturbation');
                         fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); 
+                        fielddisplay(self,'steric_rate','rate of steric ocean expansion (in mm/yr)'); 
                         fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
                         fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
@@ -159 +166 @@
                         WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean');
                         WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean');
+                        WriteData(fid,prefix,'object',self,'fieldname','steric_rate','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
                         WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double');
                         WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
@@ -192 +200 @@
                         writejsdouble(fid,[modelname '.slr.rotation'],self.rotation);
                         writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling);
+                        writejs1Darray(fid,[modelname '.slr.steric_rate'],self.steric_rate);
                         writejsdouble(fid,[modelname '.slr.degacc'],self.degacc);
                         writejscellstring(fid,[modelname '.slr.requested_outputs'],self.requested_outputs);