Index: /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 13520) +++ /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 13521) @@ -197,4 +197,6 @@ SurfaceforcingsMonthlytemperaturesEnum, SurfaceforcingsHcEnum, + SurfaceforcingsHrefEnum, + SurfaceforcingsSmbrefEnum, SurfaceforcingsSmbPosMaxEnum, SurfaceforcingsSmbPosMinEnum, Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 13520) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 13521) @@ -2657,4 +2657,6 @@ IssmDouble b_neg[NUMVERTICES]; // Hs-SMB relation paremeter IssmDouble Hc[NUMVERTICES]; // elevation of transition between accumulation regime and ablation regime + IssmDouble Href[NUMVERTICES]; // reference elevation from which deviations are used to calculate the SMB adjustment + IssmDouble Smbref[NUMVERTICES]; // reference SMB to which deviations are added IssmDouble smb_pos_max[NUMVERTICES]; // maximum SMB value in the accumulation regime IssmDouble smb_pos_min[NUMVERTICES]; // minimum SMB value in the accumulation regime @@ -2667,4 +2669,6 @@ /*Recover SmbGradients*/ GetInputListOnVertices(&Hc[0],SurfaceforcingsHcEnum); + GetInputListOnVertices(&Href[0],SurfaceforcingsHrefEnum); + GetInputListOnVertices(&Smbref[0],SurfaceforcingsSmbrefEnum); GetInputListOnVertices(&smb_pos_max[0],SurfaceforcingsSmbPosMaxEnum); GetInputListOnVertices(&smb_pos_min[0],SurfaceforcingsSmbPosMinEnum); @@ -2683,15 +2687,16 @@ // loop over all vertices - for(i=0;iHc[i]){ - smb[i]=a_pos[i]+b_pos[i]*s[i]; - if(smb[i]>smb_pos_max[i]){smb[i]=smb_pos_max[i];} - if(smb[i]Hc[i]){smb[i]=Smbref[i]+b_pos[i]*(s[i]-Href[i]);} + if(Href[i]<=Hc[i]){smb[i]=a_pos[i]+b_pos[i]*s[i];} + if(smb[i]>smb_pos_max[i]){smb[i]=smb_pos_max[i];} + if(smb[i]Hc[i]){smb[i]=a_neg[i]+b_neg[i]*s[i];} + if(Href[i]<=Hc[i]){smb[i]=Smbref[i]+b_neg[i]*(s[i]-Href[i]);} } smb[i]=smb[i]/rho_ice; // SMB in m/y ice - } //end of the loop over the vertices /*Update inputs*/ @@ -3654,8 +3659,8 @@ vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14; h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.)); - - K[0][0]=h/(2*vel)*fabs(vx*vx); K[0][1]=h/(2*vel)*fabs(vx*vy); K[0][2]=h/(2*vel)*fabs(vx*vz); - K[1][0]=h/(2*vel)*fabs(vy*vx); K[1][1]=h/(2*vel)*fabs(vy*vy); K[1][2]=h/(2*vel)*fabs(vy*vz); - K[2][0]=h/(2*vel)*fabs(vz*vx); K[2][1]=h/(2*vel)*fabs(vz*vy); K[2][2]=h/(2*vel)*fabs(vz*vz); + vz=0; + K[0][0]=h/(2*vel)*(vx*vx); K[0][1]=h/(2*vel)*(vx*vy); K[0][2]=h/(2*vel)*(vx*vz); + K[1][0]=h/(2*vel)*(vy*vx); K[1][1]=h/(2*vel)*(vy*vy); K[1][2]=h/(2*vel)*(vy*vz); + K[2][0]=h/(2*vel)*(vz*vx); K[2][1]=h/(2*vel)*(vz*vy); K[2][2]=h/(2*vel)*(vz*vz); D_scalar_stab=gauss->weight*Jdet; @@ -4314,5 +4319,5 @@ IssmDouble B_average,s_average; int* doflist=NULL; - //IssmDouble pressure[numdof]; + IssmDouble pressure[numdof]; /*Get dof list: */ @@ -4322,7 +4327,7 @@ for(i=0;imatpar->TMeltingPoint(pressure[i])) values[i]=matpar->TMeltingPoint(pressure[i]); - //if(values[i]TMeltingPoint(pressure[i])-50) values[i]=matpar->TMeltingPoint(pressure[i])-50; + GetInputListOnVertices(&pressure[0],PressureEnum); + if(values[i]>matpar->TMeltingPoint(pressure[i])) values[i]=matpar->TMeltingPoint(pressure[i]); + if(values[i]TMeltingPoint(pressure[i])-50) values[i]=matpar->TMeltingPoint(pressure[i])-50; /*Check solution*/ Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 13520) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 13521) @@ -2263,4 +2263,6 @@ IssmDouble b_neg[NUMVERTICES]; // Hs-SMB relation paremeter IssmDouble Hc[NUMVERTICES]; // elevation of transition between accumulation regime and ablation regime + IssmDouble Href[NUMVERTICES]; // reference elevation from which deviations are used to calculate the SMB adjustment + IssmDouble Smbref[NUMVERTICES]; // reference SMB to which deviations are added IssmDouble smb_pos_max[NUMVERTICES]; // maximum SMB value in the accumulation regime IssmDouble smb_pos_min[NUMVERTICES]; // minimum SMB value in the accumulation regime @@ -2273,4 +2275,6 @@ /*Recover SmbGradients*/ GetInputListOnVertices(&Hc[0],SurfaceforcingsHcEnum); + GetInputListOnVertices(&Href[0],SurfaceforcingsHrefEnum); + GetInputListOnVertices(&Smbref[0],SurfaceforcingsSmbrefEnum); GetInputListOnVertices(&smb_pos_max[0],SurfaceforcingsSmbPosMaxEnum); GetInputListOnVertices(&smb_pos_min[0],SurfaceforcingsSmbPosMinEnum); @@ -2291,13 +2295,22 @@ for(i=0;iHc[i]){ - smb[i]=a_pos[i]+b_pos[i]*s[i]; - if(smb[i]>smb_pos_max[i]){smb[i]=smb_pos_max[i];} - if(smb[i]Hc[i]){smb[i]=Smbref[i]+b_pos[i]*(s[i]-Href[i]);} + if(Href[i]<=Hc[i]){smb[i]=a_pos[i]+b_pos[i]*s[i];} + if(smb[i]>smb_pos_max[i]){smb[i]=smb_pos_max[i];} + if(smb[i]Hc[i]){smb[i]=a_neg[i]+b_neg[i]*s[i];} + if(Href[i]<=Hc[i]){smb[i]=Smbref[i]+b_neg[i]*(s[i]-Href[i]);} } smb[i]=smb[i]/rho_ice; // SMB in m/y ice - } //end of the loop over the vertices + /* printf("s %e \n",s[i]); + printf("Hsref %e \n",Href[i]); + printf("Hc %e \n",Hc[i]); + printf("Smbref %e \n",Smbref[i]); + printf("b_neg %e \n",b_neg[i]); + printf("smb %e \n",smb[i]); + _error_("stop-in-code"); */ + } //end of the loop over the vertices /*Update inputs*/ this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&smb[0])); Index: /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 13520) +++ /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 13521) @@ -202,4 +202,6 @@ case SurfaceforcingsMonthlytemperaturesEnum : return "SurfaceforcingsMonthlytemperatures"; case SurfaceforcingsHcEnum : return "SurfaceforcingsHc"; + case SurfaceforcingsHrefEnum : return "SurfaceforcingsHref"; + case SurfaceforcingsSmbrefEnum : return "SurfaceforcingsSmbref"; case SurfaceforcingsSmbPosMaxEnum : return "SurfaceforcingsSmbPosMax"; case SurfaceforcingsSmbPosMinEnum : return "SurfaceforcingsSmbPosMin"; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp (revision 13520) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp (revision 13521) @@ -73,4 +73,6 @@ if(issmbgradients){ iomodel->FetchDataToInput(elements,SurfaceforcingsHcEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsHrefEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsSmbrefEnum); iomodel->FetchDataToInput(elements,SurfaceforcingsSmbPosMaxEnum); iomodel->FetchDataToInput(elements,SurfaceforcingsSmbPosMinEnum); Index: /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 13520) +++ /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 13521) @@ -206,4 +206,6 @@ else if (strcmp(name,"SurfaceforcingsMonthlytemperatures")==0) return SurfaceforcingsMonthlytemperaturesEnum; else if (strcmp(name,"SurfaceforcingsHc")==0) return SurfaceforcingsHcEnum; + else if (strcmp(name,"SurfaceforcingsHref")==0) return SurfaceforcingsHrefEnum; + else if (strcmp(name,"SurfaceforcingsSmbref")==0) return SurfaceforcingsSmbrefEnum; else if (strcmp(name,"SurfaceforcingsSmbPosMax")==0) return SurfaceforcingsSmbPosMaxEnum; else if (strcmp(name,"SurfaceforcingsSmbPosMin")==0) return SurfaceforcingsSmbPosMinEnum; @@ -259,10 +261,10 @@ else if (strcmp(name,"PrognosticSolution")==0) return PrognosticSolutionEnum; else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum; - else if (strcmp(name,"SurfaceSlopeAnalysis")==0) return SurfaceSlopeAnalysisEnum; - else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"SurfaceSlopeXAnalysis")==0) return SurfaceSlopeXAnalysisEnum; + if (strcmp(name,"SurfaceSlopeAnalysis")==0) return SurfaceSlopeAnalysisEnum; + else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; + else if (strcmp(name,"SurfaceSlopeXAnalysis")==0) return SurfaceSlopeXAnalysisEnum; else if (strcmp(name,"SurfaceSlopeYAnalysis")==0) return SurfaceSlopeYAnalysisEnum; else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum; @@ -382,10 +384,10 @@ else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum; else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; - else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; - else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; + if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; + else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; + else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum; else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum; @@ -505,10 +507,10 @@ else if (strcmp(name,"Separate")==0) return SeparateEnum; else if (strcmp(name,"Sset")==0) return SsetEnum; - else if (strcmp(name,"Verbose")==0) return VerboseEnum; - else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum; + if (strcmp(name,"Verbose")==0) return VerboseEnum; + else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum; + else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum; else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum; else if (strcmp(name,"XY")==0) return XYEnum; Index: /issm/trunk-jpl/src/c/shared/Numerics/UnitConversion.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Numerics/UnitConversion.cpp (revision 13520) +++ /issm/trunk-jpl/src/c/shared/Numerics/UnitConversion.cpp (revision 13521) @@ -72,4 +72,5 @@ case SurfaceforcingsANegEnum: scale=yts;break; //m/yr case SurfaceforcingsBNegEnum: scale=yts;break; //m/yr + case SurfaceforcingsSmbrefEnum: scale=yts;break; //m/yr case MisfitEnum: scale=pow(yts,2);break; //(m/yr)^2 case MassFluxEnum: scale=pow((IssmDouble)10,-12)*yts;break; // (GigaTon/year) Index: /issm/trunk-jpl/src/m/classes/surfaceforcings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/surfaceforcings.m (revision 13520) +++ /issm/trunk-jpl/src/m/classes/surfaceforcings.m (revision 13521) @@ -12,4 +12,6 @@ isdelta18o = 0; hc = NaN; + href = NaN; + smbref = NaN; smb_pos_max = NaN; smb_pos_min = NaN; @@ -60,4 +62,6 @@ elseif(obj.issmbgradients) md = checkfield(md,'surfaceforcings.hc','forcing',1,'NaN',1); + md = checkfield(md,'surfaceforcings.href','forcing',1,'NaN',1); + md = checkfield(md,'surfaceforcings.smbref','forcing',1,'NaN',1); md = checkfield(md,'surfaceforcings.smb_pos_max','forcing',1,'NaN',1); md = checkfield(md,'surfaceforcings.smb_pos_min','forcing',1,'NaN',1); @@ -89,4 +93,6 @@ fielddisplay(obj,'issmbgradients','is smb gradients method activated (0 or 1, default is 0)'); fielddisplay(obj,'hc',' elevation of intersection between accumulation and ablation regime required if smb gradients is activated'); + fielddisplay(obj,'href',' reference elevation from which deviation is used to calculate SMB adjustment in smb gradients method'); + fielddisplay(obj,'smbref',' reference smb from which deviation is calculated in smb gradients method'); fielddisplay(obj,'smb_pos_max',' maximum value of positive smb required if smb gradients is activated'); fielddisplay(obj,'smb_pos_min',' minimum value of positive smb required if smb gradients is activated'); @@ -117,4 +123,6 @@ if obj.issmbgradients, WriteData(fid,'object',obj,'fieldname','hc','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','href','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','smbref','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'fieldname','smb_pos_max','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'fieldname','smb_pos_min','format','DoubleMat','mattype',1);