Changeset 25958


Ignore:
Timestamp:
01/27/21 17:49:18 (4 years ago)
Author:
Eric.Larour
Message:

CHG: switched to rates on solidearth surafce loads.

Location:
issm/trunk-jpl/src
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

    r25947 r25958  
    795795void           MasstransportAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    796796
     797
    797798        /*Only update if on base*/
    798799        if(!element->IsOnBase()) return;
    799 
     800       
    800801        /*Fetch dof list and allocate solution vector*/
    801802        int *doflist = NULL;
     
    805806        IssmDouble* newthickness = xNew<IssmDouble>(numnodes);
    806807        IssmDouble* thicknessresidual = xNew<IssmDouble>(numnodes);
     808       
     809        /*recover time step:*/
     810        IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum);
    807811
    808812        /*Use the dof list to index into the solution vector: */
     
    840844        IssmDouble* oldthickness      = xNew<IssmDouble>(numvertices);
    841845        IssmDouble* cumdeltathickness = xNew<IssmDouble>(numvertices);
    842         IssmDouble* deltathickness    = xNew<IssmDouble>(numvertices);
     846        IssmDouble* icethicknessrate    = xNew<IssmDouble>(numvertices);
    843847        IssmDouble* newbase           = xNew<IssmDouble>(numvertices);
    844848        IssmDouble* bed               = xNew<IssmDouble>(numvertices);
     
    866870        for(int i=0;i<numvertices;i++){
    867871                cumdeltathickness[i] += newthickness[i]-oldthickness[i];
    868                 deltathickness[i]     = newthickness[i]-oldthickness[i];
     872                icethicknessrate[i]     = (newthickness[i]-oldthickness[i])/dt;
    869873        }
    870874
     
    903907        element->AddBasalInput(BaseEnum,newbase,P1Enum);
    904908        element->AddBasalInput(SealevelchangeCumDeltathicknessEnum,cumdeltathickness,P1Enum);
    905         element->AddBasalInput(SurfaceloadIceThicknessChangeEnum,deltathickness,P1Enum);
     909        element->AddBasalInput(SurfaceloadIceThicknessRateEnum,icethicknessrate,P1Enum);
    906910
    907911        /*Free ressources:*/
     
    911915        xDelete<IssmDouble>(newsurface);
    912916        xDelete<IssmDouble>(oldthickness);
    913         xDelete<IssmDouble>(deltathickness);
     917        xDelete<IssmDouble>(icethicknessrate);
    914918        xDelete<IssmDouble>(oldbase);
    915919        xDelete<IssmDouble>(oldsurface);
  • issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp

    r25956 r25958  
    4040        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    4141        iomodel->FetchData(&geodetic,"md.solidearth.settings.computesealevelchange");
    42         iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.icethicknesschange",SurfaceloadIceThicknessChangeEnum);
    4342        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
    4443        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);
    4647               
    4748        /*dynamic sea level: */
  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r25956 r25958  
    55235523/*}}}*/
    55245524void    Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
     5525
     5526
    55255527        /*early return if we are not on an ice cap OR ocean:*/
    55265528        if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){
     
    55995601        }
    56005602        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: */
    56045607                rho_ice=FindParam(MaterialsRhoIceEnum);
     5608                dt=FindParam(TimesteppingTimeStepEnum);
    56055609
    56065610                /*Compute ice thickness change: */
    5607                 Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessChangeEnum);
     5611                Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessRateEnum);
    56085612                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;
    56145619        }
    56155620
     
    58195824        IssmDouble area;
    58205825        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)
    58225827        bool notfullygrounded=false;
    58235828        bool scaleoceanarea= false;
    58245829        bool computerigid= false;
    58255830        int  glfraction=1;
     5831        IssmDouble dt=0;
    58265832
    58275833        /*output: */
     
    58685874        rho_ice=FindParam(MaterialsRhoIceEnum);
    58695875        rho_water=FindParam(MaterialsRhoSeawaterEnum);
     5876        dt=FindParam(TimesteppingTimeStepEnum);
    58705877        this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    58715878        this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
     
    58925899
    58935900        /*Retrieve ice thickness at vertices: */
    5894         Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessChangeEnum);
     5901        Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessRateEnum);
    58955902        if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!");
    58965903
    58975904        /*/Average ice thickness over grounded area of the element only: {{{*/
    5898         if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
     5905        if(!notfullygrounded)deltathickness_input->GetInputAverage(&Idot);
    58995906        else{
    59005907                IssmDouble total_weight=0;
     
    59095916                /* Start  looping on the number of gaussian points and average over these gaussian points: */
    59105917                total_weight=0;
    5911                 I=0;
     5918                Idot=0;
    59125919                while(gauss->next()){
    5913                         IssmDouble Ig=0;
    5914                         deltathickness_input->GetInputValue(&Ig,gauss);
    5915                         I+=Ig*gauss->weight;
     5920                        IssmDouble Idotg=0;
     5921                        deltathickness_input->GetInputValue(&Idotg,gauss);
     5922                        Idot+=Idotg*gauss->weight;
    59165923                        total_weight+=gauss->weight;
    59175924                }
    5918                 I=I/total_weight;
     5925                Idot=Idot/total_weight;
    59195926                delete gauss;
    59205927        }
     
    59245931        _assert_(oceanarea>0.);
    59255932        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);
    59275934        _assert_(!xIsNan<IssmDouble>(bslcice));
    59285935
    59295936        if(computerigid){
    59305937                /*convert from m to kg/m^2:*/
    5931                 I=I*rho_ice*phi;
     5938                Idot=Idot*rho_ice*phi;
    59325939
    59335940                /*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;
    59355942        }
    59365943
     
    59505957        IssmDouble area;
    59515958        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;
    59535961        bool notfullygrounded=false;
    59545962        bool scaleoceanarea= false;
     
    59835991        rho_water=FindParam(MaterialsRhoSeawaterEnum);
    59845992        rho_freshwater=FindParam(MaterialsRhoFreshwaterEnum);
     5993        dt=FindParam(TimesteppingTimeStepEnum);
    59855994        this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
    59865995        this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
     
    59936002
    59946003        /*Retrieve water height at vertices: */
    5995         Input* deltathickness_input=this->GetInput(SurfaceloadWaterHeightChangeEnum);
    5996         if (!deltathickness_input)_error_("SurfaceloadWaterHeightChangeEnum 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);
    59986007
    59996008        /*Compute barystatic component:*/
    60006009        _assert_(oceanarea>0.);
    60016010        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);
    60036012        _assert_(!xIsNan<IssmDouble>(bslchydro));
    60046013
    60056014        if(computeelastic){
    60066015                /*convert from m to kg/m^2:*/
    6007                 W=W*rho_freshwater*phi;
     6016                Wdot=Wdot*rho_freshwater*phi;
    60086017
    60096018                /*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;
    60116020        }
    60126021
     
    61186127        /*diverse:*/
    61196128        int gsize;
    6120         IssmDouble I, S, BP;            //change in relative ice thickness and sea level
     6129        IssmDouble Idot, S, BP;         //change in relative ice thickness and sea level
    61216130        IssmDouble rho_ice,rho_water;
    61226131        int horiz;
     
    61826191        }
    61836192        else if (masks->isiceonly[this->lid]){
     6193                IssmDouble dt=0;
     6194                dt=FindParam(TimesteppingTimeStepEnum);
    61846195
    61856196                /*Compute ice thickness change: */
    6186                 Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessChangeEnum);
     6197                Input* deltathickness_input=this->GetInput(SurfaceloadIceThicknessRateEnum);
    61876198                if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!");
    6188                 deltathickness_input->GetInputAverage(&I);
     6199                deltathickness_input->GetInputAverage(&Idot);
    61896200
    61906201                /*convert to kg/m^2*/
    6191                 I=I*rho_ice;
     6202                Idot=Idot*rho_ice;
    61926203
    61936204                for(int i=0;i<gsize;i++){
    6194                         Up[i]+=I*GU[i];
     6205                        Up[i]+=Idot*dt*GU[i];
    61956206                        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];
    61986209                        }
    61996210                }
  • issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp

    r25957 r25958  
    255255
    256256                /*we have accumulated thicknesses, dump them in deltathcikness: */
    257                 if(modelid==earthid)InputDuplicatex(femmodel,SealevelchangeCumDeltathicknessEnum,SurfaceloadIceThicknessChangeEnum);
     257                if(modelid==earthid)InputDuplicatex(femmodel,SealevelchangeCumDeltathicknessEnum,SurfaceloadIceThicknessRateEnum);
    258258        }
    259259
  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r25956 r25958  
    734734syn keyword cConstant SealevelchangeCumDeltathicknessEnum
    735735syn keyword cConstant SealevelchangeCumDeltathicknessOldEnum
    736 syn keyword cConstant SurfaceloadOtherEnum
    737 syn keyword cConstant SurfaceloadIceThicknessChangeEnum
    738 syn keyword cConstant SurfaceloadWaterHeightChangeEnum
     736syn keyword cConstant SurfaceloadRateEnum
     737syn keyword cConstant SurfaceloadIceThicknessRateEnum
     738syn keyword cConstant SurfaceloadWaterHeightRateEnum
     739syn keyword cConstant SurfaceloadOtherRateEnum
    739740syn keyword cConstant SealevelchangeIndicesEnum
    740741syn keyword cConstant SealevelchangeGEnum
  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r25956 r25958  
    730730        SealevelchangeCumDeltathicknessEnum,
    731731        SealevelchangeCumDeltathicknessOldEnum,
    732         SurfaceloadOtherEnum,
    733         SurfaceloadIceThicknessChangeEnum,
    734         SurfaceloadWaterHeightChangeEnum,
     732        SurfaceloadRateEnum,
     733        SurfaceloadIceThicknessRateEnum,
     734        SurfaceloadWaterHeightRateEnum,
     735        SurfaceloadOtherRateEnum,
    735736        SealevelchangeIndicesEnum,
    736737        SealevelchangeGEnum,
  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r25956 r25958  
    736736                case SealevelchangeCumDeltathicknessEnum : return "SealevelchangeCumDeltathickness";
    737737                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";
    741742                case SealevelchangeIndicesEnum : return "SealevelchangeIndices";
    742743                case SealevelchangeGEnum : return "SealevelchangeG";
  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r25956 r25958  
    751751              else if (strcmp(name,"SealevelchangeCumDeltathickness")==0) return SealevelchangeCumDeltathicknessEnum;
    752752              else if (strcmp(name,"SealevelchangeCumDeltathicknessOld")==0) return SealevelchangeCumDeltathicknessOldEnum;
    753               else if (strcmp(name,"SurfaceloadOther")==0) return SurfaceloadOtherEnum;
     753              else if (strcmp(name,"SurfaceloadRate")==0) return SurfaceloadRateEnum;
    754754         else stage=7;
    755755   }
    756756   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;
    759760              else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
    760761              else if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum;
     
    874875              else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
    875876              else if (strcmp(name,"SolidearthExternalBarystaticSeaLevelRate")==0) return SolidearthExternalBarystaticSeaLevelRateEnum;
    876               else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
    877877         else stage=8;
    878878   }
    879879   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;
    881882              else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum;
    882883              else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum;
     
    997998              else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum;
    998999              else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum;
    999               else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
    10001000         else stage=9;
    10011001   }
    10021002   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;
    10041005              else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
    10051006              else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum;
     
    11201121              else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
    11211122              else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
    1122               else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
    11231123         else stage=10;
    11241124   }
    11251125   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;
    11271128              else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
    11281129              else if (strcmp(name,"Element")==0) return ElementEnum;
     
    12431244              else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum;
    12441245              else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum;
    1245               else if (strcmp(name,"Materials")==0) return MaterialsEnum;
    12461246         else stage=11;
    12471247   }
    12481248   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;
    12501251              else if (strcmp(name,"Matice")==0) return MaticeEnum;
    12511252              else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
     
    13661367              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
    13671368              else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
    1368               else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum;
    13691369         else stage=12;
    13701370   }
    13711371   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;
    13731374              else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;
    13741375              else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum;
  • issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp

    r25784 r25958  
    163163                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    164164        }
    165         else if(strcmp(string_in,"SurfaceloadIceThicknessChange")==0){
     165        else if(strcmp(string_in,"SurfaceloadIceThicknessRate")==0){
    166166                const char* field = "md.solidearth.surfaceload.icethicknesschange";
    167                 input_enum        = SurfaceloadIceThicknessChangeEnum;
     167                input_enum        = SurfaceloadIceThicknessRateEnum;
    168168                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    169169        }
     
    178178                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    179179        }
    180         else if(strcmp(string_in,"SurfaceloadWaterHeightChange")==0){
     180        else if(strcmp(string_in,"SurfaceloadWaterHeightRate")==0){
    181181                const char* field = "md.solidearth.surfaceload.waterheightchange";
    182                 input_enum        = SurfaceloadWaterHeightChangeEnum;
     182                input_enum        = SurfaceloadWaterHeightRateEnum;
    183183                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    184184        }
  • issm/trunk-jpl/src/m/classes/mmeadditionalsolidearthsolution.m

    r25956 r25958  
    6969                function marshall(self,prefix,md,fid) % {{{
    7070                        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');
    7273                        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);
    78104
    79105                end % }}}
     106
    80107                function savemodeljs(self,fid,modelname) % {{{
    81108                        error('mmeadditionalsolidearthsolution error message: not implemented yet');
  • issm/trunk-jpl/src/m/classes/mmeofflinesolidearthsolution.m

    r25956 r25958  
    7070                        WriteData(fid,prefix,'object',self,'data',4,'name','md.solidearth.external.nature','format','Integer'); %code 4 for mmeofflinesolidearthsolution  class
    7171                        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);
    78104
    79105                end % }}}
  • issm/trunk-jpl/src/m/classes/surfaceload.m

    r25956 r25958  
    6565                end % }}}
    6666                function marshall(self,prefix,md,fid) % {{{
    67 
     67                       
     68                        %deal with ice thickness change: {{{
    6869                        if isempty(self.icethicknesschange),
    6970                                self.icethicknesschange=zeros(md.mesh.numberofelements+1,1);
    7071                        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: {{{
    7197                        if isempty(self.waterheightchange),
    7298                                self.waterheightchange=zeros(md.mesh.numberofelements+1,1);
    7399                        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: {{{
    74125                        if isempty(self.other),
    75126                                self.other=zeros(md.mesh.numberofelements+1,1);
    76127                        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);
    80141                        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);
    83150                        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                        %}}}
    89152
    90153                end % }}}
Note: See TracChangeset for help on using the changeset viewer.