Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26603) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26604) @@ -12,37 +12,8 @@ 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"); - - /*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); + IoModelToDynamicConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,finiteelement); + //IoModelToConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,finiteelement); } /*}}}*/ @@ -83,9 +54,15 @@ /*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 26603) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26604) @@ -99,22 +99,75 @@ xDelete(phi_basin); }/*}}}*/ -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; +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; IssmDouble* phi_basin = xNew(arorder); IssmDouble* varlist = xNew(numvertices); IssmDouble* valuesautoregression = NULL; - - /*Get Basin ID and Basin coefficients*/ - if(enum_type==SMBautoregressionEnum) this->GetInputValue(&basinid,SmbBasinsIdEnum); - if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum); + 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*/ 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*/ @@ -131,24 +184,22 @@ /*Stochastic variable value*/ - varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin; + varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noiseterm; } - /*Update autoregression TF values*/ + /*Update autoregression values*/ IssmDouble* temparray = xNew(numvertices*arorder); /*Assign newest values and shift older values*/ for(int i=0;iinputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); - if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); + this->inputs->SetArrayInput(arenum_type,this->lid,temparray,numvertices*arorder); xDelete(temparray); } /*Add input to element*/ - if(enum_type==SMBautoregressionEnum) this->AddInput(SmbMassBalanceEnum,varlist,P1Enum); - if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->AddInput(FrontalForcingsThermalForcingEnum,varlist,P1Enum); + this->AddInput(outenum_type,varlist,P1Enum); /*Cleanup*/ xDelete(phi_basin); xDelete(varlist); - xDelete(valuesautoregression); + xDelete(valuesautoregression); }/*}}}*/ void Element::ComputeLambdaS(){/*{{{*/ @@ -3778,5 +3829,4 @@ IssmDouble teValue=1.0; IssmDouble aValue=0.0; - IssmDouble dulwrfValue=0.0; IssmDouble szaValue=0.0; IssmDouble cotValue=0.0; @@ -3785,5 +3835,4 @@ IssmDouble dt,time,smb_dt; int aIdx=0; - int eIdx=0; int denIdx=0; int dsnowIdx=0; @@ -3815,9 +3864,7 @@ 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:{{{ */ @@ -3869,5 +3916,4 @@ 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); @@ -3886,9 +3932,7 @@ 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: {{{*/ @@ -4080,5 +4124,4 @@ 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); @@ -4096,5 +4139,4 @@ 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] @@ -4122,5 +4164,5 @@ } /*Thermal profile computation:*/ - 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); + 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); /*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 26603) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26604) @@ -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,IssmDouble* noiseterms,int enum_type); + void Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,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 26603) +++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26604) @@ -80,5 +80,6 @@ /*Load parameters*/ bool isstochastic; - int M,N,Nphi,arorder,numbasins,my_rank; + bool istfstochastic = false; + int M,N,Nphi,arorder,numbasins,my_rank; femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum); @@ -87,5 +88,4 @@ IssmDouble* beta1 = NULL; IssmDouble* phi = NULL; - IssmDouble* noiseterms = xNewZeroInit(numbasins); femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum); @@ -94,5 +94,4 @@ 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){ @@ -102,8 +101,7 @@ femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields); for(int i=0;iparameters->FindParam(&noiseterms,&M,ThermalforcingAutoregressionNoiseEnum); _assert_(M==numbasins); - } + if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum) istfstochastic = true; } + xDelete(stochasticfields); } /*Time elapsed with respect to AR model initial time*/ @@ -113,5 +111,5 @@ for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,FrontalForcingsRignotAutoregressionEnum); + element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,istfstochastic,FrontalForcingsRignotAutoregressionEnum); } @@ -120,4 +118,3 @@ 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 26603) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26604) @@ -502,7 +502,8 @@ parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum); if(isstochasticforcing){ - int num_fields; + int num_fields,stochastic_dim; 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"); @@ -515,6 +516,5 @@ 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 26603) +++ /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp (revision 26604) @@ -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,15 +46,66 @@ int i=0; for(int j=0;j(dimensions[j]); for(int k=0;kparameters->SetParam(noisefield,dimensions[j],enum_type); - i=i+dimensions[j]; + + int dimensionid; + + /*Deal with the autoregressive models*/ + if(fields[j]==SMBautoregressionEnum || fields[j]==FrontalForcingsRignotAutoregressionEnum){ + switch(fields[j]){ + case SMBautoregressionEnum: + dimenum_type = SmbBasinsIdEnum; + noiseenum_type = SmbAutoregressionNoiseEnum; + break; + case FrontalForcingsRignotAutoregressionEnum: + dimenum_type = FrontalForcingsBasinIdEnum; + noiseenum_type = ThermalforcingAutoregressionNoiseEnum; + break; + } + for(Object* &object:femmodel->elements->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 "<(noisefield); } Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26603) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26604) @@ -197,4 +197,5 @@ /*Load parameters*/ bool isstochastic; + bool issmbstochastic = false; int M,N,Nphi,arorder,numbasins,my_rank; femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum); @@ -210,6 +211,4 @@ 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){ @@ -219,7 +218,5 @@ femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields); for(int i=0;iparameters->FindParam(&noiseterms,&M,SmbAutoregressionNoiseEnum); _assert_(M==numbasins); - } + if(stochasticfields[i]==SMBautoregressionEnum) issmbstochastic = true; } xDelete(stochasticfields); @@ -231,5 +228,5 @@ for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,SMBautoregressionEnum); + element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum); } @@ -238,5 +235,4 @@ 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 26603) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26604) @@ -398,4 +398,5 @@ syn keyword cConstant SolidearthSettingsOceanAreaScalingEnum syn keyword cConstant StochasticForcingCovarianceEnum +syn keyword cConstant StochasticForcingDefaultDimensionEnum syn keyword cConstant StochasticForcingDimensionsEnum syn keyword cConstant StochasticForcingFieldsEnum @@ -428,5 +429,4 @@ syn keyword cConstant SmbAdThreshEnum syn keyword cConstant SmbAutoregressionInitialTimeEnum -syn keyword cConstant SmbAutoregressionNoiseEnum syn keyword cConstant SmbAutoregressionTimestepEnum syn keyword cConstant SmbAutoregressiveOrderEnum @@ -504,5 +504,4 @@ syn keyword cConstant StressbalanceRiftPenaltyThresholdEnum syn keyword cConstant StressbalanceShelfDampeningEnum -syn keyword cConstant ThermalforcingAutoregressionNoiseEnum syn keyword cConstant ThermalIsdrainicecolumnEnum syn keyword cConstant ThermalIsdynamicbasalspcEnum @@ -599,4 +598,5 @@ syn keyword cConstant BasalStressEnum syn keyword cConstant BaseEnum +syn keyword cConstant BaselineCalvingCalvingrateEnum syn keyword cConstant BaseOldEnum syn keyword cConstant BaseSlopeXEnum @@ -613,4 +613,5 @@ syn keyword cConstant BottomPressureOldEnum syn keyword cConstant CalvingCalvingrateEnum +syn keyword cConstant CalvingCalvingrateNoiseEnum syn keyword cConstant CalvingHabFractionEnum syn keyword cConstant CalvingMeltingrateEnum @@ -895,4 +896,5 @@ syn keyword cConstant SmbAdiffiniEnum syn keyword cConstant SmbAiniEnum +syn keyword cConstant SmbAutoregressionNoiseEnum syn keyword cConstant SmbBasinsIdEnum syn keyword cConstant SmbBMaxEnum @@ -997,4 +999,5 @@ syn keyword cConstant SolidearthExternalDisplacementUpRateEnum syn keyword cConstant SolidearthExternalGeoidRateEnum +syn keyword cConstant StochasticForcingDefaultIdEnum syn keyword cConstant StrainRateeffectiveEnum syn keyword cConstant StrainRateparallelEnum @@ -1032,4 +1035,5 @@ 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 26603) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26604) @@ -392,4 +392,5 @@ SolidearthSettingsOceanAreaScalingEnum, StochasticForcingCovarianceEnum, + StochasticForcingDefaultDimensionEnum, StochasticForcingDimensionsEnum, StochasticForcingFieldsEnum, @@ -422,5 +423,4 @@ SmbAdThreshEnum, SmbAutoregressionInitialTimeEnum, - SmbAutoregressionNoiseEnum, SmbAutoregressionTimestepEnum, SmbAutoregressiveOrderEnum, @@ -498,5 +498,4 @@ StressbalanceRiftPenaltyThresholdEnum, StressbalanceShelfDampeningEnum, - ThermalforcingAutoregressionNoiseEnum, ThermalIsdrainicecolumnEnum, ThermalIsdynamicbasalspcEnum, @@ -595,4 +594,5 @@ BasalStressEnum, BaseEnum, + BaselineCalvingCalvingrateEnum, BaseOldEnum, BaseSlopeXEnum, @@ -609,4 +609,5 @@ BottomPressureOldEnum, CalvingCalvingrateEnum, + CalvingCalvingrateNoiseEnum, CalvingHabFractionEnum, CalvingMeltingrateEnum, @@ -891,4 +892,5 @@ SmbAdiffiniEnum, SmbAiniEnum, + SmbAutoregressionNoiseEnum, SmbBasinsIdEnum, SmbBMaxEnum, @@ -994,4 +996,5 @@ SolidearthExternalDisplacementUpRateEnum, SolidearthExternalGeoidRateEnum, + StochasticForcingDefaultIdEnum, StrainRateeffectiveEnum, StrainRateparallelEnum, @@ -1029,4 +1032,5 @@ 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 26603) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26604) @@ -400,4 +400,5 @@ case SolidearthSettingsOceanAreaScalingEnum : return "SolidearthSettingsOceanAreaScaling"; case StochasticForcingCovarianceEnum : return "StochasticForcingCovariance"; + case StochasticForcingDefaultDimensionEnum : return "StochasticForcingDefaultDimension"; case StochasticForcingDimensionsEnum : return "StochasticForcingDimensions"; case StochasticForcingFieldsEnum : return "StochasticForcingFields"; @@ -430,5 +431,4 @@ case SmbAdThreshEnum : return "SmbAdThresh"; case SmbAutoregressionInitialTimeEnum : return "SmbAutoregressionInitialTime"; - case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise"; case SmbAutoregressionTimestepEnum : return "SmbAutoregressionTimestep"; case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder"; @@ -506,5 +506,4 @@ case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold"; case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening"; - case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise"; case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn"; case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc"; @@ -601,4 +600,5 @@ case BasalStressEnum : return "BasalStress"; case BaseEnum : return "Base"; + case BaselineCalvingCalvingrateEnum : return "BaselineCalvingCalvingrate"; case BaseOldEnum : return "BaseOld"; case BaseSlopeXEnum : return "BaseSlopeX"; @@ -615,4 +615,5 @@ case BottomPressureOldEnum : return "BottomPressureOld"; case CalvingCalvingrateEnum : return "CalvingCalvingrate"; + case CalvingCalvingrateNoiseEnum : return "CalvingCalvingrateNoise"; case CalvingHabFractionEnum : return "CalvingHabFraction"; case CalvingMeltingrateEnum : return "CalvingMeltingrate"; @@ -897,4 +898,5 @@ case SmbAdiffiniEnum : return "SmbAdiffini"; case SmbAiniEnum : return "SmbAini"; + case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise"; case SmbBasinsIdEnum : return "SmbBasinsId"; case SmbBMaxEnum : return "SmbBMax"; @@ -999,4 +1001,5 @@ case SolidearthExternalDisplacementUpRateEnum : return "SolidearthExternalDisplacementUpRate"; case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate"; + case StochasticForcingDefaultIdEnum : return "StochasticForcingDefaultId"; case StrainRateeffectiveEnum : return "StrainRateeffective"; case StrainRateparallelEnum : return "StrainRateparallel"; @@ -1034,4 +1037,5 @@ 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 26603) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26604) @@ -409,4 +409,5 @@ 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; @@ -439,5 +440,4 @@ 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,5 +518,4 @@ 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; @@ -613,4 +612,5 @@ 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,10 +627,11 @@ 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,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; + if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum; + else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum; else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; @@ -751,9 +752,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,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum; + if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum; + else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum; else if (strcmp(name,"HydrologyTws")==0) return HydrologyTwsEnum; else if (strcmp(name,"HydrologyTwsSpc")==0) return HydrologyTwsSpcEnum; @@ -874,9 +875,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,"SealevelchangeGU")==0) return SealevelchangeGUEnum; + if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum; + else if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum; else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum; else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum; @@ -918,4 +919,5 @@ 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; @@ -996,10 +998,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,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum; + if (strcmp(name,"SmbTa")==0) return SmbTaEnum; + else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum; + else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum; else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum; else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum; @@ -1023,4 +1025,5 @@ 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; @@ -1058,4 +1061,5 @@ 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; @@ -1117,12 +1121,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 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 stage=10; - } - if(stage==10){ - if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum; + else if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum; else if (strcmp(name,"Outputdefinition30")==0) return Outputdefinition30Enum; else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum; @@ -1240,12 +1244,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 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 stage=11; - } - if(stage==11){ - if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum; + else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum; else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum; else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum; @@ -1363,12 +1367,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 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 stage=12; - } - if(stage==12){ - if (strcmp(name,"J")==0) return JEnum; + else if (strcmp(name,"J")==0) return JEnum; else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum; else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum; @@ -1486,12 +1490,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 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 stage=13; - } - if(stage==13){ - if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; + else 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 26603) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26604) @@ -8,5 +8,6 @@ isstochasticforcing = 0; fields = NaN; - dimensions = NaN; + defaultdimension = 0; + default_id = NaN; covariance = NaN; randomflag = 1; @@ -34,25 +35,59 @@ num_fields = numel(self.fields); - size_tot = sum(self.dimensions); + + %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 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]); - - %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 + 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 end % }}} function disp(self) % {{{ @@ -60,5 +95,6 @@ fielddisplay(self,'isstochasticforcing','is stochasticity activated?'); fielddisplay(self,'fields','fields with stochasticity applied, ex: {''SMBautoregression''}, or {''FrontalForcingsRignotAutoregression''}'); - fielddisplay(self,'dimensions','dimensionality of each field'); + 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,'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)'); @@ -77,10 +113,25 @@ 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 = {'SMBautoregression'}; %list of fields that need scaling *1/yts + scaledfields = {'DefaultCalving','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(self.dimensions(1:i-1)):1:sum(self.dimensions(1:i))]; + inds = [1+sum(dimensions(1:i-1)):1:sum(dimensions(1:i))]; for row=inds %scale rows corresponding to scaled field tempcovariance(row,:) = 1./yts*tempcovariance(row,:); @@ -93,5 +144,7 @@ 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,'object',self,'fieldname','dimensions','format','IntMat'); + 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,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat'); WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean'); @@ -105,6 +158,7 @@ list = {... - 'SMBautoregression',... - 'FrontalForcingsRignotAutoregression' + 'DefaultCalving',... + 'FrontalForcingsRignotAutoregression',... + 'SMBautoregression' }; end % }}} Index: /issm/trunk-jpl/test/NightlyRun/test257.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26603) +++ /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26604) @@ -49,5 +49,4 @@ md.stochasticforcing.isstochasticforcing = 1; md.stochasticforcing.fields = [{'SMBautoregression'}]; -md.stochasticforcing.dimensions = [md.smb.num_basins]; %dimension of each field 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 Index: /issm/trunk-jpl/test/NightlyRun/test543.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.m (revision 26603) +++ /issm/trunk-jpl/test/NightlyRun/test543.m (revision 26604) @@ -39,5 +39,4 @@ md.stochasticforcing.isstochasticforcing = 1; md.stochasticforcing.fields = [{'FrontalForcingsRignotAutoregression'}]; -md.stochasticforcing.dimensions = [md.frontalforcings.num_basins]; %dimension of each field md.stochasticforcing.covariance = 1e-4*[[1.5,0.5];[0.5,0.4]]; %global covariance among- and between-fields md.stochasticforcing.randomflag = 0; %determines true/false randomness