Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23905) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23906) @@ -2726,4 +2726,5 @@ IssmDouble tf,gamma0,base,delta_t_basin,mean_tf_basin,absval; IssmDouble basalmeltrate[NUMVERTICES]; + bool islocal; IssmDouble* delta_t = NULL; IssmDouble* mean_tf = NULL; @@ -2744,7 +2745,11 @@ this->parameters->FindParam(&gamma0,BasalforcingsIsmp6Gamma0Enum); this->parameters->FindParam(&delta_t,&M,BasalforcingsIsmp6DeltaTEnum); _assert_(M==num_basins); - this->parameters->FindParam(&mean_tf,&N,BasalforcingsIsmp6AverageTfEnum); _assert_(N==num_basins); + this->parameters->FindParam(&islocal,BasalforcingsIsmp6IsLocalEnum); + if(!islocal) { + this->parameters->FindParam(&mean_tf,&N,BasalforcingsIsmp6AverageTfEnum); _assert_(N==num_basins); + } Input* tf_input=this->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input); - delta_t_basin = delta_t[basinid]; mean_tf_basin = mean_tf[basinid]; + delta_t_basin = delta_t[basinid]; + if(!islocal) mean_tf_basin = mean_tf[basinid]; /* Compute melt rate */ @@ -2753,7 +2758,9 @@ gauss->GaussVertex(i); tf_input->GetInputValue(&tf,gauss); - absval = mean_tf_basin+delta_t_basin; - if (absval<0) {absval = -1.*absval;} - basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*(tf+delta_t_basin)*absval; + if(!islocal) { + absval = mean_tf_basin+delta_t_basin; + if (absval<0) {absval = -1.*absval;} + basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*(tf+delta_t_basin)*absval;} + else {basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*pow(max(tf+delta_t_basin,0.),2);} } @@ -3489,5 +3496,5 @@ /*Compute ice shelf base slope (radians)*/ - alpha = atan(sqrt(slopex*slopex + slopey*slopey)); + alpha = atan(sqrt(slopex*slopex + slopey*slopey)); if(alpha>=M_PI) alpha = M_PI - 0.001; //ensure sin(alpha) > 0 for meltrate calculations Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 23905) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 23906) @@ -75,8 +75,10 @@ int num_basins, basinid,num_depths; IssmDouble area, tf, base, time; + bool islocal; IssmDouble* tf_depths = NULL; femmodel->parameters->FindParam(&num_basins,BasalforcingsIsmp6NumBasinsEnum); femmodel->parameters->FindParam(&tf_depths,&num_depths,BasalforcingsIsmp6TfDepthsEnum); _assert_(tf_depths); + femmodel->parameters->FindParam(&islocal,BasalforcingsIsmp6IsLocalEnum); /*Binary search works for vectors that are sorted in increasing order only, make depths positive*/ @@ -148,25 +150,27 @@ } - /*Compute sums of tf*area and shelf-area per cpu*/ - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - if(!element->IsIceInElement() || !element->IsFloating()) continue; - Input* tf_input=element->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input); - element->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum); - Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0); - tf_input->GetInputValue(&tf,gauss); - delete gauss; - area=element->GetHorizontalSurfaceArea(); - tf_weighted_avg[basinid]+=tf*area; - areas_summed[basinid] +=area; + if(!islocal) { + /*Compute sums of tf*area and shelf-area per cpu*/ + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + if(!element->IsIceInElement() || !element->IsFloating()) continue; + Input* tf_input=element->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input); + element->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum); + Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0); + tf_input->GetInputValue(&tf,gauss); + delete gauss; + area=element->GetHorizontalSurfaceArea(); + tf_weighted_avg[basinid]+=tf*area; + areas_summed[basinid] +=area; + } + + /*Syncronize across cpus*/ + ISSM_MPI_Allreduce(tf_weighted_avg,tf_weighted_avg_cpu,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + ISSM_MPI_Allreduce(areas_summed,areas_summed_cpu,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + + /*Compute weighted means and save*/ + for(int k=0;kparameters->AddObject(new DoubleVecParam(BasalforcingsIsmp6AverageTfEnum,tf_weighted_avg_cpu,num_basins)); } - - /*Syncronize across cpus*/ - ISSM_MPI_Allreduce(tf_weighted_avg,tf_weighted_avg_cpu,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); - ISSM_MPI_Allreduce(areas_summed,areas_summed_cpu,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); - - /*Compute weighted means and save*/ - for(int k=0;kparameters->AddObject(new DoubleVecParam(BasalforcingsIsmp6AverageTfEnum,tf_weighted_avg_cpu,num_basins)); /*Cleanup and return */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23905) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23906) @@ -238,4 +238,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_basins",BasalforcingsIsmp6NumBasinsEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_0",BasalforcingsIsmp6Gamma0Enum)); + parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.islocal",BasalforcingsIsmp6IsLocalEnum)); iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.delta_t"); parameters->AddObject(new DoubleVecParam(BasalforcingsIsmp6DeltaTEnum,transparam,N)); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23905) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23906) @@ -67,4 +67,5 @@ BasalforcingsIsmp6DeltaTEnum, BasalforcingsIsmp6Gamma0Enum, + BasalforcingsIsmp6IsLocalEnum, BasalforcingsIsmp6NumBasinsEnum, BasalforcingsIsmp6TfDepthsEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23905) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23906) @@ -75,4 +75,5 @@ case BasalforcingsIsmp6DeltaTEnum : return "BasalforcingsIsmp6DeltaT"; case BasalforcingsIsmp6Gamma0Enum : return "BasalforcingsIsmp6Gamma0"; + case BasalforcingsIsmp6IsLocalEnum : return "BasalforcingsIsmp6IsLocal"; case BasalforcingsIsmp6NumBasinsEnum : return "BasalforcingsIsmp6NumBasins"; case BasalforcingsIsmp6TfDepthsEnum : return "BasalforcingsIsmp6TfDepths"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23905) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23906) @@ -75,4 +75,5 @@ else if (strcmp(name,"BasalforcingsIsmp6DeltaT")==0) return BasalforcingsIsmp6DeltaTEnum; else if (strcmp(name,"BasalforcingsIsmp6Gamma0")==0) return BasalforcingsIsmp6Gamma0Enum; + else if (strcmp(name,"BasalforcingsIsmp6IsLocal")==0) return BasalforcingsIsmp6IsLocalEnum; else if (strcmp(name,"BasalforcingsIsmp6NumBasins")==0) return BasalforcingsIsmp6NumBasinsEnum; else if (strcmp(name,"BasalforcingsIsmp6TfDepths")==0) return BasalforcingsIsmp6TfDepthsEnum; @@ -136,9 +137,9 @@ else if (strcmp(name,"FemModelComm")==0) return FemModelCommEnum; else if (strcmp(name,"FlowequationFeFS")==0) return FlowequationFeFSEnum; - else if (strcmp(name,"FlowequationIsFS")==0) return FlowequationIsFSEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"FlowequationIsHO")==0) return FlowequationIsHOEnum; + if (strcmp(name,"FlowequationIsFS")==0) return FlowequationIsFSEnum; + else if (strcmp(name,"FlowequationIsHO")==0) return FlowequationIsHOEnum; else if (strcmp(name,"FlowequationIsL1L2")==0) return FlowequationIsL1L2Enum; else if (strcmp(name,"FlowequationIsSIA")==0) return FlowequationIsSIAEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"MaterialsTemperateiceconductivity")==0) return MaterialsTemperateiceconductivityEnum; else if (strcmp(name,"MaterialsEffectiveconductivityAveraging")==0) return MaterialsEffectiveconductivityAveragingEnum; - else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum; + if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum; + else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum; else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum; else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum; else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum; - else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"Step")==0) return StepEnum; + if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; + else if (strcmp(name,"Step")==0) return StepEnum; else if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum; else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"EffectivePressureHydrostep")==0) return EffectivePressureHydrostepEnum; else if (strcmp(name,"EffectivePressureStacked")==0) return EffectivePressureStackedEnum; - else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; + if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum; + else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; else if (strcmp(name,"EplHeadStacked")==0) return EplHeadStackedEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum; else if (strcmp(name,"SealevelUEsa")==0) return SealevelUEsaEnum; - else if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum; + if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum; + else if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum; else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum; else if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; - else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; + if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; + else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum; else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum; else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum; - else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum; + if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum; + else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum; else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum; else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum; @@ -997,9 +998,9 @@ else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum; - else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum; + if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum; + else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum; else if (strcmp(name,"Fset")==0) return FsetEnum; else if (strcmp(name,"FSpressure")==0) return FSpressureEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"Mumps")==0) return MumpsEnum; else if (strcmp(name,"Nodal")==0) return NodalEnum; - else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; + if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum; + else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; else if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum; else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"VertexLId")==0) return VertexLIdEnum; else if (strcmp(name,"Vertices")==0) return VerticesEnum; - else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"Water")==0) return WaterEnum; + if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; + else if (strcmp(name,"Water")==0) return WaterEnum; else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum; else if (strcmp(name,"XY")==0) return XYEnum; Index: /issm/trunk-jpl/src/m/classes/basalforcingsismip6.m =================================================================== --- /issm/trunk-jpl/src/m/classes/basalforcingsismip6.m (revision 23905) +++ /issm/trunk-jpl/src/m/classes/basalforcingsismip6.m (revision 23906) @@ -12,4 +12,5 @@ tf_depths = NaN; delta_t = NaN; + islocal = 0; geothermalflux = NaN; groundedice_melting_rate = NaN; @@ -47,4 +48,5 @@ function self = setdefaultparameters(self) % {{{ self.gamma_0 = 14477; %m/yr + self.islocal = false; end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -55,4 +57,5 @@ md = checkfield(md,'fieldname','basalforcings.tf_depths','NaN',1,'Inf',1); md = checkfield(md,'fieldname','basalforcings.delta_t','NaN',1,'Inf',1,'numel',md.basalforcings.num_basins); + md = checkfield(md,'fieldname','basalforcings.islocal','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); @@ -72,4 +75,5 @@ fielddisplay(self,'tf','thermal forcing (ocean temperature minus freezing point) (degrees C)'); fielddisplay(self,'delta_t','Ocean temperature correction per basin (degrees C)'); + fielddisplay(self,'islocal','boolean to use the local version of the ISMIP6 melt rate parameterization (default false)'); fielddisplay(self,'geothermalflux','geothermal heat flux (W/m^2)'); fielddisplay(self,'groundedice_melting_rate','basal melting rate (positive if melting) (m/yr)'); @@ -87,4 +91,5 @@ WriteData(fid,prefix,'object',self,'fieldname','tf','format','MatArray','name','md.basalforcings.tf','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','delta_t','format','DoubleMat','name','md.basalforcings.delta_t','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'fieldname','islocal','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/test472.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test472.m (revision 23905) +++ /issm/trunk-jpl/test/NightlyRun/test472.m (revision 23906) @@ -24,4 +24,5 @@ md.basalforcings.tf_depths = [0 -1000 -2000]; md.basalforcings.gamma_0 = 14477; +md.basalforcings.islocal = 0; %Build an artificial tf field (for times 0 and 1, 3 depth layers)