Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 27296) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 27297) @@ -419,4 +419,7 @@ /*Nothing to add to parameters*/ break; + case SMBdebrisMLEnum: + /*Nothing to add to parameters*/ + break; default: _error_("Surface mass balance model "<AddObject(iomodel->CopyConstantObject("md.smb.rdl",SmbRdlEnum)); break; + case SMBdebrisMLEnum: + break; default: _error_("Surface mass balance model "<parameters->FindParam(&Rhoi,MaterialsRhoIceEnum); // Note Carlo's model used as benchmark was run with different densities for debris and FS + femmodel->parameters->FindParam(&Lm,MaterialsLatentheatEnum); + femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); + PhiD=0.; + //if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum); + + /* Loop over all the elements of this partition */ + for(Object* & object : femmodel->elements->objects){ + Element* element=xDynamicCast(object); + + /* Allocate all arrays */ + int numvertices=element->GetNumberOfVertices(); + IssmDouble* surfacelist=xNew(numvertices); + IssmDouble* smb=xNew(numvertices); + IssmDouble* debriscover=xNew(numvertices); + element->GetInputListOnVertices(surfacelist,SurfaceEnum); + + /* Get inputs */ + element->GetInputListOnVertices(debriscover,DebrisThicknessEnum); + + /*Loop over all vertices of element and calculate SMB as function of Debris Cover and z */ + for(int v=0;vnmax){ + Alphaeff=Alphad; + } else { + Alphaeff=Alphai+n*(Alphad-Alphai)/nmax; + }*/ + + // M&L + IssmDouble Alphaeff=Alphad; + + /* compute smb */ + for (int ismb=0;ismb<2;ismb++){ + if(ismb==0){ + // calc a reference smb to identify accum and melt region; debris only develops in ablation area + debris=0.; + }else{ + // only in the meltregime debris develops + if(-MassBalanceCmDayDebris<0) debris=debriscover[v]; + } + MassBalanceCmDayDebris=(((In-(z-ReferenceElevation)*LW/100.)-(Eps*Sigma*(Tf*Tf*Tf*Tf))+ + (Q+(z-ReferenceElevation)*SW/100.)*(1.-Alphaeff)+ + (Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* + WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))*(Tm-(z- + ReferenceElevation)*AirTemp/100.))/((1-PhiD)*Rhoi*Lm)/(1.+ + ((Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* + WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/ + K*debris)-(Lv*Ustar*Ustar*(Qh-(Qh*(Humidity-(z- + ReferenceElevation)*HumidityG/100.)))*(exp(-Gamma*Xr)))/((1.-PhiD)* + Rhoi*Lm*Ustar)/((((Um-(z-ReferenceElevation)*WindSpeed/100.) + -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)))*100.*24.*60.*60.; + } + + /* account form ablation days, and convert to m/s */ + MassBalanceMYearDebris=-MassBalanceCmDayDebris/100.*AblationDays/yts; + + /*Update array accordingly*/ + smb[v]=MassBalanceMYearDebris; + } + + /*Add input to element and Free memory*/ + element->AddInput(SmbMassBalanceEnum,smb,P1Enum); + xDelete(surfacelist); + xDelete(smb); + xDelete(debriscover); + } +}/*}}}*/ void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 27296) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 27297) @@ -23,4 +23,5 @@ void SmbMeltComponentsx(FemModel* femmodel); void SmbGradientsComponentsx(FemModel* femmodel); +void SmbDebrisMLx(FemModel* femmodel); /* SEMIC: */ void SmbSemicx(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27296) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27297) @@ -157,4 +157,5 @@ DamageStressUBoundEnum, DebugProfilingEnum, + DebrisThicknessEnum, DomainDimensionEnum, DomainTypeEnum, @@ -1541,4 +1542,5 @@ SMBarmaEnum, SMBcomponentsEnum, + SMBdebrisMLEnum, SMBd18opddEnum, SMBforcingEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27296) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27297) @@ -165,4 +165,5 @@ case DamageStressUBoundEnum : return "DamageStressUBound"; case DebugProfilingEnum : return "DebugProfiling"; + case DebrisThicknessEnum : return "DebrisThickness"; case DomainDimensionEnum : return "DomainDimension"; case DomainTypeEnum : return "DomainType"; @@ -1544,4 +1545,5 @@ case SMBarmaEnum : return "SMBarma"; case SMBcomponentsEnum : return "SMBcomponents"; + case SMBdebrisMLEnum : return "SMBdebrisML"; case SMBd18opddEnum : return "SMBd18opdd"; case SMBforcingEnum : return "SMBforcing"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27296) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27297) @@ -168,4 +168,5 @@ else if (strcmp(name,"DamageStressUBound")==0) return DamageStressUBoundEnum; else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum; + else if (strcmp(name,"DebrisThickness")==0) return DebrisThicknessEnum; else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum; else if (strcmp(name,"DomainType")==0) return DomainTypeEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"HydrologydcSedimentThickness")==0) return HydrologydcSedimentThicknessEnum; else if (strcmp(name,"HydrologyStepAdapt")==0) return HydrologyStepAdaptEnum; - else if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"HydrologydcUnconfinedFlag")==0) return HydrologydcUnconfinedFlagEnum; + if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum; + else if (strcmp(name,"HydrologydcUnconfinedFlag")==0) return HydrologydcUnconfinedFlagEnum; else if (strcmp(name,"HydrologyshreveStabilization")==0) return HydrologyshreveStabilizationEnum; else if (strcmp(name,"IcecapToEarthComm")==0) return IcecapToEarthCommEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"QmuResponsedescriptors")==0) return QmuResponsedescriptorsEnum; else if (strcmp(name,"QmuVariableDescriptors")==0) return QmuVariableDescriptorsEnum; - else if (strcmp(name,"QmuVariablePartitions")==0) return QmuVariablePartitionsEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"QmuVariablePartitionsNpart")==0) return QmuVariablePartitionsNpartEnum; + if (strcmp(name,"QmuVariablePartitions")==0) return QmuVariablePartitionsEnum; + else if (strcmp(name,"QmuVariablePartitionsNpart")==0) return QmuVariablePartitionsNpartEnum; else if (strcmp(name,"QmuVariablePartitionsNt")==0) return QmuVariablePartitionsNtEnum; else if (strcmp(name,"QmuResponsePartitions")==0) return QmuResponsePartitionsEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"SmbElevationBins")==0) return SmbElevationBinsEnum; else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum; - else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"SmbDelta18oSurface")==0) return SmbDelta18oSurfaceEnum; + if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum; + else if (strcmp(name,"SmbDelta18oSurface")==0) return SmbDelta18oSurfaceEnum; else if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum; else if (strcmp(name,"SmbDt")==0) return SmbDtEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum; else if (strcmp(name,"Velocity")==0) return VelocityEnum; - else if (strcmp(name,"Xxe")==0) return XxeEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"Yye")==0) return YyeEnum; + if (strcmp(name,"Xxe")==0) return XxeEnum; + else if (strcmp(name,"Yye")==0) return YyeEnum; else if (strcmp(name,"Zze")==0) return ZzeEnum; else if (strcmp(name,"Areae")==0) return AreaeEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"Domain3D")==0) return Domain3DEnum; else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum; - else if (strcmp(name,"DrivingStressX")==0) return DrivingStressXEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"DrivingStressY")==0) return DrivingStressYEnum; + if (strcmp(name,"DrivingStressX")==0) return DrivingStressXEnum; + else if (strcmp(name,"DrivingStressY")==0) return DrivingStressYEnum; else if (strcmp(name,"Dummy")==0) return DummyEnum; else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum; else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum; - else if (strcmp(name,"Misfit")==0) return MisfitEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum; + if (strcmp(name,"Misfit")==0) return MisfitEnum; + else if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum; else if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum; else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; @@ -997,9 +998,9 @@ else if (strcmp(name,"SmbC")==0) return SmbCEnum; else if (strcmp(name,"SmbCcsnowValue")==0) return SmbCcsnowValueEnum; - else if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum; + if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum; + else if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum; else if (strcmp(name,"SmbD")==0) return SmbDEnum; else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum; - else if (strcmp(name,"Surface")==0) return SurfaceEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum; + if (strcmp(name,"Surface")==0) return SurfaceEnum; + else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum; else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"Outputdefinition66")==0) return Outputdefinition66Enum; else if (strcmp(name,"Outputdefinition67")==0) return Outputdefinition67Enum; - else if (strcmp(name,"Outputdefinition68")==0) return Outputdefinition68Enum; else stage=11; } if(stage==11){ - if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum; + if (strcmp(name,"Outputdefinition68")==0) return Outputdefinition68Enum; + else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum; else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum; else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; @@ -1366,9 +1367,9 @@ else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum; else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; - else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"Element")==0) return ElementEnum; + if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; + else if (strcmp(name,"Element")==0) return ElementEnum; else if (strcmp(name,"ElementHook")==0) return ElementHookEnum; else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum; @@ -1489,9 +1490,9 @@ else if (strcmp(name,"MantlePlumeGeothermalFlux")==0) return MantlePlumeGeothermalFluxEnum; else if (strcmp(name,"MassFlux")==0) return MassFluxEnum; - else if (strcmp(name,"Masscon")==0) return MassconEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; + if (strcmp(name,"Masscon")==0) return MassconEnum; + else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; @@ -1580,4 +1581,5 @@ else if (strcmp(name,"SMBarma")==0) return SMBarmaEnum; else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; + else if (strcmp(name,"SMBdebrisML")==0) return SMBdebrisMLEnum; else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; @@ -1611,10 +1613,10 @@ else if (strcmp(name,"Separate")==0) return SeparateEnum; else if (strcmp(name,"Seq")==0) return SeqEnum; - else if (strcmp(name,"SmbAnalysis")==0) return SmbAnalysisEnum; - else if (strcmp(name,"SmbSolution")==0) return SmbSolutionEnum; else stage=14; } if(stage==14){ - if (strcmp(name,"SmoothAnalysis")==0) return SmoothAnalysisEnum; + if (strcmp(name,"SmbAnalysis")==0) return SmbAnalysisEnum; + else if (strcmp(name,"SmbSolution")==0) return SmbSolutionEnum; + else if (strcmp(name,"SmoothAnalysis")==0) return SmoothAnalysisEnum; else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum; else if (strcmp(name,"SpatialLinearFloatingMeltRate")==0) return SpatialLinearFloatingMeltRateEnum; Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 27296) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 27297) @@ -240,6 +240,7 @@ case 10: return SMBpddSicopolisEnum; case 11: return SMBgradientscomponentsEnum; - case 12: return SMBsemicEnum; + case 12: return SMBsemicEnum; case 13: return SMBarmaEnum; + case 14: return SMBdebrisMLEnum; default: _error_("Marshalled SMB code \""<=',1,'numel',[1]); + md = checkfield(md,'fieldname','smb.averaging','numel',[1],'values',[0 1 2]); + md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1); + end % }}} + function disp(self) % {{{ + disp(sprintf(' surface forcings parameters:')); + fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'); + fielddisplay(self, 'averaging', 'averaging methods from short to long steps'); + disp(sprintf('%51s 0: Arithmetic (default)',' ')); + disp(sprintf('%51s 1: Geometric',' ')); + disp(sprintf('%51s 2: Harmonic',' ')); + fielddisplay(self,'requested_outputs','additional outputs requested'); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + yts=md.constants.yts; + + WriteData(fid,prefix,'name','md.smb.model','data',14,'format','Integer'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer'); + + %process requested outputs + outputs = self.requested_outputs; + pos = find(ismember(outputs,'default')); + if ~isempty(pos), + outputs(pos) = []; %remove 'default' from outputs + outputs = [outputs defaultoutputs(self,md)]; %add defaults + end + WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray'); + + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/initialization.m =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.m (revision 27296) +++ /issm/trunk-jpl/src/m/classes/initialization.m (revision 27297) @@ -25,4 +25,5 @@ str = NaN; sample = NaN; + debris = NaN; end methods @@ -121,4 +122,9 @@ if ~isnan(md.initialization.sample) md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + end + end + if ismember('DebrisAnalysis',analyses), + if ~isnan(md.initialization.debris) + md = checkfield(md,'fieldname','initialization.debris','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); end end @@ -145,4 +151,5 @@ fielddisplay(self,'dsl','Dynamic sea level.'); fielddisplay(self,'str','Steric sea level.'); + fielddisplay(self,'debris','Surface debris layer [m]'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -167,4 +174,5 @@ WriteData(fid,prefix,'object',self,'fieldname','hydraulic_potential','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','debris','format','DoubleMat','mattype',1); if md.thermal.isenthalpy, @@ -197,4 +205,5 @@ self.dsl=project3d(md,'vector',self.dsl,'type','node','layer',1); self.str=project3d(md,'vector',self.str,'type','node','layer',1); + self.str=project3d(md,'vector',self.debris,'type','node','layer',1); %Lithostatic pressure by default @@ -218,4 +227,5 @@ writejs1Darray(fid,[modelname '.initialization.channel'],self.channelarea); writejs1Darray(fid,[modelname '.initialization.sample'],self.sample); + writejs1Darray(fid,[modelname '.initialization.debris'],self.debris); end % }}}