Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 24144) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 24145) @@ -110,4 +110,8 @@ case BasalforcingsIsmip6Enum: iomodel->FetchDataToInput(elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum); + break; + case BeckmannGoosseFloatingMeltRateEnum: + iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum); + iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum); break; default: Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 24144) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 24145) @@ -198,4 +198,8 @@ xDelete(array3d); } + break; + case BeckmannGoosseFloatingMeltRateEnum: + iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum); + iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum); break; default: Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 24144) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 24145) @@ -786,4 +786,8 @@ iomodel->FetchDataToInput(elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum); break; + case BeckmannGoosseFloatingMeltRateEnum: + iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum); + iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum); + break; default: _error_("Basal forcing model "<FetchDataToInput(elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum); break; - + case BeckmannGoosseFloatingMeltRateEnum: + iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum); + iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum); + break; default: _error_("Basal forcing model "<(values); +}/*}}}*/ +void Element::BeckmannGoosseFloatingiceMeltingRate(){/*{{{*/ + + int numvertices = this->GetNumberOfVertices(); + IssmDouble meltratefactor,T_f,ocean_heat_flux; + IssmDouble rho_water = this->FindParam(MaterialsRhoSeawaterEnum); + IssmDouble rho_ice = this->FindParam(MaterialsRhoIceEnum); + IssmDouble latentheat = this->FindParam(MaterialsLatentheatEnum); + IssmDouble mixed_layer_capacity = this->FindParam(MaterialsMixedLayerCapacityEnum); + IssmDouble thermal_exchange_vel = this->FindParam(MaterialsThermalExchangeVelocityEnum); + + IssmDouble* base = xNew(numvertices); + IssmDouble* values = xNew(numvertices); + IssmDouble* oceansalinity = xNew(numvertices); + IssmDouble* oceantemp = xNew(numvertices); + + this->GetInputListOnVertices(base,BaseEnum); + this->GetInputListOnVertices(oceansalinity,BasalforcingsOceanSalinityEnum); + this->GetInputListOnVertices(oceantemp,BasalforcingsOceanTempEnum); + parameters->FindParam(&meltratefactor,BasalforcingsMeltrateFactorEnum); + + Gauss* gauss=this->NewGauss(); + for(int i=0;i T_ice ==> heat flux FROM ocean TO ice + ocean_heat_flux = meltratefactor * rho_water * mixed_layer_capacity * thermal_exchange_vel * (oceantemp[i] - T_f); // in W/m^2 + + // shelfbmassflux is positive if ice is freezing on; here it is always negative: + // same sign as ocean_heat_flux (positive if massflux FROM ice TO ocean) + values[i] = ocean_heat_flux / (latentheat * rho_ice); // m s-1 + } + + this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum); + xDelete(base); + xDelete(values); + xDelete(oceantemp); + xDelete(oceansalinity); + delete gauss; }/*}}}*/ void Element::MungsmtpParameterization(void){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24144) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24145) @@ -144,9 +144,10 @@ void MigrateGroundingLine(IssmDouble* sheet_ungrounding); void MismipFloatingiceMeltingRate(); + void BeckmannGoosseFloatingiceMeltingRate(); void MungsmtpParameterization(void); ElementMatrix* NewElementMatrix(int approximation_enum=NoneApproximationEnum); ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum); ElementVector* NewElementVector(int approximation_enum=NoneApproximationEnum); - void PicoUpdateBoxid(int* pmax_boxid_basin); + void PicoUpdateBoxid(int* pmax_boxid_basin); void PicoUpdateBox(int loopboxid); void PicoComputeBasalMelt(); Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 24144) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 24145) @@ -40,4 +40,8 @@ if(VerboseSolution())_printf0_(" call ISMIP 6 Floating melting rate module\n"); FloatingiceMeltingRateIsmip6x(femmodel); + break; + case BeckmannGoosseFloatingMeltRateEnum: + if(VerboseSolution())_printf0_(" call BeckmannGoosse Floating melting rate module\n"); + BeckmannGoosseFloatingiceMeltingRatex(femmodel); break; default: @@ -193,2 +197,10 @@ } /*}}}*/ +void BeckmannGoosseFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/ + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->BeckmannGoosseFloatingiceMeltingRate(); + } +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h (revision 24144) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h (revision 24145) @@ -15,4 +15,5 @@ void MismipFloatingiceMeltingRatex(FemModel* femmodel); void FloatingiceMeltingRateIsmip6x(FemModel* femmodel); +void BeckmannGoosseFloatingiceMeltingRatex(FemModel* femmodel); #endif /* _FloatingiceMeltingRatex_H*/ Index: /issm/trunk-jpl/src/c/modules/GeothermalFluxx/GeothermalFluxx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GeothermalFluxx/GeothermalFluxx.cpp (revision 24144) +++ /issm/trunk-jpl/src/c/modules/GeothermalFluxx/GeothermalFluxx.cpp (revision 24145) @@ -28,4 +28,7 @@ MantlePlumeGeothermalFluxx(femmodel); break; + case BeckmannGoosseFloatingMeltRateEnum: + /*Nothing to be done*/ + break; default: _error_("Basal forcing model "<AddObject(iomodel->CopyConstantObject("md.basalforcings.num_basins",BasalforcingsIsmip6NumBasinsEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_0",BasalforcingsIsmip6Gamma0Enum)); - parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.islocal",BasalforcingsIsmip6IsLocalEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.islocal",BasalforcingsIsmip6IsLocalEnum)); iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.delta_t"); parameters->AddObject(new DoubleVecParam(BasalforcingsIsmip6DeltaTEnum,transparam,N)); @@ -245,4 +245,7 @@ parameters->AddObject(new DoubleVecParam(BasalforcingsIsmip6TfDepthsEnum,transparam,N)); xDelete(transparam); + break; + case BeckmannGoosseFloatingMeltRateEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.meltrate_factor",BasalforcingsMeltrateFactorEnum)); break; default: Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24144) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24145) @@ -457,4 +457,6 @@ BasalforcingsIsmip6TfShelfEnum, BasalforcingsIsmip6MeltAnomalyEnum, + BasalforcingsOceanSalinityEnum, + BasalforcingsOceanTempEnum, BasalforcingsPicoBasinIdEnum, BasalforcingsPicoBoxIdEnum, @@ -929,4 +931,5 @@ BasalforcingsIsmip6Enum, BasalforcingsPicoEnum, + BeckmannGoosseFloatingMeltRateEnum, BedSlopeSolutionEnum, BoolExternalResultEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24144) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24145) @@ -463,4 +463,6 @@ case BasalforcingsIsmip6TfShelfEnum : return "BasalforcingsIsmip6TfShelf"; case BasalforcingsIsmip6MeltAnomalyEnum : return "BasalforcingsIsmip6MeltAnomaly"; + case BasalforcingsOceanSalinityEnum : return "BasalforcingsOceanSalinity"; + case BasalforcingsOceanTempEnum : return "BasalforcingsOceanTemp"; case BasalforcingsPicoBasinIdEnum : return "BasalforcingsPicoBasinId"; case BasalforcingsPicoBoxIdEnum : return "BasalforcingsPicoBoxId"; @@ -933,4 +935,5 @@ case BasalforcingsIsmip6Enum : return "BasalforcingsIsmip6"; case BasalforcingsPicoEnum : return "BasalforcingsPico"; + case BeckmannGoosseFloatingMeltRateEnum : return "BeckmannGoosseFloatingMeltRate"; case BedSlopeSolutionEnum : return "BedSlopeSolution"; case BoolExternalResultEnum : return "BoolExternalResult"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24144) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24145) @@ -472,4 +472,6 @@ else if (strcmp(name,"BasalforcingsIsmip6TfShelf")==0) return BasalforcingsIsmip6TfShelfEnum; else if (strcmp(name,"BasalforcingsIsmip6MeltAnomaly")==0) return BasalforcingsIsmip6MeltAnomalyEnum; + else if (strcmp(name,"BasalforcingsOceanSalinity")==0) return BasalforcingsOceanSalinityEnum; + else if (strcmp(name,"BasalforcingsOceanTemp")==0) return BasalforcingsOceanTempEnum; else if (strcmp(name,"BasalforcingsPicoBasinId")==0) return BasalforcingsPicoBasinIdEnum; else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum; @@ -504,10 +506,10 @@ else if (strcmp(name,"DamageF")==0) return DamageFEnum; else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum; - else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; - else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; + if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; + else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; + else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"MaterialsRheologyEs")==0) return MaterialsRheologyEsEnum; else if (strcmp(name,"MaterialsRheologyEsbar")==0) return MaterialsRheologyEsbarEnum; - else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum; - else if (strcmp(name,"MeshScaleFactor")==0) return MeshScaleFactorEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum; + if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum; + else if (strcmp(name,"MeshScaleFactor")==0) return MeshScaleFactorEnum; + else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum; else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum; else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum; @@ -750,10 +752,10 @@ else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum; else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum; - else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; - else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbTz")==0) return SmbTzEnum; + if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; + else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum; + else if (strcmp(name,"SmbTz")==0) return SmbTzEnum; else if (strcmp(name,"SmbV")==0) return SmbVEnum; else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; - else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum; - else if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum; + if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum; + else if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum; + else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum; else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum; else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum; @@ -954,4 +956,5 @@ else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum; else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum; + else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum; else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; @@ -995,11 +998,11 @@ else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum; - else if (strcmp(name,"Divergence")==0) return DivergenceEnum; - else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum; - else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; + if (strcmp(name,"Divergence")==0) return DivergenceEnum; + else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum; + else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum; + else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum; else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; @@ -1118,11 +1121,11 @@ else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; - else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; - else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; - else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"Materials")==0) return MaterialsEnum; + if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; + else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; + else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum; + else if (strcmp(name,"Materials")==0) return MaterialsEnum; else if (strcmp(name,"Matestar")==0) return MatestarEnum; else if (strcmp(name,"Matice")==0) return MaticeEnum; @@ -1241,11 +1244,11 @@ else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum; else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum; - else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum; - else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum; - else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum; + if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum; + else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum; + else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum; + else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum; else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum; else if (strcmp(name,"StringParam")==0) return StringParamEnum; Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 24144) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 24145) @@ -198,4 +198,5 @@ case 6: return SpatialLinearFloatingMeltRateEnum; case 7: return BasalforcingsIsmip6Enum; + case 8: return BeckmannGoosseFloatingMeltRateEnum; default: _error_("Marshalled Basal Forcings code \""<'); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','', '
',' geothermal heat flux [W/m^2] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' basal melting rate (positive if melting) [m/yr] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' ocean_temp [degC] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' ocean_salinity [psu] ',''); + end % }}} + function self = extrude(self,md) % {{{ + self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1); + self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux + end % }}} + function self = basalforcingsbeckmanngoosse(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + case 1 + self=structtoobj(basalforcingsbeckmanngoosse(),varargin{1}); + otherwise + error('constructor not supported'); + end + end % }}} + function self = initialize(self,md) % {{{ + + if isnan(self.groundedice_melting_rate), + self.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); + disp(' no basalforcings.groundedice_melting_rate specified: values set as zero'); + end + if isnan(self.ocean_temp), + self.ocean_temp=-1.7*ones(md.mesh.numberofvertices,1); + disp(' no basalforcings.ocean_temp specified: values set as -1.7degC'); + end + if isnan(self.ocean_salinity), + self.ocean_salinity=35.0*ones(md.mesh.numberofvertices,1); + disp(' no basalforcings.ocean_salinity specified: values set as 35 psu'); + end + + + end % }}} + function self = setdefaultparameters(self) % {{{ + + %default values for melting parameterization + self.meltrate_factor = 0.5; + + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + if ismember('MasstransportAnalysis',analyses) & ~(solution=='TransientSolution' & md.transient.ismasstransport==0), + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','basalforcings.ocean_temp','NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','basalforcings.ocean_salinity','NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'numel',1); + end + if ismember('BalancethicknessAnalysis',analyses), + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','basalforcings.ocean_temp','NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','basalforcings.ocean_salinity','NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'numel',1); + end + if ismember('ThermalAnalysis',analyses) & ~(solution=='TransientSolution' & md.transient.isthermal==0), + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'numel',1); + md = checkfield(md,'fieldname','basalforcings.geothermalflux','NaN',1,'Inf',1,'timeseries',1,'>=',0); + end + end % }}} + function disp(self) % {{{ + disp(sprintf(' Beckmann & Goosse (2003) basal melt parameterization:')); + + fielddisplay(self,'groundedice_melting_rate','basal melting rate (positive if melting) (m/yr)'); + fielddisplay(self,'geothermalflux','geothermal heat flux (W/m^2)'); + fielddisplay(self,'meltrate_factor','Melt-rate rate factor'); + fielddisplay(self,'ocean_temp','ocean temperature (degC)'); + fielddisplay(self,'ocean_salinity','ocean ocean_salinity (psu)'); + + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + yts=365.2422*24.0*3600.0; + + floatingice_melting_rate=zeros(md.mesh.numberofvertices,1); + T_f=(0.0939 - 0.057 * md.basalforcings.ocean_salinity + 7.64e-4 * md.geometry.base); + ocean_heat_flux = 1.5 * md.materials.rho_water * md.materials.mixed_layer_capacity * md.materials.thermal_exchange_velocity * (md.basalforcings.ocean_temp - T_f); + floatingice_melting_rate=ocean_heat_flux/(md.materials.latentheat*md.materials.rho_ice); + + +WriteData(fid,prefix,'name','md.basalforcings.model','data',8,'format','Integer'); + WriteData(fid,prefix,'data',floatingice_melting_rate,'format','DoubleMat','name','md.basalforcings.floatingice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1) + WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','name','md.basalforcings.groundedice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1) + WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','name','md.basalforcings.geothermalflux','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1); + WriteData(fid,prefix,'object',self,'fieldname','meltrate_factor','format','Double','name','md.basalforcings.meltrate_factor'); + WriteData(fid,prefix,'object',self,'fieldname','ocean_temp','format','DoubleMat','name','md.basalforcings.ocean_temp','mattype',1,'timeserieslength',md.mesh.numberofvertices+1); + WriteData(fid,prefix,'object',self,'fieldname','ocean_salinity','format','DoubleMat','name','md.basalforcings.ocean_salinity','mattype',1,'timeserieslength',md.mesh.numberofvertices+1); + + end % }}} + end +end Index: /issm/trunk-jpl/test/NightlyRun/test476.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test476.m (revision 24145) +++ /issm/trunk-jpl/test/NightlyRun/test476.m (revision 24145) @@ -0,0 +1,74 @@ +%Test Name: BeckmannGoosseMeltRate_HO +md=triangle(model(),'../Exp/Square.exp',90000.); +md=setmask(md,'../Exp/SquareShelf.exp',''); +md=parameterize(md,'../Par/SquareSheetShelf.par'); +md.initialization.vx(:)=1.; +md.initialization.vy(:)=1.; +md.geometry.thickness(:)=500-md.mesh.x/10000; +md.geometry.bed =-100-md.mesh.x/1000; +md.geometry.base=-md.geometry.thickness*md.materials.rho_ice/md.materials.rho_water; +md.mask.groundedice_levelset=md.geometry.thickness+md.materials.rho_water/md.materials.rho_ice*md.geometry.bed; +pos=find(md.mask.groundedice_levelset>=0); +md.geometry.base(pos)=md.geometry.bed(pos); +md.geometry.surface=md.geometry.base+md.geometry.thickness; +md=extrude(md,3,1.1); +md=setflowequation(md,'HO','all'); + +%Set Pico Parameters +md.basalforcings=basalforcingsbeckmanngoosse(md.basalforcings); +md.basalforcings.ocean_temp=-1.7*ones(md.mesh.numberofvertices,1); +md.basalforcings.ocean_salinity=35.0*ones(md.mesh.numberofvertices,1); +md.basalforcings.meltrate_factor=1; + +%Boundary conditions: +md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1); +md.mask.ice_levelset(find(md.mesh.x==max(md.mesh.x)))=0; + +%Model conditions +md.transient.isthermal=0; +md.transient.isstressbalance=1; +md.transient.isgroundingline=1; +md.transient.ismasstransport=1; +md.transient.issmb=1; +md.transient.requested_outputs={'default','BasalforcingsFloatingiceMeltingRate'}; +md.groundingline.migration='SubelementMigration'; +md.groundingline.friction_interpolation='SubelementFriction1'; +md.groundingline.melt_interpolation='SubelementMelt1'; +md.timestepping.final_time=1.5; +md.timestepping.time_step=0.5; + +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,'Transient'); + +field_names ={'Bed1','Surface1','Thickness1','Floatingice1','Vx1','Vy1','Pressure1','FloatingiceMeltingrate1',... + 'Bed2','Surface2','Thickness2','Floatingice2','Vx2','Vy2','Pressure2','FloatingiceMeltingrate2',... + 'Bed3','Surface3','Thickness3','Floatingice3','Vx3','Vy3','Pressure3','FloatingiceMeltingrate3'}; +field_tolerances={7e-09,8e-09,8e-09,7e-09,6e-08,7e-08,6e-09,8e-7,... + 7e-09,8e-09,8e-09,7e-09,6e-08,7e-08,6e-09,8e-7,... + 7e-09,8e-09,8e-09,7e-09,6e-08,7e-08,6e-09,8e-7}; +field_values={... + (md.results.TransientSolution(1).Base),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).MaskGroundediceLevelset),... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Pressure),... + (md.results.TransientSolution(1).BasalforcingsFloatingiceMeltingRate),... + (md.results.TransientSolution(2).Base),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(2).MaskGroundediceLevelset),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Pressure),... + (md.results.TransientSolution(2).BasalforcingsFloatingiceMeltingRate),... + (md.results.TransientSolution(3).Base),... + (md.results.TransientSolution(3).Surface),... + (md.results.TransientSolution(3).Thickness),... + (md.results.TransientSolution(3).MaskGroundediceLevelset),... + (md.results.TransientSolution(3).Vx),... + (md.results.TransientSolution(3).Vy),... + (md.results.TransientSolution(3).Pressure),... + (md.results.TransientSolution(3).BasalforcingsFloatingiceMeltingRate),... + };