Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26607) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26608) @@ -12,8 +12,37 @@ void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + + /*intermediary: */ int finiteelement; + int code,vector_layout; + IssmDouble *spcdata = NULL; + int M,N; + + /*Get finite element type for this analysis*/ iomodel->FindConstant(&finiteelement,"md.levelset.fe"); - IoModelToDynamicConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,finiteelement); - //IoModelToConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,finiteelement); + + /*First of, find the record for the enum, and get code of data type: */ + iomodel->SetFilePointerToData(&code, &vector_layout,"md.levelset.spclevelset"); + if(code!=7)_error_("expecting a IssmDouble vector for constraints md.levelset.spclevelset"); + if(vector_layout!=1)_error_("expecting a nodal vector for constraints md.levelset.spclevelset"); + + /*Fetch vector:*/ + iomodel->FetchData(&spcdata,&M,&N,"md.levelset.spclevelset"); + + /*Call IoModelToConstraintsx*/ + if(N>1){ + /*If it is a time series, most likely we are forcing the ice front position and do not want to have a Dynamic Constraint*/ + _assert_(M==iomodel->numberofvertices+1); + IoModelToConstraintsx(constraints,iomodel,spcdata,M,N,LevelsetAnalysisEnum,finiteelement); + } + else{ + /*This is not a time series, we probably have calving on, we need the levelset constraints to update as the levelset moves*/ + _assert_(N==1); + _assert_(M==iomodel->numberofvertices); + IoModelToDynamicConstraintsx(constraints,iomodel,spcdata,M,N,LevelsetAnalysisEnum,finiteelement); + } + + /*Clean up*/ + xDelete(spcdata); } /*}}}*/ @@ -54,15 +83,9 @@ /*Get moving front parameters*/ - bool isstochastic; int calvinglaw; iomodel->FindConstant(&calvinglaw,"md.calving.law"); - iomodel->FindConstant(&isstochastic,"md.stochasticforcing.isstochasticforcing"); switch(calvinglaw){ case DefaultCalvingEnum: iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",CalvingCalvingrateEnum); - if(isstochastic){ - iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum); - iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",BaselineCalvingCalvingrateEnum); - } break; case CalvingLevermannEnum: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26607) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26608) @@ -99,75 +99,22 @@ xDelete(phi_basin); }/*}}}*/ -void Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type){/*{{{*/ - - const int numvertices = this->GetNumberOfVertices(); - int basinid,M,N,arenum_type,basinenum_type,noiseenum_type,outenum_type; - IssmDouble beta0_basin,beta1_basin,noiseterm; +void Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type){/*{{{*/ + + const int numvertices = this->GetNumberOfVertices(); + int basinid,M,N; + IssmDouble beta0_basin,beta1_basin,noise_basin; IssmDouble* phi_basin = xNew(arorder); IssmDouble* varlist = xNew(numvertices); IssmDouble* valuesautoregression = NULL; - Input* noiseterm_input = NULL; - - /*Get field-specific enums*/ - switch(enum_type){ - case(SMBautoregressionEnum): - arenum_type = SmbValuesAutoregressionEnum; - basinenum_type = SmbBasinsIdEnum; - noiseenum_type = SmbAutoregressionNoiseEnum; - outenum_type = SmbMassBalanceEnum; - break; - case(FrontalForcingsRignotAutoregressionEnum): - arenum_type = ThermalforcingValuesAutoregressionEnum; - basinenum_type = FrontalForcingsBasinIdEnum; - noiseenum_type = ThermalforcingAutoregressionNoiseEnum; - outenum_type = FrontalForcingsThermalForcingEnum; - break; - } - - /*Get noise and autoregressive terms*/ - this->GetInputValue(&basinid,basinenum_type); - if(isfieldstochastic){ - noiseterm_input = this->GetInput(noiseenum_type); - Gauss* gauss = this->NewGauss(); - noiseterm_input->GetInputValue(&noiseterm,gauss); - delete gauss; - } - else noiseterm = 0.0; - this->inputs->GetArray(arenum_type,this->lid,&valuesautoregression,&M); - - - /* - this->GetInputValue(&basinid,SmbBasinsIdEnum); - if(isfieldstochastic){ - noiseterm_input = this->GetInput(SmbAutoregressionNoiseEnum); - Gauss* gauss = this->NewGauss(); - noiseterm_input->GetInputValue(&noiseterm,gauss); - delete gauss; - } - else noiseterm = 0.0; - this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&valuesautoregression,&M); - arenum_type = SmbValuesAutoregressionEnum; - outenum_type = SmbMassBalanceEnum; - break; - case(FrontalForcingsRignotAutoregressionEnum): - this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum); - if(isfieldstochastic){ - noiseterm_input = this->GetInput(ThermalforcingAutoregressionNoiseEnum); - Gauss* gauss = this->NewGauss(); - noiseterm_input->GetInputValue(&noiseterm,gauss); - delete gauss; - } - else noiseterm = 0.0; - this->inputs->GetArray(ThermalforcingValuesAutoregressionEnum,this->lid,&valuesautoregression,&M); - arenum_type = ThermalforcingValuesAutoregressionEnum; - outenum_type = FrontalForcingsThermalForcingEnum; - break; - } - */ - - /*Get basin coefficients*/ + + /*Get Basin ID and Basin coefficients*/ + if(enum_type==SMBautoregressionEnum) this->GetInputValue(&basinid,SmbBasinsIdEnum); + if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum); for(int ii=0;iiinputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&valuesautoregression,&M); + if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->GetArray(ThermalforcingValuesAutoregressionEnum,this->lid,&valuesautoregression,&M); /*If not AR model timestep: take the old values of variable*/ @@ -184,22 +131,24 @@ /*Stochastic variable value*/ - varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noiseterm; + varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin; } - /*Update autoregression values*/ + /*Update autoregression TF values*/ IssmDouble* temparray = xNew(numvertices*arorder); /*Assign newest values and shift older values*/ for(int i=0;iinputs->SetArrayInput(arenum_type,this->lid,temparray,numvertices*arorder); + if(enum_type==SMBautoregressionEnum) this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); + if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); xDelete(temparray); } /*Add input to element*/ - this->AddInput(outenum_type,varlist,P1Enum); + if(enum_type==SMBautoregressionEnum) this->AddInput(SmbMassBalanceEnum,varlist,P1Enum); + if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->AddInput(FrontalForcingsThermalForcingEnum,varlist,P1Enum); /*Cleanup*/ xDelete(phi_basin); xDelete(varlist); - xDelete(valuesautoregression); + xDelete(valuesautoregression); }/*}}}*/ void Element::ComputeLambdaS(){/*{{{*/ @@ -3829,4 +3778,5 @@ IssmDouble teValue=1.0; IssmDouble aValue=0.0; + IssmDouble dulwrfValue=0.0; IssmDouble szaValue=0.0; IssmDouble cotValue=0.0; @@ -3835,4 +3785,5 @@ IssmDouble dt,time,smb_dt; int aIdx=0; + int eIdx=0; int denIdx=0; int dsnowIdx=0; @@ -3864,7 +3815,9 @@ bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux; bool isconstrainsurfaceT=false; + bool isdeltaLWup=false; IssmDouble init_scaling=0.0; IssmDouble thermo_scaling=1.0; IssmDouble adThresh=1023.0; + IssmDouble teThresh=10; /*}}}*/ /*Output variables:{{{ */ @@ -3916,4 +3869,5 @@ parameters->FindParam(&smb_dt,SmbDtEnum); /*time period for the smb solution, usually smaller than the glaciological dt*/ parameters->FindParam(&aIdx,SmbAIdxEnum); + parameters->FindParam(&eIdx,SmbEIdxEnum); parameters->FindParam(&denIdx,SmbDenIdxEnum); parameters->FindParam(&swIdx,SmbSwIdxEnum); @@ -3932,7 +3886,9 @@ parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum); parameters->FindParam(&isconstrainsurfaceT,SmbIsconstrainsurfaceTEnum); + parameters->FindParam(&isdeltaLWup,SmbIsdeltaLWupEnum); parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum); parameters->FindParam(&thermo_scaling,SmbThermoDeltaTScalingEnum); parameters->FindParam(&adThresh,SmbAdThreshEnum); + parameters->FindParam(&teThresh,SmbTeThreshEnum); /*}}}*/ /*Retrieve inputs: {{{*/ @@ -4124,4 +4080,5 @@ Input *teValue_input= this->GetInput(SmbTeValueEnum,timeinputs); _assert_(teValue_input); Input *aValue_input= this->GetInput(SmbAValueEnum,timeinputs); _assert_(aValue_input); + Input *dulwrfValue_input= this->GetInput(SmbDulwrfValueEnum,timeinputs); _assert_(dulwrfValue_input); Input *szaValue_input= this->GetInput(SmbSzaValueEnum,timeinputs); _assert_(szaValue_input); Input *cotValue_input= this->GetInput(SmbCotValueEnum,timeinputs); _assert_(cotValue_input); @@ -4139,4 +4096,5 @@ pAir_input->GetInputValue(&pAir,gauss); // screen level air pressure [Pa] teValue_input->GetInputValue(&teValue,gauss); // Emissivity [0-1] + dulwrfValue_input->GetInputValue(&dulwrfValue,gauss); // LWup perturbation [W m-2] aValue_input->GetInputValue(&aValue,gauss); // Albedo [0 1] szaValue_input->GetInputValue(&szaValue,gauss); // Solar Zenith Angle [degree] @@ -4164,5 +4122,5 @@ } /*Thermal profile computation:*/ - if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid(),isconstrainsurfaceT); + if(isthermal)thermo(&EC, &T, &ulw, re, dz, d, swf, dlw, Ta, V, eAir, pAir, eIdx, teValue, dulwrfValue, teThresh, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid(),isconstrainsurfaceT,isdeltaLWup); /*Change in thickness of top cell due to evaporation/condensation assuming same density as top cell. Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26607) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26608) @@ -69,5 +69,5 @@ bool AnyFSet(void); void AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,int enum_type); - void Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type); + void Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type); void ComputeLambdaS(void); void ComputeNewDamage(); Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26607) +++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26608) @@ -80,6 +80,5 @@ /*Load parameters*/ bool isstochastic; - bool istfstochastic = false; - int M,N,Nphi,arorder,numbasins,my_rank; + int M,N,Nphi,arorder,numbasins,my_rank; femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum); @@ -88,4 +87,5 @@ IssmDouble* beta1 = NULL; IssmDouble* phi = NULL; + IssmDouble* noiseterms = xNewZeroInit(numbasins); femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum); @@ -94,4 +94,5 @@ femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); + /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/ femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); if(isstochastic){ @@ -101,7 +102,8 @@ femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields); for(int i=0;iparameters->FindParam(&noiseterms,&M,ThermalforcingAutoregressionNoiseEnum); _assert_(M==numbasins); + } } - xDelete(stochasticfields); } /*Time elapsed with respect to AR model initial time*/ @@ -111,5 +113,5 @@ for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,istfstochastic,FrontalForcingsRignotAutoregressionEnum); + element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,FrontalForcingsRignotAutoregressionEnum); } @@ -118,3 +120,4 @@ xDelete(beta1); xDelete(phi); + xDelete(noiseterms); }/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26607) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26608) @@ -502,8 +502,7 @@ parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum); if(isstochasticforcing){ - int num_fields,stochastic_dim; + int num_fields; char** fields; parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.defaultdimension",StochasticForcingDefaultDimensionEnum)); iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields"); if(num_fields<1) _error_("no stochasticforcing fields found"); @@ -516,5 +515,6 @@ parameters->AddObject(new IntVecParam(StochasticForcingFieldsEnum,stochasticforcing_enums,num_fields)); xDelete(stochasticforcing_enums); - parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.randomflag",StochasticForcingRandomflagEnum)); + + parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.randomflag",StochasticForcingRandomflagEnum)); iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.dimensions"); parameters->AddObject(new IntVecParam(StochasticForcingDimensionsEnum,transparam,N)); Index: /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp (revision 26607) +++ /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp (revision 26608) @@ -36,5 +36,5 @@ if(randomflag) fixedseed=-1; else fixedseed = reCast((time-starttime)/dt); - /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/ + /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/ IssmDouble* temparray = NULL; multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed); @@ -46,66 +46,15 @@ int i=0; for(int j=0;j(dimensions[j]); for(int k=0;kelements->objects){ - Element* element = xDynamicCast(object); - int numvertices = element->GetNumberOfVertices(); - IssmDouble* noise_element = xNew(numvertices); - element->GetInputValue(&dimensionid,dimenum_type); - for(int i=0;iAddInput(noiseenum_type,noise_element,P0Enum); - xDelete(noise_element); - } - } - else{ - switch(fields[j]){ - case SMBautoregressionEnum: - case FrontalForcingsRignotAutoregressionEnum: - /*Already done above*/ - break; - case DefaultCalvingEnum: - /*Delete CalvingCalvingrateEnum at previous time step (required if it is transient)*/ - femmodel->inputs->DeleteInput(CalvingCalvingrateEnum); - for(Object* &object:femmodel->elements->objects){ - Element* element = xDynamicCast(object); - int numvertices = element->GetNumberOfVertices(); - IssmDouble baselinecalvingrate; - IssmDouble calvingrate_tot[numvertices]; - Input* baselinecalvingrate_input = NULL; - baselinecalvingrate_input = element->GetInput(BaselineCalvingCalvingrateEnum); _assert_(baselinecalvingrate_input); - element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum); - Gauss* gauss = element->NewGauss(); - for(int i=0;iGaussVertex(i); - baselinecalvingrate_input->GetInputValue(&baselinecalvingrate,gauss); - calvingrate_tot[i] = max(0.0,baselinecalvingrate+noisefield[dimensionid]); - } - element->AddInput(CalvingCalvingrateEnum,&calvingrate_tot[0],P1DGEnum); - delete gauss; - } - break; - default: - _error_("Field "<parameters->SetParam(noisefield,dimensions[j],enum_type); + i=i+dimensions[j]; xDelete(noisefield); } Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26607) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26608) @@ -197,5 +197,4 @@ /*Load parameters*/ bool isstochastic; - bool issmbstochastic = false; int M,N,Nphi,arorder,numbasins,my_rank; femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum); @@ -211,4 +210,6 @@ femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); + /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/ + IssmDouble* noiseterms = xNewZeroInit(numbasins); femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); if(isstochastic){ @@ -218,5 +219,7 @@ femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields); for(int i=0;iparameters->FindParam(&noiseterms,&M,SmbAutoregressionNoiseEnum); _assert_(M==numbasins); + } } xDelete(stochasticfields); @@ -228,5 +231,5 @@ for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum); + element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,SMBautoregressionEnum); } @@ -235,4 +238,5 @@ xDelete(beta1); xDelete(phi); + xDelete(noiseterms); }/*}}}*/ void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26607) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26608) @@ -398,5 +398,4 @@ syn keyword cConstant SolidearthSettingsOceanAreaScalingEnum syn keyword cConstant StochasticForcingCovarianceEnum -syn keyword cConstant StochasticForcingDefaultDimensionEnum syn keyword cConstant StochasticForcingDimensionsEnum syn keyword cConstant StochasticForcingFieldsEnum @@ -429,4 +428,5 @@ syn keyword cConstant SmbAdThreshEnum syn keyword cConstant SmbAutoregressionInitialTimeEnum +syn keyword cConstant SmbAutoregressionNoiseEnum syn keyword cConstant SmbAutoregressionTimestepEnum syn keyword cConstant SmbAutoregressiveOrderEnum @@ -504,4 +504,5 @@ syn keyword cConstant StressbalanceRiftPenaltyThresholdEnum syn keyword cConstant StressbalanceShelfDampeningEnum +syn keyword cConstant ThermalforcingAutoregressionNoiseEnum syn keyword cConstant ThermalIsdrainicecolumnEnum syn keyword cConstant ThermalIsdynamicbasalspcEnum @@ -598,5 +599,4 @@ syn keyword cConstant BasalStressEnum syn keyword cConstant BaseEnum -syn keyword cConstant BaselineCalvingCalvingrateEnum syn keyword cConstant BaseOldEnum syn keyword cConstant BaseSlopeXEnum @@ -613,5 +613,4 @@ syn keyword cConstant BottomPressureOldEnum syn keyword cConstant CalvingCalvingrateEnum -syn keyword cConstant CalvingCalvingrateNoiseEnum syn keyword cConstant CalvingHabFractionEnum syn keyword cConstant CalvingMeltingrateEnum @@ -896,5 +895,4 @@ syn keyword cConstant SmbAdiffiniEnum syn keyword cConstant SmbAiniEnum -syn keyword cConstant SmbAutoregressionNoiseEnum syn keyword cConstant SmbBasinsIdEnum syn keyword cConstant SmbBMaxEnum @@ -999,5 +997,4 @@ syn keyword cConstant SolidearthExternalDisplacementUpRateEnum syn keyword cConstant SolidearthExternalGeoidRateEnum -syn keyword cConstant StochasticForcingDefaultIdEnum syn keyword cConstant StrainRateeffectiveEnum syn keyword cConstant StrainRateparallelEnum @@ -1035,5 +1032,4 @@ syn keyword cConstant TemperaturePicardEnum syn keyword cConstant TemperatureSEMICEnum -syn keyword cConstant ThermalforcingAutoregressionNoiseEnum syn keyword cConstant ThermalforcingValuesAutoregressionEnum syn keyword cConstant ThermalSpctemperatureEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26607) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26608) @@ -392,5 +392,4 @@ SolidearthSettingsOceanAreaScalingEnum, StochasticForcingCovarianceEnum, - StochasticForcingDefaultDimensionEnum, StochasticForcingDimensionsEnum, StochasticForcingFieldsEnum, @@ -423,4 +422,5 @@ SmbAdThreshEnum, SmbAutoregressionInitialTimeEnum, + SmbAutoregressionNoiseEnum, SmbAutoregressionTimestepEnum, SmbAutoregressiveOrderEnum, @@ -498,4 +498,5 @@ StressbalanceRiftPenaltyThresholdEnum, StressbalanceShelfDampeningEnum, + ThermalforcingAutoregressionNoiseEnum, ThermalIsdrainicecolumnEnum, ThermalIsdynamicbasalspcEnum, @@ -594,5 +595,4 @@ BasalStressEnum, BaseEnum, - BaselineCalvingCalvingrateEnum, BaseOldEnum, BaseSlopeXEnum, @@ -609,5 +609,4 @@ BottomPressureOldEnum, CalvingCalvingrateEnum, - CalvingCalvingrateNoiseEnum, CalvingHabFractionEnum, CalvingMeltingrateEnum, @@ -892,5 +891,4 @@ SmbAdiffiniEnum, SmbAiniEnum, - SmbAutoregressionNoiseEnum, SmbBasinsIdEnum, SmbBMaxEnum, @@ -996,5 +994,4 @@ SolidearthExternalDisplacementUpRateEnum, SolidearthExternalGeoidRateEnum, - StochasticForcingDefaultIdEnum, StrainRateeffectiveEnum, StrainRateparallelEnum, @@ -1032,5 +1029,4 @@ TemperaturePicardEnum, TemperatureSEMICEnum, - ThermalforcingAutoregressionNoiseEnum, ThermalforcingValuesAutoregressionEnum, ThermalSpctemperatureEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26607) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26608) @@ -400,5 +400,4 @@ case SolidearthSettingsOceanAreaScalingEnum : return "SolidearthSettingsOceanAreaScaling"; case StochasticForcingCovarianceEnum : return "StochasticForcingCovariance"; - case StochasticForcingDefaultDimensionEnum : return "StochasticForcingDefaultDimension"; case StochasticForcingDimensionsEnum : return "StochasticForcingDimensions"; case StochasticForcingFieldsEnum : return "StochasticForcingFields"; @@ -431,4 +430,5 @@ case SmbAdThreshEnum : return "SmbAdThresh"; case SmbAutoregressionInitialTimeEnum : return "SmbAutoregressionInitialTime"; + case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise"; case SmbAutoregressionTimestepEnum : return "SmbAutoregressionTimestep"; case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder"; @@ -506,4 +506,5 @@ case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold"; case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening"; + case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise"; case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn"; case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc"; @@ -600,5 +601,4 @@ case BasalStressEnum : return "BasalStress"; case BaseEnum : return "Base"; - case BaselineCalvingCalvingrateEnum : return "BaselineCalvingCalvingrate"; case BaseOldEnum : return "BaseOld"; case BaseSlopeXEnum : return "BaseSlopeX"; @@ -615,5 +615,4 @@ case BottomPressureOldEnum : return "BottomPressureOld"; case CalvingCalvingrateEnum : return "CalvingCalvingrate"; - case CalvingCalvingrateNoiseEnum : return "CalvingCalvingrateNoise"; case CalvingHabFractionEnum : return "CalvingHabFraction"; case CalvingMeltingrateEnum : return "CalvingMeltingrate"; @@ -898,5 +897,4 @@ case SmbAdiffiniEnum : return "SmbAdiffini"; case SmbAiniEnum : return "SmbAini"; - case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise"; case SmbBasinsIdEnum : return "SmbBasinsId"; case SmbBMaxEnum : return "SmbBMax"; @@ -1001,5 +999,4 @@ case SolidearthExternalDisplacementUpRateEnum : return "SolidearthExternalDisplacementUpRate"; case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate"; - case StochasticForcingDefaultIdEnum : return "StochasticForcingDefaultId"; case StrainRateeffectiveEnum : return "StrainRateeffective"; case StrainRateparallelEnum : return "StrainRateparallel"; @@ -1037,5 +1034,4 @@ case TemperaturePicardEnum : return "TemperaturePicard"; case TemperatureSEMICEnum : return "TemperatureSEMIC"; - case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise"; case ThermalforcingValuesAutoregressionEnum : return "ThermalforcingValuesAutoregression"; case ThermalSpctemperatureEnum : return "ThermalSpctemperature"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26607) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26608) @@ -409,5 +409,4 @@ else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum; else if (strcmp(name,"StochasticForcingCovariance")==0) return StochasticForcingCovarianceEnum; - else if (strcmp(name,"StochasticForcingDefaultDimension")==0) return StochasticForcingDefaultDimensionEnum; else if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum; else if (strcmp(name,"StochasticForcingFields")==0) return StochasticForcingFieldsEnum; @@ -440,4 +439,5 @@ else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum; else if (strcmp(name,"SmbAutoregressionInitialTime")==0) return SmbAutoregressionInitialTimeEnum; + else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum; else if (strcmp(name,"SmbAutoregressionTimestep")==0) return SmbAutoregressionTimestepEnum; else if (strcmp(name,"SmbAutoregressiveOrder")==0) return SmbAutoregressiveOrderEnum; @@ -518,4 +518,5 @@ else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum; else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum; + else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum; else if (strcmp(name,"ThermalIsdrainicecolumn")==0) return ThermalIsdrainicecolumnEnum; else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum; @@ -612,5 +613,4 @@ else if (strcmp(name,"BasalStress")==0) return BasalStressEnum; else if (strcmp(name,"Base")==0) return BaseEnum; - else if (strcmp(name,"BaselineCalvingCalvingrate")==0) return BaselineCalvingCalvingrateEnum; else if (strcmp(name,"BaseOld")==0) return BaseOldEnum; else if (strcmp(name,"BaseSlopeX")==0) return BaseSlopeXEnum; @@ -627,11 +627,10 @@ else if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum; else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum; - else if (strcmp(name,"CalvingCalvingrateNoise")==0) return CalvingCalvingrateNoiseEnum; else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum; + else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum; - else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; + if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum; else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; @@ -752,9 +751,9 @@ else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum; else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum; + else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum; - else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum; + if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum; else if (strcmp(name,"HydrologyTws")==0) return HydrologyTwsEnum; else if (strcmp(name,"HydrologyTwsSpc")==0) return HydrologyTwsSpcEnum; @@ -875,9 +874,9 @@ else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum; else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum; + else if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum; - else if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum; + if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum; else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum; else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum; @@ -919,5 +918,4 @@ else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum; else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum; - else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum; else if (strcmp(name,"SmbBasinsId")==0) return SmbBasinsIdEnum; else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum; @@ -998,10 +996,10 @@ else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum; else if (strcmp(name,"SmbT")==0) return SmbTEnum; + else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; + else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SmbTa")==0) return SmbTaEnum; - else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum; - else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum; + if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum; else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum; else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum; @@ -1025,5 +1023,4 @@ else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum; else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum; - else if (strcmp(name,"StochasticForcingDefaultId")==0) return StochasticForcingDefaultIdEnum; else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum; else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; @@ -1061,5 +1058,4 @@ else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; else if (strcmp(name,"TemperatureSEMIC")==0) return TemperatureSEMICEnum; - else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum; else if (strcmp(name,"ThermalforcingValuesAutoregression")==0) return ThermalforcingValuesAutoregressionEnum; else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum; @@ -1121,12 +1117,12 @@ else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum; else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum; + else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; + else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; + else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum; + else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum; else stage=10; } if(stage==10){ - if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; - else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; - else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum; - else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum; - else if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum; + if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum; else if (strcmp(name,"Outputdefinition30")==0) return Outputdefinition30Enum; else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum; @@ -1244,12 +1240,12 @@ else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum; else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum; + else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum; + else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum; + else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum; + else if (strcmp(name,"Channel")==0) return ChannelEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum; - else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum; - else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum; - else if (strcmp(name,"Channel")==0) return ChannelEnum; - else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum; + if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum; else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum; else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum; @@ -1367,12 +1363,12 @@ else if (strcmp(name,"IntParam")==0) return IntParamEnum; else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum; + else if (strcmp(name,"Inputs")==0) return InputsEnum; + else if (strcmp(name,"Internal")==0) return InternalEnum; + else if (strcmp(name,"Intersect")==0) return IntersectEnum; + else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"Inputs")==0) return InputsEnum; - else if (strcmp(name,"Internal")==0) return InternalEnum; - else if (strcmp(name,"Intersect")==0) return IntersectEnum; - else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum; - else if (strcmp(name,"J")==0) return JEnum; + if (strcmp(name,"J")==0) return JEnum; else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum; else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum; @@ -1490,12 +1486,12 @@ else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; + else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; + else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; + else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; + else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; - else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; - else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; - else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; - else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; + if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.m =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26607) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26608) @@ -8,6 +8,5 @@ isstochasticforcing = 0; fields = NaN; - defaultdimension = 0; - default_id = NaN; + dimensions = NaN; covariance = NaN; randomflag = 1; @@ -35,59 +34,25 @@ num_fields = numel(self.fields); - - %Check that covariance matrix is positive definite - try - chol(self.covariance); - catch - error('md.stochasticforcing.covariance is not positive definite'); - end - - %Check that all fields agree with the corresponding md class and if any field needs the default params - checkdefaults = false; %need to check defaults only if one of the field does not have its own dimensionality - for field=self.fields - %Checking agreement of classes - if(contains(field,'SMB')) - if~(isequal(class(md.smb),char(field))) - error('md.smb does not agree with stochasticforcing field %s', char(field)); - end - end - if(contains(field,'frontalforcings')) - if~(isequal(class(md.frontalforcings),char(field))) - error('md.frontalforcings does not agree with stochasticforcing field %s', char(field)); - end - end - %Checking for specific dimensions - if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression')) - checkdefaults = true; %field with non-specific dimensionality - end - end - - %Retrieve sum of all the field dimensionalities - size_tot = self.defaultdimension*num_fields; - indSMBar = -1; %about to check for index of SMBautoregression - indTFar = -1; %about to check for index of FrontalForcingsRignotAutoregression - if any(contains(self.fields,'SMBautoregression')) - size_tot = size_tot-self.defaultdimension+md.smb.num_basins; - indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with TFar timestep (08Nov2021) - end - if any(contains(self.fields,'FrontalForcingsRignotAutoregression')) - size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins; - indTFar = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with SMBar timestep (08Nov2021) - end - - if(indSMBar~=-1 && indTFar~=-1) %both autoregressive models are used: check autoregressive time step consistency - if((md.smb.ar_timestep~=md.frontalforcings.ar_timestep) && any(self.covariance(1+sum(self.dimensions(1:indSMBar-1)):sum(self.dimensions(1:indSMBar)),1+sum(self.dimensions(1:indTFar-1)):sum(self.dimensions(1:indTFar))))) - error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance'); - end - end + size_tot = sum(self.dimensions); md = checkfield(md,'fieldname','stochasticforcing.isstochasticforcing','values',[0 1]); md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings()); %VV check here 'cell' (19Oct2021) + md = checkfield(md,'fieldname','stochasticforcing.dimensions','NaN',1,'Inf',1,'>',0,'size',[1,num_fields]); %specific dimension for each field md = checkfield(md,'fieldname','stochasticforcing.covariance','NaN',1,'Inf',1,'size',[size_tot,size_tot]); %global covariance matrix md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]); - if(checkdefaults) %need to check the defaults - md = checkfield(md,'fieldname','stochasticforcing.defaultdimension','numel',1,'NaN',1,'Inf',1,'>',0); - md = checkfield(md,'fieldname','stochasticforcing.default_id','Inf',1,'>=',0,'<=',self.defaultdimension,'size',[md.mesh.numberofelements,1]); - end + + %Check that all fields agree with the corresponding md class + for field=self.fields + if(contains(field,'SMB')) + if~(isequal(class(md.smb),char(field))) + error('md.smb does not agree with stochasticforcing field %s', char(field)); + end + end + if(contains(field,'frontalforcings')) + if~(isequal(class(md.frontalforcings),char(field))) + error('md.frontalforcings does not agree with stochasticforcing field %s', char(field)); + end + end + end end % }}} function disp(self) % {{{ @@ -95,6 +60,5 @@ fielddisplay(self,'isstochasticforcing','is stochasticity activated?'); fielddisplay(self,'fields','fields with stochasticity applied, ex: {''SMBautoregression''}, or {''FrontalForcingsRignotAutoregression''}'); - fielddisplay(self,'defaultdimension','dimensionality of the noise terms (does not apply to fields with their specific dimension)'); - fielddisplay(self,'default_id','id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)'); + fielddisplay(self,'dimensions','dimensionality of each field'); fielddisplay(self,'covariance','covariance matrix for within- and between-fields covariance (units must be squared field units)'); fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'); @@ -113,25 +77,10 @@ return else - - %Retrieve dimensionality of each field - dimensions = self.defaultdimension*ones(1,num_fields); - ind = 1; - for field=self.fields - %Checking for specific dimensions - if(strcmp(field,'SMBautoregression')) - dimensions(ind) = md.smb.num_basins; - end - if(strcmp(field,'FrontalForcingsRignotAutoregression')) - dimensions(ind) = md.frontalforcings.num_basins; - end - ind = ind+1; - end - %Scaling covariance matrix (scale column-by-column and row-by-row) - scaledfields = {'DefaultCalving','SMBautoregression'}; %list of fields that need scaling *1/yts + scaledfields = {'SMBautoregression'}; %list of fields that need scaling *1/yts tempcovariance = self.covariance; %copy of covariance to avoid writing back in member variable for i=1:num_fields if any(strcmp(scaledfields,self.fields(i))) - inds = [1+sum(dimensions(1:i-1)):1:sum(dimensions(1:i))]; + inds = [1+sum(self.dimensions(1:i-1)):1:sum(self.dimensions(1:i))]; for row=inds %scale rows corresponding to scaled field tempcovariance(row,:) = 1./yts*tempcovariance(row,:); @@ -144,7 +93,5 @@ WriteData(fid,prefix,'data',num_fields,'name','md.stochasticforcing.num_fields','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','fields','format','StringArray'); - WriteData(fid,prefix,'data',dimensions,'name','md.stochasticforcing.dimensions','format','IntMat'); - WriteData(fid,prefix,'object',self,'fieldname','default_id','data',self.default_id-1,'format','IntMat','mattype',2); %0-indexed - WriteData(fid,prefix,'object',self,'fieldname','defaultdimension','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','dimensions','format','IntMat'); WriteData(fid,prefix,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat'); WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean'); @@ -158,7 +105,6 @@ list = {... - 'DefaultCalving',... - 'FrontalForcingsRignotAutoregression',... - 'SMBautoregression' + 'SMBautoregression',... + 'FrontalForcingsRignotAutoregression' }; end % }}}