Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 26858)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 26859)
@@ -545,4 +545,5 @@
       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));
       iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields");
       if(num_fields<1) _error_("no stochasticforcing fields found");
Index: /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp	(revision 26858)
+++ /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp	(revision 26859)
@@ -26,24 +26,43 @@
    femmodel->parameters->FindParam(&covariance,&M,&N,StochasticForcingCovarianceEnum); _assert_(M==dimtot); _assert_(N==dimtot);
 
+	/*Check if this is a timestep for new noiseterms computation*/
+	bool isstepforstoch = false;
+	IssmDouble time,dt,starttime,tstep_stoch;
+   femmodel->parameters->FindParam(&time,TimeEnum);
+   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+   femmodel->parameters->FindParam(&tstep_stoch,StochasticForcingTimestepEnum);
+		
+	#ifndef _HAVE_AD_
+   if((fmod(time,tstep_stoch)<fmod((time-dt),tstep_stoch)) || (time<=starttime+dt) || tstep_stoch==dt) isstepforstoch = true;
+   #else
+   _error_("not implemented yet");
+   #endif
+
    /*Compute noise terms*/
-   IssmDouble* noiseterms = xNew<IssmDouble>(dimtot);
-   my_rank=IssmComm::GetRank();
-   if(my_rank==0){
-      int fixedseed;
-      IssmDouble time,dt,starttime;
-      femmodel->parameters->FindParam(&time,TimeEnum);
-      femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-      femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
-		/*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
-      if(randomflag) fixedseed=-1;
-      else fixedseed = reCast<int,IssmDouble>((time-starttime)/dt);
-		/*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
-      IssmDouble* temparray = NULL;
-      multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed);
-      for(int i=0;i<dimtot;i++) noiseterms[i]=temparray[i];
+	IssmDouble* noiseterms = xNew<IssmDouble>(dimtot);
+   if(isstepforstoch){
+		my_rank=IssmComm::GetRank();
+   	if(my_rank==0){
+   	   int fixedseed;
+			/*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
+   	   if(randomflag) fixedseed=-1;
+   	   else fixedseed = reCast<int,IssmDouble>((time-starttime)/dt);
+			/*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
+   	   IssmDouble* temparray = NULL;
+   	   multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed);
+   	   for(int i=0;i<dimtot;i++) noiseterms[i]=temparray[i];
+			xDelete<IssmDouble>(temparray);
+   	}
+   	ISSM_MPI_Bcast(noiseterms,dimtot,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+		femmodel->parameters->SetParam(noiseterms,dimtot,StochasticForcingNoisetermsEnum);
+	}
+	else{
+		IssmDouble* temparray = NULL;
+		femmodel->parameters->FindParam(&temparray,&N,StochasticForcingNoisetermsEnum); _assert_(N==dimtot);
+		for(int i=0;i<dimtot;i++) noiseterms[i] = temparray[i];
 		xDelete<IssmDouble>(temparray);
-   }
-   ISSM_MPI_Bcast(noiseterms,dimtot,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-   
+	}
+
 	int i=0;
    for(int j=0;j<numfields;j++){
Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 26858)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 26859)
@@ -113,4 +113,9 @@
 syn keyword cConstant CalvingMinthicknessEnum
 syn keyword cConstant CalvingTestSpeedfactorEnum
+syn keyword cConstant CalvingUseParamEnum
+syn keyword cConstant CalvingScaleThetaEnum
+syn keyword cConstant CalvingAmpAlphaEnum
+syn keyword cConstant CalvingMidpointEnum
+syn keyword cConstant CalvingNonlinearLawEnum
 syn keyword cConstant ConfigurationTypeEnum
 syn keyword cConstant ConstantsGEnum
@@ -428,6 +433,8 @@
 syn keyword cConstant StochasticForcingIsStochasticForcingEnum
 syn keyword cConstant StochasticForcingIsWaterPressureEnum
+syn keyword cConstant StochasticForcingNoisetermsEnum
 syn keyword cConstant StochasticForcingNumFieldsEnum
 syn keyword cConstant StochasticForcingRandomflagEnum
+syn keyword cConstant StochasticForcingTimestepEnum
 syn keyword cConstant SolidearthSettingsReltolEnum
 syn keyword cConstant SolidearthSettingsSelfAttractionEnum
@@ -1248,4 +1255,5 @@
 syn keyword cConstant CalvingLevermannEnum
 syn keyword cConstant CalvingTestEnum
+syn keyword cConstant CalvingParameterizationEnum
 syn keyword cConstant CalvingVonmisesEnum
 syn keyword cConstant CfdragcoeffabsgradEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26858)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26859)
@@ -427,6 +427,8 @@
 	StochasticForcingIsStochasticForcingEnum,
 	StochasticForcingIsWaterPressureEnum,
+	StochasticForcingNoisetermsEnum,
 	StochasticForcingNumFieldsEnum,
 	StochasticForcingRandomflagEnum,
+	StochasticForcingTimestepEnum,
 	SolidearthSettingsReltolEnum,
 	SolidearthSettingsSelfAttractionEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 26858)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 26859)
@@ -115,4 +115,9 @@
 		case CalvingMinthicknessEnum : return "CalvingMinthickness";
 		case CalvingTestSpeedfactorEnum : return "CalvingTestSpeedfactor";
+		case CalvingUseParamEnum : return "CalvingUseParam";
+		case CalvingScaleThetaEnum : return "CalvingScaleTheta";
+		case CalvingAmpAlphaEnum : return "CalvingAmpAlpha";
+		case CalvingMidpointEnum : return "CalvingMidpoint";
+		case CalvingNonlinearLawEnum : return "CalvingNonlinearLaw";
 		case ConfigurationTypeEnum : return "ConfigurationType";
 		case ConstantsGEnum : return "ConstantsG";
@@ -430,6 +435,8 @@
 		case StochasticForcingIsStochasticForcingEnum : return "StochasticForcingIsStochasticForcing";
 		case StochasticForcingIsWaterPressureEnum : return "StochasticForcingIsWaterPressure";
+		case StochasticForcingNoisetermsEnum : return "StochasticForcingNoiseterms";
 		case StochasticForcingNumFieldsEnum : return "StochasticForcingNumFields";
 		case StochasticForcingRandomflagEnum : return "StochasticForcingRandomflag";
+		case StochasticForcingTimestepEnum : return "StochasticForcingTimestep";
 		case SolidearthSettingsReltolEnum : return "SolidearthSettingsReltol";
 		case SolidearthSettingsSelfAttractionEnum : return "SolidearthSettingsSelfAttraction";
@@ -1250,4 +1257,5 @@
 		case CalvingLevermannEnum : return "CalvingLevermann";
 		case CalvingTestEnum : return "CalvingTest";
+		case CalvingParameterizationEnum : return "CalvingParameterization";
 		case CalvingVonmisesEnum : return "CalvingVonmises";
 		case CfdragcoeffabsgradEnum : return "Cfdragcoeffabsgrad";
Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 26858)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 26859)
@@ -106,4 +106,9 @@
 syn keyword juliaConstC CalvingMinthicknessEnum
 syn keyword juliaConstC CalvingTestSpeedfactorEnum
+syn keyword juliaConstC CalvingUseParamEnum
+syn keyword juliaConstC CalvingScaleThetaEnum
+syn keyword juliaConstC CalvingAmpAlphaEnum
+syn keyword juliaConstC CalvingMidpointEnum
+syn keyword juliaConstC CalvingNonlinearLawEnum
 syn keyword juliaConstC ConfigurationTypeEnum
 syn keyword juliaConstC ConstantsGEnum
@@ -421,6 +426,8 @@
 syn keyword juliaConstC StochasticForcingIsStochasticForcingEnum
 syn keyword juliaConstC StochasticForcingIsWaterPressureEnum
+syn keyword juliaConstC StochasticForcingNoisetermsEnum
 syn keyword juliaConstC StochasticForcingNumFieldsEnum
 syn keyword juliaConstC StochasticForcingRandomflagEnum
+syn keyword juliaConstC StochasticForcingTimestepEnum
 syn keyword juliaConstC SolidearthSettingsReltolEnum
 syn keyword juliaConstC SolidearthSettingsSelfAttractionEnum
@@ -1241,4 +1248,5 @@
 syn keyword juliaConstC CalvingLevermannEnum
 syn keyword juliaConstC CalvingTestEnum
+syn keyword juliaConstC CalvingParameterizationEnum
 syn keyword juliaConstC CalvingVonmisesEnum
 syn keyword juliaConstC CfdragcoeffabsgradEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26858)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26859)
@@ -115,4 +115,9 @@
 	      else if (strcmp(name,"CalvingMinthickness")==0) return CalvingMinthicknessEnum;
 	      else if (strcmp(name,"CalvingTestSpeedfactor")==0) return CalvingTestSpeedfactorEnum;
+	      else if (strcmp(name,"CalvingUseParam")==0) return CalvingUseParamEnum;
+	      else if (strcmp(name,"CalvingScaleTheta")==0) return CalvingScaleThetaEnum;
+	      else if (strcmp(name,"CalvingAmpAlpha")==0) return CalvingAmpAlphaEnum;
+	      else if (strcmp(name,"CalvingMidpoint")==0) return CalvingMidpointEnum;
+	      else if (strcmp(name,"CalvingNonlinearLaw")==0) return CalvingNonlinearLawEnum;
 	      else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
 	      else if (strcmp(name,"ConstantsG")==0) return ConstantsGEnum;
@@ -132,13 +137,13 @@
 	      else if (strcmp(name,"CumGmtslc")==0) return CumGmtslcEnum;
 	      else if (strcmp(name,"CumGmslc")==0) return CumGmslcEnum;
-	      else if (strcmp(name,"DamageC1")==0) return DamageC1Enum;
+         else stage=2;
+   }
+   if(stage==2){
+	      if (strcmp(name,"DamageC1")==0) return DamageC1Enum;
 	      else if (strcmp(name,"DamageC2")==0) return DamageC2Enum;
 	      else if (strcmp(name,"DamageC3")==0) return DamageC3Enum;
 	      else if (strcmp(name,"DamageC4")==0) return DamageC4Enum;
 	      else if (strcmp(name,"Damage")==0) return DamageEnum;
-         else stage=2;
-   }
-   if(stage==2){
-	      if (strcmp(name,"DamageEquivStress")==0) return DamageEquivStressEnum;
+	      else if (strcmp(name,"DamageEquivStress")==0) return DamageEquivStressEnum;
 	      else if (strcmp(name,"DamageEvolutionNumRequestedOutputs")==0) return DamageEvolutionNumRequestedOutputsEnum;
 	      else if (strcmp(name,"DamageEvolutionRequestedOutputs")==0) return DamageEvolutionRequestedOutputsEnum;
@@ -255,13 +260,13 @@
 	      else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
 	      else if (strcmp(name,"InversionControlScalingFactors")==0) return InversionControlScalingFactorsEnum;
-	      else if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum;
+         else stage=3;
+   }
+   if(stage==3){
+	      if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum;
 	      else if (strcmp(name,"InversionDxmin")==0) return InversionDxminEnum;
 	      else if (strcmp(name,"InversionGatol")==0) return InversionGatolEnum;
 	      else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum;
 	      else if (strcmp(name,"InversionGrtol")==0) return InversionGrtolEnum;
-         else stage=3;
-   }
-   if(stage==3){
-	      if (strcmp(name,"InversionGttol")==0) return InversionGttolEnum;
+	      else if (strcmp(name,"InversionGttol")==0) return InversionGttolEnum;
 	      else if (strcmp(name,"InversionIncompleteAdjoint")==0) return InversionIncompleteAdjointEnum;
 	      else if (strcmp(name,"InversionIscontrol")==0) return InversionIscontrolEnum;
@@ -378,13 +383,13 @@
 	      else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
 	      else if (strcmp(name,"SolidearthPartitionIce")==0) return SolidearthPartitionIceEnum;
-	      else if (strcmp(name,"SolidearthPartitionHydro")==0) return SolidearthPartitionHydroEnum;
+         else stage=4;
+   }
+   if(stage==4){
+	      if (strcmp(name,"SolidearthPartitionHydro")==0) return SolidearthPartitionHydroEnum;
 	      else if (strcmp(name,"SolidearthPartitionOcean")==0) return SolidearthPartitionOceanEnum;
 	      else if (strcmp(name,"SolidearthNpartIce")==0) return SolidearthNpartIceEnum;
 	      else if (strcmp(name,"SolidearthNpartOcean")==0) return SolidearthNpartOceanEnum;
 	      else if (strcmp(name,"SolidearthNpartHydro")==0) return SolidearthNpartHydroEnum;
-         else stage=4;
-   }
-   if(stage==4){
-	      if (strcmp(name,"SolidearthPlanetRadius")==0) return SolidearthPlanetRadiusEnum;
+	      else if (strcmp(name,"SolidearthPlanetRadius")==0) return SolidearthPlanetRadiusEnum;
 	      else if (strcmp(name,"SolidearthPlanetArea")==0) return SolidearthPlanetAreaEnum;
 	      else if (strcmp(name,"SolidearthSettingsAbstol")==0) return SolidearthSettingsAbstolEnum;
@@ -439,6 +444,8 @@
 	      else if (strcmp(name,"StochasticForcingIsStochasticForcing")==0) return StochasticForcingIsStochasticForcingEnum;
 	      else if (strcmp(name,"StochasticForcingIsWaterPressure")==0) return StochasticForcingIsWaterPressureEnum;
+	      else if (strcmp(name,"StochasticForcingNoiseterms")==0) return StochasticForcingNoisetermsEnum;
 	      else if (strcmp(name,"StochasticForcingNumFields")==0) return StochasticForcingNumFieldsEnum;
 	      else if (strcmp(name,"StochasticForcingRandomflag")==0) return StochasticForcingRandomflagEnum;
+	      else if (strcmp(name,"StochasticForcingTimestep")==0) return StochasticForcingTimestepEnum;
 	      else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum;
 	      else if (strcmp(name,"SolidearthSettingsSelfAttraction")==0) return SolidearthSettingsSelfAttractionEnum;
@@ -499,5 +506,8 @@
 	      else if (strcmp(name,"SmbLapseRatePos")==0) return SmbLapseRatePosEnum;
 	      else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum;
-	      else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
+         else stage=5;
+   }
+   if(stage==5){
+	      if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
 	      else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum;
 	      else if (strcmp(name,"SmbPhi")==0) return SmbPhiEnum;
@@ -506,8 +516,5 @@
 	      else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
 	      else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum;
-         else stage=5;
-   }
-   if(stage==5){
-	      if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum;
+	      else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum;
 	      else if (strcmp(name,"SmbRunoffalti")==0) return SmbRunoffaltiEnum;
 	      else if (strcmp(name,"SmbRunoffgrad")==0) return SmbRunoffgradEnum;
@@ -622,5 +629,8 @@
 	      else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum;
 	      else if (strcmp(name,"BasalforcingsLinearBasinId")==0) return BasalforcingsLinearBasinIdEnum;
-	      else if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum;
+         else stage=6;
+   }
+   if(stage==6){
+	      if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum;
 	      else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum;
 	      else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum;
@@ -629,8 +639,5 @@
 	      else if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum;
 	      else if (strcmp(name,"BasalforcingsIsmip6Tf")==0) return BasalforcingsIsmip6TfEnum;
-         else stage=6;
-   }
-   if(stage==6){
-	      if (strcmp(name,"BasalforcingsIsmip6TfShelf")==0) return BasalforcingsIsmip6TfShelfEnum;
+	      else if (strcmp(name,"BasalforcingsIsmip6TfShelf")==0) return BasalforcingsIsmip6TfShelfEnum;
 	      else if (strcmp(name,"BasalforcingsIsmip6MeltAnomaly")==0) return BasalforcingsIsmip6MeltAnomalyEnum;
 	      else if (strcmp(name,"BasalforcingsMeltrateFactor")==0) return BasalforcingsMeltrateFactorEnum;
@@ -745,5 +752,8 @@
 	      else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum;
 	      else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
-	      else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
+         else stage=7;
+   }
+   if(stage==7){
+	      if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
 	      else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
 	      else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
@@ -752,8 +762,5 @@
 	      else if (strcmp(name,"FrontalForcingsSubglacialDischarge")==0) return FrontalForcingsSubglacialDischargeEnum;
 	      else if (strcmp(name,"FrontalForcingsThermalForcing")==0) return FrontalForcingsThermalForcingEnum;
-         else stage=7;
-   }
-   if(stage==7){
-	      if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
+	      else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
 	      else if (strcmp(name,"NGia")==0) return NGiaEnum;
 	      else if (strcmp(name,"NGiaRate")==0) return NGiaRateEnum;
@@ -868,5 +875,8 @@
 	      else if (strcmp(name,"SealevelBarystaticIceArea")==0) return SealevelBarystaticIceAreaEnum;
 	      else if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum;
-	      else if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum;
+         else stage=8;
+   }
+   if(stage==8){
+	      if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum;
 	      else if (strcmp(name,"SealevelBarystaticIceLoad")==0) return SealevelBarystaticIceLoadEnum;
 	      else if (strcmp(name,"SealevelBarystaticHydroMask")==0) return SealevelBarystaticHydroMaskEnum;
@@ -875,8 +885,5 @@
 	      else if (strcmp(name,"SealevelBarystaticHydroLatbar")==0) return SealevelBarystaticHydroLatbarEnum;
 	      else if (strcmp(name,"SealevelBarystaticHydroLongbar")==0) return SealevelBarystaticHydroLongbarEnum;
-         else stage=8;
-   }
-   if(stage==8){
-	      if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum;
+	      else if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum;
 	      else if (strcmp(name,"SealevelBarystaticBpMask")==0) return SealevelBarystaticBpMaskEnum;
 	      else if (strcmp(name,"SealevelBarystaticBpWeights")==0) return SealevelBarystaticBpWeightsEnum;
@@ -991,5 +998,8 @@
 	      else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum;
 	      else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
-	      else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
+         else stage=9;
+   }
+   if(stage==9){
+	      if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
 	      else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
 	      else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
@@ -998,8 +1008,5 @@
 	      else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum;
 	      else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum;
-         else stage=9;
-   }
-   if(stage==9){
-	      if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
+	      else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
 	      else if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum;
 	      else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum;
@@ -1114,5 +1121,8 @@
 	      else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
 	      else if (strcmp(name,"VyBase")==0) return VyBaseEnum;
-	      else if (strcmp(name,"Vy")==0) return VyEnum;
+         else stage=10;
+   }
+   if(stage==10){
+	      if (strcmp(name,"Vy")==0) return VyEnum;
 	      else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
 	      else if (strcmp(name,"VyObs")==0) return VyObsEnum;
@@ -1121,8 +1131,5 @@
 	      else if (strcmp(name,"Vz")==0) return VzEnum;
 	      else if (strcmp(name,"VzFS")==0) return VzFSEnum;
-         else stage=10;
-   }
-   if(stage==10){
-	      if (strcmp(name,"VzHO")==0) return VzHOEnum;
+	      else if (strcmp(name,"VzHO")==0) return VzHOEnum;
 	      else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
 	      else if (strcmp(name,"VzSSA")==0) return VzSSAEnum;
@@ -1237,5 +1244,8 @@
 	      else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
 	      else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
-	      else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum;
+         else stage=11;
+   }
+   if(stage==11){
+	      if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum;
 	      else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum;
 	      else if (strcmp(name,"Outputdefinition100")==0) return Outputdefinition100Enum;
@@ -1244,8 +1254,5 @@
 	      else if (strcmp(name,"AdaptiveTimestepping")==0) return AdaptiveTimesteppingEnum;
 	      else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
-         else stage=11;
-   }
-   if(stage==11){
-	      if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
+	      else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
 	      else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
 	      else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum;
@@ -1280,4 +1287,5 @@
 	      else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
 	      else if (strcmp(name,"CalvingTest")==0) return CalvingTestEnum;
+	      else if (strcmp(name,"CalvingParameterization")==0) return CalvingParameterizationEnum;
 	      else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
 	      else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
@@ -1359,5 +1367,8 @@
 	      else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
 	      else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
-	      else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
+         else stage=12;
+   }
+   if(stage==12){
+	      if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
 	      else if (strcmp(name,"GenericExternalResult")==0) return GenericExternalResultEnum;
 	      else if (strcmp(name,"Gradient1")==0) return Gradient1Enum;
@@ -1367,8 +1378,5 @@
 	      else if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum;
 	      else if (strcmp(name,"GroundedAreaScaled")==0) return GroundedAreaScaledEnum;
-         else stage=12;
-   }
-   if(stage==12){
-	      if (strcmp(name,"GroundingOnly")==0) return GroundingOnlyEnum;
+	      else if (strcmp(name,"GroundingOnly")==0) return GroundingOnlyEnum;
 	      else if (strcmp(name,"GroundinglineMassFlux")==0) return GroundinglineMassFluxEnum;
 	      else if (strcmp(name,"Gset")==0) return GsetEnum;
@@ -1482,5 +1490,8 @@
 	      else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
 	      else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
-	      else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
+         else stage=13;
+   }
+   if(stage==13){
+	      if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
 	      else if (strcmp(name,"None")==0) return NoneEnum;
 	      else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
@@ -1490,8 +1501,5 @@
 	      else if (strcmp(name,"OceantransportAnalysis")==0) return OceantransportAnalysisEnum;
 	      else if (strcmp(name,"OceantransportSolution")==0) return OceantransportSolutionEnum;
-         else stage=13;
-   }
-   if(stage==13){
-	      if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
+	      else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
 	      else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
 	      else if (strcmp(name,"Open")==0) return OpenEnum;
@@ -1605,5 +1613,8 @@
 	      else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
 	      else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
-	      else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
+         else stage=14;
+   }
+   if(stage==14){
+	      if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
 	      else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
 	      else if (strcmp(name,"Tria")==0) return TriaEnum;
@@ -1613,8 +1624,5 @@
 	      else if (strcmp(name,"Vertex")==0) return VertexEnum;
 	      else if (strcmp(name,"VertexLId")==0) return VertexLIdEnum;
-         else stage=14;
-   }
-   if(stage==14){
-	      if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
+	      else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
 	      else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
 	      else if (strcmp(name,"Vertices")==0) return VerticesEnum;
Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/stochasticforcing.m	(revision 26858)
+++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m	(revision 26859)
@@ -11,4 +11,5 @@
 		default_id				= NaN;
 		covariance				= NaN;
+		stochastictimestep   = 0;
 		randomflag				= 1;
 	end
@@ -35,4 +36,7 @@
 
 			num_fields = numel(self.fields);
+			if(self.stochastictimestep==0)
+				md.stochasticforcing.stochastictimestep = md.timestepping.time_step; %by default: stochastictimestep set to ISSM time step
+			end
 
 			%Check that covariance matrix is positive definite
@@ -116,12 +120,21 @@
 				indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with other ar timesteps 
 				dimensions(indSMBar) = md.smb.num_basins;
+				if(md.smb.ar_timestep<self.stochastictimestep)
+					error('SMBautoregression cannot have a timestep shorter than stochastictimestep');
+				end
 			end
 			if any(contains(self.fields,'FrontalForcingsRignotAutoregression'))
 				indTFar	= find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with other ar timesteps 
 				dimensions(indTFar) = md.frontalforcings.num_basins;
+				if(md.frontalforcings.ar_timestep<self.stochastictimestep)
+					error('FrontalForcingsRignotAutoregression cannot have a timestep shorter than stochastictimestep');
+				end
 			end
 			if any(contains(self.fields,'BasalforcingsDeepwaterMeltingRateAutoregression'))
 				indBDWar	= find(contains(self.fields,'BasalforcingsDeepwaterMeltingRateAutoregression')); %index of BDWar, now check for consistency with other ar timesteps 
 				dimensions(indBDWar) = md.basalforcings.num_basins;
+				if(md.basalforcings.ar_timestep<self.stochastictimestep)
+					error('BasalforcingsDeepwaterMeltingRateAutoregression cannot have a timestep shorter than stochastictimestep');
+				end
 			end
 			size_tot = sum(dimensions);
@@ -154,4 +167,5 @@
 			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.stochastictimestep','NaN',1,'Inf',1,'>=',md.timestepping.time_step);
 			md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]);
 			if(checkdefaults) %need to check the defaults
@@ -167,4 +181,5 @@
 			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,'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)');
 			disp('Available fields:');
@@ -182,4 +197,9 @@
 				return
 			else
+
+				if(self.stochastictimestep==0)
+					disp('      stochasticforcing.stocahstictimestep not specified: set to md.timestepping.time_step');
+					self.stochastictimestep = md.timestepping.time_step; %by default: stochastictimestep set to ISSM time step
+				end
 
 				%Retrieve dimensionality of each field
@@ -225,4 +245,5 @@
 				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','stochastictimestep','format','Double','scale',yts);
 				WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean');
 			end
Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/stochasticforcing.py	(revision 26858)
+++ /issm/trunk-jpl/src/m/classes/stochasticforcing.py	(revision 26859)
@@ -18,4 +18,5 @@
         self.default_id = np.nan
         self.covariance = np.nan
+        self.stochastictimestep = 0
         self.randomflag = 1
 
@@ -32,4 +33,5 @@
         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, '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)'))
         s += 'Available fields:\n'
@@ -58,4 +60,7 @@
 
         num_fields = len(self.fields)
+        if(self.stochastictimestep==0):
+            md.stochasticforcing.stochastictimestep = md.timestepping.time_step #by default: stochastictimestep set to ISSM time step
+            print('      stochasticforcing.stocahstictimestep not specified: set to md.timestepping.time_step')
 
         # Check that covariance matrix is positive definite (this is done internally by linalg)
@@ -119,10 +124,16 @@
             indSMBar = self.fields.index('SMBautoregression')  # Index of SMBar, now check for consistency with other timesteps
             dimensions[indSMBar] = md.smb.num_basins
+            if(md.smb.ar_timestep<self.stochastictimestep):
+                raise TypeError('SMBautoregression cannot have a timestep shorter than stochastictimestep')
         if ('FrontalForcingsRignotAutoregression' in self.fields):
             indTFar = self.fields.index('FrontalForcingsRignotAutoregression')  # Index of TFar, now check for consistency with other timesteps
             dimensions[indTFar] = md.frontalforcings.num_basins
+            if(md.frontalforcings.ar_timestep<self.stochastictimestep):
+                raise TypeError('FrontalForcingsRignotAutoregression cannot have a timestep shorter than stochastictimestep')
         if ('BasalforcingsDeepwaterMeltingRateAutoregression' in self.fields):
             indBDWar = self.fields.index('BasalforcingsDeepwaterMeltingRateAutoregression')  # Index of BDWar, now check for consistency with other timesteps
             dimensions[indTFar] = md.basalforcings.num_basins
+            if(md.basalforcings.ar_timestep<self.stochastictimestep):
+                raise TypeError('BasalforcingsDeepwaterMeltingRateAutoregression cannot have a timestep shorter than stochastictimestep')
         size_tot = np.sum(dimensions)
 
@@ -143,4 +154,5 @@
         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.stochastictimestep', 'NaN', 1,'Inf', 1, '>=', md.timestepping.time_step)
         md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1])
         if (checkdefaults):
@@ -164,4 +176,6 @@
         else:
             num_fields = len(self.fields)
+            if(self.stochastictimestep==0):
+                md.stochasticforcing.stochastictimestep = md.timestepping.time_step #by default: stochastictimestep set to ISSM time step
             # Retrieve dimensionality of each field
             dimensions = self.defaultdimension * np.ones((num_fields))
@@ -197,4 +211,5 @@
             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', 'stochastictimestep', 'format', 'Double', 'scale', yts)
             WriteData(fid, prefix, 'object', self, 'fieldname', 'randomflag', 'format', 'Boolean')
     # }}}
Index: /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m	(revision 26858)
+++ /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m	(revision 26859)
@@ -43,4 +43,5 @@
 		['/Users/larour/ModelData/BedMachine/' basename '-' ncdate '.nc'],...
 		['./' basename '-' ncdate '.nc'],...
+		'/media/vincent/TOSH4TB/GeorgiaTech/DataSearch/BedMachine/BedMachineGreenland-2021-04-20.nc',...
 		};
 
@@ -87,7 +88,9 @@
 if strcmp(string,'mask') | strcmp(string,'source'),
 	%Need nearest neighbor to avoid interpolation between 0 and 2
-	output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest');
+	%output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest');
+	output = InterpFromGridToMesh(xdata,flipud(ydata),flipud(data),double(X),double(Y),NaN); %VV
 else
-	output = InterpFromGrid(xdata,ydata,data,double(X),double(Y));
+	%output = InterpFromGrid(xdata,ydata,data,double(X),double(Y));
+	output = InterpFromGridToMesh(xdata,flipud(ydata),flipud(data),double(X),double(Y),NaN); %VV
 end