Index: /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 22917) +++ /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 22918) @@ -41,5 +41,6 @@ iomodel->FetchDataToInput(elements,"md.geometry.surface",SurfaceEnum); - iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum); + iomodel->FetchDataToInput(elements,"md.geometry.bed",BedEnum); + iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum); iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0); iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); @@ -165,6 +166,8 @@ case SurfaceSlopeXEnum: input2 = basalelement->GetInput(SurfaceEnum); _assert_(input2); break; case SurfaceSlopeYEnum: input2 = basalelement->GetInput(SurfaceEnum); _assert_(input2); break; - case BedSlopeXEnum: input2 = basalelement->GetInput(BaseEnum); _assert_(input2); break; - case BedSlopeYEnum: input2 = basalelement->GetInput(BaseEnum); _assert_(input2); break; + case BedSlopeXEnum: input2 = basalelement->GetInput(BedEnum); _assert_(input2); break; + case BedSlopeYEnum: input2 = basalelement->GetInput(BedEnum); _assert_(input2); break; + case BaseSlopeXEnum: input2 = basalelement->GetInput(BaseEnum); _assert_(input2); break; + case BaseSlopeYEnum: input2 = basalelement->GetInput(BaseEnum); _assert_(input2); break; case LevelsetfunctionSlopeXEnum: input2 = basalelement->GetInput(MaskIceLevelsetEnum); _assert_(input2); break; case LevelsetfunctionSlopeYEnum: input2 = basalelement->GetInput(MaskIceLevelsetEnum); _assert_(input2); break; @@ -182,6 +185,6 @@ if(input2) input2->GetInputDerivativeValue(&slopes[0],xyz_list,gauss); switch(input_enum){ - case SurfaceSlopeXEnum: case BedSlopeXEnum: case LevelsetfunctionSlopeXEnum: value = slopes[0]; break; - case SurfaceSlopeYEnum: case BedSlopeYEnum: case LevelsetfunctionSlopeYEnum: value = slopes[1]; break; + case SurfaceSlopeXEnum: case BedSlopeXEnum: case BaseSlopeXEnum: case LevelsetfunctionSlopeXEnum: value = slopes[0]; break; + case SurfaceSlopeYEnum: case BedSlopeYEnum: case BaseSlopeYEnum: case LevelsetfunctionSlopeYEnum: value = slopes[1]; break; default: input->GetInputValue(&value,gauss); } Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22917) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22918) @@ -281,4 +281,7 @@ virtual void PicoUpdateBoxid(int* pmax_boxid_basin)=0; virtual void PicoUpdateBox(int loopboxid)=0; + virtual void PicoComputeBaseSlope(void)=0; + virtual void PicoComputeGroundingLineDepth(void)=0; + virtual void PicoComputeBasalMelt(void)=0; virtual void PotentialUngrounding(Vector* potential_sheet_ungrounding)=0; virtual int PressureInterpolation()=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22917) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22918) @@ -143,4 +143,7 @@ void PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");}; void PicoUpdateBox(int loopboxid){_error_("not implemented yet");}; + void PicoComputeBaseSlope(void){_error_("not implemented yet");}; + void PicoComputeGroundingLineDepth(void){_error_("not implemented yet");}; + void PicoComputeBasalMelt(void){_error_("not implemented yet");}; void PotentialUngrounding(Vector* potential_sheet_ungrounding); int PressureInterpolation(); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 22917) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 22918) @@ -132,4 +132,7 @@ void PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");}; void PicoUpdateBox(int loopboxid){_error_("not implemented yet");}; + void PicoComputeBaseSlope(void){_error_("not implemented yet");}; + void PicoComputeGroundingLineDepth(void){_error_("not implemented yet");}; + void PicoComputeBasalMelt(void){_error_("not implemented yet");}; void PotentialUngrounding(Vector* potential_sheet_ungrounding){_error_("not implemented yet");}; int PressureInterpolation(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 22917) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 22918) @@ -140,4 +140,7 @@ void PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");}; void PicoUpdateBox(int loopboxid){_error_("not implemented yet");}; + void PicoComputeBaseSlope(void){_error_("not implemented yet");}; + void PicoComputeGroundingLineDepth(void){_error_("not implemented yet");}; + void PicoComputeBasalMelt(void){_error_("not implemented yet");}; void PotentialUngrounding(Vector* potential_sheet_ungrounding){_error_("not implemented yet");}; int PressureInterpolation(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22917) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22918) @@ -2937,4 +2937,5 @@ IssmDouble gamma_T, overturning_coeff, thickness; IssmDouble pressure, T_star,p_coeff, q_coeff; + bool isplume; /*Get variables*/ @@ -2959,4 +2960,5 @@ this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); + this->parameters->FindParam(&isplume, BasalforcingsPicoIsplumeEnum); Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); _assert_(basinid<=num_basins); @@ -3001,9 +3003,9 @@ Socs[i] = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Tocs[i]); potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure; - basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]); + if(!isplume){basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);} overturnings[i] = overturning_coeff*rho_star*(Beta*(soc_farocean-Socs[i])-alpha*(toc_farocean-Tocs[i])); } - this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); + if(!isplume){this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);} this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); @@ -3038,8 +3040,8 @@ Socs[i] = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda)); potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure; - basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]); - } - - this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); + if(!isplume){basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);} + } + + if(!isplume){this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);} this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); @@ -3056,4 +3058,142 @@ } /*}}}*/ +void Tria::PicoComputeBaseSlope(void){/*{{{*/ + + //if(!this->IsIceInElement() || !this->IsFloating()) return; + + /*Define variables and get inputs */ + IssmDouble slopex,slopey; + IssmDouble slope[NUMVERTICES]; + Input* baseslopex_input = this->GetInput(BaseSlopeXEnum); _assert_(baseslopex_input); + Input* baseslopey_input = this->GetInput(BaseSlopeYEnum); _assert_(baseslopey_input); + + Gauss* gauss=this->NewGauss(); + for(int i=0;iGaussVertex(i); + baseslopex_input->GetInputValue(&slopex,gauss); + baseslopey_input->GetInputValue(&slopey,gauss); + slope[i] = sqrt(slopex*slopex + slopey*slopey); + } + + this->AddInput(BasalforcingsPicoShelfSlopeEnum,slope,P1Enum); + + /*Cleanup and return*/ + delete gauss; +}/*}}}*/ +void Tria::PicoComputeGroundingLineDepth(void){/*{{{*/ + + if(!this->IsIceInElement() || !this->IsFloating()) return; + + /*Define variables and get inputs */ + IssmDouble bed; + IssmDouble zgl[NUMVERTICES]; + Input* bed_input = this->GetInput(BedEnum); _assert_(bed_input); + + Gauss* gauss=this->NewGauss(); + for(int i=0;iGaussVertex(i); + bed_input->GetInputValue(&bed,gauss); + zgl[i] = bed -50.; + } + + this->AddInput(BasalforcingsPicoGroundingLineDepthEnum,zgl,P1Enum); + + /*Cleanup and return*/ + delete gauss; +}/*}}}*/ +void Tria::PicoComputeBasalMelt(void){/*{{{*/ + + if(!this->IsIceInElement() || !this->IsFloating()) return; + + IssmDouble E0, Cd, CdT, YT, lam1, lam2, lam3, M0, CdTS0, y1, y2, x0; + IssmDouble Tf_gl, YTS, CdTS, G1, G2, G3, g_alpha, M, l, X_hat, M_hat; + IssmDouble alpha, zgl, Toc, Soc, z_base, yts; + + /*Get variables*/ + E0 = 3.6e-2; //Entrainment coefficient + Cd = 2.5e-3; //Drag coefficient + CdT = 1.1e-3; //Turbulent heat exchange coefficient + YT = CdT/sqrt(Cd); //Heat exchange coefficient + lam1 = -5.73e-2; //Freezing point-salinity coefficient (degrees C) + lam2 = 8.32e-2; //Freezing point offset (degrees C) + lam3 = 7.61e-4; //Freezing point-depth coefficient (K m-1) + M0 = 10.; //Melt-rate parameter (m yr-1 C-2) + CdTS0 = 6e-4; //Heat exchange parameter + y1 = 0.545; //Heat exchange parameter + y2 = 3.5e-5; //Heat exchange parameter + x0 = 0.56; //Dimentionless scaling factor + + /*Define arrays*/ + IssmDouble basalmeltrates_shelf[NUMVERTICES]; //Basal melt-rate + + /*Polynomial coefficients*/ + IssmDouble p[12]; + p[0] = 0.1371330075095435; + p[1] = 5.527656234709359E1; + p[2] = -8.951812433987858E2; + p[3] = 8.927093637594877E3; + p[4] = -5.563863123811898E4; + p[5] = 2.218596970948727E5; + p[6] = -5.820015295669482E5; + p[7] = 1.015475347943186E6; + p[8] = -1.166290429178556E6; + p[9] = 8.466870335320488E5; + p[10] = -3.520598035764990E5; + p[11] = 6.387953795485420E4; + + /*Get inputs*/ + Input* alpha_input=this->GetInput(BasalforcingsPicoShelfSlopeEnum); _assert_(alpha_input); + Input* zgl_input=this->GetInput(BasalforcingsPicoGroundingLineDepthEnum); _assert_(zgl_input); + Input* toc_input=this->GetInput(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(toc_input); + Input* soc_input=this->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(soc_input); + Input* base_input=this->GetInput(BaseEnum); _assert_(base_input); + this->FindParam(&yts, ConstantsYtsEnum); + + /*Loop over the number of vertices in this element*/ + Gauss* gauss=this->NewGauss(); + for(int i=0;iGaussVertex(i); + + /*Get inputs*/ + alpha_input->GetInputValue(&alpha,gauss); + zgl_input->GetInputValue(&zgl,gauss); + toc_input->GetInputValue(&Toc,gauss); //(K) + soc_input->GetInputValue(&Soc,gauss); + base_input->GetInputValue(&z_base,gauss); + + /*Make necessary conversions*/ + Toc = Toc-273.15; // (Degrees C) + alpha = atan(alpha); // (Degrees) + + /*Low bound for Toc to ensure X_hat is between 0 and 1*/ + if(TocAddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); + + /*Cleanup and return*/ + delete gauss; +}/*}}}*/ void Tria::PotentialUngrounding(Vector* potential_ungrounding){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22917) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22918) @@ -114,4 +114,7 @@ void PicoUpdateBoxid(int* pmax_boxid_basin); void PicoUpdateBox(int loopboxid); + void PicoComputeBaseSlope(void); + void PicoComputeGroundingLineDepth(void); + void PicoComputeBasalMelt(void); void PotentialUngrounding(Vector* potential_sheet_ungrounding); int PressureInterpolation(); Index: /issm/trunk-jpl/src/c/cores/bmb_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/bmb_core.cpp (revision 22917) +++ /issm/trunk-jpl/src/c/cores/bmb_core.cpp (revision 22918) @@ -20,12 +20,10 @@ /*In some cases we need to run additional analyses to get the required input data*/ if(basalforcing_model==BasalforcingsPicoEnum){ - //femmodel->parameters->FindParam(&isplume,BasalforcingsIsplumeEnum); + femmodel->parameters->FindParam(&isplume,BasalforcingsPicoIsplumeEnum); if(isplume){ - femmodel->SetCurrentConfiguration(L2ProjectionBaseAnalysisEnum); - _error_("STOP"); - //femmodel->parameters->SetParam(BaseSlopeXEnum,InputToL2ProjectEnum); + femmodel->parameters->SetParam(BaseSlopeXEnum,InputToL2ProjectEnum); solutionsequence_linear(femmodel); - //femmodel->parameters->SetParam(BaseSlopeYEnum,InputToL2ProjectEnum); + femmodel->parameters->SetParam(BaseSlopeYEnum,InputToL2ProjectEnum); solutionsequence_linear(femmodel); } Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp (revision 22917) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp (revision 22918) @@ -10,4 +10,5 @@ int maxbox; + bool isplume; /*First, reset all melt to 0 */ @@ -20,4 +21,5 @@ } + /*PICO nelt rate parameterization (Reese et al., 2018)*/ femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); UpdateBoxIdsPico(femmodel); @@ -26,4 +28,12 @@ UpdateBoxPico(femmodel,i); ComputeAverageOceanvarsPico(femmodel,i); + } + + /*Optional buoyant plume melt rate parameterization (Lazeroms et al., 2018) */ + femmodel->parameters->FindParam(&isplume,BasalforcingsPicoIsplumeEnum); + if(isplume){ + ComputeSubshelfSlopePlume(femmodel); + ComputeGroundingLineDepthPlume(femmodel); + ComputeBasalMeltPlume(femmodel); } }/*}}}*/ @@ -222,2 +232,20 @@ xDelete(boxareas); }/*}}}*/ +void ComputeSubshelfSlopePlume(FemModel* femmodel){/*{{{*/ + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->PicoComputeBaseSlope(); + } +}/*}}}*/ +void ComputeGroundingLineDepthPlume(FemModel* femmodel){/*{{{*/ + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->PicoComputeGroundingLineDepth(); + } +}/*}}}*/ +void ComputeBasalMeltPlume(FemModel* femmodel){/*{{{*/ + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->PicoComputeBasalMelt(); + } +}/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.h =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.h (revision 22917) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.h (revision 22918) @@ -15,4 +15,7 @@ void UpdateBoxPico(FemModel* femmodel, int loopboxid); void ComputeAverageOceanvarsPico(FemModel* femmodel, int boxid); +void ComputeSubshelfSlopePlume(FemModel* femmodel); +void ComputeGroundingLineDepthPlume(FemModel* femmodel); +void ComputeBasalMeltPlume(FemModel* femmodel); #endif /* _FloatingiceMeltingRatePicox_H*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22917) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22918) @@ -233,4 +233,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.overturning_coeff",BasalforcingsPicoOverturningCoeffEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_T",BasalforcingsPicoGammaTEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.isplume",BasalforcingsPicoIsplumeEnum)); iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_temperature"); _assert_(M>=1 && N>=1); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22917) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22918) @@ -73,4 +73,5 @@ BasalforcingsPicoFarOceantemperatureEnum, BasalforcingsPicoGammaTEnum, + BasalforcingsPicoIsplumeEnum, BasalforcingsPicoMaxboxcountEnum, BasalforcingsPicoNumBasinsEnum, @@ -372,8 +373,12 @@ BasalforcingsPicoBasinIdEnum, BasalforcingsPicoBoxIdEnum, + BasalforcingsPicoGroundingLineDepthEnum, + BasalforcingsPicoShelfSlopeEnum, BasalforcingsPicoSubShelfOceanOverturningEnum, BasalforcingsPicoSubShelfOceanSalinityEnum, BasalforcingsPicoSubShelfOceanTempEnum, BaseEnum, + BaseSlopeXEnum, + BaseSlopeYEnum, BedEnum, BedSlopeXEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22917) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22918) @@ -81,4 +81,5 @@ case BasalforcingsPicoFarOceantemperatureEnum : return "BasalforcingsPicoFarOceantemperature"; case BasalforcingsPicoGammaTEnum : return "BasalforcingsPicoGammaT"; + case BasalforcingsPicoIsplumeEnum : return "BasalforcingsPicoIsplume"; case BasalforcingsPicoMaxboxcountEnum : return "BasalforcingsPicoMaxboxcount"; case BasalforcingsPicoNumBasinsEnum : return "BasalforcingsPicoNumBasins"; @@ -378,8 +379,12 @@ case BasalforcingsPicoBasinIdEnum : return "BasalforcingsPicoBasinId"; case BasalforcingsPicoBoxIdEnum : return "BasalforcingsPicoBoxId"; + case BasalforcingsPicoGroundingLineDepthEnum : return "BasalforcingsPicoGroundingLineDepth"; + case BasalforcingsPicoShelfSlopeEnum : return "BasalforcingsPicoShelfSlope"; case BasalforcingsPicoSubShelfOceanOverturningEnum : return "BasalforcingsPicoSubShelfOceanOverturning"; case BasalforcingsPicoSubShelfOceanSalinityEnum : return "BasalforcingsPicoSubShelfOceanSalinity"; case BasalforcingsPicoSubShelfOceanTempEnum : return "BasalforcingsPicoSubShelfOceanTemp"; case BaseEnum : return "Base"; + case BaseSlopeXEnum : return "BaseSlopeX"; + case BaseSlopeYEnum : return "BaseSlopeY"; case BedEnum : return "Bed"; case BedSlopeXEnum : return "BedSlopeX"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22917) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22918) @@ -81,4 +81,5 @@ else if (strcmp(name,"BasalforcingsPicoFarOceantemperature")==0) return BasalforcingsPicoFarOceantemperatureEnum; else if (strcmp(name,"BasalforcingsPicoGammaT")==0) return BasalforcingsPicoGammaTEnum; + else if (strcmp(name,"BasalforcingsPicoIsplume")==0) return BasalforcingsPicoIsplumeEnum; else if (strcmp(name,"BasalforcingsPicoMaxboxcount")==0) return BasalforcingsPicoMaxboxcountEnum; else if (strcmp(name,"BasalforcingsPicoNumBasins")==0) return BasalforcingsPicoNumBasinsEnum; @@ -136,9 +137,9 @@ else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum; else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum; - else if (strcmp(name,"HydrologydcEplflipLock")==0) return HydrologydcEplflipLockEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"HydrologydcEplThickComp")==0) return HydrologydcEplThickCompEnum; + if (strcmp(name,"HydrologydcEplflipLock")==0) return HydrologydcEplflipLockEnum; + else if (strcmp(name,"HydrologydcEplThickComp")==0) return HydrologydcEplThickCompEnum; else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum; else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"SealevelriseTidalLoveH")==0) return SealevelriseTidalLoveHEnum; else if (strcmp(name,"SealevelriseTidalLoveK")==0) return SealevelriseTidalLoveKEnum; - else if (strcmp(name,"SealevelriseTransitions")==0) return SealevelriseTransitionsEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum; + if (strcmp(name,"SealevelriseTransitions")==0) return SealevelriseTransitionsEnum; + else if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum; else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum; else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum; @@ -382,13 +383,17 @@ else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum; else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum; - else if (strcmp(name,"BasalforcingsPicoBasinId")==0) return BasalforcingsPicoBasinIdEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum; + if (strcmp(name,"BasalforcingsPicoBasinId")==0) return BasalforcingsPicoBasinIdEnum; + else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum; + else if (strcmp(name,"BasalforcingsPicoGroundingLineDepth")==0) return BasalforcingsPicoGroundingLineDepthEnum; + else if (strcmp(name,"BasalforcingsPicoShelfSlope")==0) return BasalforcingsPicoShelfSlopeEnum; else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum; else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum; else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum; else if (strcmp(name,"Base")==0) return BaseEnum; + else if (strcmp(name,"BaseSlopeX")==0) return BaseSlopeXEnum; + else if (strcmp(name,"BaseSlopeY")==0) return BaseSlopeYEnum; else if (strcmp(name,"Bed")==0) return BedEnum; else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum; @@ -501,13 +506,13 @@ else if (strcmp(name,"Misfit")==0) return MisfitEnum; else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; - else if (strcmp(name,"Node")==0) return NodeEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"Node")==0) return NodeEnum; else if (strcmp(name,"OmegaAbsGradient")==0) return OmegaAbsGradientEnum; else if (strcmp(name,"P0")==0) return P0Enum; else if (strcmp(name,"P1")==0) return P1Enum; else if (strcmp(name,"Pressure")==0) return PressureEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum; + else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum; else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; @@ -624,13 +629,13 @@ else if (strcmp(name,"Vy")==0) return VyEnum; else if (strcmp(name,"VyMesh")==0) return VyMeshEnum; - else if (strcmp(name,"VyObs")==0) return VyObsEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"VyObs")==0) return VyObsEnum; else if (strcmp(name,"Vz")==0) return VzEnum; else if (strcmp(name,"VzFS")==0) return VzFSEnum; else if (strcmp(name,"VzHO")==0) return VzHOEnum; else if (strcmp(name,"VzMesh")==0) return VzMeshEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"VzSSA")==0) return VzSSAEnum; + else if (strcmp(name,"VzSSA")==0) return VzSSAEnum; else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum; else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum; @@ -747,13 +752,13 @@ else if (strcmp(name,"FlowequationBorderFS")==0) return FlowequationBorderFSEnum; else if (strcmp(name,"Free")==0) return FreeEnum; - else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum; else if (strcmp(name,"Fset")==0) return FsetEnum; else if (strcmp(name,"FSpressure")==0) return FSpressureEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"FSSolver")==0) return FSSolverEnum; + else if (strcmp(name,"FSSolver")==0) return FSSolverEnum; else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum; @@ -870,13 +875,13 @@ else if (strcmp(name,"Matestar")==0) return MatestarEnum; else if (strcmp(name,"Matice")==0) return MaticeEnum; - else if (strcmp(name,"Matlitho")==0) return MatlithoEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"Matlitho")==0) return MatlithoEnum; else if (strcmp(name,"Matpar")==0) return MatparEnum; else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum; else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum; + else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum; else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum; else if (strcmp(name,"MaxVel")==0) return MaxVelEnum; @@ -993,13 +998,13 @@ else if (strcmp(name,"Outputdefinition81")==0) return Outputdefinition81Enum; else if (strcmp(name,"Outputdefinition82")==0) return Outputdefinition82Enum; - else if (strcmp(name,"Outputdefinition83")==0) return Outputdefinition83Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition83")==0) return Outputdefinition83Enum; else if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum; else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum; else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum; else if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum; + else if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum; else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum; else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; @@ -1116,13 +1121,13 @@ else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum; else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum; - else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum; else if (strcmp(name,"TotalFloatingBmbScaled")==0) return TotalFloatingBmbScaledEnum; else if (strcmp(name,"TotalGroundedBmb")==0) return TotalGroundedBmbEnum; else if (strcmp(name,"TotalGroundedBmbScaled")==0) return TotalGroundedBmbScaledEnum; else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; + else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; Index: /issm/trunk-jpl/src/m/classes/basalforcingspico.m =================================================================== --- /issm/trunk-jpl/src/m/classes/basalforcingspico.m (revision 22917) +++ /issm/trunk-jpl/src/m/classes/basalforcingspico.m (revision 22918) @@ -13,4 +13,5 @@ farocean_temperature = NaN; farocean_salinity = NaN; + isplume = NaN; geothermalflux = NaN; groundedice_melting_rate = NaN; @@ -38,10 +39,8 @@ disp(' no maximum number of boxes set, setting value to 5'); end - if isnan(self.overturning_coeff) self.overturning_coeff = 1e6; %m^3/s disp(' no overturning strength set, setting value to 1e6'); end - if isnan(self.gamma_T) self.gamma_T = 2e-5; %m/s @@ -56,7 +55,8 @@ function self = setdefaultparameters(self) % {{{ - self.maxboxcount = 5; + self.maxboxcount = 5; self.overturning_coeff = 1e6; %m^3/s - self.gamma_T = 2e-5; %m/s + self.gamma_T = 2e-5; %m/s + self.isplume = false; end % }}} @@ -70,4 +70,5 @@ md = checkfield(md,'fieldname','basalforcings.farocean_temperature','NaN',1,'Inf',1,'>',0,'size',[md.basalforcings.num_basins+1 NaN]); md = checkfield(md,'fieldname','basalforcings.farocean_salinity','NaN',1,'Inf',1,'>',0,'size',[md.basalforcings.num_basins+1 NaN]); + md = checkfield(md,'fieldname','basalforcings.isplume','values',[0 1]); md = checkfield(md,'fieldname','basalforcings.geothermalflux','NaN',1,'Inf',1,'>=',0,'timeseries',1); md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1); @@ -79,8 +80,9 @@ fielddisplay(self,'basin_id','basin number assigned to each node [unitless]'); fielddisplay(self,'maxboxcount','maximum number of boxes initialized under all ice shelves'); - fielddisplay(self,'overturning_coeff','Overturning strength [m^3/s]'); - fielddisplay(self,'gamma_T','Turbulent temperature exchange velocity [m/s]'); + fielddisplay(self,'overturning_coeff','overturning strength [m^3/s]'); + fielddisplay(self,'gamma_T','turbulent temperature exchange velocity [m/s]'); fielddisplay(self,'farocean_temperature','depth averaged ocean temperature in front of the ice shelf for basin i [K]'); fielddisplay(self,'farocean_salinity','depth averaged ocean salinity in front of the ice shelf for basin i [psu]'); + fielddisplay(self,'isplume','boolean to use buoyant plume melt rate parameterization from Lazeroms et al., 2018 (default false)'); fielddisplay(self,'geothermalflux','geothermal heat flux [W/m^2]'); fielddisplay(self,'groundedice_melting_rate','basal melting rate (positive if melting) [m/yr]'); @@ -98,6 +100,6 @@ WriteData(fid,prefix,'object',self,'fieldname','farocean_temperature','format','DoubleMat','name','md.basalforcings.farocean_temperature','timeserieslength',md.basalforcings.num_basins+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','farocean_salinity','format','DoubleMat','name','md.basalforcings.farocean_salinity','timeserieslength',md.basalforcings.num_basins+1,'yts',md.constants.yts); - %WriteData(fid,prefix,'object',self,'fieldname','basin_id','format','DoubleMat','name','md.basalforcings.basin_id','mattype',2); WriteData(fid,prefix,'object',self,'fieldname','basin_id','data',self.basin_id-1,'name','md.basalforcings.basin_id','format','IntMat','mattype',2); %Change to 0-indexing + WriteData(fid,prefix,'object',self,'fieldname','isplume','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','format','DoubleMat','name','md.basalforcings.geothermalflux','mattype',1,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); Index: /issm/trunk-jpl/test/NightlyRun/test470.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test470.m (revision 22917) +++ /issm/trunk-jpl/test/NightlyRun/test470.m (revision 22918) @@ -24,4 +24,5 @@ md.basalforcings.farocean_salinity = [31 32 33; 34 35 36; 0.5 1 1.5]; %PSU md.basalforcings.maxboxcount=5; +md.basalforcings.isplume = 0; %Boundary conditions: