Changeset 25958
- Timestamp:
- 01/27/21 17:49:18 (4 years ago)
- Location:
- issm/trunk-jpl/src
- Files:
-
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
r25947 r25958 795 795 void MasstransportAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ 796 796 797 797 798 /*Only update if on base*/ 798 799 if(!element->IsOnBase()) return; 799 800 800 801 /*Fetch dof list and allocate solution vector*/ 801 802 int *doflist = NULL; … … 805 806 IssmDouble* newthickness = xNew<IssmDouble>(numnodes); 806 807 IssmDouble* thicknessresidual = xNew<IssmDouble>(numnodes); 808 809 /*recover time step:*/ 810 IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum); 807 811 808 812 /*Use the dof list to index into the solution vector: */ … … 840 844 IssmDouble* oldthickness = xNew<IssmDouble>(numvertices); 841 845 IssmDouble* cumdeltathickness = xNew<IssmDouble>(numvertices); 842 IssmDouble* deltathickness= xNew<IssmDouble>(numvertices);846 IssmDouble* icethicknessrate = xNew<IssmDouble>(numvertices); 843 847 IssmDouble* newbase = xNew<IssmDouble>(numvertices); 844 848 IssmDouble* bed = xNew<IssmDouble>(numvertices); … … 866 870 for(int i=0;i<numvertices;i++){ 867 871 cumdeltathickness[i] += newthickness[i]-oldthickness[i]; 868 deltathickness[i] = newthickness[i]-oldthickness[i];872 icethicknessrate[i] = (newthickness[i]-oldthickness[i])/dt; 869 873 } 870 874 … … 903 907 element->AddBasalInput(BaseEnum,newbase,P1Enum); 904 908 element->AddBasalInput(SealevelchangeCumDeltathicknessEnum,cumdeltathickness,P1Enum); 905 element->AddBasalInput(SurfaceloadIceThickness ChangeEnum,deltathickness,P1Enum);909 element->AddBasalInput(SurfaceloadIceThicknessRateEnum,icethicknessrate,P1Enum); 906 910 907 911 /*Free ressources:*/ … … 911 915 xDelete<IssmDouble>(newsurface); 912 916 xDelete<IssmDouble>(oldthickness); 913 xDelete<IssmDouble>( deltathickness);917 xDelete<IssmDouble>(icethicknessrate); 914 918 xDelete<IssmDouble>(oldbase); 915 919 xDelete<IssmDouble>(oldsurface); -
issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp
r25956 r25958 40 40 iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); 41 41 iomodel->FetchData(&geodetic,"md.solidearth.settings.computesealevelchange"); 42 iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.icethicknesschange",SurfaceloadIceThicknessChangeEnum);43 42 iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0); 44 43 iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); 45 iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.waterheightchange",SurfaceloadWaterHeightChangeEnum); 44 iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.icethicknesschange",SurfaceloadIceThicknessRateEnum); 45 iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.waterheightchange",SurfaceloadWaterHeightRateEnum); 46 iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.other",SurfaceloadOtherRateEnum); 46 47 47 48 /*dynamic sea level: */ -
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r25956 r25958 5523 5523 /*}}}*/ 5524 5524 void Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/ 5525 5526 5525 5527 /*early return if we are not on an ice cap OR ocean:*/ 5526 5528 if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){ … … 5599 5601 } 5600 5602 else if(masks->isiceonly[this->lid]){ 5601 IssmDouble rho_ice, I; 5602 5603 /*recover material parameters: */ 5603 IssmDouble rho_ice, dIdt, dt; 5604 5605 5606 /*recover parameters: */ 5604 5607 rho_ice=FindParam(MaterialsRhoIceEnum); 5608 dt=FindParam(TimesteppingTimeStepEnum); 5605 5609 5606 5610 /*Compute ice thickness change: */ 5607 Input* deltathickness_input=this->GetInput(SurfaceloadIceThickness ChangeEnum);5611 Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessRateEnum); 5608 5612 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!"); 5609 deltathickness_input->GetInputAverage(&I); 5610 5611 dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; 5612 dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea; 5613 dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea; 5613 deltathickness_input->GetInputAverage(&dIdt); 5614 5615 5616 dI_list[0] = -4*PI*(rho_ice*dIdt*dt*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; 5617 dI_list[1] = -4*PI*(rho_ice*dIdt*dt*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea; 5618 dI_list[2] = +4*PI*(rho_ice*dIdt*dt*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea; 5614 5619 } 5615 5620 … … 5819 5824 IssmDouble area; 5820 5825 IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic 5821 IssmDouble I ; //change in ice thickness or water level(Farrel and Clarke, Equ. 4)5826 IssmDouble Idot; //ice thickness rate (Farrel and Clarke, Equ. 4) 5822 5827 bool notfullygrounded=false; 5823 5828 bool scaleoceanarea= false; 5824 5829 bool computerigid= false; 5825 5830 int glfraction=1; 5831 IssmDouble dt=0; 5826 5832 5827 5833 /*output: */ … … 5868 5874 rho_ice=FindParam(MaterialsRhoIceEnum); 5869 5875 rho_water=FindParam(MaterialsRhoSeawaterEnum); 5876 dt=FindParam(TimesteppingTimeStepEnum); 5870 5877 this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); 5871 5878 this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); … … 5892 5899 5893 5900 /*Retrieve ice thickness at vertices: */ 5894 Input* deltathickness_input=this->GetInput(SurfaceloadIceThickness ChangeEnum);5901 Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessRateEnum); 5895 5902 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!"); 5896 5903 5897 5904 /*/Average ice thickness over grounded area of the element only: {{{*/ 5898 if(!notfullygrounded)deltathickness_input->GetInputAverage(&I );5905 if(!notfullygrounded)deltathickness_input->GetInputAverage(&Idot); 5899 5906 else{ 5900 5907 IssmDouble total_weight=0; … … 5909 5916 /* Start looping on the number of gaussian points and average over these gaussian points: */ 5910 5917 total_weight=0; 5911 I =0;5918 Idot=0; 5912 5919 while(gauss->next()){ 5913 IssmDouble I g=0;5914 deltathickness_input->GetInputValue(&I g,gauss);5915 I +=Ig*gauss->weight;5920 IssmDouble Idotg=0; 5921 deltathickness_input->GetInputValue(&Idotg,gauss); 5922 Idot+=Idotg*gauss->weight; 5916 5923 total_weight+=gauss->weight; 5917 5924 } 5918 I =I/total_weight;5925 Idot=Idot/total_weight; 5919 5926 delete gauss; 5920 5927 } … … 5924 5931 _assert_(oceanarea>0.); 5925 5932 if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 5926 bslcice = rho_ice*area*phi*I /(oceanarea*rho_water);5933 bslcice = rho_ice*area*phi*Idot*dt/(oceanarea*rho_water); 5927 5934 _assert_(!xIsNan<IssmDouble>(bslcice)); 5928 5935 5929 5936 if(computerigid){ 5930 5937 /*convert from m to kg/m^2:*/ 5931 I =I*rho_ice*phi;5938 Idot=Idot*rho_ice*phi; 5932 5939 5933 5940 /*convolve:*/ 5934 for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*I ;5941 for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*Idot*dt; 5935 5942 } 5936 5943 … … 5950 5957 IssmDouble area; 5951 5958 IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic 5952 IssmDouble W; //change in water height thickness (Farrel and Clarke, Equ. 4) 5959 IssmDouble Wdot; //change in water height thickness (Farrel and Clarke, Equ. 4) 5960 IssmDouble dt=0; 5953 5961 bool notfullygrounded=false; 5954 5962 bool scaleoceanarea= false; … … 5983 5991 rho_water=FindParam(MaterialsRhoSeawaterEnum); 5984 5992 rho_freshwater=FindParam(MaterialsRhoFreshwaterEnum); 5993 dt=FindParam(TimesteppingTimeStepEnum); 5985 5994 this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); 5986 5995 this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); … … 5993 6002 5994 6003 /*Retrieve water height at vertices: */ 5995 Input* deltathickness_input=this->GetInput(SurfaceloadWaterHeight ChangeEnum);5996 if (!deltathickness_input)_error_("SurfaceloadWaterHeight ChangeEnum input needed to compute sea level change!");5997 deltathickness_input->GetInputAverage(&W );6004 Input* deltathickness_input=this->GetInput(SurfaceloadWaterHeightRateEnum); 6005 if (!deltathickness_input)_error_("SurfaceloadWaterHeightRateEnum input needed to compute sea level change!"); 6006 deltathickness_input->GetInputAverage(&Wdot); 5998 6007 5999 6008 /*Compute barystatic component:*/ 6000 6009 _assert_(oceanarea>0.); 6001 6010 if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 6002 bslchydro = rho_freshwater*area*phi*W /(oceanarea*rho_water);6011 bslchydro = rho_freshwater*area*phi*Wdot*dt/(oceanarea*rho_water); 6003 6012 _assert_(!xIsNan<IssmDouble>(bslchydro)); 6004 6013 6005 6014 if(computeelastic){ 6006 6015 /*convert from m to kg/m^2:*/ 6007 W =W*rho_freshwater*phi;6016 Wdot=Wdot*rho_freshwater*phi; 6008 6017 6009 6018 /*convolve:*/ 6010 for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*W ;6019 for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*Wdot*dt; 6011 6020 } 6012 6021 … … 6118 6127 /*diverse:*/ 6119 6128 int gsize; 6120 IssmDouble I , S, BP; //change in relative ice thickness and sea level6129 IssmDouble Idot, S, BP; //change in relative ice thickness and sea level 6121 6130 IssmDouble rho_ice,rho_water; 6122 6131 int horiz; … … 6182 6191 } 6183 6192 else if (masks->isiceonly[this->lid]){ 6193 IssmDouble dt=0; 6194 dt=FindParam(TimesteppingTimeStepEnum); 6184 6195 6185 6196 /*Compute ice thickness change: */ 6186 Input* deltathickness_input=this->GetInput(SurfaceloadIceThickness ChangeEnum);6197 Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessRateEnum); 6187 6198 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!"); 6188 deltathickness_input->GetInputAverage(&I );6199 deltathickness_input->GetInputAverage(&Idot); 6189 6200 6190 6201 /*convert to kg/m^2*/ 6191 I =I*rho_ice;6202 Idot=Idot*rho_ice; 6192 6203 6193 6204 for(int i=0;i<gsize;i++){ 6194 Up[i]+=I *GU[i];6205 Up[i]+=Idot*dt*GU[i]; 6195 6206 if(horiz){ 6196 North[i]+=I *GN[i];6197 East[i]+=I *GE[i];6207 North[i]+=Idot*dt*GN[i]; 6208 East[i]+=Idot*dt*GE[i]; 6198 6209 } 6199 6210 } -
issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp
r25957 r25958 255 255 256 256 /*we have accumulated thicknesses, dump them in deltathcikness: */ 257 if(modelid==earthid)InputDuplicatex(femmodel,SealevelchangeCumDeltathicknessEnum,SurfaceloadIceThickness ChangeEnum);257 if(modelid==earthid)InputDuplicatex(femmodel,SealevelchangeCumDeltathicknessEnum,SurfaceloadIceThicknessRateEnum); 258 258 } 259 259 -
issm/trunk-jpl/src/c/shared/Enum/Enum.vim
r25956 r25958 734 734 syn keyword cConstant SealevelchangeCumDeltathicknessEnum 735 735 syn keyword cConstant SealevelchangeCumDeltathicknessOldEnum 736 syn keyword cConstant SurfaceloadOtherEnum 737 syn keyword cConstant SurfaceloadIceThicknessChangeEnum 738 syn keyword cConstant SurfaceloadWaterHeightChangeEnum 736 syn keyword cConstant SurfaceloadRateEnum 737 syn keyword cConstant SurfaceloadIceThicknessRateEnum 738 syn keyword cConstant SurfaceloadWaterHeightRateEnum 739 syn keyword cConstant SurfaceloadOtherRateEnum 739 740 syn keyword cConstant SealevelchangeIndicesEnum 740 741 syn keyword cConstant SealevelchangeGEnum -
issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
r25956 r25958 730 730 SealevelchangeCumDeltathicknessEnum, 731 731 SealevelchangeCumDeltathicknessOldEnum, 732 SurfaceloadOtherEnum, 733 SurfaceloadIceThicknessChangeEnum, 734 SurfaceloadWaterHeightChangeEnum, 732 SurfaceloadRateEnum, 733 SurfaceloadIceThicknessRateEnum, 734 SurfaceloadWaterHeightRateEnum, 735 SurfaceloadOtherRateEnum, 735 736 SealevelchangeIndicesEnum, 736 737 SealevelchangeGEnum, -
issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
r25956 r25958 736 736 case SealevelchangeCumDeltathicknessEnum : return "SealevelchangeCumDeltathickness"; 737 737 case SealevelchangeCumDeltathicknessOldEnum : return "SealevelchangeCumDeltathicknessOld"; 738 case SurfaceloadOtherEnum : return "SurfaceloadOther"; 739 case SurfaceloadIceThicknessChangeEnum : return "SurfaceloadIceThicknessChange"; 740 case SurfaceloadWaterHeightChangeEnum : return "SurfaceloadWaterHeightChange"; 738 case SurfaceloadRateEnum : return "SurfaceloadRate"; 739 case SurfaceloadIceThicknessRateEnum : return "SurfaceloadIceThicknessRate"; 740 case SurfaceloadWaterHeightRateEnum : return "SurfaceloadWaterHeightRate"; 741 case SurfaceloadOtherRateEnum : return "SurfaceloadOtherRate"; 741 742 case SealevelchangeIndicesEnum : return "SealevelchangeIndices"; 742 743 case SealevelchangeGEnum : return "SealevelchangeG"; -
issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
r25956 r25958 751 751 else if (strcmp(name,"SealevelchangeCumDeltathickness")==0) return SealevelchangeCumDeltathicknessEnum; 752 752 else if (strcmp(name,"SealevelchangeCumDeltathicknessOld")==0) return SealevelchangeCumDeltathicknessOldEnum; 753 else if (strcmp(name,"Surfaceload Other")==0) return SurfaceloadOtherEnum;753 else if (strcmp(name,"SurfaceloadRate")==0) return SurfaceloadRateEnum; 754 754 else stage=7; 755 755 } 756 756 if(stage==7){ 757 if (strcmp(name,"SurfaceloadIceThicknessChange")==0) return SurfaceloadIceThicknessChangeEnum; 758 else if (strcmp(name,"SurfaceloadWaterHeightChange")==0) return SurfaceloadWaterHeightChangeEnum; 757 if (strcmp(name,"SurfaceloadIceThicknessRate")==0) return SurfaceloadIceThicknessRateEnum; 758 else if (strcmp(name,"SurfaceloadWaterHeightRate")==0) return SurfaceloadWaterHeightRateEnum; 759 else if (strcmp(name,"SurfaceloadOtherRate")==0) return SurfaceloadOtherRateEnum; 759 760 else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum; 760 761 else if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum; … … 874 875 else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum; 875 876 else if (strcmp(name,"SolidearthExternalBarystaticSeaLevelRate")==0) return SolidearthExternalBarystaticSeaLevelRateEnum; 876 else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;877 877 else stage=8; 878 878 } 879 879 if(stage==8){ 880 if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; 880 if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum; 881 else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; 881 882 else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; 882 883 else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; … … 997 998 else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum; 998 999 else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum; 999 else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;1000 1000 else stage=9; 1001 1001 } 1002 1002 if(stage==9){ 1003 if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum; 1003 if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum; 1004 else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum; 1004 1005 else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum; 1005 1006 else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum; … … 1120 1121 else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; 1121 1122 else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; 1122 else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;1123 1123 else stage=10; 1124 1124 } 1125 1125 if(stage==10){ 1126 if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; 1126 if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum; 1127 else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; 1127 1128 else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; 1128 1129 else if (strcmp(name,"Element")==0) return ElementEnum; … … 1243 1244 else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; 1244 1245 else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum; 1245 else if (strcmp(name,"Materials")==0) return MaterialsEnum;1246 1246 else stage=11; 1247 1247 } 1248 1248 if(stage==11){ 1249 if (strcmp(name,"Matestar")==0) return MatestarEnum; 1249 if (strcmp(name,"Materials")==0) return MaterialsEnum; 1250 else if (strcmp(name,"Matestar")==0) return MatestarEnum; 1250 1251 else if (strcmp(name,"Matice")==0) return MaticeEnum; 1251 1252 else if (strcmp(name,"Matlitho")==0) return MatlithoEnum; … … 1366 1367 else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum; 1367 1368 else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum; 1368 else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum;1369 1369 else stage=12; 1370 1370 } 1371 1371 if(stage==12){ 1372 if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum; 1372 if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum; 1373 else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum; 1373 1374 else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum; 1374 1375 else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum; -
issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp
r25784 r25958 163 163 fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1)); 164 164 } 165 else if(strcmp(string_in,"SurfaceloadIceThickness Change")==0){165 else if(strcmp(string_in,"SurfaceloadIceThicknessRate")==0){ 166 166 const char* field = "md.solidearth.surfaceload.icethicknesschange"; 167 input_enum = SurfaceloadIceThickness ChangeEnum;167 input_enum = SurfaceloadIceThicknessRateEnum; 168 168 fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1)); 169 169 } … … 178 178 fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1)); 179 179 } 180 else if(strcmp(string_in,"SurfaceloadWaterHeight Change")==0){180 else if(strcmp(string_in,"SurfaceloadWaterHeightRate")==0){ 181 181 const char* field = "md.solidearth.surfaceload.waterheightchange"; 182 input_enum = SurfaceloadWaterHeight ChangeEnum;182 input_enum = SurfaceloadWaterHeightRateEnum; 183 183 fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1)); 184 184 } -
issm/trunk-jpl/src/m/classes/mmeadditionalsolidearthsolution.m
r25956 r25958 69 69 function marshall(self,prefix,md,fid) % {{{ 70 70 WriteData(fid,prefix,'object',self,'data',3,'name','md.solidearth.external.nature','format','Integer'); %code 3 for mmeadditionalsolidearthsolution class 71 WriteData(fid,prefix,'name','md.solidearth.external.nummodels','data',length(self.displacementeast),'format','Integer'); 71 nummodels=length(self.displacementeast); 72 WriteData(fid,prefix,'name','md.solidearth.external.nummodels','data',nummodels,'format','Integer'); 72 73 WriteData(fid,prefix,'object',self,'fieldname','modelid','format','Double'); 73 WriteData(fid,prefix,'object',self,'fieldname','displacementeast','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); 74 WriteData(fid,prefix,'object',self,'fieldname','displacementup','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); 75 WriteData(fid,prefix,'object',self,'fieldname','displacementnorth','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); 76 WriteData(fid,prefix,'object',self,'fieldname','geoid','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); 77 WriteData(fid,prefix,'object',self,'fieldname','barystaticsealevel','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); 74 75 %transform our cell array of time series into cell array of time series of rates 76 for i=1:nummodels, 77 displacementeast=self.displacementeast{i}; 78 displacementnorth=self.displacementnorth{i}; 79 displacementup=self.displacementup{i}; 80 geoid=self.geoid{i}; 81 barystaticsealevel=self.barystaticsealevel{i}; 82 83 time=displacementeast(end,:); 84 dt=diff(time,1,2); 85 86 displacementeast_rate=diff(displacementeast(1:end-1,:),1,2)./dt; 87 displacementnorth_rate=diff(displacementnorth(1:end-1,:),1,2)./dt; 88 displacementup_rate=diff(displacementup(1:end-1,:),1,2)./dt; 89 geoid_rate=diff(geoid(1:end-1,:),1,2)./dt; 90 barystaticsealevel_rate=diff(barystaticsealevel(1:end-1,:),1,2)./dt; 91 92 self.displacementeast{i}=displacementeast_rate; 93 self.displacementnorth{i}=displacementnorth_rate; 94 self.displacementup{i}=displacementup_rate; 95 self.geoid{i}=geoid_rate; 96 self.barystaticsealevel{i}=barystaticsealevel_rate; 97 end 98 99 WriteData(fid,prefix,'object',self,'fieldname','displacementeast','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1/yts); 100 WriteData(fid,prefix,'object',self,'fieldname','displacementup','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1/yts); 101 WriteData(fid,prefix,'object',self,'fieldname','displacementnorth','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1/yts); 102 WriteData(fid,prefix,'object',self,'fieldname','geoid','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1/yts); 103 WriteData(fid,prefix,'object',self,'fieldname','barystaticsealevel','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1/yts); 78 104 79 105 end % }}} 106 80 107 function savemodeljs(self,fid,modelname) % {{{ 81 108 error('mmeadditionalsolidearthsolution error message: not implemented yet'); -
issm/trunk-jpl/src/m/classes/mmeofflinesolidearthsolution.m
r25956 r25958 70 70 WriteData(fid,prefix,'object',self,'data',4,'name','md.solidearth.external.nature','format','Integer'); %code 4 for mmeofflinesolidearthsolution class 71 71 WriteData(fid,prefix,'object',self,'fieldname','modelid','format','Double'); 72 WriteData(fid,prefix,'name','md.solidearth.external.nummodels','data',length(self.displacementeast),'format','Integer'); 73 WriteData(fid,prefix,'object',self,'fieldname','displacementeast','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); 74 WriteData(fid,prefix,'object',self,'fieldname','displacementup','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); 75 WriteData(fid,prefix,'object',self,'fieldname','displacementnorth','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); 76 WriteData(fid,prefix,'object',self,'fieldname','geoid','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); 77 WriteData(fid,prefix,'object',self,'fieldname','barystaticsealevel','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3); 72 nummodels=length(self.displacementeast); 73 WriteData(fid,prefix,'name','md.solidearth.external.nummodels','data',nummodels,'format','Integer'); 74 75 %transform our cell array of time series into cell array of time series of rates 76 for i=1:nummodels, 77 displacementeast=self.displacementeast{i}; 78 displacementnorth=self.displacementnorth{i}; 79 displacementup=self.displacementup{i}; 80 geoid=self.geoid{i}; 81 barystaticsealevel=self.barystaticsealevel{i}; 82 83 time=displacementeast(end,:); 84 dt=diff(time,1,2); 85 86 displacementeast_rate=diff(displacementeast(1:end-1,:),1,2)./dt; 87 displacementnorth_rate=diff(displacementnorth(1:end-1,:),1,2)./dt; 88 displacementup_rate=diff(displacementup(1:end-1,:),1,2)./dt; 89 geoid_rate=diff(geoid(1:end-1,:),1,2)./dt; 90 barystaticsealevel_rate=diff(barystaticsealevel(1:end-1,:),1,2)./dt; 91 92 self.displacementeast{i}=displacementeast_rate; 93 self.displacementnorth{i}=displacementnorth_rate; 94 self.displacementup{i}=displacementup_rate; 95 self.geoid{i}=geoid_rate; 96 self.barystaticsealevel{i}=barystaticsealevel_rate; 97 end 98 99 WriteData(fid,prefix,'object',self,'fieldname','displacementeast','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1/yts); 100 WriteData(fid,prefix,'object',self,'fieldname','displacementup','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1/yts); 101 WriteData(fid,prefix,'object',self,'fieldname','displacementnorth','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1/yts); 102 WriteData(fid,prefix,'object',self,'fieldname','geoid','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1/yts); 103 WriteData(fid,prefix,'object',self,'fieldname','barystaticsealevel','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1/yts); 78 104 79 105 end % }}} -
issm/trunk-jpl/src/m/classes/surfaceload.m
r25956 r25958 65 65 end % }}} 66 66 function marshall(self,prefix,md,fid) % {{{ 67 67 68 %deal with ice thickness change: {{{ 68 69 if isempty(self.icethicknesschange), 69 70 self.icethicknesschange=zeros(md.mesh.numberofelements+1,1); 70 71 end 72 73 if isa(self.icethicknesschange,'cell'), 74 %transform our cell array of time series into cell array of time series of rates 75 nummodels=length(self.icethicknesschange); 76 for i=1:nummodels, 77 icethicknesschange=self.icethicknesschange{i}; 78 time=icethicknesschange(end,:); 79 dt=diff(time,1,2); 80 icethicknesschange_rate=diff(icethicknesschange(1:end-1,:),1,2)./dt; 81 self.icethicknesschange{i}=icethicknesschange_rate; 82 end 83 WriteData(fid,prefix,'object',self,'fieldname','icethicknesschange','name','md.solidearth.surfaceload.icethicknesschange',... 84 'format','MatArray','timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts,'scale',1/yts); 85 else 86 icethicknesschange=self.icethicknesschange; 87 time=icethicknesschange(end,:); 88 dt=diff(time,1,2); 89 icethicknesschange_rate=diff(icethicknesschange(1:end-1,:),1,2)./dt; 90 self.icethicknesschange=icethicknesschange_rate; 91 92 WriteData(fid,prefix,'object',self,'fieldname','icethicknesschange','name','md.solidearth.surfaceload.icethicknesschange',... 93 'format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts,'scale',1/yts); 94 end 95 %}}} 96 %deal with water height change: {{{ 71 97 if isempty(self.waterheightchange), 72 98 self.waterheightchange=zeros(md.mesh.numberofelements+1,1); 73 99 end 100 101 if isa(self.waterheightchange,'cell'), 102 %transform our cell array of time series into cell array of time series of rates 103 nummodels=length(self.waterheightchange); 104 for i=1:nummodels, 105 waterheightchange=self.waterheightchange{i}; 106 time=waterheightchange(end,:); 107 dt=diff(time,1,2); 108 waterheightchange_rate=diff(waterheightchange(1:end-1,:),1,2)./dt; 109 self.waterheightchange{i}=waterheightchange_rate; 110 end 111 WriteData(fid,prefix,'object',self,'fieldname','waterheightchange','name','md.solidearth.surfaceload.waterheightchange',... 112 'format','MatArray','timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts,'scale',1/yts); 113 else 114 waterheightchange=self.waterheightchange; 115 time=waterheightchange(end,:); 116 dt=diff(time,1,2); 117 waterheightchange_rate=diff(waterheightchange(1:end-1,:),1,2)./dt; 118 self.waterheightchange=waterheightchange_rate; 119 120 WriteData(fid,prefix,'object',self,'fieldname','waterheightchange','name','md.solidearth.surfaceload.waterheightchange',... 121 'format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts,'scale',1/yts); 122 end 123 %}}} 124 %deal with other: {{{ 74 125 if isempty(self.other), 75 126 self.other=zeros(md.mesh.numberofelements+1,1); 76 127 end 77 if isa(self.icethicknesschange,'cell'), 78 WriteData(fid,prefix,'object',self,'fieldname','icethicknesschange','name','md.solidearth.surfaceload.icethicknesschange',... 79 'format','MatArray','timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); 128 129 if isa(self.otherchange,'cell'), 130 %transform our cell array of time series into cell array of time series of rates 131 nummodels=length(self.otherchange); 132 for i=1:nummodels, 133 otherchange=self.otherchange{i}; 134 time=otherchange(end,:); 135 dt=diff(time,1,2); 136 otherchange_rate=diff(otherchange(1:end-1,:),1,2)./dt; 137 self.otherchange{i}=otherchange_rate; 138 end 139 WriteData(fid,prefix,'object',self,'fieldname','otherchange','name','md.solidearth.surfaceload.otherchange',... 140 'format','MatArray','timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts,'scale',1/yts); 80 141 else 81 WriteData(fid,prefix,'object',self,'fieldname','icethicknesschange','name','md.solidearth.surfaceload.icethicknesschange',... 82 'format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); 142 otherchange=self.otherchange; 143 time=otherchange(end,:); 144 dt=diff(time,1,2); 145 otherchange_rate=diff(otherchange(1:end-1,:),1,2)./dt; 146 self.otherchange=otherchange_rate; 147 148 WriteData(fid,prefix,'object',self,'fieldname','otherchange','name','md.solidearth.surfaceload.otherchange',... 149 'format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts,'scale',1/yts); 83 150 end 84 85 WriteData(fid,prefix,'object',self,'fieldname','waterheightchange','name','md.solidearth.surfaceload.waterheightchange',... 86 'format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); 87 WriteData(fid,prefix,'object',self,'fieldname','other','name','md.solidearth.surfaceload.other',... 88 'format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); 151 %}}} 89 152 90 153 end % }}}
Note:
See TracChangeset
for help on using the changeset viewer.