Changeset 26947


Ignore:
Timestamp:
04/09/22 06:47:21 (3 years ago)
Author:
vverjans
Message:

CHG: changed use of lapse rate coefficients in SMBautoregression

Location:
issm/trunk-jpl
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r26850 r26947  
    22962296
    22972297}/*}}}*/
    2298 void       Element::LapseRateBasinSMB(IssmDouble* lapseratepos, IssmDouble* lapserateneg,IssmDouble* refelevation){/*{{{*/
     2298void       Element::LapseRateBasinSMB(int numelevbins, IssmDouble* lapserates, IssmDouble* elevbins,IssmDouble* refelevation){/*{{{*/
    22992299
    23002300        /*Variable declaration*/
    23012301   const int numvertices = this->GetNumberOfVertices();
    2302    bool isadjustsmb = true;
    2303         int basinid;
    2304    IssmDouble ela,refelevation_b,lapseratepos_b,lapserateneg_b;
    2305    IssmDouble* surfacelist = xNew<IssmDouble>(numvertices);
    2306    IssmDouble* smblist     = xNew<IssmDouble>(numvertices);
     2302   bool isadjustsmb = false;
     2303        int basinid,bb1,bb2;
     2304   IssmDouble ela,refelevation_b;
     2305   IssmDouble* surfacelist  = xNew<IssmDouble>(numvertices);
     2306   IssmDouble* smblist      = xNew<IssmDouble>(numvertices);
     2307   /* numelevbins values of lapse rates */
     2308        IssmDouble* lapserates_b = xNew<IssmDouble>(numelevbins);
     2309   /* (numelevbins-1) limits between elevation bins (be cautious with indexing) */
     2310        IssmDouble* elevbins_b   = xNew<IssmDouble>(numelevbins-1);
    23072311
    23082312   /*Retrieve SMB values non-adjusted for SMB lapse rate*/
    23092313   this->GetInputListOnVertices(smblist,SmbMassBalanceEnum);
     2314        /*Get surface elevation on vertices*/
     2315        this->GetInputListOnVertices(surfacelist,SurfaceEnum);
    23102316   /*Get basin-specific parameters of the element*/
    23112317   this->GetInputValue(&basinid,SmbBasinsIdEnum);
    23122318   refelevation_b = refelevation[basinid];
    2313    lapseratepos_b = lapseratepos[basinid];
    2314    lapserateneg_b = lapserateneg[basinid];
    2315    /*Get surface elevation on vertices*/
    2316    this->GetInputListOnVertices(surfacelist,SurfaceEnum);
    2317 
    2318    for(int v=0;v<numvertices;v++){
    2319       /*Find ELA*/
    2320                 if(smblist[v]>=0 && lapseratepos_b!=0)     ela = refelevation_b-smblist[v]/lapseratepos_b;
    2321       else if(smblist[v]<0 && lapserateneg_b!=0) ela = refelevation_b-smblist[v]/lapserateneg_b;
    2322                 /*Lapse rate is zero: SMB not adjusted*/
    2323                 else                                       isadjustsmb = false;
    2324                 /*Adjust SMB value*/
    2325       if(isadjustsmb) smblist[v] = lapseratepos_b*max(surfacelist[v]-ela,0.0)+lapserateneg_b*min(surfacelist[v]-ela,0.0);
     2319        for(int ii=0;ii<(numelevbins-1);ii++) elevbins_b[ii] = elevbins[basinid*(numelevbins-1)+ii];
     2320        for(int ii=0;ii<numelevbins;ii++){
     2321                lapserates_b[ii] = lapserates[basinid*numelevbins+ii];
     2322                if(lapserates_b[ii]!=0) isadjustsmb=true;
     2323        }
     2324        /*Adjust SMB if any lapse rate value is non-zero*/
     2325        if(isadjustsmb){
     2326       
     2327           for(int v=0;v<numvertices;v++){
     2328              /*Find elevation bin of Reference elevation and of Vertex*/
     2329                        for(int ii=0;ii<(numelevbins-1);ii++){
     2330                                if(surfacelist[v]<=elevbins_b[ii]) bb1 = ii;   
     2331                                if(refelevation_b<=elevbins_b[ii]) bb2 = ii;
     2332                        }
     2333                        /*Check for elevations above highest bin limit */
     2334                        if(surfacelist[v]>elevbins_b[numelevbins-1-1]) bb1 = numelevbins-1;
     2335                        if(refelevation_b>elevbins_b[numelevbins-1-1]) bb2 = numelevbins-1;
     2336                        /*Vertex and Reference elevation in same elevation bin*/
     2337                        if(bb1==bb2){
     2338                                smblist[v] = smblist[v]+(surfacelist[v]-refelevation_b)*lapserates_b[bb2];
     2339                        }
     2340                        /*Vertex in lower elevation bin than Reference elevation*/
     2341                        if(bb1<bb2){
     2342                                smblist[v] = smblist[v]+(elevbins_b[bb2-1]-refelevation_b)*lapserates_b[bb2];
     2343                                for(int ii=bb2-1;ii>bb1;ii--) smblist[v] = smblist[v]+(elevbins_b[ii-1]-elevbins_b[ii])*lapserates_b[ii];
     2344                                smblist[v] = smblist[v]+(surfacelist[v]-elevbins_b[bb1])*lapserates_b[bb1];
     2345                        }
     2346                        /*Vertex in higher elevation bin than Reference elevation*/
     2347                        if(bb1>bb2){
     2348                                smblist[v] = smblist[v]+(elevbins_b[bb2]-refelevation_b)*lapserates_b[bb2];
     2349                                for(int ii=bb2+1;ii<bb1;ii++) smblist[v] = smblist[v]+(elevbins_b[ii]-elevbins_b[ii-1])*lapserates_b[ii];
     2350                                smblist[v] = smblist[v]+(surfacelist[v]-elevbins_b[bb1-1])*lapserates_b[bb1];
     2351                        }
     2352                }
    23262353        }
    23272354
    … …  
    23302357
    23312358   /*Cleanup*/
     2359   xDelete<IssmDouble>(lapserates_b);
     2360   xDelete<IssmDouble>(elevbins_b);
    23322361   xDelete<IssmDouble>(surfacelist);
    23332362   xDelete<IssmDouble>(smblist);

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r26850 r26947  
    154154                bool               IsLandInElement();
    155155                void               Ismip6FloatingiceMeltingRate();
    156                 void               LapseRateBasinSMB(IssmDouble* lapseratepos, IssmDouble* lapserateneg,IssmDouble* refelevation);
     156                void               LapseRateBasinSMB(int numelevbins, IssmDouble* lapserates, IssmDouble* elevbins,IssmDouble* refelevation);
    157157                void               LinearFloatingiceMeltingRate();
    158158                void               SpatialLinearFloatingiceMeltingRate();
    … …  
    229229                virtual void       AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;
    230230                virtual void             BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){_error_("not implemented yet");};
     231                virtual void       CalvingProjectionXY(void){_error_("not implemented yet");};
    231232                virtual void       CalvingRateParameterization(void){_error_("not implemented yet");};
    232233                virtual void       CalvingRateVonmises(void){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r26859 r26947  
    428428         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_initialtime",SmbAutoregressionInitialTimeEnum));
    429429         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_timestep",SmbAutoregressionTimestepEnum));
     430         parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_bins",SmbNumElevationBinsEnum));
    430431         iomodel->FetchData(&transparam,&M,&N,"md.smb.beta0");
    431432         parameters->AddObject(new DoubleVecParam(SmbBeta0Enum,transparam,N));
    … …  
    437438         parameters->AddObject(new DoubleMatParam(SmbPhiEnum,transparam,M,N));
    438439         xDelete<IssmDouble>(transparam);
    439          iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserate_pos");
    440          parameters->AddObject(new DoubleVecParam(SmbLapseRatePosEnum,transparam,N));
    441          xDelete<IssmDouble>(transparam);
    442          iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserate_neg");
    443          parameters->AddObject(new DoubleVecParam(SmbLapseRateNegEnum,transparam,N));
     440         iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserates");
     441         parameters->AddObject(new DoubleMatParam(SmbLapseRatesEnum,transparam,M,N));
     442         xDelete<IssmDouble>(transparam);
     443         iomodel->FetchData(&transparam,&M,&N,"md.smb.elevationbins");
     444         parameters->AddObject(new DoubleMatParam(SmbElevationBinsEnum,transparam,M,N));
    444445         xDelete<IssmDouble>(transparam);
    445446                        iomodel->FetchData(&transparam,&M,&N,"md.smb.refelevation");

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

    r26810 r26947  
    198198   bool isstochastic;
    199199   bool issmbstochastic = false;
    200    int M,N,Nphi,arorder,numbasins,my_rank;
     200   int M,N,Nphi,arorder,numbasins,numelevbins,my_rank;
    201201   femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
    202202   femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
     203   femmodel->parameters->FindParam(&numelevbins,SmbNumElevationBinsEnum);
    203204   IssmDouble tinit_ar;
    204    IssmDouble* beta0        = NULL;
    205    IssmDouble* beta1        = NULL;
    206    IssmDouble* phi          = NULL;
    207    IssmDouble* lapseratepos = NULL;
    208    IssmDouble* lapserateneg = NULL;
    209    IssmDouble* refelevation = NULL;
     205   IssmDouble* beta0         = NULL;
     206   IssmDouble* beta1         = NULL;
     207   IssmDouble* phi           = NULL;
     208   IssmDouble* lapserates    = NULL;
     209   IssmDouble* elevbins      = NULL;
     210   IssmDouble* refelevation  = NULL;
    210211
    211212   femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
    212    femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);               _assert_(M==numbasins);
    213    femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);               _assert_(M==numbasins);
    214    femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);             _assert_(M==numbasins); _assert_(Nphi==arorder);
    215    femmodel->parameters->FindParam(&lapseratepos,&M,SmbLapseRatePosEnum); _assert_(M==numbasins);
    216    femmodel->parameters->FindParam(&lapserateneg,&M,SmbLapseRateNegEnum); _assert_(M==numbasins);
    217    femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum); _assert_(M==numbasins);
     213   femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);                  _assert_(M==numbasins);
     214   femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);                  _assert_(M==numbasins);
     215   femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);                _assert_(M==numbasins); _assert_(Nphi==arorder);
     216   femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum);     _assert_(M==numbasins); _assert_(N==numelevbins);
     217   femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum);    _assert_(M==numbasins); _assert_(N==numelevbins-1);
     218   femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum);    _assert_(M==numbasins);
    218219
    219220   femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
    … …  
    237238                element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum);
    238239                /*Compute lapse rate adjustment*/
    239                 element->LapseRateBasinSMB(lapseratepos,lapserateneg,refelevation);
     240                element->LapseRateBasinSMB(numelevbins,lapserates,elevbins,refelevation);
    240241        }
    241242
    … …  
    244245   xDelete<IssmDouble>(beta1);
    245246   xDelete<IssmDouble>(phi);
    246    xDelete<IssmDouble>(lapseratepos);
    247    xDelete<IssmDouble>(lapserateneg);
     247   xDelete<IssmDouble>(lapserates);
     248   xDelete<IssmDouble>(elevbins);
    248249   xDelete<IssmDouble>(refelevation);
    249250}/*}}}*/

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r26874 r26947  
    465465syn keyword cConstant SmbDpermilEnum
    466466syn keyword cConstant SmbDsnowIdxEnum
     467syn keyword cConstant SmbElevationBinsEnum
    467468syn keyword cConstant SmbCldFracEnum
    468469syn keyword cConstant SmbDelta18oEnum
    … …  
    492493syn keyword cConstant SmbIsturbulentfluxEnum
    493494syn keyword cConstant SmbKEnum
    494 syn keyword cConstant SmbLapseRateNegEnum
    495 syn keyword cConstant SmbLapseRatePosEnum
     495syn keyword cConstant SmbLapseRatesEnum
    496496syn keyword cConstant SmbNumBasinsEnum
     497syn keyword cConstant SmbNumElevationBinsEnum
    497498syn keyword cConstant SmbNumRequestedOutputsEnum
    498499syn keyword cConstant SmbPfacEnum
    … …  
    10861087syn keyword cConstant TransientAccumulatedDeltaIceThicknessEnum
    10871088syn keyword cConstant VelEnum
     1089syn keyword cConstant VVmismipglposinputEnum
    10881090syn keyword cConstant VxAverageEnum
    10891091syn keyword cConstant VxBaseEnum
    … …  
    16231625syn keyword cType Cfsurfacesquare
    16241626syn keyword cType Channel
     1627syn keyword cType classes
    16251628syn keyword cType Constraint
    16261629syn keyword cType Constraints
    … …  
    16291632syn keyword cType ControlInput
    16301633syn keyword cType Covertree
     1634syn keyword cType DatasetInput
    16311635syn keyword cType DataSetParam
    1632 syn keyword cType DatasetInput
    16331636syn keyword cType Definition
    16341637syn keyword cType DependentObject
    … …  
    16431646syn keyword cType ElementInput
    16441647syn keyword cType ElementMatrix
     1648syn keyword cType Elements
    16451649syn keyword cType ElementVector
    1646 syn keyword cType Elements
    16471650syn keyword cType ExponentialVariogram
    16481651syn keyword cType ExternalResult
    … …  
    16511654syn keyword cType Friction
    16521655syn keyword cType Gauss
     1656syn keyword cType GaussianVariogram
     1657syn keyword cType gaussobjects
    16531658syn keyword cType GaussPenta
    16541659syn keyword cType GaussSeg
    16551660syn keyword cType GaussTetra
    16561661syn keyword cType GaussTria
    1657 syn keyword cType GaussianVariogram
    16581662syn keyword cType GenericExternalResult
    16591663syn keyword cType GenericOption
    … …  
    16711675syn keyword cType IssmDirectApplicInterface
    16721676syn keyword cType IssmParallelDirectApplicInterface
     1677syn keyword cType krigingobjects
    16731678syn keyword cType Load
    16741679syn keyword cType Loads
    … …  
    16811686syn keyword cType Matice
    16821687syn keyword cType Matlitho
     1688syn keyword cType matrixobjects
    16831689syn keyword cType MatrixParam
    16841690syn keyword cType Misfit
    … …  
    16931699syn keyword cType Observations
    16941700syn keyword cType Option
     1701syn keyword cType Options
    16951702syn keyword cType OptionUtilities
    1696 syn keyword cType Options
    16971703syn keyword cType Param
    16981704syn keyword cType Parameters
    … …  
    17081714syn keyword cType Regionaloutput
    17091715syn keyword cType Results
     1716syn keyword cType Riftfront
    17101717syn keyword cType RiftStruct
    1711 syn keyword cType Riftfront
    17121718syn keyword cType SealevelGeometry
    17131719syn keyword cType Seg
    17141720syn keyword cType SegInput
     1721syn keyword cType Segment
    17151722syn keyword cType SegRef
    1716 syn keyword cType Segment
    17171723syn keyword cType SpcDynamic
    17181724syn keyword cType SpcStatic
    … …  
    17331739syn keyword cType Vertex
    17341740syn keyword cType Vertices
    1735 syn keyword cType classes
    1736 syn keyword cType gaussobjects
    1737 syn keyword cType krigingobjects
    1738 syn keyword cType matrixobjects
    17391741syn keyword cType AdjointBalancethickness2Analysis
    17401742syn keyword cType AdjointBalancethicknessAnalysis

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

    r26859 r26947  
    459459        SmbDpermilEnum,
    460460        SmbDsnowIdxEnum,
     461   SmbElevationBinsEnum,
    461462        SmbCldFracEnum,
    462463        SmbDelta18oEnum,
    … …  
    486487        SmbIsturbulentfluxEnum,
    487488        SmbKEnum,
    488         SmbLapseRateNegEnum,
    489    SmbLapseRatePosEnum,
     489   SmbLapseRatesEnum,
    490490        SmbNumBasinsEnum,
     491        SmbNumElevationBinsEnum,
    491492        SmbNumRequestedOutputsEnum,
    492493        SmbPfacEnum,

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

    r26859 r26947  
    467467                case SmbDpermilEnum : return "SmbDpermil";
    468468                case SmbDsnowIdxEnum : return "SmbDsnowIdx";
     469                case SmbElevationBinsEnum : return "SmbElevationBins";
    469470                case SmbCldFracEnum : return "SmbCldFrac";
    470471                case SmbDelta18oEnum : return "SmbDelta18o";
    … …  
    494495                case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux";
    495496                case SmbKEnum : return "SmbK";
    496                 case SmbLapseRateNegEnum : return "SmbLapseRateNeg";
    497                 case SmbLapseRatePosEnum : return "SmbLapseRatePos";
     497                case SmbLapseRatesEnum : return "SmbLapseRates";
    498498                case SmbNumBasinsEnum : return "SmbNumBasins";
     499                case SmbNumElevationBinsEnum : return "SmbNumElevationBins";
    499500                case SmbNumRequestedOutputsEnum : return "SmbNumRequestedOutputs";
    500501                case SmbPfacEnum : return "SmbPfac";
    … …  
    10881089                case TransientAccumulatedDeltaIceThicknessEnum : return "TransientAccumulatedDeltaIceThickness";
    10891090                case VelEnum : return "Vel";
     1091                case VVmismipglposinputEnum : return "VVmismipglposinput";
    10901092                case VxAverageEnum : return "VxAverage";
    10911093                case VxBaseEnum : return "VxBase";

  • issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim

    r26859 r26947  
    458458syn keyword juliaConstC SmbDpermilEnum
    459459syn keyword juliaConstC SmbDsnowIdxEnum
     460syn keyword juliaConstC SmbElevationBinsEnum
    460461syn keyword juliaConstC SmbCldFracEnum
    461462syn keyword juliaConstC SmbDelta18oEnum
    … …  
    485486syn keyword juliaConstC SmbIsturbulentfluxEnum
    486487syn keyword juliaConstC SmbKEnum
    487 syn keyword juliaConstC SmbLapseRateNegEnum
    488 syn keyword juliaConstC SmbLapseRatePosEnum
     488syn keyword juliaConstC SmbLapseRatesEnum
    489489syn keyword juliaConstC SmbNumBasinsEnum
     490syn keyword juliaConstC SmbNumElevationBinsEnum
    490491syn keyword juliaConstC SmbNumRequestedOutputsEnum
    491492syn keyword juliaConstC SmbPfacEnum
    … …  
    10791080syn keyword juliaConstC TransientAccumulatedDeltaIceThicknessEnum
    10801081syn keyword juliaConstC VelEnum
     1082syn keyword juliaConstC VVmismipglposinputEnum
    10811083syn keyword juliaConstC VxAverageEnum
    10821084syn keyword juliaConstC VxBaseEnum

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

    r26859 r26947  
    476476              else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum;
    477477              else if (strcmp(name,"SmbDsnowIdx")==0) return SmbDsnowIdxEnum;
     478              else if (strcmp(name,"SmbElevationBins")==0) return SmbElevationBinsEnum;
    478479              else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum;
    479480              else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum;
    … …  
    503504              else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum;
    504505              else if (strcmp(name,"SmbK")==0) return SmbKEnum;
    505               else if (strcmp(name,"SmbLapseRateNeg")==0) return SmbLapseRateNegEnum;
    506               else if (strcmp(name,"SmbLapseRatePos")==0) return SmbLapseRatePosEnum;
     506              else if (strcmp(name,"SmbLapseRates")==0) return SmbLapseRatesEnum;
    507507              else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
     511              if (strcmp(name,"SmbNumElevationBins")==0) return SmbNumElevationBinsEnum;
     512              else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
    512513              else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum;
    513514              else if (strcmp(name,"SmbPhi")==0) return SmbPhiEnum;
    … …  
    628629              else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum;
    629630              else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum;
    630               else if (strcmp(name,"BasalforcingsLinearBasinId")==0) return BasalforcingsLinearBasinIdEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum;
     634              if (strcmp(name,"BasalforcingsLinearBasinId")==0) return BasalforcingsLinearBasinIdEnum;
     635              else if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum;
    635636              else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum;
    636637              else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum;
    … …  
    751752              else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
    752753              else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum;
    753               else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
     757              if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
     758              else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
    758759              else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
    759760              else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
    … …  
    874875              else if (strcmp(name,"SealevelBarystaticIceWeights")==0) return SealevelBarystaticIceWeightsEnum;
    875876              else if (strcmp(name,"SealevelBarystaticIceArea")==0) return SealevelBarystaticIceAreaEnum;
    876               else if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum;
     880              if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum;
     881              else if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum;
    881882              else if (strcmp(name,"SealevelBarystaticIceLoad")==0) return SealevelBarystaticIceLoadEnum;
    882883              else if (strcmp(name,"SealevelBarystaticHydroMask")==0) return SealevelBarystaticHydroMaskEnum;
    … …  
    997998              else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum;
    998999              else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum;
    999               else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
     1003              if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
     1004              else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
    10041005              else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
    10051006              else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
    … …  
    11121113              else if (strcmp(name,"TransientAccumulatedDeltaIceThickness")==0) return TransientAccumulatedDeltaIceThicknessEnum;
    11131114              else if (strcmp(name,"Vel")==0) return VelEnum;
     1115              else if (strcmp(name,"VVmismipglposinput")==0) return VVmismipglposinputEnum;
    11141116              else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
    11151117              else if (strcmp(name,"VxBase")==0) return VxBaseEnum;
    … …  
    11191121              else if (strcmp(name,"VxShear")==0) return VxShearEnum;
    11201122              else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum;
    1121               else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
    1122               else if (strcmp(name,"VyBase")==0) return VyBaseEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"Vy")==0) return VyEnum;
     1126              if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
     1127              else if (strcmp(name,"VyBase")==0) return VyBaseEnum;
     1128              else if (strcmp(name,"Vy")==0) return VyEnum;
    11271129              else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
    11281130              else if (strcmp(name,"VyObs")==0) return VyObsEnum;
    … …  
    12421244              else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum;
    12431245              else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum;
    1244               else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
    1245               else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum;
     1249              if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
     1250              else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
     1251              else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum;
    12501252              else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum;
    12511253              else if (strcmp(name,"Outputdefinition100")==0) return Outputdefinition100Enum;
    … …  
    13651367              else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
    13661368              else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
    1367               else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
    1368               else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
     1372              if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
     1373              else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
     1374              else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
    13731375              else if (strcmp(name,"GenericExternalResult")==0) return GenericExternalResultEnum;
    13741376              else if (strcmp(name,"Gradient1")==0) return Gradient1Enum;
    … …  
    14881490              else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum;
    14891491              else if (strcmp(name,"Nodal")==0) return NodalEnum;
    1490               else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
    1491               else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    14921492         else stage=13;
    14931493   }
    14941494   if(stage==13){
    1495               if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
     1495              if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
     1496              else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
     1497              else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
    14961498              else if (strcmp(name,"None")==0) return NoneEnum;
    14971499              else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
    … …  
    16111613              else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
    16121614              else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum;
    1613               else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
    1614               else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
    16151615         else stage=14;
    16161616   }
    16171617   if(stage==14){
    1618               if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
     1618              if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
     1619              else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
     1620              else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
    16191621              else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
    16201622              else if (strcmp(name,"Tria")==0) return TriaEnum;

  • issm/trunk-jpl/src/m/classes/SMBautoregression.m

    r26832 r26947  
    1414                phi               = NaN;
    1515                basin_id          = NaN;
    16                 lapserate_pos     = NaN;
    17                 lapserate_neg     = NaN;
     16                lapserates        = NaN;
     17                elevationbins     = NaN;
    1818                refelevation      = NaN;
    1919                steps_per_step    = 1;
    … …  
    6262                        self.ar_order    = 0.0; %autoregression model of order 0
    6363                end % }}}
    64                 function md = checkconsistency(self,md,solution,analyses)
     64                function md = checkconsistency(self,md,solution,analyses) % {{{
    6565
    6666                        if ismember('MasstransportAnalysis',analyses),
    … …  
    7373                                md = checkfield(md,'fieldname','smb.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep
    7474                                md = checkfield(md,'fieldname','smb.phi','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]);
    75                        
    76                                 if (any(isnan(md.smb.lapserate_pos)==0) || numel(md.smb.lapserate_pos)>1)
    77                md = checkfield(md,'fieldname','smb.lapserate_pos','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins);
     75
     76                                if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1)
     77               md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins);
    7878            end
    79                                 if (any(isnan(md.smb.lapserate_neg)==0) || numel(md.smb.lapserate_neg)>1)
    80                md = checkfield(md,'fieldname','smb.lapserate_neg','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins);
    81             end
    82                                 if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1)
    83                                         md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins);
     79                                [nbas,nbins] = size(md.smb.lapserates);
     80                                if (any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins]))==0) || numel(md.smb.lapserates)>1)
     81                                        md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins],'numel',md.smb.num_basins*nbins);
     82                                        md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins-1],'numel',md.smb.num_basins*(nbins-1));
     83                                        if(issorted(md.smb.elevationbins,2)==0)
     84                                                error('md.smb.elevationbins should have rows in order of increasing elevation');
     85                                        end
     86                                elseif (isnan(md.smb.elevationbins(1,1))==0 || numel(md.smb.elevationbins)>1)
     87                                        %elevationbins specified but not lapserates: this will inevitably lead to inconsistencies
     88                                        [nbas,nbins] = size(md.smb.elevationbins);
     89                                        nbins        = nbins+1;
     90                                        md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins],'numel',md.smb.num_basins*nbins);
     91                                        md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins-1],'numel',md.smb.num_basins*(nbins-1));
    8492                                end
    85 
    8693                        end
    8794                        md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
    8895                        md = checkfield(md,'fieldname','smb.averaging','numel',[1],'values',[0 1 2]);
    8996                        md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
    90                 end
     97                end % }}}
    9198                function disp(self) % {{{
    9299                        disp(sprintf('   surface forcings parameters:'));
    … …  
    99106                        fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]');
    100107                        fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]');
    101                         fielddisplay(self,'lapserate_pos','basin-specific SMB lapse rates applied in range of SMB>=0 [m ice eq yr^-1 m^-1] (default: no lapse rate)');
    102                         fielddisplay(self,'lapserate_neg','basin-specific SMB lapse rates applied in range of SMB<0 [m ice eq yr^-1 m^-1] (default: no lapse rate');
    103                         fielddisplay(self,'refelevation','basin-specific reference elevations on which SMB lapse rates are applied (default: basin mean elevation) [m]');
     108                        fielddisplay(self,'lapserates','basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin [m ice eq yr^-1 m^-1] (default: no lapse rate)');
     109                        fielddisplay(self,'elevationbins','basin-specific separations between elevation bins, 1 row per basin, 1 column per limit between bins [m] (default: no basin separation)');
     110                        fielddisplay(self,'refelevation','basin-specific reference elevations at which SMB is calculated, and from which SMB is downscaled using lapserates (default: basin mean elevation) [m]');
    104111                        fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
    105112                        fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
    … …  
    114121                        yts=md.constants.yts;
    115122
    116                         templapserate_pos = md.smb.lapserate_pos;
    117                         templapserate_neg = md.smb.lapserate_neg;
     123                        templapserates    = md.smb.lapserates;
     124                        tempelevationbins = md.smb.elevationbins;
    118125                        temprefelevation  = md.smb.refelevation;
    119                         if(any(isnan(md.smb.lapserate_pos)))
    120                                 templapserate_pos = zeros(1,md.smb.num_basins);
    121                                 disp('      smb.lapserate_pos not specified: set to 0');
     126                        [nbas,nbins]      = size(md.smb.lapserates);
     127                        if(any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins]))))
     128                                templapserates = zeros(md.smb.num_basins,2);
     129                                disp('      smb.lapserates not specified: set to 0');
     130                           tempelevationbins = zeros(md.smb.num_basins,1); %dummy elevation bins
    122131                        end
    123                         if(any(isnan(md.smb.lapserate_neg)))
    124             templapserate_neg = zeros(1,md.smb.num_basins);
    125             disp('      smb.lapserate_neg not specified: set to 0');
    126          end
    127132                        if(any(isnan(md.smb.refelevation)))
    128133                                temprefelevation = zeros(1,md.smb.num_basins);
    … …  
    136141                                        temprefelevation(ii) = sum(areas(indices).*elemsh)/sum(areas(indices));
    137142                                end
    138                                 if(any([templapserate_pos,templapserate_neg])~=0)
     143                                if(any(reshape(md.smb.lapserates,[1,nbas*nbins])~=0))
    139144                                        disp('      smb.refelevation not specified: Reference elevations set to mean surface elevation of basins');
    140145                                end
    141146                        end
     147                        [nbas,nbins] = size(templapserates);
    142148
    143149                        WriteData(fid,prefix,'name','md.smb.model','data',13,'format','Integer');
    … …  
    150156                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta1','format','DoubleMat','name','md.smb.beta1','scale',1./(yts^2),'yts',yts);
    151157                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',yts);
    152                         WriteData(fid,prefix,'data',templapserate_pos,'format','DoubleMat','name','md.smb.lapserate_pos','scale',1./yts,'yts',yts);
    153                         WriteData(fid,prefix,'data',templapserate_neg,'format','DoubleMat','name','md.smb.lapserate_neg','scale',1./yts,'yts',yts);
     158                        WriteData(fid,prefix,'data',templapserates,'format','DoubleMat','name','md.smb.lapserates','scale',1./yts,'yts',yts);
     159                        WriteData(fid,prefix,'data',tempelevationbins,'format','DoubleMat','name','md.smb.elevationbins');
    154160                        WriteData(fid,prefix,'data',temprefelevation,'format','DoubleMat','name','md.smb.refelevation');
     161                        WriteData(fid,prefix,'data',nbins,'format','Integer','name','md.smb.num_bins');
    155162                        WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer');
    156163                        WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer');

  • issm/trunk-jpl/src/m/classes/SMBautoregression.py

    r26905 r26947  
    2323        self.phi = np.nan
    2424        self.basin_id = np.nan
    25         self.lapserate_pos = np.nan
    26         self.lapserate_neg = np.nan
     25        self.lapserates = np.nan
     26        self.elevationbins = np.nan
    2727        self.refelevation = np.nan
    2828        self.steps_per_step = 1
    … …  
    4747        s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]'))
    4848        s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]'))
    49         s += '{}\n'.format(fielddisplay(self, 'lapserate_pos', 'basin-specific SMB lapse rates applied in range of SMB>=0 [m ice eq yr^-1 m^-1] (default: no lapse rate)'))
    50         s += '{}\n'.format(fielddisplay(self, 'lapserate_neg', 'basin-specific SMB lapse rates applied in range of SMB<0 [m ice eq yr^-1 m^-1] (default: no lapse rate)'))
    51         s += '{}\n'.format(fielddisplay(self, 'refelevation', 'basin-specific reference elevations on which SMB lapse rates are applied (default: basin mean elevation) [m]'))
     49        s += '{}\n'.format(fielddisplay(self, 'lapserates', 'basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin [m ice eq yr^-1 m^-1] (default: no lapse rate)'))
     50        s += '{}\n'.format(fielddisplay(self, 'elevationbins', 'basin-specific SMB lapse rates applied in range of SMB<0 [m ice eq yr^-1 m^-1] (default: no lapse rate)'))
     51        s += '{}\n'.format(fielddisplay(self, 'refelevation', 'basin-specific reference elevations at which SMB is calculated, and from which SMB is downscaled using lapserates (default: basin mean elevation) [m]'))
    5252        s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
    5353        s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
    … …  
    104104            md = checkfield(md, 'fieldname', 'smb.ar_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep
    105105            md = checkfield(md, 'fieldname', 'smb.phi', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ar_order])
    106 
    107             if(np.any(np.isnan(self.lapserate_pos) is False) or np.size(self.lapserate_pos) > 1):
    108                 if len(np.shape(self.lapserate_pos)) == 1:
    109                     self.lapserate_pos = np.array([self.lapserate_pos])
    110                 md = checkfield(md, 'fieldname', 'smb.lapserate_pos', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins)
    111 
    112             if(np.any(np.isnan(self.lapserate_neg) is False) or np.size(self.lapserate_neg) > 1):
    113                 if len(np.shape(self.lapserate_neg)) == 1:
    114                     self.lapserate_neg = np.array([self.lapserate_neg])
    115                 md = checkfield(md, 'fieldname', 'smb.lapserate_neg', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins)
    116 
     106           
    117107            if(np.any(np.isnan(self.refelevation) is False) or np.size(self.refelevation) > 1):
    118108                if len(np.shape(self.refelevation)) == 1:
    119109                    self.refelevation = np.array([self.refelevation])
    120110                md = checkfield(md, 'fieldname', 'smb.refelevation', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins)
     111
     112            nbins = np.shape(self.lapserates)[1]
     113            if(np.any(np.isnan(self.lapserates) is False) or np.size(self.lapserates) > 1):
     114                if len(np.shape(self.lapserates)) == 1:
     115                    self.lapserates = np.array([self.lapserates])
     116                if len(np.shape(self.elevationbins)) == 1:
     117                    self.elevationbins = np.array([self.elevationbins])
     118                md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins], 'numel', md.smb.num_basins*nbins)
     119                md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins-1], 'numel', md.smb.num_basins*(nbins-1))
     120                for rr in range(md.smb.num_basins):
     121                    if(np.all(self.elevationbins[rr,0:-1]<=self.elevationbins[rr,1:])==False):
     122                        raise TypeError('md.smb.elevationbins should have rows in order of increasing elevation')
     123            elif(np.any(np.isnan(self.elevationbins) is False) or np.size(self.elevationbins) > 1):
     124                #elevationbins specified but not lapserates: this will inevitably lead to inconsistencies
     125                nbins = np.shape(self.elevationbins)[1]
     126                nbins += 1
     127                md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins], 'numel', md.smb.num_basins*nbins)
     128                md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins-1], 'numel', md.smb.num_basins*(nbins-1))
    121129
    122130        md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
    … …  
    129137        yts = md.constants.yts
    130138
    131         templapserate_pos = np.copy(md.smb.lapserate_pos)
    132         templapserate_neg = np.copy(md.smb.lapserate_neg)
     139        templapserates    = np.copy(md.smb.lapserates)
     140        tempelevationbins = np.copy(md.smb.elevationbins)
    133141        temprefelevation  = np.copy(md.smb.refelevation)
    134         if(np.any(np.isnan(md.smb.lapserate_pos))):
    135             templapserate_pos = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins)
    136             print('      smb.lapserate_pos not specified: set to 0')
    137         if(np.any(np.isnan(md.smb.lapserate_neg))):
    138             templapserate_neg = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins)
    139             print('      smb.lapserate_neg not specified: set to 0')
     142        if(np.any(np.isnan(md.smb.lapserates))):
     143            templapserates = np.zeros((md.smb.num_basins,2))
     144            print('      smb.lapserates not specified: set to 0')
     145            tempelevationbins = np.zeros((md.smb.num_basins,1)) #dummy elevation bins
    140146        if(np.any(np.isnan(md.smb.refelevation))):
    141147            temprefelevation = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins)
    … …  
    147153                    elemsh[jj] = np.mean(md.geometry.surface[md.mesh.elements[indices[jj], :] - 1])
    148154                temprefelevation[0, ii] = np.sum(areas[indices] * elemsh) / np.sum(areas[indices])
    149             if(np.any(templapserate_pos != 0) or np.any(templapserate_neg != 0)):
     155            if(np.any(templapserates != 0)):
    150156                print('      smb.refelevation not specified: Reference elevations set to mean surface elevation of basins')
     157        nbins = np.shape(templapserates)[1]
    151158
    152159        WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 13, 'format', 'Integer')
    … …  
    159166        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.smb.beta1', 'scale', 1 / (yts ** 2), 'yts', yts)
    160167        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.smb.phi', 'yts', yts)
    161         WriteData(fid, prefix, 'data', templapserate_pos, 'name', 'md.smb.lapserate_pos', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts)
    162         WriteData(fid, prefix, 'data', templapserate_neg, 'name', 'md.smb.lapserate_neg', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts)
     168        WriteData(fid, prefix, 'data', templapserates, 'name', 'md.smb.lapserates', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts)
     169        WriteData(fid, prefix, 'data', tempelevationbins, 'name', 'md.smb.elevationbins', 'format', 'DoubleMat')
    163170        WriteData(fid, prefix, 'data', temprefelevation, 'name', 'md.smb.refelevation', 'format', 'DoubleMat')
     171        WriteData(fid, prefix, 'data', nbins, 'name', 'md.smb.num_bins', 'format', 'Integer')
    164172        WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
    165173        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')

  • issm/trunk-jpl/test/NightlyRun/test257.m

    r26810 r26947  
    4545md.smb.ar_timestep    = 2.0; %timestep of the autoregressive model [yr]
    4646md.smb.phi            = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1];[0.4,-0.4,0.1,-0.1]];
    47 md.smb.lapserate_pos  = [0.01,0.0,0.0];
    48 md.smb.lapserate_neg  = [0.01,0.0,0.0];
     47md.smb.lapserates     = [0.01,0.0;0.01,-0.01;0.0,-0.01];
     48md.smb.elevationbins  = [100;150;100];
    4949
    5050%Stochastic forcing
    … …  
    5353md.stochasticforcing.covariance          = [[0.15 0.08 -0.02];[0.08 0.12 -0.05];[-0.02 -0.05 0.1]]; %global covariance among- and between-fields
    5454md.stochasticforcing.randomflag          = 0; %fixed random seeds
     55md.stochasticforcing.stochastictimestep  = 1.0;
    5556
    5657md=solve(md,'Transient');

  • issm/trunk-jpl/test/NightlyRun/test257.py

    r26907 r26947  
    5050md.smb.ar_timestep = 2.0  #timestep of the autoregressive model [yr]
    5151md.smb.phi = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1], [0.4, -0.4, 0.1, -0.1]])
    52 md.smb.lapserate_pos = np.array([[0.01,0,0]])
    53 md.smb.lapserate_neg = np.array([[0.01,0,0]])
     52md.smb.lapserates        = np.array([[0.01,0.0],[0.01,-0.01],[0.0,-0.01]])
     53md.smb.elevationbins  = np.array([100,150,100]).reshape(md.smb.num_basins,1)
    5454
    5555# Stochastic forcing
    … …  
    5858md.stochasticforcing.covariance = np.array([[0.15, 0.08, -0.02], [0.08, 0.12, -0.05], [-0.02, -0.05, 0.1]])  # global covariance among- and between-fields
    5959md.stochasticforcing.randomflag = 0  # fixed random seeds
    60 
     60md.stochasticforcing.stochastictimestep  = 1.0
    6161
    6262md = solve(md, 'Transient')
Note: See TracChangeset for help on using the changeset viewer.