Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26946) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26947) @@ -2296,32 +2296,59 @@ }/*}}}*/ -void Element::LapseRateBasinSMB(IssmDouble* lapseratepos, IssmDouble* lapserateneg,IssmDouble* refelevation){/*{{{*/ +void Element::LapseRateBasinSMB(int numelevbins, IssmDouble* lapserates, IssmDouble* elevbins,IssmDouble* refelevation){/*{{{*/ /*Variable declaration*/ const int numvertices = this->GetNumberOfVertices(); - bool isadjustsmb = true; - int basinid; - IssmDouble ela,refelevation_b,lapseratepos_b,lapserateneg_b; - IssmDouble* surfacelist = xNew(numvertices); - IssmDouble* smblist = xNew(numvertices); + bool isadjustsmb = false; + int basinid,bb1,bb2; + IssmDouble ela,refelevation_b; + IssmDouble* surfacelist = xNew(numvertices); + IssmDouble* smblist = xNew(numvertices); + /* numelevbins values of lapse rates */ + IssmDouble* lapserates_b = xNew(numelevbins); + /* (numelevbins-1) limits between elevation bins (be cautious with indexing) */ + IssmDouble* elevbins_b = xNew(numelevbins-1); /*Retrieve SMB values non-adjusted for SMB lapse rate*/ this->GetInputListOnVertices(smblist,SmbMassBalanceEnum); + /*Get surface elevation on vertices*/ + this->GetInputListOnVertices(surfacelist,SurfaceEnum); /*Get basin-specific parameters of the element*/ this->GetInputValue(&basinid,SmbBasinsIdEnum); refelevation_b = refelevation[basinid]; - lapseratepos_b = lapseratepos[basinid]; - lapserateneg_b = lapserateneg[basinid]; - /*Get surface elevation on vertices*/ - this->GetInputListOnVertices(surfacelist,SurfaceEnum); - - for(int v=0;v=0 && lapseratepos_b!=0) ela = refelevation_b-smblist[v]/lapseratepos_b; - else if(smblist[v]<0 && lapserateneg_b!=0) ela = refelevation_b-smblist[v]/lapserateneg_b; - /*Lapse rate is zero: SMB not adjusted*/ - else isadjustsmb = false; - /*Adjust SMB value*/ - if(isadjustsmb) smblist[v] = lapseratepos_b*max(surfacelist[v]-ela,0.0)+lapserateneg_b*min(surfacelist[v]-ela,0.0); + for(int ii=0;ii<(numelevbins-1);ii++) elevbins_b[ii] = elevbins[basinid*(numelevbins-1)+ii]; + for(int ii=0;iielevbins_b[numelevbins-1-1]) bb1 = numelevbins-1; + if(refelevation_b>elevbins_b[numelevbins-1-1]) bb2 = numelevbins-1; + /*Vertex and Reference elevation in same elevation bin*/ + if(bb1==bb2){ + smblist[v] = smblist[v]+(surfacelist[v]-refelevation_b)*lapserates_b[bb2]; + } + /*Vertex in lower elevation bin than Reference elevation*/ + if(bb1bb1;ii--) smblist[v] = smblist[v]+(elevbins_b[ii-1]-elevbins_b[ii])*lapserates_b[ii]; + smblist[v] = smblist[v]+(surfacelist[v]-elevbins_b[bb1])*lapserates_b[bb1]; + } + /*Vertex in higher elevation bin than Reference elevation*/ + if(bb1>bb2){ + smblist[v] = smblist[v]+(elevbins_b[bb2]-refelevation_b)*lapserates_b[bb2]; + for(int ii=bb2+1;ii(lapserates_b); + xDelete(elevbins_b); xDelete(surfacelist); xDelete(smblist); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26946) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26947) @@ -154,5 +154,5 @@ bool IsLandInElement(); void Ismip6FloatingiceMeltingRate(); - void LapseRateBasinSMB(IssmDouble* lapseratepos, IssmDouble* lapserateneg,IssmDouble* refelevation); + void LapseRateBasinSMB(int numelevbins, IssmDouble* lapserates, IssmDouble* elevbins,IssmDouble* refelevation); void LinearFloatingiceMeltingRate(); void SpatialLinearFloatingiceMeltingRate(); @@ -229,4 +229,5 @@ virtual void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0; virtual void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){_error_("not implemented yet");}; + virtual void CalvingProjectionXY(void){_error_("not implemented yet");}; virtual void CalvingRateParameterization(void){_error_("not implemented yet");}; virtual void CalvingRateVonmises(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26946) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26947) @@ -428,4 +428,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_initialtime",SmbAutoregressionInitialTimeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_timestep",SmbAutoregressionTimestepEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_bins",SmbNumElevationBinsEnum)); iomodel->FetchData(&transparam,&M,&N,"md.smb.beta0"); parameters->AddObject(new DoubleVecParam(SmbBeta0Enum,transparam,N)); @@ -437,9 +438,9 @@ parameters->AddObject(new DoubleMatParam(SmbPhiEnum,transparam,M,N)); xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserate_pos"); - parameters->AddObject(new DoubleVecParam(SmbLapseRatePosEnum,transparam,N)); - xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserate_neg"); - parameters->AddObject(new DoubleVecParam(SmbLapseRateNegEnum,transparam,N)); + iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserates"); + parameters->AddObject(new DoubleMatParam(SmbLapseRatesEnum,transparam,M,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.smb.elevationbins"); + parameters->AddObject(new DoubleMatParam(SmbElevationBinsEnum,transparam,M,N)); xDelete(transparam); iomodel->FetchData(&transparam,&M,&N,"md.smb.refelevation"); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26946) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26947) @@ -198,22 +198,23 @@ bool isstochastic; bool issmbstochastic = false; - int M,N,Nphi,arorder,numbasins,my_rank; + int M,N,Nphi,arorder,numbasins,numelevbins,my_rank; femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum); femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum); + femmodel->parameters->FindParam(&numelevbins,SmbNumElevationBinsEnum); IssmDouble tinit_ar; - IssmDouble* beta0 = NULL; - IssmDouble* beta1 = NULL; - IssmDouble* phi = NULL; - IssmDouble* lapseratepos = NULL; - IssmDouble* lapserateneg = NULL; - IssmDouble* refelevation = NULL; + IssmDouble* beta0 = NULL; + IssmDouble* beta1 = NULL; + IssmDouble* phi = NULL; + IssmDouble* lapserates = NULL; + IssmDouble* elevbins = NULL; + IssmDouble* refelevation = NULL; femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum); - femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); - femmodel->parameters->FindParam(&lapseratepos,&M,SmbLapseRatePosEnum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&lapserateneg,&M,SmbLapseRateNegEnum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); + femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum); _assert_(M==numbasins); _assert_(N==numelevbins); + femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum); _assert_(M==numbasins); _assert_(N==numelevbins-1); + femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum); _assert_(M==numbasins); femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); @@ -237,5 +238,5 @@ element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum); /*Compute lapse rate adjustment*/ - element->LapseRateBasinSMB(lapseratepos,lapserateneg,refelevation); + element->LapseRateBasinSMB(numelevbins,lapserates,elevbins,refelevation); } @@ -244,6 +245,6 @@ xDelete(beta1); xDelete(phi); - xDelete(lapseratepos); - xDelete(lapserateneg); + xDelete(lapserates); + xDelete(elevbins); xDelete(refelevation); }/*}}}*/ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26946) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26947) @@ -465,4 +465,5 @@ syn keyword cConstant SmbDpermilEnum syn keyword cConstant SmbDsnowIdxEnum +syn keyword cConstant SmbElevationBinsEnum syn keyword cConstant SmbCldFracEnum syn keyword cConstant SmbDelta18oEnum @@ -492,7 +493,7 @@ syn keyword cConstant SmbIsturbulentfluxEnum syn keyword cConstant SmbKEnum -syn keyword cConstant SmbLapseRateNegEnum -syn keyword cConstant SmbLapseRatePosEnum +syn keyword cConstant SmbLapseRatesEnum syn keyword cConstant SmbNumBasinsEnum +syn keyword cConstant SmbNumElevationBinsEnum syn keyword cConstant SmbNumRequestedOutputsEnum syn keyword cConstant SmbPfacEnum @@ -1086,4 +1087,5 @@ syn keyword cConstant TransientAccumulatedDeltaIceThicknessEnum syn keyword cConstant VelEnum +syn keyword cConstant VVmismipglposinputEnum syn keyword cConstant VxAverageEnum syn keyword cConstant VxBaseEnum @@ -1623,4 +1625,5 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel +syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1629,6 +1632,6 @@ syn keyword cType ControlInput syn keyword cType Covertree +syn keyword cType DatasetInput syn keyword cType DataSetParam -syn keyword cType DatasetInput syn keyword cType Definition syn keyword cType DependentObject @@ -1643,6 +1646,6 @@ syn keyword cType ElementInput syn keyword cType ElementMatrix +syn keyword cType Elements syn keyword cType ElementVector -syn keyword cType Elements syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1651,9 +1654,10 @@ syn keyword cType Friction syn keyword cType Gauss +syn keyword cType GaussianVariogram +syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria -syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1671,4 +1675,5 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface +syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1681,4 +1686,5 @@ syn keyword cType Matice syn keyword cType Matlitho +syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1693,6 +1699,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType Options syn keyword cType OptionUtilities -syn keyword cType Options syn keyword cType Param syn keyword cType Parameters @@ -1708,11 +1714,11 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType Riftfront syn keyword cType RiftStruct -syn keyword cType Riftfront syn keyword cType SealevelGeometry syn keyword cType Seg syn keyword cType SegInput +syn keyword cType Segment syn keyword cType SegRef -syn keyword cType Segment syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1733,8 +1739,4 @@ syn keyword cType Vertex syn keyword cType Vertices -syn keyword cType classes -syn keyword cType gaussobjects -syn keyword cType krigingobjects -syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26946) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26947) @@ -459,4 +459,5 @@ SmbDpermilEnum, SmbDsnowIdxEnum, + SmbElevationBinsEnum, SmbCldFracEnum, SmbDelta18oEnum, @@ -486,7 +487,7 @@ SmbIsturbulentfluxEnum, SmbKEnum, - SmbLapseRateNegEnum, - SmbLapseRatePosEnum, + SmbLapseRatesEnum, SmbNumBasinsEnum, + SmbNumElevationBinsEnum, SmbNumRequestedOutputsEnum, SmbPfacEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26946) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26947) @@ -467,4 +467,5 @@ case SmbDpermilEnum : return "SmbDpermil"; case SmbDsnowIdxEnum : return "SmbDsnowIdx"; + case SmbElevationBinsEnum : return "SmbElevationBins"; case SmbCldFracEnum : return "SmbCldFrac"; case SmbDelta18oEnum : return "SmbDelta18o"; @@ -494,7 +495,7 @@ case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux"; case SmbKEnum : return "SmbK"; - case SmbLapseRateNegEnum : return "SmbLapseRateNeg"; - case SmbLapseRatePosEnum : return "SmbLapseRatePos"; + case SmbLapseRatesEnum : return "SmbLapseRates"; case SmbNumBasinsEnum : return "SmbNumBasins"; + case SmbNumElevationBinsEnum : return "SmbNumElevationBins"; case SmbNumRequestedOutputsEnum : return "SmbNumRequestedOutputs"; case SmbPfacEnum : return "SmbPfac"; @@ -1088,4 +1089,5 @@ case TransientAccumulatedDeltaIceThicknessEnum : return "TransientAccumulatedDeltaIceThickness"; case VelEnum : return "Vel"; + case VVmismipglposinputEnum : return "VVmismipglposinput"; case VxAverageEnum : return "VxAverage"; case VxBaseEnum : return "VxBase"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26946) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26947) @@ -458,4 +458,5 @@ syn keyword juliaConstC SmbDpermilEnum syn keyword juliaConstC SmbDsnowIdxEnum +syn keyword juliaConstC SmbElevationBinsEnum syn keyword juliaConstC SmbCldFracEnum syn keyword juliaConstC SmbDelta18oEnum @@ -485,7 +486,7 @@ syn keyword juliaConstC SmbIsturbulentfluxEnum syn keyword juliaConstC SmbKEnum -syn keyword juliaConstC SmbLapseRateNegEnum -syn keyword juliaConstC SmbLapseRatePosEnum +syn keyword juliaConstC SmbLapseRatesEnum syn keyword juliaConstC SmbNumBasinsEnum +syn keyword juliaConstC SmbNumElevationBinsEnum syn keyword juliaConstC SmbNumRequestedOutputsEnum syn keyword juliaConstC SmbPfacEnum @@ -1079,4 +1080,5 @@ syn keyword juliaConstC TransientAccumulatedDeltaIceThicknessEnum syn keyword juliaConstC VelEnum +syn keyword juliaConstC VVmismipglposinputEnum syn keyword juliaConstC VxAverageEnum syn keyword juliaConstC VxBaseEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26946) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26947) @@ -476,4 +476,5 @@ else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; else if (strcmp(name,"SmbDsnowIdx")==0) return SmbDsnowIdxEnum; + else if (strcmp(name,"SmbElevationBins")==0) return SmbElevationBinsEnum; else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum; else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum; @@ -503,11 +504,11 @@ else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum; else if (strcmp(name,"SmbK")==0) return SmbKEnum; - else if (strcmp(name,"SmbLapseRateNeg")==0) return SmbLapseRateNegEnum; - else if (strcmp(name,"SmbLapseRatePos")==0) return SmbLapseRatePosEnum; + else if (strcmp(name,"SmbLapseRates")==0) return SmbLapseRatesEnum; else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum; + if (strcmp(name,"SmbNumElevationBins")==0) return SmbNumElevationBinsEnum; + else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum; else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum; else if (strcmp(name,"SmbPhi")==0) return SmbPhiEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum; else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum; - else if (strcmp(name,"BasalforcingsLinearBasinId")==0) return BasalforcingsLinearBasinIdEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum; + if (strcmp(name,"BasalforcingsLinearBasinId")==0) return BasalforcingsLinearBasinIdEnum; + else if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum; else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum; else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"FrictionP")==0) return FrictionPEnum; else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum; - else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; + if (strcmp(name,"FrictionQ")==0) return FrictionQEnum; + else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"SealevelBarystaticIceWeights")==0) return SealevelBarystaticIceWeightsEnum; else if (strcmp(name,"SealevelBarystaticIceArea")==0) return SealevelBarystaticIceAreaEnum; - else if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum; + if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum; + else if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum; else if (strcmp(name,"SealevelBarystaticIceLoad")==0) return SealevelBarystaticIceLoadEnum; else if (strcmp(name,"SealevelBarystaticHydroMask")==0) return SealevelBarystaticHydroMaskEnum; @@ -997,9 +998,9 @@ else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum; else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum; - else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SmbHref")==0) return SmbHrefEnum; + if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum; + else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum; else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum; else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; @@ -1112,4 +1113,5 @@ else if (strcmp(name,"TransientAccumulatedDeltaIceThickness")==0) return TransientAccumulatedDeltaIceThicknessEnum; else if (strcmp(name,"Vel")==0) return VelEnum; + else if (strcmp(name,"VVmismipglposinput")==0) return VVmismipglposinputEnum; else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; else if (strcmp(name,"VxBase")==0) return VxBaseEnum; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"VxShear")==0) return VxShearEnum; else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum; - else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; - else if (strcmp(name,"VyBase")==0) return VyBaseEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"Vy")==0) return VyEnum; + if (strcmp(name,"VyAverage")==0) return VyAverageEnum; + else if (strcmp(name,"VyBase")==0) return VyBaseEnum; + else if (strcmp(name,"Vy")==0) return VyEnum; else if (strcmp(name,"VyMesh")==0) return VyMeshEnum; else if (strcmp(name,"VyObs")==0) return VyObsEnum; @@ -1242,10 +1244,10 @@ else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum; else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum; - else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum; - else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum; else stage=11; } if(stage==11){ - if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum; + if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum; + else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum; + else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum; else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum; else if (strcmp(name,"Outputdefinition100")==0) return Outputdefinition100Enum; @@ -1365,10 +1367,10 @@ else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; - else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; - else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"GenericParam")==0) return GenericParamEnum; + if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; + else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; + else if (strcmp(name,"GenericParam")==0) return GenericParamEnum; else if (strcmp(name,"GenericExternalResult")==0) return GenericExternalResultEnum; else if (strcmp(name,"Gradient1")==0) return Gradient1Enum; @@ -1488,10 +1490,10 @@ else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum; else if (strcmp(name,"Nodal")==0) return NodalEnum; - else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum; - else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum; + if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum; + else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; + else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum; else if (strcmp(name,"None")==0) return NoneEnum; else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum; @@ -1611,10 +1613,10 @@ else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; - else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; - else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; else stage=14; } if(stage==14){ - if (strcmp(name,"TransientParam")==0) return TransientParamEnum; + if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; + else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; + else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; else if (strcmp(name,"Tria")==0) return TriaEnum; Index: /issm/trunk-jpl/src/m/classes/SMBautoregression.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBautoregression.m (revision 26946) +++ /issm/trunk-jpl/src/m/classes/SMBautoregression.m (revision 26947) @@ -14,6 +14,6 @@ phi = NaN; basin_id = NaN; - lapserate_pos = NaN; - lapserate_neg = NaN; + lapserates = NaN; + elevationbins = NaN; refelevation = NaN; steps_per_step = 1; @@ -62,5 +62,5 @@ self.ar_order = 0.0; %autoregression model of order 0 end % }}} - function md = checkconsistency(self,md,solution,analyses) + function md = checkconsistency(self,md,solution,analyses) % {{{ if ismember('MasstransportAnalysis',analyses), @@ -73,20 +73,27 @@ 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 md = checkfield(md,'fieldname','smb.phi','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]); - - if (any(isnan(md.smb.lapserate_pos)==0) || numel(md.smb.lapserate_pos)>1) - md = checkfield(md,'fieldname','smb.lapserate_pos','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); + + if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1) + md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); end - if (any(isnan(md.smb.lapserate_neg)==0) || numel(md.smb.lapserate_neg)>1) - md = checkfield(md,'fieldname','smb.lapserate_neg','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); - end - if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1) - md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); + [nbas,nbins] = size(md.smb.lapserates); + if (any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins]))==0) || numel(md.smb.lapserates)>1) + md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins],'numel',md.smb.num_basins*nbins); + md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins-1],'numel',md.smb.num_basins*(nbins-1)); + if(issorted(md.smb.elevationbins,2)==0) + error('md.smb.elevationbins should have rows in order of increasing elevation'); + end + elseif (isnan(md.smb.elevationbins(1,1))==0 || numel(md.smb.elevationbins)>1) + %elevationbins specified but not lapserates: this will inevitably lead to inconsistencies + [nbas,nbins] = size(md.smb.elevationbins); + nbins = nbins+1; + md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins],'numel',md.smb.num_basins*nbins); + md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins-1],'numel',md.smb.num_basins*(nbins-1)); end - end md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]); md = checkfield(md,'fieldname','smb.averaging','numel',[1],'values',[0 1 2]); md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1); - end + end % }}} function disp(self) % {{{ disp(sprintf(' surface forcings parameters:')); @@ -99,7 +106,7 @@ fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]'); fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]'); - 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)'); - 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'); - fielddisplay(self,'refelevation','basin-specific reference elevations on which SMB lapse rates are applied (default: basin mean elevation) [m]'); + 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)'); + fielddisplay(self,'elevationbins','basin-specific separations between elevation bins, 1 row per basin, 1 column per limit between bins [m] (default: no basin separation)'); + 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]'); fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'); fielddisplay(self, 'averaging', 'averaging methods from short to long steps'); @@ -114,15 +121,13 @@ yts=md.constants.yts; - templapserate_pos = md.smb.lapserate_pos; - templapserate_neg = md.smb.lapserate_neg; + templapserates = md.smb.lapserates; + tempelevationbins = md.smb.elevationbins; temprefelevation = md.smb.refelevation; - if(any(isnan(md.smb.lapserate_pos))) - templapserate_pos = zeros(1,md.smb.num_basins); - disp(' smb.lapserate_pos not specified: set to 0'); + [nbas,nbins] = size(md.smb.lapserates); + if(any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins])))) + templapserates = zeros(md.smb.num_basins,2); + disp(' smb.lapserates not specified: set to 0'); + tempelevationbins = zeros(md.smb.num_basins,1); %dummy elevation bins end - if(any(isnan(md.smb.lapserate_neg))) - templapserate_neg = zeros(1,md.smb.num_basins); - disp(' smb.lapserate_neg not specified: set to 0'); - end if(any(isnan(md.smb.refelevation))) temprefelevation = zeros(1,md.smb.num_basins); @@ -136,8 +141,9 @@ temprefelevation(ii) = sum(areas(indices).*elemsh)/sum(areas(indices)); end - if(any([templapserate_pos,templapserate_neg])~=0) + if(any(reshape(md.smb.lapserates,[1,nbas*nbins])~=0)) disp(' smb.refelevation not specified: Reference elevations set to mean surface elevation of basins'); end end + [nbas,nbins] = size(templapserates); WriteData(fid,prefix,'name','md.smb.model','data',13,'format','Integer'); @@ -150,7 +156,8 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta1','format','DoubleMat','name','md.smb.beta1','scale',1./(yts^2),'yts',yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',yts); - WriteData(fid,prefix,'data',templapserate_pos,'format','DoubleMat','name','md.smb.lapserate_pos','scale',1./yts,'yts',yts); - WriteData(fid,prefix,'data',templapserate_neg,'format','DoubleMat','name','md.smb.lapserate_neg','scale',1./yts,'yts',yts); + WriteData(fid,prefix,'data',templapserates,'format','DoubleMat','name','md.smb.lapserates','scale',1./yts,'yts',yts); + WriteData(fid,prefix,'data',tempelevationbins,'format','DoubleMat','name','md.smb.elevationbins'); WriteData(fid,prefix,'data',temprefelevation,'format','DoubleMat','name','md.smb.refelevation'); + WriteData(fid,prefix,'data',nbins,'format','Integer','name','md.smb.num_bins'); WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer'); Index: /issm/trunk-jpl/src/m/classes/SMBautoregression.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBautoregression.py (revision 26946) +++ /issm/trunk-jpl/src/m/classes/SMBautoregression.py (revision 26947) @@ -23,6 +23,6 @@ self.phi = np.nan self.basin_id = np.nan - self.lapserate_pos = np.nan - self.lapserate_neg = np.nan + self.lapserates = np.nan + self.elevationbins = np.nan self.refelevation = np.nan self.steps_per_step = 1 @@ -47,7 +47,7 @@ s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]')) s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]')) - 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)')) - 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)')) - s += '{}\n'.format(fielddisplay(self, 'refelevation', 'basin-specific reference elevations on which SMB lapse rates are applied (default: basin mean elevation) [m]')) + 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)')) + 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)')) + 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]')) s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step')) s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) @@ -104,19 +104,27 @@ 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 md = checkfield(md, 'fieldname', 'smb.phi', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ar_order]) - - if(np.any(np.isnan(self.lapserate_pos) is False) or np.size(self.lapserate_pos) > 1): - if len(np.shape(self.lapserate_pos)) == 1: - self.lapserate_pos = np.array([self.lapserate_pos]) - md = checkfield(md, 'fieldname', 'smb.lapserate_pos', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) - - if(np.any(np.isnan(self.lapserate_neg) is False) or np.size(self.lapserate_neg) > 1): - if len(np.shape(self.lapserate_neg)) == 1: - self.lapserate_neg = np.array([self.lapserate_neg]) - md = checkfield(md, 'fieldname', 'smb.lapserate_neg', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) - + if(np.any(np.isnan(self.refelevation) is False) or np.size(self.refelevation) > 1): if len(np.shape(self.refelevation)) == 1: self.refelevation = np.array([self.refelevation]) md = checkfield(md, 'fieldname', 'smb.refelevation', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) + + nbins = np.shape(self.lapserates)[1] + if(np.any(np.isnan(self.lapserates) is False) or np.size(self.lapserates) > 1): + if len(np.shape(self.lapserates)) == 1: + self.lapserates = np.array([self.lapserates]) + if len(np.shape(self.elevationbins)) == 1: + self.elevationbins = np.array([self.elevationbins]) + md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins], 'numel', md.smb.num_basins*nbins) + md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins-1], 'numel', md.smb.num_basins*(nbins-1)) + for rr in range(md.smb.num_basins): + if(np.all(self.elevationbins[rr,0:-1]<=self.elevationbins[rr,1:])==False): + raise TypeError('md.smb.elevationbins should have rows in order of increasing elevation') + elif(np.any(np.isnan(self.elevationbins) is False) or np.size(self.elevationbins) > 1): + #elevationbins specified but not lapserates: this will inevitably lead to inconsistencies + nbins = np.shape(self.elevationbins)[1] + nbins += 1 + md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins], 'numel', md.smb.num_basins*nbins) + md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins-1], 'numel', md.smb.num_basins*(nbins-1)) md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1]) @@ -129,13 +137,11 @@ yts = md.constants.yts - templapserate_pos = np.copy(md.smb.lapserate_pos) - templapserate_neg = np.copy(md.smb.lapserate_neg) + templapserates = np.copy(md.smb.lapserates) + tempelevationbins = np.copy(md.smb.elevationbins) temprefelevation = np.copy(md.smb.refelevation) - if(np.any(np.isnan(md.smb.lapserate_pos))): - templapserate_pos = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins) - print(' smb.lapserate_pos not specified: set to 0') - if(np.any(np.isnan(md.smb.lapserate_neg))): - templapserate_neg = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins) - print(' smb.lapserate_neg not specified: set to 0') + if(np.any(np.isnan(md.smb.lapserates))): + templapserates = np.zeros((md.smb.num_basins,2)) + print(' smb.lapserates not specified: set to 0') + tempelevationbins = np.zeros((md.smb.num_basins,1)) #dummy elevation bins if(np.any(np.isnan(md.smb.refelevation))): temprefelevation = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins) @@ -147,6 +153,7 @@ elemsh[jj] = np.mean(md.geometry.surface[md.mesh.elements[indices[jj], :] - 1]) temprefelevation[0, ii] = np.sum(areas[indices] * elemsh) / np.sum(areas[indices]) - if(np.any(templapserate_pos != 0) or np.any(templapserate_neg != 0)): + if(np.any(templapserates != 0)): print(' smb.refelevation not specified: Reference elevations set to mean surface elevation of basins') + nbins = np.shape(templapserates)[1] WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 13, 'format', 'Integer') @@ -159,7 +166,8 @@ WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.smb.beta1', 'scale', 1 / (yts ** 2), 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.smb.phi', 'yts', yts) - WriteData(fid, prefix, 'data', templapserate_pos, 'name', 'md.smb.lapserate_pos', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts) - WriteData(fid, prefix, 'data', templapserate_neg, 'name', 'md.smb.lapserate_neg', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts) + WriteData(fid, prefix, 'data', templapserates, 'name', 'md.smb.lapserates', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts) + WriteData(fid, prefix, 'data', tempelevationbins, 'name', 'md.smb.elevationbins', 'format', 'DoubleMat') WriteData(fid, prefix, 'data', temprefelevation, 'name', 'md.smb.refelevation', 'format', 'DoubleMat') + WriteData(fid, prefix, 'data', nbins, 'name', 'md.smb.num_bins', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer') Index: /issm/trunk-jpl/test/NightlyRun/test257.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26946) +++ /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26947) @@ -45,6 +45,6 @@ md.smb.ar_timestep = 2.0; %timestep of the autoregressive model [yr] md.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]]; -md.smb.lapserate_pos = [0.01,0.0,0.0]; -md.smb.lapserate_neg = [0.01,0.0,0.0]; +md.smb.lapserates = [0.01,0.0;0.01,-0.01;0.0,-0.01]; +md.smb.elevationbins = [100;150;100]; %Stochastic forcing @@ -53,4 +53,5 @@ md.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 md.stochasticforcing.randomflag = 0; %fixed random seeds +md.stochasticforcing.stochastictimestep = 1.0; md=solve(md,'Transient'); Index: /issm/trunk-jpl/test/NightlyRun/test257.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.py (revision 26946) +++ /issm/trunk-jpl/test/NightlyRun/test257.py (revision 26947) @@ -50,6 +50,6 @@ md.smb.ar_timestep = 2.0 #timestep of the autoregressive model [yr] md.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]]) -md.smb.lapserate_pos = np.array([[0.01,0,0]]) -md.smb.lapserate_neg = np.array([[0.01,0,0]]) +md.smb.lapserates = np.array([[0.01,0.0],[0.01,-0.01],[0.0,-0.01]]) +md.smb.elevationbins = np.array([100,150,100]).reshape(md.smb.num_basins,1) # Stochastic forcing @@ -58,5 +58,5 @@ md.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 md.stochasticforcing.randomflag = 0 # fixed random seeds - +md.stochasticforcing.stochastictimestep = 1.0 md = solve(md, 'Transient')