Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 27779) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 27780) @@ -575,9 +575,10 @@ parameters->AddObject(new BoolParam(StochasticForcingIsWaterPressureEnum,false)); if(isstochasticforcing){ - int num_fields,stochastic_dim; + int num_fields,num_tcov,stochastic_dim; char** fields; parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.defaultdimension",StochasticForcingDefaultDimensionEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.stochastictimestep",StochasticForcingTimestepEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_timescovariance",StochasticForcingNumTimesCovarianceEnum)); iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields"); if(num_fields<1) _error_("no stochasticforcing fields found"); @@ -594,4 +595,7 @@ parameters->AddObject(new IntVecParam(StochasticForcingDimensionsEnum,transparam,N)); xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.timecovariance"); + parameters->AddObject(new DoubleVecParam(StochasticForcingTimeCovarianceEnum,transparam,N)); + xDelete(transparam); iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.covariance"); parameters->AddObject(new DoubleMatParam(StochasticForcingCovarianceEnum,transparam,M,N)); Index: /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp (revision 27779) +++ /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp (revision 27780) @@ -14,15 +14,18 @@ /*Retrieve parameters*/ bool randomflag; - int M,N,numfields,my_rank; - int* fields = NULL; - int* dimensions = NULL; - IssmDouble* covariance = NULL; + int M,N,numfields,numtcov,my_rank; + int* fields = NULL; + int* dimensions = NULL; + IssmDouble* timecovariance = NULL; + IssmDouble* covariance = NULL; femmodel->parameters->FindParam(&randomflag,StochasticForcingRandomflagEnum); femmodel->parameters->FindParam(&numfields,StochasticForcingNumFieldsEnum); + femmodel->parameters->FindParam(&numtcov,StochasticForcingNumTimesCovarianceEnum); femmodel->parameters->FindParam(&fields,&N,StochasticForcingFieldsEnum); _assert_(N==numfields); femmodel->parameters->FindParam(&dimensions,&N,StochasticForcingDimensionsEnum); _assert_(N==numfields); + femmodel->parameters->FindParam(&timecovariance,&N,StochasticForcingTimeCovarianceEnum); _assert_(N==numtcov); int dimtot=0; for(int i=0;iparameters->FindParam(&covariance,&M,&N,StochasticForcingCovarianceEnum); _assert_(M==dimtot); _assert_(N==dimtot); + femmodel->parameters->FindParam(&covariance,&M,&N,StochasticForcingCovarianceEnum); _assert_(M==numtcov); _assert_(N==dimtot*dimtot); /*Check if this is a timestep for new noiseterms computation*/ @@ -45,6 +48,16 @@ /*Compute noise terms*/ - IssmDouble* noiseterms = xNew(dimtot); + IssmDouble* timestepcovariance = xNew(dimtot*dimtot); + IssmDouble* noiseterms = xNew(dimtot); if(isstepforstoch){ + /*Find covariance to be applied at current time step*/ + int itime; + if(numtcov>1){ + for(int i=0;i=timecovariance[i]) itime=i; + } + } + else itime=0; + for(int i=0;i(temparray); @@ -233,4 +246,6 @@ xDelete(dimensions); xDelete(covariance); + xDelete(timecovariance); + xDelete(timestepcovariance); xDelete(noiseterms); }/*}}}*/ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27779) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27780) @@ -507,5 +507,7 @@ syn keyword cConstant StochasticForcingNoisetermsEnum syn keyword cConstant StochasticForcingNumFieldsEnum +syn keyword cConstant StochasticForcingNumTimesCovarianceEnum syn keyword cConstant StochasticForcingRandomflagEnum +syn keyword cConstant StochasticForcingTimeCovarianceEnum syn keyword cConstant StochasticForcingTimestepEnum syn keyword cConstant SolidearthSettingsReltolEnum @@ -1775,4 +1777,5 @@ syn keyword cType Cfsurfacesquaretransient syn keyword cType Channel +syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1782,6 +1785,6 @@ syn keyword cType ControlParam syn keyword cType Covertree +syn keyword cType DatasetInput syn keyword cType DataSetParam -syn keyword cType DatasetInput syn keyword cType Definition syn keyword cType DependentObject @@ -1796,6 +1799,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 @@ -1804,9 +1807,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 @@ -1825,4 +1829,5 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface +syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1835,4 +1840,5 @@ syn keyword cType Matice syn keyword cType Matlitho +syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1847,6 +1853,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 @@ -1862,11 +1868,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 @@ -1887,8 +1893,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 27779) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27780) @@ -501,5 +501,7 @@ StochasticForcingNoisetermsEnum, StochasticForcingNumFieldsEnum, + StochasticForcingNumTimesCovarianceEnum, StochasticForcingRandomflagEnum, + StochasticForcingTimeCovarianceEnum, StochasticForcingTimestepEnum, SolidearthSettingsReltolEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27779) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27780) @@ -509,5 +509,7 @@ case StochasticForcingNoisetermsEnum : return "StochasticForcingNoiseterms"; case StochasticForcingNumFieldsEnum : return "StochasticForcingNumFields"; + case StochasticForcingNumTimesCovarianceEnum : return "StochasticForcingNumTimesCovariance"; case StochasticForcingRandomflagEnum : return "StochasticForcingRandomflag"; + case StochasticForcingTimeCovarianceEnum : return "StochasticForcingTimeCovariance"; case StochasticForcingTimestepEnum : return "StochasticForcingTimestep"; case SolidearthSettingsReltolEnum : return "SolidearthSettingsReltol"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27779) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27780) @@ -500,5 +500,7 @@ syn keyword juliaConstC StochasticForcingNoisetermsEnum syn keyword juliaConstC StochasticForcingNumFieldsEnum +syn keyword juliaConstC StochasticForcingNumTimesCovarianceEnum syn keyword juliaConstC StochasticForcingRandomflagEnum +syn keyword juliaConstC StochasticForcingTimeCovarianceEnum syn keyword juliaConstC StochasticForcingTimestepEnum syn keyword juliaConstC SolidearthSettingsReltolEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27779) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27780) @@ -521,5 +521,7 @@ else if (strcmp(name,"StochasticForcingNoiseterms")==0) return StochasticForcingNoisetermsEnum; else if (strcmp(name,"StochasticForcingNumFields")==0) return StochasticForcingNumFieldsEnum; + else if (strcmp(name,"StochasticForcingNumTimesCovariance")==0) return StochasticForcingNumTimesCovarianceEnum; else if (strcmp(name,"StochasticForcingRandomflag")==0) return StochasticForcingRandomflagEnum; + else if (strcmp(name,"StochasticForcingTimeCovariance")==0) return StochasticForcingTimeCovarianceEnum; else if (strcmp(name,"StochasticForcingTimestep")==0) return StochasticForcingTimestepEnum; else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum; else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum; - else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum; - else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; + if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum; + else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum; + else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; else if (strcmp(name,"Step")==0) return StepEnum; else if (strcmp(name,"Steps")==0) return StepsEnum; @@ -750,10 +752,10 @@ else if (strcmp(name,"BasalforcingsPicoOverturningCoeff")==0) return BasalforcingsPicoOverturningCoeffEnum; else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum; - else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum; - else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"BasalStressx")==0) return BasalStressxEnum; + if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum; + else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum; + else if (strcmp(name,"BasalStressx")==0) return BasalStressxEnum; else if (strcmp(name,"BasalStressy")==0) return BasalStressyEnum; else if (strcmp(name,"BasalStress")==0) return BasalStressEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"UGia")==0) return UGiaEnum; else if (strcmp(name,"UGiaRate")==0) return UGiaRateEnum; - else if (strcmp(name,"Gradient")==0) return GradientEnum; - else if (strcmp(name,"GroundinglineHeight")==0) return GroundinglineHeightEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"HydraulicPotential")==0) return HydraulicPotentialEnum; + if (strcmp(name,"Gradient")==0) return GradientEnum; + else if (strcmp(name,"GroundinglineHeight")==0) return GroundinglineHeightEnum; + else if (strcmp(name,"HydraulicPotential")==0) return HydraulicPotentialEnum; else if (strcmp(name,"HydraulicPotentialOld")==0) return HydraulicPotentialOldEnum; else if (strcmp(name,"HydrologyBasalFlux")==0) return HydrologyBasalFluxEnum; @@ -996,10 +998,10 @@ else if (strcmp(name,"SealevelBarystaticHydroLatbar")==0) return SealevelBarystaticHydroLatbarEnum; else if (strcmp(name,"SealevelBarystaticHydroLongbar")==0) return SealevelBarystaticHydroLongbarEnum; - else if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum; - else if (strcmp(name,"SealevelBarystaticBpMask")==0) return SealevelBarystaticBpMaskEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SealevelBarystaticBpWeights")==0) return SealevelBarystaticBpWeightsEnum; + if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum; + else if (strcmp(name,"SealevelBarystaticBpMask")==0) return SealevelBarystaticBpMaskEnum; + else if (strcmp(name,"SealevelBarystaticBpWeights")==0) return SealevelBarystaticBpWeightsEnum; else if (strcmp(name,"SealevelBarystaticBpArea")==0) return SealevelBarystaticBpAreaEnum; else if (strcmp(name,"SealevelBarystaticBpLoad")==0) return SealevelBarystaticBpLoadEnum; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum; else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; - else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; - else if (strcmp(name,"SmbMassBalanceSnow")==0) return SmbMassBalanceSnowEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"SmbMassBalanceIce")==0) return SmbMassBalanceIceEnum; + if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; + else if (strcmp(name,"SmbMassBalanceSnow")==0) return SmbMassBalanceSnowEnum; + else if (strcmp(name,"SmbMassBalanceIce")==0) return SmbMassBalanceIceEnum; else if (strcmp(name,"SmbMassBalanceSemic")==0) return SmbMassBalanceSemicEnum; else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum; @@ -1242,10 +1244,10 @@ else if (strcmp(name,"ThicknessResidual")==0) return ThicknessResidualEnum; else if (strcmp(name,"TransientAccumulatedDeltaIceThickness")==0) return TransientAccumulatedDeltaIceThicknessEnum; - else if (strcmp(name,"Vel")==0) return VelEnum; - else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"VxBase")==0) return VxBaseEnum; + if (strcmp(name,"Vel")==0) return VelEnum; + else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; + else if (strcmp(name,"VxBase")==0) return VxBaseEnum; else if (strcmp(name,"VxDebris")==0) return VxDebrisEnum; else if (strcmp(name,"Vx")==0) return VxEnum; @@ -1365,10 +1367,10 @@ else if (strcmp(name,"Outputdefinition82")==0) return Outputdefinition82Enum; else if (strcmp(name,"Outputdefinition83")==0) return Outputdefinition83Enum; - else if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum; - else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum; else stage=12; } if(stage==12){ - if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum; + if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum; + else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum; + else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum; else if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum; else if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum; @@ -1488,10 +1490,10 @@ else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; - else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum; - else if (strcmp(name,"FSSolver")==0) return FSSolverEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"FSpressure")==0) return FSpressureEnum; + if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum; + else if (strcmp(name,"FSSolver")==0) return FSSolverEnum; + else if (strcmp(name,"FSpressure")==0) return FSpressureEnum; else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; else if (strcmp(name,"FemModel")==0) return FemModelEnum; @@ -1611,10 +1613,10 @@ else if (strcmp(name,"Matestar")==0) return MatestarEnum; else if (strcmp(name,"Matice")==0) return MaticeEnum; - else if (strcmp(name,"Matlitho")==0) return MatlithoEnum; - else if (strcmp(name,"Mathydro")==0) return MathydroEnum; else stage=14; } if(stage==14){ - if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; + if (strcmp(name,"Matlitho")==0) return MatlithoEnum; + else if (strcmp(name,"Mathydro")==0) return MathydroEnum; + else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum; else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum; @@ -1734,10 +1736,10 @@ else if (strcmp(name,"Sset")==0) return SsetEnum; else if (strcmp(name,"StatisticsSolution")==0) return StatisticsSolutionEnum; - else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum; - else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum; else stage=15; } if(stage==15){ - if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum; + if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum; + else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum; + else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum; else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum; else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum; Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.m =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 27779) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 27780) @@ -11,4 +11,5 @@ default_id = NaN; covariance = NaN; + timecovariance = NaN; stochastictimestep = 0; randomflag = 1; @@ -40,9 +41,26 @@ end - %Check that covariance matrix is positive definite - try - chol(self.covariance); - catch - error('md.stochasticforcing.covariance is not positive definite'); + + if(numel(size(self.covariance)==3)) + numtcovmat = numel(self.covariance(1,1,:)); %number of covariance matrices in time + lsCovmats = {}; + for ii=[1:numtcovmat] %loop over 3rd dimension + lsCovmats{ii} = self.covariance(:,:,ii); + %Check that covariance matrix is positive definite + try + chol(self.covariance(:,:,ii)); + catch + error('an entry in md.stochasticforcing.covariance is not positive definite'); + end + end + elseif(numel(size(self.covariance)==2)) + numtcovmat = 1; %number of covariance matrices in time + lsCovmats = {self.covariance}; + %Check that covariance matrix is positive definite + try + chol(self.covariance); + catch + error('md.stochasticforcing.covariance is not positive definite'); + end end @@ -176,28 +194,40 @@ if(indPwarma~=-1) if(md.basalforcings.arma_timestep~=md.hydrology.arma_timestep) - crossentries = reshape(self.covariance(1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]); - if any(crossentries~=0) - error('BasalforcingsDeepwaterMeltingRatearma and hydrologyarmapw have different arma_timestep and non-zero covariance'); + for(ii=[1:numel(lsCovmats)]) + covm = lsCovmats{ii}; + crossentries = reshape(covm(1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]); + if any(crossentries~=0) + error('BasalforcingsDeepwaterMeltingRatearma and hydrologyarmapw have different arma_timestep and non-zero covariance'); + end end end elseif(indSdarma~=-1) if(md.frontalforcings.sd_arma_timestep~=md.basalforcings.arma_timestep) - crossentries = reshape(self.covariance(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]); - if any(crossentries~=0) - error('FrontalForcingsSubglacialDischargearma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance'); + for(ii=[1:numel(lsCovmats)]) + covm = lsCovmats{ii}; + crossentries = reshape(covm(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]); + if any(crossentries~=0) + error('FrontalForcingsSubglacialDischargearma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance'); + end end end elseif(indSMBarma~=-1) if(md.smb.arma_timestep~=md.basalforcings.arma_timestep) - crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]); - if any(crossentries~=0) - error('SMBarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance'); + for(ii=[1:numel(lsCovmats)]) + covm = lsCovmats{ii}; + crossentries = reshape(covm(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]); + if any(crossentries~=0) + error('SMBarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance'); + end end end elseif(indTFarma~=-1) if(md.frontalforcings.arma_timestep~=md.basalforcings.arma_timestep) - crossentries = reshape(self.covariance(1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]); - if any(crossentries~=0) - error('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance'); + for(ii=[1:numel(lsCovmats)]) + covm = lsCovmats{ii}; + crossentries = reshape(covm(1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]); + if any(crossentries~=0) + error('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance'); + end end end @@ -206,21 +236,30 @@ if(indSdarma~=-1) if(md.frontalforcings.sd_arma_timestep~=md.hydrology.arma_timestep) - crossentries = reshape(self.covariance(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]); - if any(crossentries~=0) - error('FrontalForcingsSubglacialDischargearma and hydrologyarmapw have different arma_timestep and non-zero covariance'); + for(ii=[1:numel(lsCovmats)]) + covm = lsCovmats{ii}; + crossentries = reshape(covm(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]); + if any(crossentries~=0) + error('FrontalForcingsSubglacialDischargearma and hydrologyarmapw have different arma_timestep and non-zero covariance'); + end end end elseif(indSMBarma~=-1) if(md.smb.arma_timestep~=md.hydrology.arma_timestep) - crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]); - if any(crossentries~=0) - error('SMBarma and hydrologyarmapw have different arma_timestep and non-zero covariance'); + for(ii=[1:numel(lsCovmats)]) + covm = lsCovmats{ii}; + crossentries = reshape(covm(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]); + if any(crossentries~=0) + error('SMBarma and hydrologyarmapw have different arma_timestep and non-zero covariance'); + end end end elseif(indTFarma~=-1) if(md.frontalforcings.arma_timestep~=md.hydrology.arma_timestep) - crossentries = reshape(self.covariance(1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]); - if any(crossentries~=0) - error('FrontalForcingsRignotarma and hydrologyarmapw have different arma_timestep and non-zero covariance'); + for(ii=[1:numel(lsCovmats)]) + covm = lsCovmats{ii}; + crossentries = reshape(covm(1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]); + if any(crossentries~=0) + error('FrontalForcingsRignotarma and hydrologyarmapw have different arma_timestep and non-zero covariance'); + end end end @@ -229,14 +268,20 @@ if(indSMBarma~=-1) if(md.smb.arma_timestep~=md.frontalforcings.sd_arma_timestep) - crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma))),1,[]); - if any(crossentries~=0) - error('SMBarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance'); + for(ii=[1:numel(lsCovmats)]) + covm = lsCovmats{ii}; + crossentries = reshape(covm(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma))),1,[]); + if any(crossentries~=0) + error('SMBarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance'); + end end end elseif(indTFarma~=-1) if(md.frontalforcings.sd_arma_timestep~=md.frontalforcings.arma_timestep) - crossentries = reshape(self.covariance(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma))),1,[]); - if any(crossentries~=0) - error('FrontalForcingsRignotarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance'); + for(ii=[1:numel(lsCovmats)]) + covm = lsCovmats{ii}; + crossentries = reshape(covm(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma))),1,[]); + if any(crossentries~=0) + error('FrontalForcingsRignotarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance'); + end end end @@ -245,7 +290,10 @@ if(indTFarma~=-1) if(md.smb.arma_timestep~=md.frontalforcings.arma_timestep) - crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma))),1,[]); - if any(crossentries~=0) - error('SMBarma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance'); + for(ii=[1:numel(lsCovmats)]) + covm = lsCovmats{ii}; + crossentries = reshape(covm(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma))),1,[]); + if any(crossentries~=0) + error('SMBarma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance'); + end end end @@ -256,7 +304,10 @@ md = checkfield(md,'fieldname','stochasticforcing.isstochasticforcing','values',[0 1]); md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings()); - md = checkfield(md,'fieldname','stochasticforcing.covariance','NaN',1,'Inf',1,'size',[size_tot,size_tot]); %global covariance matrix + md = checkfield(md,'fieldname','stochasticforcing.covariance','NaN',1,'Inf',1,'size',[size_tot,size_tot,numtcovmat]); %global covariance matrix md = checkfield(md,'fieldname','stochasticforcing.stochastictimestep','NaN',1,'Inf',1,'>=',md.timestepping.time_step); md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]); + if(numtcovmat>1) %check the time steps at which each covariance matrix starts to be applied + md = checkfield(md,'fieldname','stochasticforcing.timecovariance','NaN',1,'Inf',1,'>=',md.timestepping.start_time,'<=',md.timestepping.final_time,'size',[1,numtcovmat]); + end if(checkdefaults) %need to check the defaults md = checkfield(md,'fieldname','stochasticforcing.defaultdimension','numel',1,'NaN',1,'Inf',1,'>',0); @@ -270,5 +321,6 @@ 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,'covariance',{'covariance matrix for within- and between-fields covariance (units must be squared field units)','multiple matrices can be concatenated along 3rd dimension to apply different covariances in time'}); + fielddisplay(self,'timecovariance','starting dates at which covariances apply (only applicabe if multiple covariance matrices are prescribed)'); fielddisplay(self,'stochastictimestep','timestep at which new stochastic noise terms are generated (default: md.timestepping.time_step)'); fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'); @@ -325,16 +377,42 @@ end + if(numel(size(self.covariance)==3)) + [nrow,ncol,numtcovmat] = size(self.covariance); + lsCovmats = {}; + for ii=[1:numtcovmat] %loop over 3rd dimension + lsCovmats{ii} = self.covariance(:,:,ii); + end + if(md.timestepping.interp_forcing==1) + disp('WARNING: md.timestepping.interp_forcing is 1, but be aware that there is no interpolation between covariance matrices'); + disp(' the changes between covariance matrices occur at the time steps specified in md.stochasticforcing.timecovariance'); + end + elseif(numel(size(self.covariance)==2)) + [nrow,ncol] = size(self.covariance); + numtcovmat = 1; %number of covariance matrices in time + lsCovmats = {self.covariance}; + end + %Scaling covariance matrix (scale column-by-column and row-by-row) scaledfields = {'BasalforcingsDeepwaterMeltingRatearma','BasalforcingsSpatialDeepwaterMeltingRate','DefaultCalving','FloatingMeltRate','SMBarma','SMBforcing'}; %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))]; - for row=inds %scale rows corresponding to scaled field - tempcovariance(row,:) = 1./yts*tempcovariance(row,:); - end - for col=inds %scale columns corresponding to scaled field - tempcovariance(:,col) = 1./yts*tempcovariance(:,col); - end + tempcovariance2d = zeros(numtcovmat,sum(nrow*ncol)); %covariance matrices in 2d array + % Loop over covariance matrices % + for kk=[1:numtcovmat] + kkcov = self.covariance(:,:,kk); %extract covariance at index kk + % Loop over the fields % + for i=1:num_fields + if any(strcmp(scaledfields,self.fields(i))) + inds = [1+sum(dimensions(1:i-1)):1:sum(dimensions(1:i))]; + for row=inds %scale rows corresponding to scaled field + kkcov(row,:) = 1./yts*kkcov(row,:); + end + for col=inds %scale columns corresponding to scaled field + kkcov(:,col) = 1./yts*kkcov(:,col); + end + end + end + % Save scaled covariance % + for rr=[1:nrow] + ind0 = 1+(rr-1)*ncol; + tempcovariance2d(kk,ind0:ind0+ncol-1) = kkcov(rr,:); end end @@ -342,4 +420,8 @@ if isnan(self.default_id) self.default_id = zeros(md.mesh.numberofelements,1); + end + %Set dummy timecovariance vector if a single covariance matrix is used + if(numtcovmat==1) + self.timecovariance = [md.timestepping.start_time]; end @@ -349,5 +431,7 @@ 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,'data',numtcovmat,'name','md.stochasticforcing.num_timescovariance','format','Integer'); + WriteData(fid,prefix,'data',tempcovariance2d,'name','md.stochasticforcing.covariance','format','DoubleMat'); + WriteData(fid,prefix,'object',self,'fieldname','timecovariance','format','DoubleMat','scale',yts); WriteData(fid,prefix,'object',self,'fieldname','stochastictimestep','format','Double','scale',yts); WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean'); Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.py =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 27779) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 27780) @@ -19,4 +19,5 @@ self.default_id = np.nan self.covariance = np.nan + self.timecovariance = np.nan self.stochastictimestep = 0 self.randomflag = 1 @@ -33,5 +34,6 @@ s += '{}\n'.format(fielddisplay(self, 'defaultdimension', 'dimensionality of the noise terms (does not apply to fields with their specific dimension)')) s += '{}\n'.format(fielddisplay(self, 'default_id', 'id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)')) - s += '{}\n'.format(fielddisplay(self, 'covariance', 'covariance matrix for within- and between-fields covariance (units must be squared field units)')) + s += '{}\n'.format(fielddisplay(self, 'covariance', 'covariance matrix for within- and between-fields covariance (units must be squared field units),multiple matrices can be concatenated along 3rd dimension to apply different covariances in time')) + s += '{}\n'.format(fielddisplay(self, 'timecovariance', 'starting dates at which covariances apply (only applicabe if multiple covariance matrices are prescribed)')) s += '{}\n'.format(fielddisplay(self, 'stochastictimestep', 'timestep at which new stochastic noise terms are generated (default: md.timestepping.time_step)')) s += '{}\n'.format(fielddisplay(self, 'randomflag', 'whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)')) @@ -69,9 +71,21 @@ print(' stochasticforcing.stocahstictimestep not specified: set to md.timestepping.time_step') - # Check that covariance matrix is positive definite (this is done internally by linalg) - try: - np.linalg.cholesky(self.covariance) - except: - raise TypeError('md.stochasticforcing.covariance is not positive definite') + if(len(np.shape(self.covariance))==3): + numtcovmat = np.shape(self.covariance)[2] #number of covariance matrices in time + lsCovmats = [] + for ii in range(numtcovmat): + lsCovmats.append(self.covariance[:,:,ii]) + try: + np.linalg.cholesky(self.covariance[:,:,ii]) + except: + raise TypeError('an entry in md.stochasticforcing.covariance is not positive definite') + elif(len(np.shape(self.covariance))==2): + numtcovmat = 1 + lsCovmats = [self.covariance] + # Check that covariance matrix is positive definite (this is done internally by linalg) + try: + np.linalg.cholesky(self.covariance) + except: + raise TypeError('md.stochasticforcing.covariance is not positive definite') # Check that all fields agree with the corresponding md class and if any field needs the default params @@ -171,52 +185,74 @@ if indBDWarma != -1: if indPwarma != -1: - covsum = self.covariance[np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)] - if md.smb.arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0): - raise IOError('BasalforcingsDeepwaterMeltingRatarma and hydrologyarmapw have different arma_timestep and non-zero covariance') + for ii in range(len(lsCovmats)): + covm = lsCovmats[ii] + covsum = covm[np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)] + if md.smb.arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0): + raise IOError('BasalforcingsDeepwaterMeltingRatarma and hydrologyarmapw have different arma_timestep and non-zero covariance') elif indSdarma != -1: - covsum = self.covariance[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)] - if md.frontalforcings.sd_arma_timestep != md.basalforcings.arma_timestep and np.any(covsum != 0): - raise IOError('FrontalForcingsSubglacialDischargearma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance') + for ii in range(len(lsCovmats)): + covm = lsCovmats[ii] + covsum = covm[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)] + if md.frontalforcings.sd_arma_timestep != md.basalforcings.arma_timestep and np.any(covsum != 0): + raise IOError('FrontalForcingsSubglacialDischargearma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance') elif indSMBarma != -1: - covsum = self.covariance[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)] - if md.smb.arma_timestep != md.basalforcings.arma_timestep and np.any(covsum != 0): - raise IOError('SMBarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance') + for ii in range(len(lsCovmats)): + covm = lsCovmats[ii] + covsum = covm[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)] + if md.smb.arma_timestep != md.basalforcings.arma_timestep and np.any(covsum != 0): + raise IOError('SMBarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance') elif indTFarma != -1: - covsum = self.covariance[np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)] - if md.frontalforcings.arma_timestep != md.basalforcings.arma_timestep and np.any(covsum != 0): - raise IOError('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance') + for ii in range(len(lsCovmats)): + covm = lsCovmats[ii] + covsum = covm[np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)] + if md.frontalforcings.arma_timestep != md.basalforcings.arma_timestep and np.any(covsum != 0): + raise IOError('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance') elif indPwarma != -1: if indSdarma != -1: - covsum = self.covariance[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)] - if md.frontalforcings.sd_arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0): - raise IOError('FrontalForingsSubglacialDischargearma and hydrologyarmapw have different arma_timestep and non-zero covariance') + for ii in range(len(lsCovmats)): + covm = lsCovmats[ii] + covsum = covm[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)] + if md.frontalforcings.sd_arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0): + raise IOError('FrontalForingsSubglacialDischargearma and hydrologyarmapw have different arma_timestep and non-zero covariance') elif indSMBarma != -1: - covsum = self.covariance[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)] - if md.smb.arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0): - raise IOError('SMBarma and hydrologyarmapw have different arma_timestep and non-zero covariance') + for ii in range(len(lsCovmats)): + covm = lsCovmats[ii] + covsum = covm[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)] + if md.smb.arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0): + raise IOError('SMBarma and hydrologyarmapw have different arma_timestep and non-zero covariance') elif indTFarma != -1: - covsum = self.covariance[np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)] - if md.frontalforcings.arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0): - raise IOError('FrontalForcingsRignotarma and hydrologyarmapw have different arma_timestep and non-zero covariance') + for ii in range(len(lsCovmats)): + covm = lsCovmats[ii] + covsum = covm[np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)] + if md.frontalforcings.arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0): + raise IOError('FrontalForcingsRignotarma and hydrologyarmapw have different arma_timestep and non-zero covariance') elif indSdarma != -1: if indSMBarma != -1: - covsum = self.covariance[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int)] - if md.smb.arma_timestep != md.frontalforcings.sd_arma_timestep and np.any(covsum != 0): - raise IOError('SMBarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance') + for ii in range(len(lsCovmats)): + covm = lsCovmats[ii] + covsum = covm[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int)] + if md.smb.arma_timestep != md.frontalforcings.sd_arma_timestep and np.any(covsum != 0): + raise IOError('SMBarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance') elif indTFarma != -1: - covsum = self.covariance[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int)] - if md.frontalforcings.sd_arma_timestep != md.frontalforcings.arma_timestep and np.any(covsum != 0): - raise IOError('FrontalForcingsSubglacialDischargearma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance') + for ii in range(len(lsCovmats)): + covm = lsCovmats[ii] + covsum = covm[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int)] + if md.frontalforcings.sd_arma_timestep != md.frontalforcings.arma_timestep and np.any(covsum != 0): + raise IOError('FrontalForcingsSubglacialDischargearma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance') elif indSMBarma != -1: if indTFarma != -1: - covsum = self.covariance[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int)] - if md.smb.arma_timestep != md.frontalforcings.arma_timestep and np.any(covsum != 0): - raise IOError('SMBarma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance') + for ii in range(len(lsCovmats)): + covm = lsCovmats[ii] + covsum = covm[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int)] + if md.smb.arma_timestep != md.frontalforcings.arma_timestep and np.any(covsum != 0): + raise IOError('SMBarma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance') md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1]) md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', self.supportedstochforcings()) - md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix + md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot, numtcovmat]) # global covariance matrix md = checkfield(md, 'fieldname', 'stochasticforcing.stochastictimestep', 'NaN', 1,'Inf', 1, '>=', md.timestepping.time_step) md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1]) + if(numtcovmat>1): #check the time steps at which each covariance matrix starts to be applie + md = checkfield(md, 'fieldname', 'stochasticforcing.timecovariance', 'NaN', 1, 'Inf', 1, '>=',md.timestepping.start_time,'<=',md.timestepping.final_time,'size',[1,numtcovmat]) # global covariance matrix if checkdefaults: md = checkfield(md, 'fieldname', 'stochasticforcing.defaultdimension', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) @@ -260,18 +296,43 @@ elif field == 'FrictionWaterPressure' and ispwHydroarma: dimensions[ind] = md.hydrology.num_basins + + if(len(np.shape(self.covariance))==3): + nrow,ncol,numtcovmat = np.shape(self.covariance) + lsCovmats = [] + for ii in range(numtcovmat): + lsCovmats.append(self.covariance[:,:,ii]) + if(md.timestepping.interp_forcing==1): + print('WARNING: md.timestepping.interp_forcing is 1, but be aware that there is no interpolation between covariance matrices') + print(' the changes between covariance matrices occur at the time steps specified in md.stochasticforcing.timecovariance') + elif(len(np.shape(self.covariance))==2): + nrow,ncol = np.shape(self.covariance) + numtcovmat = 1 + lsCovmats = [self.covariance] # Scaling covariance matrix (scale column-by-column and row-by-row) scaledfields = ['BasalforcingsDeepwaterMeltingRatearma','BasalforcingsSpatialDeepwaterMeltingRate','DefaultCalving', 'FloatingMeltRate', 'SMBarma', 'SMBforcing'] # list of fields that need scaling * 1 / yts - tempcovariance = np.copy(self.covariance) - for i in range(num_fields): - if self.fields[i] in scaledfields: - inds = range(int(np.sum(dimensions[0:i])), int(np.sum(dimensions[0:i + 1]))) - for row in inds: # scale rows corresponding to scaled field - tempcovariance[row, :] = 1 / yts * tempcovariance[row, :] - for col in inds: # scale columns corresponding to scaled field - tempcovariance[:, col] = 1 / yts * tempcovariance[:, col] + tempcovariance2d = np.zeros((numtcovmat,nrow*ncol)) + # Loop over covariance matrices # + for kk in range(numtcovmat): + kkcov = lsCovmats[kk] + # Loop over the fields # + for i in range(num_fields): + if self.fields[i] in scaledfields: + inds = range(int(np.sum(dimensions[0:i])), int(np.sum(dimensions[0:i + 1]))) + for row in inds: # scale rows corresponding to scaled field + kkcov[row, :] = 1 / yts * kkcov[row, :] + for col in inds: # scale columns corresponding to scaled field + kkcov[:, col] = 1 / yts * kkcov[:, col] + # Save scaled covariance # + for rr in range(nrow): + ind0 = rr*ncol + tempcovariance2d[kk,ind0:ind0+ncol] = np.copy(kkcov[rr,:]) + # Set dummy default_id vector if defaults not used if np.any(np.isnan(self.default_id)): self.default_id = np.zeros(md.mesh.numberofelements) + # Set dummy timecovariance vector if a single covariance matrix is used + if(numtcovmat==1): + self.timecovariance = np.array([md.timestepping.start_time]) # Reshape dimensions as column array for marshalling dimensions = dimensions.reshape(1, len(dimensions)) @@ -282,5 +343,7 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'default_id', 'data', self.default_id - 1, 'format', 'IntMat', 'mattype', 2) #12Nov2021 make sure this is zero-indexed! WriteData(fid, prefix, 'object', self, 'fieldname', 'defaultdimension', 'format', 'Integer') - WriteData(fid, prefix, 'data', tempcovariance, 'name', 'md.stochasticforcing.covariance', 'format', 'DoubleMat') + WriteData(fid, prefix, 'data', numtcovmat, 'name', 'md.stochasticforcing.num_timescovariance', 'format', 'Integer') + WriteData(fid, prefix, 'data', tempcovariance2d, 'name', 'md.stochasticforcing.covariance', 'format', 'DoubleMat') + WriteData(fid, prefix, 'object', self, 'fieldname', 'timecovariance', 'format', 'DoubleMat', 'scale', yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'stochastictimestep', 'format', 'Double', 'scale', yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'randomflag', 'format', 'Boolean') Index: /issm/trunk-jpl/test/NightlyRun/test546.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test546.m (revision 27780) +++ /issm/trunk-jpl/test/NightlyRun/test546.m (revision 27780) @@ -0,0 +1,144 @@ +%Test Name: PigTranMultcovStochasticforcings +md=triangle(model(),'../Exp/Pig.exp',8000.); +md=setmask(md,'../Exp/PigShelves.exp','../Exp/PigIslands.exp'); +md=parameterize(md,'../Par/Pig.par'); +md=setflowequation(md,'SSA','all'); +md.timestepping.start_time = 0; +md.timestepping.time_step = 1; +md.timestepping.final_time = 10; + +%Basin separation +idb = zeros(md.mesh.numberofelements,1); +iid1 = find(md.mesh.x>=-1.6e6); +for ii=1:md.mesh.numberofelements + for vertex=1:3 + if any(iid1==md.mesh.elements(ii,vertex)) %one vertex in basin 1 + idb(ii) = 1; + end + end + if idb(ii)==0 %no vertex was found in basin 1 + idb(ii) = 2; + end +end +nb_bas = 2; + +%SMB +numparams = 1; +numbreaks = 0; +intercept = [0.5;0.01]; +polynomialparams = intercept; +datebreaks = NaN; +md.smb = SMBarma(); +md.smb.num_basins = nb_bas; %number of basins +md.smb.basin_id = idb; %prescribe basin ID number to elements +md.smb.num_params = numparams; %number of parameters in the polynomial +md.smb.num_breaks = numbreaks; %number of breakpoints +md.smb.polynomialparams = polynomialparams; +md.smb.datebreaks = datebreaks; +md.smb.ar_order = 4; +md.smb.ma_order = 4; +md.smb.arma_timestep = 2.0; %timestep of the ARMA model [yr] +md.smb.arlag_coefs = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1]]; +md.smb.malag_coefs = [[0.1,0.1,0.2,0.3];[0.5,0.8,1.3,2.4]]; + +%Calving +md.mask.ice_levelset = 1e4*(md.mask.ice_levelset + 0.5); +md.calving.calvingrate = 0.1*ones(md.mesh.numberofvertices,1); +md.levelset.spclevelset = NaN(md.mesh.numberofvertices,1); +md.levelset.migration_max = 10.0; +md.frontalforcings.meltingrate = zeros(md.mesh.numberofvertices,1); + +% Basal forcing implementation +numparams = 2; +numbreaks = 1; +intercept = [3.0,4.0;1.0,0.5]; +trendlin = [0.0,0.1;0,0]; +polynomialparams = cat(3,intercept,trendlin); +datebreaks = [6;7]; +md.basalforcings = linearbasalforcingsarma(); +md.basalforcings.num_basins = nb_bas; %number of basins +md.basalforcings.basin_id = idb; %prescribe basin ID number to elements +md.basalforcings.polynomialparams = polynomialparams; +md.basalforcings.datebreaks = datebreaks; +md.basalforcings.num_params = numparams; %number of parameters in the polynomial +md.basalforcings.num_breaks = numbreaks; %number of breakpoints +md.basalforcings.ar_order = 1; +md.basalforcings.ma_order = 1; +md.basalforcings.arma_timestep = 1.0; %timestep of the ARMA model [yr] +md.basalforcings.arlag_coefs = [0.0;0.1]; +md.basalforcings.malag_coefs = [0.55;0.34]; +md.basalforcings.deepwater_elevation = [-1000,-1520]; +md.basalforcings.upperwater_elevation = [0,-50]; +md.basalforcings.upperwater_melting_rate = [0.0,0.0]; +md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1); + +% Covariance matrices +sdvsmb = [1,1]; +sdvclv = [0.1,0.01]; +sdvdwm = [300,300]; +vecsdv = [sdvsmb,sdvclv,sdvdwm]; +corrmat = [1,0,0,0,0,0; + 0,1,0,0,0,0; + 0,0,1,0.4,0.1,0.1; + 0,0,0.4,1,0.1,0.1; + 0,0,0.1,0.1,1,0.3; + 0,0,0.1,0.1,0.3,1]; +covglob0 = diag(vecsdv)*corrmat*diag(vecsdv); +covglob1 = 2*covglob0; +multcov = cat(3,covglob0,covglob1); +tmcov = [0,5]; + +% Stochastic forcing +md.stochasticforcing.isstochasticforcing = 1; +md.stochasticforcing.fields = [{'SMBarma'},{'DefaultCalving'},{'BasalforcingsDeepwaterMeltingRatearma'}]; +md.stochasticforcing.defaultdimension = 2; +md.stochasticforcing.default_id = idb; +md.stochasticforcing.covariance = multcov; %global covariance among- and between-fields +md.stochasticforcing.timecovariance = tmcov; %simulation times when covariance matrix switches +md.stochasticforcing.randomflag = 0; %determines true/false randomness + +md.transient.ismovingfront = 1; +md.transient.requested_outputs = {'default','SmbMassBalance','BasalforcingsFloatingiceMeltingRate','BasalforcingsSpatialDeepwaterMeltingRate'}; +md.transient.isstressbalance = 1; +md.transient.ismasstransport = 1; +md.transient.issmb = 1; +md.transient.isthermal = 0; +md.transient.isgroundingline = 1; + +md.cluster=generic('name',oshostname(),'np',2); +md=solve(md,'Transient'); + +%Fields and tolerances to track changes +field_names ={... + 'Vx1' ,'Vy1' ,'Vel1' ,'Thickness1', 'SmbMassBalance1', 'BasalforcingsFloatingiceMeltingRate1', 'BasalforcingsSpatialDeepwaterMeltingRate1',... + 'Vx2' ,'Vy2' ,'Vel2' ,'Thickness2', 'SmbMassBalance2' ,'BasalforcingsFloatingiceMeltingRate2', 'BasalforcingsSpatialDeepwaterMeltingRate2',... + 'Vx3' ,'Vy3' ,'Vel3' ,'Thickness3', 'SmbMassBalance3' ,'BasalforcingsFloatingiceMeltingRate3', 'BasalforcingsSpatialDeepwaterMeltingRate3',... + }; +field_tolerances={... + 1e-11,1e-11,2e-11,1e-11,1e-10,1e-9,1e-10,... + 1e-11,1e-11,2e-11,9e-11,1e-10,1e-9,1e-10,... + 2e-10,2e-10,2e-10,1e-10,1e-10,1e-9,1e-10,... + }; +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).SmbMassBalance),... + (md.results.TransientSolution(1).BasalforcingsFloatingiceMeltingRate),... + (md.results.TransientSolution(1).BasalforcingsSpatialDeepwaterMeltingRate),... + (md.results.TransientSolution(5).Vx),... + (md.results.TransientSolution(5).Vy),... + (md.results.TransientSolution(5).Vel),... + (md.results.TransientSolution(5).Thickness),... + (md.results.TransientSolution(5).SmbMassBalance),... + (md.results.TransientSolution(5).BasalforcingsFloatingiceMeltingRate),... + (md.results.TransientSolution(5).BasalforcingsSpatialDeepwaterMeltingRate),... + (md.results.TransientSolution(10).Vx),... + (md.results.TransientSolution(10).Vy),... + (md.results.TransientSolution(10).Vel),... + (md.results.TransientSolution(10).Thickness),... + (md.results.TransientSolution(10).SmbMassBalance),... + (md.results.TransientSolution(10).BasalforcingsFloatingiceMeltingRate),... + (md.results.TransientSolution(10).BasalforcingsSpatialDeepwaterMeltingRate),... + }; Index: /issm/trunk-jpl/test/NightlyRun/test546.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test546.py (revision 27780) +++ /issm/trunk-jpl/test/NightlyRun/test546.py (revision 27780) @@ -0,0 +1,156 @@ +#Test Name: PigTranARMAandStochasticforcings +import numpy as np +from linearbasalforcingsarma import linearbasalforcingsarma +from SMBarma import SMBarma +from stochasticforcing import stochasticforcing +from socket import gethostname +from model import * +from parameterize import parameterize +from setflowequation import setflowequation +from setmask import setmask +from solve import solve +from triangle import triangle + + +md = triangle(model(), '../Exp/Pig.exp', 8000) +md = setmask(md, '../Exp/PigShelves.exp', '../Exp/PigIslands.exp') +md = parameterize(md, '../Par/Pig.py') +md = setflowequation(md, 'SSA', 'all') +md.timestepping.start_time = 0 +md.timestepping.time_step = 1 +md.timestepping.final_time = 10 + +# Basin separation +idb = np.zeros((md.mesh.numberofelements,)) +iid1 = np.where(md.mesh.x >= -1.6e6)[0] +for ii in range(md.mesh.numberofelements): + for vertex in range(3): + if md.mesh.elements[ii][vertex] - 1 in iid1: # one vertex in basin 1; NOTE: offset because of 1-based vertex indexing + idb[ii] = 1 + if idb[ii] == 0: # no vertex was found in basin 1 + for vertex in range(3): + idb[ii] = 2 +nb_bas = 2 + +#SMB +numparams = 1; +numbreaks = 0; +intercept = np.array([[0.5],[0.01]]) +polynomialparams = np.copy(intercept) +datebreaks = np.nan +md.smb = SMBarma() +md.smb.num_basins = nb_bas # number of basins +md.smb.basin_id = idb # prescribe basin ID number to elements; +md.smb.num_params = 1*numparams +md.smb.num_breaks = 1*numbreaks +md.smb.polynomialparams = 1*polynomialparams +md.smb.datebreaks = 1*datebreaks +md.smb.ar_order = 4 +md.smb.ma_order = 4 +md.smb.arma_timestep = 2.0 #timestep of the ARMA model [yr] +md.smb.arlag_coefs = np.array([[0.2,0.1,0.05,0.01],[0.4,0.2,-0.2,0.1]]) +md.smb.malag_coefs = np.array([[0.1,0.1,0.2,0.3],[0.5,0.8,1.3,2.4]]) + +#Calving +md.mask.ice_levelset = 1e4*(md.mask.ice_levelset + 0.5) +md.calving.calvingrate = 0.1*np.ones((md.mesh.numberofvertices,)) +md.levelset.spclevelset = np.full((md.mesh.numberofvertices,), np.nan) +md.levelset.migration_max = 10.0 +md.frontalforcings.meltingrate = np.zeros((md.mesh.numberofvertices,)) + +#Basal forcing implementation +numparams = 2 +numbreaks = 1 +intercept = np.array([[3.0,4.0],[1.0,0.5]]) +trendlin = np.array([[0.0,0.1],[0,0]]) +polynomialparams = np.transpose(np.stack((intercept,trendlin)),(1,2,0)) +datebreaks = np.array([[6],[7]]) + +md.basalforcings = linearbasalforcingsarma() +md.basalforcings.num_basins = nb_bas +md.basalforcings.basin_id = idb +md.basalforcings.const = np.array([[1.0, 2.50]]) # intercept values of DeepwaterMelt in basins [m/yr] +md.basalforcings.trend = np.array([[0.2, 0.01]]) # trend values of DeepwaterMelt in basins [m/yr^2] +md.basalforcings.arma_initialtime = md.timestepping.start_time # initial time in the ARMA model parameterization [yr] +md.basalforcings.ar_order = 1 +md.basalforcings.ma_order = 1 +md.basalforcings.polynomialparams = 1*polynomialparams; +md.basalforcings.datebreaks = 1*datebreaks; +md.basalforcings.num_params = 1*numparams +md.basalforcings.num_breaks = 1*numbreaks +md.basalforcings.arma_timestep = 1.0 # timestep of the ARMA model [yr] +md.basalforcings.arlag_coefs = np.array([[0.0], [0.1]]) # autoregressive parameters +md.basalforcings.malag_coefs = np.array([[0.55], [0.34]]) # moving-average parameters +md.basalforcings.deepwater_elevation = np.array([[-1000, -1520]]) +md.basalforcings.upperwater_elevation = np.array([[0, -50]]) +md.basalforcings.upperwater_melting_rate = np.array([[0,0]]) +md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices,)) + +#Covariance matrices +sdvsmb = np.array([1,1]) +sdvclv = np.array([0.1,0.01]) +sdvdwm = np.array([300,300]) +vecsdv = np.concatenate((sdvsmb,sdvclv,sdvdwm)) +corrmat = np.array([[1.0, 0., 0., 0., 0., 0.], + [0., 1.0, 0., 0., 0., 0.], + [0., 0., 1.0, 0.4, 0.1, 0.1], + [0., 0., 0.4, 1.0, 0.1, 0.1], + [0., 0., 0.1, 0.1, 1.0, 0.3], + [0., 0., 0.1, 0.1, 0.3, 1.0]]) +covglob0 = np.diag(vecsdv) @ corrmat @ np.diag(vecsdv) +covglob1 = 2*covglob0 +multcov = np.stack((covglob0,covglob1),axis=2) +tmcov = np.array([[0,5]]) + +#Stochastic forcing +md.stochasticforcing.isstochasticforcing = 1 +md.stochasticforcing.fields = ['SMBarma', 'DefaultCalving', 'BasalforcingsDeepwaterMeltingRatearma'] +md.stochasticforcing.defaultdimension = 2 +md.stochasticforcing.default_id = idb +md.stochasticforcing.covariance = multcov # global covariance among- and between-fields +md.stochasticforcing.timecovariance = tmcov #simulation times when covariance matrix switches +md.stochasticforcing.randomflag = 0 # determines true/false randomness + +md.transient.ismovingfront = 1 +md.transient.requested_outputs = ['default', 'SmbMassBalance', 'BasalforcingsFloatingiceMeltingRate', 'BasalforcingsSpatialDeepwaterMeltingRate'] +md.transient.isstressbalance = 1 +md.transient.ismasstransport = 1 +md.transient.issmb = 1 +md.transient.isthermal = 0 +md.transient.isgroundingline = 1 + +md.cluster = generic('name', gethostname(), 'np', 2) +md = solve(md, 'Transient') + +# Fields and tolerances to track changes +field_names = [ + 'Vx1', 'Vy1', 'Vel1', 'Thickness1', 'SmbMassBalance1', 'BasalforcingsFloatingiceMeltingRate1', 'BasalforcingsSpatialDeepwaterMeltingRate1', + 'Vx5', 'Vy5', 'Vel5', 'Thickness5', 'SmbMassBalance5', 'BasalforcingsFloatingiceMeltingRate5', 'BasalforcingsSpatialDeepwaterMeltingRate5', + 'Vx10', 'Vy10', 'Vel10', 'Thickness10', 'SmbMassBalance10', 'BasalforcingsFloatingiceMeltingRate10', 'BasalforcingsSpatialDeepwaterMeltingRate10'] + +field_tolerances = [ + 1e-11, 1e-11, 2e-11, 1e-11, 1e10, 1e-9, 1e-10, + 1e-11, 1e-11, 2e-11, 9e-11, 1e10, 1e-9, 1e-10, + 2e-10, 2e-10, 2e-10, 1e-10, 1e10, 1e-9, 1e-10] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].SmbMassBalance, + md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[0].BasalforcingsSpatialDeepwaterMeltingRate, + md.results.TransientSolution[4].Vx, + md.results.TransientSolution[4].Vy, + md.results.TransientSolution[4].Vel, + md.results.TransientSolution[4].Thickness, + md.results.TransientSolution[4].SmbMassBalance, + md.results.TransientSolution[4].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[4].BasalforcingsSpatialDeepwaterMeltingRate, + md.results.TransientSolution[9].Vx, + md.results.TransientSolution[9].Vy, + md.results.TransientSolution[9].Vel, + md.results.TransientSolution[9].Thickness, + md.results.TransientSolution[9].SmbMassBalance, + md.results.TransientSolution[9].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[9].BasalforcingsSpatialDeepwaterMeltingRate]