Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 28114) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 28115) @@ -354,6 +354,6 @@ /*Intermediaries*/ int domaintype,dim,stabilization; - IssmDouble Jdet,dt; - IssmDouble gmb,fmb,mb,bed,vx,vy,vz,tau; + IssmDouble Jdet,dt,intrusiondist; + IssmDouble gmb,fmb,mb,bed,vx,vy,vz,tau,gldistance; Element* basalelement = NULL; IssmDouble *xyz_list = NULL; @@ -361,4 +361,5 @@ /*Get basal element*/ element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ case Domain2DhorizontalEnum: @@ -394,4 +395,6 @@ basalelement->FindParam(&dt,TimesteppingTimeStepEnum); basalelement->FindParam(&melt_style,GroundinglineMeltInterpolationEnum); + basalelement->FindParam(&intrusiondist,GroundinglineIntrusionDistanceEnum); + Input* groundedice_input = basalelement->GetInput(MaskOceanLevelsetEnum); _assert_(groundedice_input); Input* gmb_input = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input); @@ -424,4 +427,8 @@ gauss = basalelement->NewGauss(point1,fraction1,fraction2,3); } + if(melt_style==IntrusionMeltEnum){ + basalelement->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + gauss = basalelement->NewGauss(point1,fraction1,fraction2,3); + } else{ gauss = basalelement->NewGauss(3); @@ -459,4 +466,17 @@ else mb=gmb; } + else if(melt_style==IntrusionMeltEnum){ + Input* gldistance_input = basalelement->GetInput(DistanceToGroundinglineEnum); _assert_(gldistance_input); + gldistance_input->GetInputValue(&gldistance,gauss); + if (intrusiondist==0) + if(gllevelset>0.) mb=gmb; + else mb=fmb; + else if(gldistance>intrusiondist) + mb=gmb; + else if(gldistance<=intrusiondist && gldistance>0) + mb=fmb*(1-gldistance/intrusiondist); + else + mb=fmb; + } else _error_("melt interpolation "<FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum); + if(isgroundingline) iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); /*Initialize ThicknessResidual input*/ @@ -272,4 +273,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.masstransport.min_thickness",MasstransportMinThicknessEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.masstransport.penalty_factor",MasstransportPenaltyFactorEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.groundingline.intrusion_distance",GroundinglineIntrusionDistanceEnum)); iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.masstransport.requested_outputs"); @@ -622,6 +624,6 @@ bool mainlyfloating; IssmDouble fraction1,fraction2; - IssmDouble Jdet,dt; - IssmDouble ms,mb,gmb,fmb,thickness,fmb_pert; + IssmDouble Jdet,dt,intrusiondist; + IssmDouble ms,mb,gmb,fmb,thickness,fmb_pert,gldistance; IssmDouble vx,vy,vel,dvxdx,dvydy,xi,h,tau; IssmDouble dvx[2],dvy[2]; @@ -652,4 +654,6 @@ element->FindParam(&dt,TimesteppingTimeStepEnum); element->FindParam(&stabilization,MasstransportStabilizationEnum); + element->FindParam(&intrusiondist,GroundinglineIntrusionDistanceEnum); + Input* gmb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input); Input* fmb_input = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input); @@ -660,5 +664,5 @@ Input* vxaverage_input = element->GetInput(VxAverageEnum); _assert_(vxaverage_input); Input* vyaverage_input = element->GetInput(VyAverageEnum); _assert_(vyaverage_input); - + //Input* gldistance_input = element->GetInput(DistanceToGroundinglineEnum); _assert_(gldistance_input); h=element->CharacteristicLength(); @@ -667,5 +671,9 @@ if(melt_style==SubelementMelt2Enum){ element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); - gauss = element->NewGauss(point1,fraction1,fraction2,3); + gauss = element->NewGauss(point1,fraction1,fraction2,3); + } + else if(melt_style==IntrusionMeltEnum){ + element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + gauss = element->NewGauss(point1,fraction1,fraction2,3); } else{ @@ -704,4 +712,17 @@ else mb=gmb; } + else if(melt_style==IntrusionMeltEnum){ + Input* gldistance_input = element->GetInput(DistanceToGroundinglineEnum); _assert_(gldistance_input); + gldistance_input->GetInputValue(&gldistance,gauss); + if (intrusiondist==0) + if(gllevelset>0.) mb=gmb; + else mb=fmb; + else if(gldistance>intrusiondist) + mb=gmb; + else if(gldistance<=intrusiondist && gldistance>0) + mb=fmb*(1-gldistance/intrusiondist); + else + mb=fmb; + } else _error_("melt interpolation "<FindParam(&dt,TimesteppingTimeStepEnum); element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum); + element->FindParam(&intrusiondist,GroundinglineIntrusionDistanceEnum); + Input* gmb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input); Input* fmb_input = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input); @@ -778,4 +801,8 @@ gauss = element->NewGauss(point1,fraction1,fraction2,3); } + else if(melt_style==IntrusionMeltEnum){ + element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + gauss = element->NewGauss(point1,fraction1,fraction2,3); + } else{ gauss = element->NewGauss(3); @@ -798,17 +825,30 @@ else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating } - else if(melt_style==SubelementMelt2Enum){ + else if(melt_style==SubelementMelt2Enum){ if(gllevelset>0.) mb=gmb; else mb=fmb; } - else if(melt_style==NoMeltOnPartiallyFloatingEnum){ + else if(melt_style==NoMeltOnPartiallyFloatingEnum){ if (phi<0.00000001) mb=fmb; else mb=gmb; } - else if(melt_style==FullMeltOnPartiallyFloatingEnum){ + else if(melt_style==FullMeltOnPartiallyFloatingEnum){ if (phi<0.99999999) mb=fmb; else mb=gmb; - } - else _error_("melt interpolation "<GetInput(DistanceToGroundinglineEnum); _assert_(gldistance_input); + gldistance_input->GetInputValue(&gldistance,gauss); + if (intrusiondist==0) + if(gllevelset>0.) mb=gmb; + else mb=fmb; + else if(gldistance>intrusiondist) + mb=gmb; + else if(gldistance<=intrusiondist && gldistance>0) + mb=fmb*(1-gldistance/intrusiondist); + else + mb=fmb; + } + else _error_("melt interpolation "<values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i]; Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 28114) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 28115) @@ -7,9 +7,25 @@ #include "../../toolkits/toolkits.h" #include "./../../classes/Inputs/DatasetInput.h" +#include "../InputDuplicatex/InputDuplicatex.h" + void FloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/ /*Intermediaties*/ - int basalforcing_model; + int basalforcing_model; + int melt_style; + + /*First, get melt_interpolation model from parameters*/ + femmodel->parameters->FindParam(&melt_style,GroundinglineMeltInterpolationEnum); + + switch(melt_style){ + case IntrusionMeltEnum: + if(VerboseSolution())_printf0_(" call Intrusion melting module\n"); + DistanceToGroundingline(femmodel); + break; + default: + /*Nothing to be done*/ + break; + } /*First, get BMB model from parameters*/ @@ -288,2 +304,12 @@ xDelete(perturbation); }/*}}}*/ + +void DistanceToGroundingline(FemModel* femmodel){/*{{{*/ + + int numvertices; + IssmDouble* distances=NULL; + + InputDuplicatex(femmodel,MaskOceanLevelsetEnum,DistanceToGroundinglineEnum); + femmodel->DistanceToFieldValue(MaskOceanLevelsetEnum,0.,DistanceToGroundinglineEnum); +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h (revision 28114) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h (revision 28115) @@ -11,4 +11,5 @@ /* local prototypes: */ void FloatingiceMeltingRatex(FemModel* femmodel); + void LinearFloatingiceMeltingRatex(FemModel* femmodel); void SpatialLinearFloatingiceMeltingRatex(FemModel* femmodel); @@ -17,4 +18,5 @@ void BeckmannGoosseFloatingiceMeltingRatex(FemModel* femmodel); void LinearFloatingiceMeltingRatearmax(FemModel* femmodel); +void DistanceToGroundingline(FemModel* femmodel); #endif /* _FloatingiceMeltingRatex_H*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 28114) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 28115) @@ -76,4 +76,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.groundingline.friction_interpolation",GroundinglineFrictionInterpolationEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.groundingline.melt_interpolation",GroundinglineMeltInterpolationEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.groundingline.intrusion_distance",GroundinglineIntrusionDistanceEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.transient.isstressbalance",TransientIsstressbalanceEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.transient.ismasstransport",TransientIsmasstransportEnum)); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 28114) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 28115) @@ -868,4 +868,5 @@ syn keyword cConstant FrontalForcingsSubglacialDischargeEnum syn keyword cConstant GeometryHydrostaticRatioEnum +syn keyword cConstant GroundinglineIntrusionDistanceEnum syn keyword cConstant NGiaEnum syn keyword cConstant NGiaRateEnum @@ -925,4 +926,5 @@ syn keyword cConstant InversionVxObsEnum syn keyword cConstant InversionVyObsEnum +syn keyword cConstant IntrusionMeltEnum syn keyword cConstant LevelsetfunctionSlopeXEnum syn keyword cConstant LevelsetfunctionSlopeYEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 28114) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 28115) @@ -240,4 +240,5 @@ GrdModelEnum, GroundinglineFrictionInterpolationEnum, + GroundinglineIntrusionDistanceEnum, GroundinglineMeltInterpolationEnum, GroundinglineMigrationEnum, @@ -1739,4 +1740,5 @@ SubelementFriction1Enum, SubelementFriction2Enum, + IntrusionMeltEnum, SubelementMelt1Enum, SubelementMelt2Enum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 28114) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 28115) @@ -870,4 +870,5 @@ case FrontalForcingsSubglacialDischargeEnum : return "FrontalForcingsSubglacialDischarge"; case GeometryHydrostaticRatioEnum : return "GeometryHydrostaticRatio"; + case GroundinglineIntrusionDistanceEnum : return "GroundinglineIntrusionDistance"; case NGiaEnum : return "NGia"; case NGiaRateEnum : return "NGiaRate"; @@ -927,4 +928,5 @@ case InversionVxObsEnum : return "InversionVxObs"; case InversionVyObsEnum : return "InversionVyObs"; + case IntrusionMeltEnum : return "IntrusionMelt"; case LevelsetfunctionSlopeXEnum : return "LevelsetfunctionSlopeX"; case LevelsetfunctionSlopeYEnum : return "LevelsetfunctionSlopeY"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 28114) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 28115) @@ -861,4 +861,5 @@ syn keyword juliaConstC FrontalForcingsSubglacialDischargeEnum syn keyword juliaConstC GeometryHydrostaticRatioEnum +syn keyword juliaConstC GroundinglineIntrusionDistanceEnum syn keyword juliaConstC NGiaEnum syn keyword juliaConstC NGiaRateEnum @@ -918,4 +919,5 @@ syn keyword juliaConstC InversionVxObsEnum syn keyword juliaConstC InversionVyObsEnum +syn keyword juliaConstC IntrusionMeltEnum syn keyword juliaConstC LevelsetfunctionSlopeXEnum syn keyword juliaConstC LevelsetfunctionSlopeYEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 28114) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 28115) @@ -891,4 +891,5 @@ else if (strcmp(name,"FrontalForcingsSubglacialDischarge")==0) return FrontalForcingsSubglacialDischargeEnum; else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum; + else if (strcmp(name,"GroundinglineIntrusionDistance")==0) return GroundinglineIntrusionDistanceEnum; else if (strcmp(name,"NGia")==0) return NGiaEnum; else if (strcmp(name,"NGiaRate")==0) return NGiaRateEnum; @@ -948,4 +949,5 @@ else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; + else if (strcmp(name,"IntrusionMelt")==0) return IntrusionMeltEnum; else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; @@ -996,10 +998,10 @@ else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; else if (strcmp(name,"Sample")==0) return SampleEnum; - else if (strcmp(name,"SampleOld")==0) return SampleOldEnum; - else if (strcmp(name,"SampleNoise")==0) return SampleNoiseEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SamplingBeta")==0) return SamplingBetaEnum; + if (strcmp(name,"SampleOld")==0) return SampleOldEnum; + else if (strcmp(name,"SampleNoise")==0) return SampleNoiseEnum; + else if (strcmp(name,"SamplingBeta")==0) return SamplingBetaEnum; else if (strcmp(name,"SamplingKappa")==0) return SamplingKappaEnum; else if (strcmp(name,"SamplingPhi")==0) return SamplingPhiEnum; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum; else if (strcmp(name,"SmbDz")==0) return SmbDzEnum; - else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum; - else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum; + if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum; + else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum; + else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum; else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum; else if (strcmp(name,"SmbEC")==0) return SmbECEnum; @@ -1242,10 +1244,10 @@ else if (strcmp(name,"SubglacialdischargeValuesMovingaverage")==0) return SubglacialdischargeValuesMovingaverageEnum; else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum; - else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum; - else if (strcmp(name,"Area")==0) return AreaEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"SealevelArea")==0) return SealevelAreaEnum; + if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum; + else if (strcmp(name,"Area")==0) return AreaEnum; + else if (strcmp(name,"SealevelArea")==0) return SealevelAreaEnum; else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum; else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; @@ -1365,10 +1367,10 @@ 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; else stage=12; } if(stage==12){ - if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum; + if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum; + else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum; + else if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum; else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum; else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum; @@ -1488,10 +1490,10 @@ else if (strcmp(name,"Cuffey")==0) return CuffeyEnum; else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum; - else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum; - else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"DataSet")==0) return DataSetEnum; + if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum; + else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum; + else if (strcmp(name,"DataSet")==0) return DataSetEnum; else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; @@ -1611,10 +1613,10 @@ else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; - else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum; - else if (strcmp(name,"LinearFloatingMeltRatearma")==0) return LinearFloatingMeltRatearmaEnum; else stage=14; } if(stage==14){ - if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; + if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum; + else if (strcmp(name,"LinearFloatingMeltRatearma")==0) return LinearFloatingMeltRatearmaEnum; + else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; else if (strcmp(name,"Loads")==0) return LoadsEnum; else if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; @@ -1734,10 +1736,10 @@ else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; - else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; - else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; else stage=15; } if(stage==15){ - if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; + if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; + else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; + else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; Index: /issm/trunk-jpl/src/m/classes/groundingline.m =================================================================== --- /issm/trunk-jpl/src/m/classes/groundingline.m (revision 28114) +++ /issm/trunk-jpl/src/m/classes/groundingline.m (revision 28115) @@ -9,4 +9,5 @@ friction_interpolation = ''; melt_interpolation = ''; + intrusion_distance = 0; requested_outputs = {}; end @@ -26,6 +27,7 @@ self.friction_interpolation= 'SubelementFriction1'; self.melt_interpolation = 'NoMeltOnPartiallyFloating'; + self.intrusion_distance = 0; %default output - self.requested_outputs = {'default'}; + self.requested_outputs = {'default'}; end % }}} @@ -34,5 +36,6 @@ md = checkfield(md,'fieldname','groundingline.migration','values',{'None' 'SubelementMigration' 'AggressiveMigration' 'SoftMigration' 'Contact' 'GroundingOnly'}); md = checkfield(md,'fieldname','groundingline.friction_interpolation','values',{'NoFrictionOnPartiallyFloating' 'SubelementFriction1' 'SubelementFriction2'}); - md = checkfield(md,'fieldname','groundingline.melt_interpolation','values',{'NoMeltOnPartiallyFloating' 'SubelementMelt1' 'SubelementMelt2' 'FullMeltOnPartiallyFloating'}); + md = checkfield(md,'fieldname','groundingline.melt_interpolation','values',{'NoMeltOnPartiallyFloating' 'SubelementMelt1' 'SubelementMelt2' 'IntrusionMelt' 'FullMeltOnPartiallyFloating'}); + md = checkfield(md,'fieldname','groundingline.intrusion_distance','>=',0); md = checkfield(md,'fieldname','groundingline.requested_outputs','stringrow',1); @@ -61,5 +64,6 @@ fielddisplay(self,'migration','type of grounding line migration: ''SoftMigration'',''SubelementMigration'',''AggressiveMigration'',''Contact'' or ''None'''); fielddisplay(self,'friction_interpolation','type of friction interpolation for partially floating elements: ''NoFrictionOnPartiallyFloating'',''SubelementFriction1'', or ''SubelementFriction2'''); - fielddisplay(self,'melt_interpolation','type of melt interpolation for partially floating elements: ''NoMeltOnPartiallyFloating'',''SubelementMelt1'',''SubelementMelt2'', or ''FullMeltOnPartiallyFloating'''); + fielddisplay(self,'melt_interpolation','type of melt interpolation for partially floating elements: ''NoMeltOnPartiallyFloating'',''SubelementMelt1'',''SubelementMelt2'',''IntrusionMelt'' or ''FullMeltOnPartiallyFloating'''); + fielddisplay(self,'intrusion_distance','distance of seawater intrusion from grounding line [m]'); fielddisplay(self,'requested_outputs','additional outputs requested'); @@ -69,4 +73,5 @@ WriteData(fid,prefix,'data',self.friction_interpolation,'name','md.groundingline.friction_interpolation','format','String'); WriteData(fid,prefix,'data',self.melt_interpolation,'name','md.groundingline.melt_interpolation','format','String'); + WriteData(fid,prefix,'object',self,'fieldname','intrusion_distance','format','Double'); %process requested outputs Index: /issm/trunk-jpl/test/NightlyRun/test456.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test456.m (revision 28115) +++ /issm/trunk-jpl/test/NightlyRun/test456.m (revision 28115) @@ -0,0 +1,84 @@ +%Test Name: MISMIP3D-IntrusionMelt +md=triangle(model(),'../Exp/Square.exp',100000.); +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.ocean_levelset=md.geometry.thickness+md.materials.rho_water/md.materials.rho_ice*md.geometry.bed; +pos=find(md.mask.ocean_levelset>=0); +md.geometry.base(pos)=md.geometry.bed(pos); +md.geometry.surface=md.geometry.base+md.geometry.thickness; +md=setflowequation(md,'SSA','all'); + +%Boundary conditions: +md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1); +md.mask.ice_levelset(find(md.mesh.x==max(md.mesh.x)))=0; +md.stressbalance.spcvx(:)=NaN; +md.stressbalance.spcvy(:)=NaN; +md.stressbalance.spcvz(:)=NaN; +pos=find((md.mesh.y<1000000.1 & md.mesh.y>999999.9) | (md.mesh.y<0.1 & md.mesh.y>-0.1)); +md.stressbalance.spcvy(pos)=0; +pos2=find(md.mesh.x<0.1 & md.mesh.x>-0.1); +md.stressbalance.spcvx(pos2)=0; +md.stressbalance.spcvy(pos2)=0; + +md.materials.rheology_B=1/((10^-25)^(1/3))*ones(md.mesh.numberofvertices,1); +md.materials.rheology_law='None'; +md.friction.coefficient(:)=sqrt(10^7)*ones(md.mesh.numberofvertices,1); +md.friction.p=3*ones(md.mesh.numberofelements,1); +md.smb.mass_balance(:)=1; +md.basalforcings.groundedice_melting_rate(:)=0; +md.basalforcings.floatingice_melting_rate(:)=30; +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='SubelementFriction2'; +md.groundingline.melt_interpolation='IntrusionMelt'; +md.groundingline.intrusion_distance = 5000; + +md.timestepping.final_time=30; +md.timestepping.time_step=10; + +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,'Transient'); + +%Fields and tolerances to track changes +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={2e-11,5e-12,2e-11,1e-11,5e-10,1e-08,1e-13,1e-13,... + 3e-11,3e-11,9e-10,7e-11,1e-09,5e-08,1e-10,1e-13,... + 1e-10,3e-11,1e-10,7e-11,1e-09,5e-08,1e-10,1e-13}; +field_values={... + (md.results.TransientSolution(1).Base),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).MaskOceanLevelset),... + (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).MaskOceanLevelset),... + (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).MaskOceanLevelset),... + (md.results.TransientSolution(3).Vx),... + (md.results.TransientSolution(3).Vy),... + (md.results.TransientSolution(3).Pressure),... + (md.results.TransientSolution(3).BasalforcingsFloatingiceMeltingRate),... + }; Index: /issm/trunk-jpl/test/NightlyRun/test520.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test520.m (revision 28115) +++ /issm/trunk-jpl/test/NightlyRun/test520.m (revision 28115) @@ -0,0 +1,59 @@ +%Test Name: PigTran-IntrusionMelt +md=triangle(model(),'../Exp/Pig.exp',20000.); +md=setmask(md,'../Exp/PigShelves.exp','../Exp/PigIslands.exp'); +md=parameterize(md,'../Par/Pig.par'); +md=setflowequation(md,'SSA','all'); +md.cluster=generic('name',oshostname(),'np',3); +md.geometry.bed=md.geometry.base; +[md.mesh.lat,md.mesh.long] = xy2ll(md.mesh.x,md.mesh.y,-1); +md.mesh.scale_factor=0.9*ones(md.mesh.numberofvertices,1); +md.transient.requested_outputs={'default','IceVolume','IceVolumeScaled','GroundedArea','GroundedAreaScaled','FloatingArea','FloatingAreaScaled','TotalSmb','TotalSmbScaled','TotalFloatingBmb','TotalFloatingBmbScaled'}; +md.transient.isgroundingline=1; +md.basalforcings.floatingice_melting_rate=50*ones(md.mesh.numberofvertices,1); + +md.groundingline.melt_interpolation = 'IntrusionMelt'; +md.groundingline.intrusion_distance = 5000; + +md=solve(md,'Transient'); + +%Fields and tolerances to track changes +field_names ={'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','IceVolume1','IceVolumeScaled1','GroundedArea1','GroundedAreaScaled1','FloatingArea1','FloatingAreaScaled1','TotalSmb1','TotalSmbScaled1','TotalFloatingBmb1','TotalFloatingBmbScaled1',... + 'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','IceVolume2','IceVolumeScaled2','GroundedArea2','GroundedAreaScaled2','FloatingArea2','FloatingAreaScaled2','TotalSmb2','TotalSmbScaled2','TotalFloatingBmb2','TotalFloatingBmbScaled2'}; +field_tolerances={1e-12,2e-12,2e-12,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,... + 1e-12,1e-12,1e-12,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Pressure),... + (md.results.TransientSolution(1).Base),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).IceVolume),... + (md.results.TransientSolution(1).IceVolumeScaled),... + (md.results.TransientSolution(1).GroundedArea),... + (md.results.TransientSolution(1).GroundedAreaScaled),... + (md.results.TransientSolution(1).FloatingArea),... + (md.results.TransientSolution(1).FloatingAreaScaled),... + (md.results.TransientSolution(1).TotalSmb),... + (md.results.TransientSolution(1).TotalSmbScaled),... + (md.results.TransientSolution(1).TotalFloatingBmb),... + (md.results.TransientSolution(1).TotalFloatingBmbScaled),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Vel),... + (md.results.TransientSolution(2).Pressure),... + (md.results.TransientSolution(2).Base),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(2).IceVolume),... + (md.results.TransientSolution(2).IceVolumeScaled),... + (md.results.TransientSolution(2).GroundedArea),... + (md.results.TransientSolution(2).GroundedAreaScaled),... + (md.results.TransientSolution(2).FloatingArea),... + (md.results.TransientSolution(2).FloatingAreaScaled),... + (md.results.TransientSolution(2).TotalSmb),... + (md.results.TransientSolution(2).TotalSmbScaled),... + (md.results.TransientSolution(2).TotalFloatingBmb),... + (md.results.TransientSolution(2).TotalFloatingBmbScaled),... + };