Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 21751) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 21752) @@ -36,4 +36,5 @@ iomodel->FetchDataToInput(elements,"md.slr.deltathickness",SealevelriseDeltathicknessEnum); iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(elements,"md.slr.steric_rate",SealevelriseStericRateEnum); }/*}}}*/ Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 21751) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 21752) @@ -15,4 +15,5 @@ Vector *Sg_absolute = NULL; Vector *Sg_eustatic = NULL; + Vector *steric_rate_g = NULL; Vector *U_radial = NULL; Vector *U_north = NULL; @@ -33,4 +34,5 @@ IssmDouble *yy = NULL; IssmDouble *zz = NULL; + IssmDouble dt; /*Recover some parameters: */ @@ -40,5 +42,5 @@ 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,9 +82,7 @@ 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,4 +95,9 @@ 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,5 +107,6 @@ 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,4 +125,5 @@ delete U_east; delete Sg_absolute; + delete steric_rate_g; if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} } Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21751) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21752) @@ -803,4 +803,5 @@ SealevelriseAngularVelocityEnum, SealevelriseOceanAreaScalingEnum, + SealevelriseStericRateEnum, SealevelriseGElasticEnum, SealevelriseUElasticEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21751) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21752) @@ -783,4 +783,5 @@ case SealevelriseAngularVelocityEnum : return "SealevelriseAngularVelocity"; case SealevelriseOceanAreaScalingEnum : return "SealevelriseOceanAreaScaling"; + case SealevelriseStericRateEnum : return "SealevelriseStericRate"; case SealevelriseGElasticEnum : return "SealevelriseGElastic"; case SealevelriseUElasticEnum : return "SealevelriseUElastic"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21751) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21752) @@ -801,4 +801,5 @@ 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,9 +875,9 @@ 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,9 +998,9 @@ 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; Index: /issm/trunk-jpl/src/m/classes/slr.m =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.m (revision 21751) +++ /issm/trunk-jpl/src/m/classes/slr.m (revision 21752) @@ -24,4 +24,5 @@ rotation = 0; ocean_area_scaling = 0; + steric_rate = 0; %rate of ocean expansion from steric effects. degacc = 0; requested_outputs = {}; @@ -68,5 +69,9 @@ %numerical discretization accuracy self.degacc=.01; - + + %steric: + self.steric_rate=0; + + %output default: self.requested_outputs={'default'}; @@ -93,4 +98,5 @@ 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,4 +140,5 @@ 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,4 +166,5 @@ 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,4 +200,5 @@ 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);