Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 24734) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 24735) @@ -245,4 +245,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.smb.isturbulentflux",SmbIsturbulentfluxEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.isconstrainsurfaceT",SmbIsconstrainsurfaceTEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum)); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24734) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24735) @@ -3631,4 +3631,5 @@ bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux; bool isclimatology=false; + bool isconstrainsurfaceT=false; IssmDouble init_scaling=0.0; IssmDouble thermo_scaling=1.0; @@ -3694,4 +3695,5 @@ parameters->FindParam(&isdensification,SmbIsdensificationEnum); parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum); + parameters->FindParam(&isconstrainsurfaceT,SmbIsconstrainsurfaceTEnum); parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum); parameters->FindParam(&thermo_scaling,SmbThermoDeltaTScalingEnum); @@ -3893,6 +3895,10 @@ netSW = netSW + cellsum(swf,m)*smb_dt/dt; + if(isconstrainsurfaceT){ + if (m>0) T[0]=Ta; + if (m>1) T[1]=Ta; + } /*Thermal profile computation:*/ - if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid()); + if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid(),isconstrainsurfaceT); /*Change in thickness of top cell due to evaporation/condensation assuming same density as top cell. @@ -3905,5 +3911,5 @@ /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/ - if(ismelt)melt(&M, &Msurf, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid()); + if(ismelt)melt(&M, &Msurf, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop, zY, rho_ice,this->Sid()); /*Allow non-melt densification and determine compaction [m]*/ @@ -3956,5 +3962,5 @@ for(int i=0;i 0) fac += dz[i]*(rho_ice - fmin(d[i],rho_ice)); } Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24734) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24735) @@ -615,5 +615,5 @@ } /*}}}*/ -void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* pulwrf, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid) { /*{{{*/ +void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* pulwrf, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid, bool isconstrainsurfaceT) { /*{{{*/ /* ENGLACIAL THERMODYNAMICS*/ @@ -1372,5 +1372,5 @@ *pm=m; } /*}}}*/ -void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid){ /*{{{*/ +void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY, IssmDouble dIce, int sid){ /*{{{*/ //// MELT ROUTINE @@ -1423,5 +1423,5 @@ IssmDouble* Zcum=NULL; IssmDouble* dzMin2=NULL; - IssmDouble zY2=1.025; + IssmDouble zY2=zY; bool lastCellFlag = false; int X1=0; @@ -1591,5 +1591,5 @@ depthice=0.0; if (d[i] >= dPHC-Dtol){ - for(int l=i;(l=dPHC-Dtol);l++) depthice=depthice+dz[l]; + for(int l=i;(l=dPHC-Dtol);l++) depthice=depthice+dz[l]; } @@ -1718,5 +1718,5 @@ for (int i=0;i dzMin * 2') // adjust variables as a linearly weighted function of mass Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 24734) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 24735) @@ -31,7 +31,7 @@ void albedo(IssmDouble** a,int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid); void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid); -void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* pulwrf, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid); +void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* pulwrf, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid, bool isconstrainsurfaceT); void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid); -void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid); +void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY, IssmDouble dIce, int sid); void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid); void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24734) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24735) @@ -354,4 +354,5 @@ syn keyword cConstant SmbIsalbedoEnum syn keyword cConstant SmbIsclimatologyEnum +syn keyword cConstant SmbIsconstrainsurfaceTEnum syn keyword cConstant SmbIsd18opdEnum syn keyword cConstant SmbIsdelta18oEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24734) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24735) @@ -348,4 +348,5 @@ SmbIsalbedoEnum, SmbIsclimatologyEnum, + SmbIsconstrainsurfaceTEnum, SmbIsd18opdEnum, SmbIsdelta18oEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24734) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24735) @@ -356,4 +356,5 @@ case SmbIsalbedoEnum : return "SmbIsalbedo"; case SmbIsclimatologyEnum : return "SmbIsclimatology"; + case SmbIsconstrainsurfaceTEnum : return "SmbIsconstrainsurfaceT"; case SmbIsd18opdEnum : return "SmbIsd18opd"; case SmbIsdelta18oEnum : return "SmbIsdelta18o"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24734) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24735) @@ -362,4 +362,5 @@ else if (strcmp(name,"SmbIsalbedo")==0) return SmbIsalbedoEnum; else if (strcmp(name,"SmbIsclimatology")==0) return SmbIsclimatologyEnum; + else if (strcmp(name,"SmbIsconstrainsurfaceT")==0) return SmbIsconstrainsurfaceTEnum; else if (strcmp(name,"SmbIsd18opd")==0) return SmbIsd18opdEnum; else if (strcmp(name,"SmbIsdelta18o")==0) return SmbIsdelta18oEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum; else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum; - else if (strcmp(name,"SmbRunoffalti")==0) return SmbRunoffaltiEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SmbRunoffgrad")==0) return SmbRunoffgradEnum; + if (strcmp(name,"SmbRunoffalti")==0) return SmbRunoffaltiEnum; + else if (strcmp(name,"SmbRunoffgrad")==0) return SmbRunoffgradEnum; else if (strcmp(name,"SmbRunoffref")==0) return SmbRunoffrefEnum; else if (strcmp(name,"SmbSealev")==0) return SmbSealevEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum; - else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum; + if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; + else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum; else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum; else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum; - else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; + if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; + else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum; else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum; - else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum; + if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum; + else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum; else if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum; else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum; else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum; - else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum; + if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum; + else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum; else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum; else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum; @@ -997,9 +998,9 @@ else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; else if (strcmp(name,"BoolInput")==0) return BoolInputEnum; - else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum; else stage=9; } if(stage==9){ - if (strcmp(name,"IntInput2")==0) return IntInput2Enum; + if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum; + else if (strcmp(name,"IntInput2")==0) return IntInput2Enum; else if (strcmp(name,"BoolParam")==0) return BoolParamEnum; else if (strcmp(name,"Boundary")==0) return BoundaryEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum; - else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum; + if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum; + else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum; else if (strcmp(name,"Incremental")==0) return IncrementalEnum; else if (strcmp(name,"Indexed")==0) return IndexedEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum; else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum; - else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"Regular")==0) return RegularEnum; + if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum; + else if (strcmp(name,"Regular")==0) return RegularEnum; else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum; else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 24734) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 24735) @@ -23,4 +23,5 @@ isturbulentflux; isclimatology; + isconstrainsurfaceT; %inputs: @@ -163,4 +164,5 @@ self.isturbulentflux=1; self.isclimatology=0; + self.isconstrainsurfaceT=0; self.aIdx = 1; @@ -218,4 +220,5 @@ md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]); md = checkfield(md,'fieldname','smb.isclimatology','values',[0 1]); + md = checkfield(md,'fieldname','smb.isconstrainsurfaceT','values',[0 1]); md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value @@ -298,4 +301,5 @@ fielddisplay(self,'isdensification','run densification module (default true)'); fielddisplay(self,'isturbulentflux','run turbulant heat fluxes module (default true)'); + fielddisplay(self,'isconstrainsurfaceT','constrain surface temperatures to air temperature, turn off EC and surface flux contribution to surface temperature change'); fielddisplay(self,'isclimatology','repeat all forcings when past last forcing time (default false)'); fielddisplay(self,'Ta','2 m air temperature, in Kelvin'); @@ -394,4 +398,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','isturbulentflux','format','Boolean'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','isclimatology','format','Boolean'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','isconstrainsurfaceT','format','Boolean'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); Index: /issm/trunk-jpl/src/m/classes/SMBgemb.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 24734) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 24735) @@ -102,5 +102,4 @@ #}}} - def __repr__(self): # {{{ @@ -118,4 +117,5 @@ string = "%s\n%s" % (string, fielddisplay(self, 'isdensification', 'run densification module (default true)')) string = "%s\n%s" % (string, fielddisplay(self, 'isturbulentflux', 'run turbulant heat fluxes module (default true)')) + string = "%s\n%s" % (string, fielddisplay(self, 'isconstrainsurfaceT', 'constrain surface temperatures to air temperature, turn off EC and surface flux contribution to surface temperature change')) string = "%s\n%s" % (string, fielddisplay(self, 'isclimatology', 'repeat all forcings when past last forcing time (default false)')) string = "%s\n%s" % (string, fielddisplay(self, 'Ta', '2 m air temperature, in Kelvin')) @@ -226,4 +226,5 @@ self.isturbulentflux = 1 self.isclimatology = 0 + self.isconstrainsurfaceT = 0; self.aIdx = 1 @@ -281,4 +282,5 @@ md = checkfield(md, 'fieldname', 'smb.isturbulentflux', 'values', [0, 1]) md = checkfield(md, 'fieldname', 'smb.isclimatology', 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'smb.isconstrainsurfaceT', 'values', [0, 1]) md = checkfield(md, 'fieldname', 'smb.Ta', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>', 273 - 100, '<', 273 + 100) #-100/100 celsius min/max value @@ -355,4 +357,5 @@ WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isturbulentflux', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isclimatology', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isconstrainsurfaceT', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Ta', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) Index: /issm/trunk-jpl/test/NightlyRun/test243.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test243.m (revision 24734) +++ /issm/trunk-jpl/test/NightlyRun/test243.m (revision 24735) @@ -46,16 +46,19 @@ md=solve(md,'Transient'); +nlayers=size(md.results.TransientSolution(end).SmbT,2); + %Fields and tolerances to track changes -field_names ={'SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC','SmbMeanSHF','SmbMeanLHF','SmbMeanULW','SmbNetLW','SmbNetSW'}; -field_tolerances ={8e-10,5e-12,3e-10,8e-10,1e-11,2e-7,4e-11,4e-12,1e-12,1e-12,1e-12,2e-10,2e-11, 2e-11, 1e-11, 9e-10, 2e-11}; +field_names ={'Layers','SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC','SmbMeanSHF','SmbMeanLHF','SmbMeanULW','SmbNetLW','SmbNetSW'}; +field_tolerances ={1e-12,1e-12,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-11,1e-12,1e-12,1e-12,1e-11,2e-11,2e-11,1e-11,9e-10,2e-11}; field_values={... - (md.results.TransientSolution(end).SmbDz(1,1:240)),... - (md.results.TransientSolution(end).SmbT(1,1:240)),... - (md.results.TransientSolution(end).SmbD(1,1:240)),... - (md.results.TransientSolution(end).SmbRe(1,1:240)),... - (md.results.TransientSolution(end).SmbGdn(1,1:240)),... - (md.results.TransientSolution(end).SmbGsp(1,1:240)),... - (md.results.TransientSolution(end).SmbA(1,1:240)),... + (nlayers),... + (md.results.TransientSolution(end).SmbDz(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbT(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbD(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbRe(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbGdn(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbGsp(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbA(1,1:nlayers)),... (md.results.TransientSolution(end).SmbEC(1)),... (md.results.TransientSolution(end).SmbMassBalance(1)),... Index: /issm/trunk-jpl/test/NightlyRun/test243.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test243.py (revision 24734) +++ /issm/trunk-jpl/test/NightlyRun/test243.py (revision 24735) @@ -60,15 +60,17 @@ md = solve(md, 'Transient') +nlayers=np.size(md.results.TransientSolution[-1].SmbT,1) + #Fields and tolerances to track changes -field_names = ['SmbDz', 'SmbT', 'SmbD', 'SmbRe', 'SmbGdn', 'SmbGsp', 'SmbA', 'SmbEC', 'SmbMassBalance', 'SmbMAdd', 'SmbDzAdd', 'SmbFAC', 'SmbMeanSHF', 'SmbMeanLHF', 'SmbMeanULW', 'SmbNetLW', 'SmbNetSW'] -field_tolerances = [8e-10, 5e-12, 3e-10, 8e-10, 1e-11, 2e-7, 4e-11, 4e-12, 1e-12, 1e-12, 1e-12, 2e-10, 2e-11, 2e-11, 1e-11, 9e-10, 2e-11] +field_names = ['Layers','SmbDz', 'SmbT', 'SmbD', 'SmbRe', 'SmbGdn', 'SmbGsp', 'SmbA', 'SmbEC', 'SmbMassBalance', 'SmbMAdd', 'SmbDzAdd', 'SmbFAC', 'SmbMeanSHF', 'SmbMeanLHF', 'SmbMeanULW', 'SmbNetLW', 'SmbNetSW'] +field_tolerances = [1e-12,1e-12,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-11,1e-12,1e-12,1e-12,1e-11,2e-11,2e-11,1e-11,9e-10,2e-11] #shape is different in python solution (fixed using reshape) which can cause test failure: -field_values = [md.results.TransientSolution[-1].SmbDz[0, 0:240].reshape(1, -1), - md.results.TransientSolution[-1].SmbT[0, 0:240].reshape(1, -1), - md.results.TransientSolution[-1].SmbD[0, 0:240].reshape(1, -1), - md.results.TransientSolution[-1].SmbRe[0, 0:240].reshape(1, -1), - md.results.TransientSolution[-1].SmbGdn[0, 0:240].reshape(1, -1), - md.results.TransientSolution[-1].SmbGsp[0, 0:240].reshape(1, -1), - md.results.TransientSolution[-1].SmbA[0, 0:240].reshape(1, -1), +field_values = [nlayers,md.results.TransientSolution[-1].SmbDz[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbT[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbD[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbRe[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbGdn[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbGsp[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbA[0, 0:nlayers].reshape(1, -1), md.results.TransientSolution[-1].SmbEC[0], md.results.TransientSolution[-1].SmbMassBalance[0], Index: /issm/trunk-jpl/test/NightlyRun/test252.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test252.m (revision 24734) +++ /issm/trunk-jpl/test/NightlyRun/test252.m (revision 24735) @@ -54,22 +54,25 @@ md=solve(md,'Transient'); +nlayers=size(md.results.TransientSolution(end).SmbT,2); + %Fields and tolerances to track changes -field_names ={'SmbDz1','SmbT1' ,'SmbD1' ,'SmbRe1','SmbGdn1','SmbGsp1','SmbA1' ,'SmbEC1','SmbMassBalance1','SmbMAdd1','SmbDzAdd1','SmbFAC1',... +field_names ={'Layers','SmbDz1','SmbT1' ,'SmbD1' ,'SmbRe1','SmbGdn1','SmbGsp1','SmbA1' ,'SmbEC1','SmbMassBalance1','SmbMAdd1','SmbDzAdd1','SmbFAC1',... 'SmbDz2','SmbT2' ,'SmbD2' ,'SmbRe2','SmbGdn2','SmbGsp2','SmbA2' ,'SmbEC2','SmbMassBalance2','SmbMAdd2','SmbDzAdd2','SmbFAC2',... 'SmbDz3','SmbT3' ,'SmbD3' ,'SmbRe3','SmbGdn3','SmbGsp3','SmbA3' ,'SmbEC3','SmbMassBalance3','SmbMAdd3','SmbDzAdd3','SmbFAC3',... 'SmbDz4','SmbT4' ,'SmbD4' ,'SmbRe4','SmbGdn4','SmbGsp4','SmbA4' ,'SmbEC4','SmbMassBalance4','SmbMAdd4','SmbDzAdd4','SmbFAC4'}; -field_tolerances ={1e-11,1e-12,1e-11,2e-11,1e-11,1e-11,1e-12,1e-12,1e-12,1e-12,1e-12,1e-11,... - 1e-11,1e-12,1e-11,8e-10,1e-11,4e-11,1e-12,4e-12,1e-12,1e-12,1e-12,2e-11,... - 1e-11,1e-12,1e-11,8e-10,1e-11,4e-11,1e-12,6e-12,6e-12,1e-12,1e-12,2e-11,... - 2e-9,1e-11,3e-9,1e-9,1e-11,2e-7,7e-11,2e-12,3e-9,1e-12,1e-12,9e-11}; +field_tolerances ={1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,... + 1e-12,1e-12,1e-11,1e-10,1e-11,1e-11,1e-12,4e-12,1e-12,1e-12,1e-12,1e-11,... + 1e-12,1e-12,2e-12,2e-11,1e-11,1e-11,1e-12,1e-11,1e-11,1e-12,1e-12,1e-11,... + 1e-11,1e-11,4e-11,1e-11,1e-12,3e-11,1e-12,1e-12,1e-10,1e-12,1e-12,2e-11}; field_values={... - (md.results.TransientSolution(1).SmbDz(1,1:230)),... - (md.results.TransientSolution(1).SmbT(1,1:230)),... - (md.results.TransientSolution(1).SmbD(1,1:230)),... - (md.results.TransientSolution(1).SmbRe(1,1:230)),... - (md.results.TransientSolution(1).SmbGdn(1,1:230)),... - (md.results.TransientSolution(1).SmbGsp(1,1:230)),... - (md.results.TransientSolution(1).SmbA(1,1:230)),... + (nlayers)... + (md.results.TransientSolution(1).SmbDz(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbT(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbD(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbRe(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbGdn(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbGsp(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbA(1,1:nlayers)),... (md.results.TransientSolution(1).SmbEC(1)),... (md.results.TransientSolution(1).SmbMassBalance(1)),... @@ -77,11 +80,11 @@ (md.results.TransientSolution(1).SmbDzAdd(1)),... (md.results.TransientSolution(1).SmbFAC(1)),... - (md.results.TransientSolution(146).SmbDz(1,1:230)),... - (md.results.TransientSolution(146).SmbT(1,1:230)),... - (md.results.TransientSolution(146).SmbD(1,1:230)),... - (md.results.TransientSolution(146).SmbRe(1,1:230)),... - (md.results.TransientSolution(146).SmbGdn(1,1:230)),... - (md.results.TransientSolution(146).SmbGsp(1,1:230)),... - (md.results.TransientSolution(146).SmbA(1,1:230)),... + (md.results.TransientSolution(146).SmbDz(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbT(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbD(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbRe(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbGdn(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbGsp(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbA(1,1:nlayers)),... (md.results.TransientSolution(146).SmbEC(1)),... (md.results.TransientSolution(146).SmbMassBalance(1)),... @@ -89,11 +92,11 @@ (md.results.TransientSolution(146).SmbDzAdd(1)),... (md.results.TransientSolution(146).SmbFAC(1)),... - (md.results.TransientSolution(147).SmbDz(1,1:230)),... - (md.results.TransientSolution(147).SmbT(1,1:230)),... - (md.results.TransientSolution(147).SmbD(1,1:230)),... - (md.results.TransientSolution(147).SmbRe(1,1:230)),... - (md.results.TransientSolution(147).SmbGdn(1,1:230)),... - (md.results.TransientSolution(147).SmbGsp(1,1:230)),... - (md.results.TransientSolution(147).SmbA(1,1:230)),... + (md.results.TransientSolution(147).SmbDz(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbT(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbD(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbRe(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbGdn(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbGsp(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbA(1,1:nlayers)),... (md.results.TransientSolution(147).SmbEC(1)),... (md.results.TransientSolution(147).SmbMassBalance(1)),... @@ -101,11 +104,11 @@ (md.results.TransientSolution(147).SmbDzAdd(1)),... (md.results.TransientSolution(147).SmbFAC(1)),... - (md.results.TransientSolution(end).SmbDz(1,1:230)),... - (md.results.TransientSolution(end).SmbT(1,1:230)),... - (md.results.TransientSolution(end).SmbD(1,1:230)),... - (md.results.TransientSolution(end).SmbRe(1,1:230)),... - (md.results.TransientSolution(end).SmbGdn(1,1:230)),... - (md.results.TransientSolution(end).SmbGsp(1,1:230)),... - (md.results.TransientSolution(end).SmbA(1,1:230)),... + (md.results.TransientSolution(end).SmbDz(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbT(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbD(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbRe(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbGdn(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbGsp(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbA(1,1:nlayers)),... (md.results.TransientSolution(end).SmbEC(1)),... (md.results.TransientSolution(end).SmbMassBalance(1)),... Index: /issm/trunk-jpl/test/NightlyRun/test252.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test252.py (revision 24734) +++ /issm/trunk-jpl/test/NightlyRun/test252.py (revision 24735) @@ -69,23 +69,26 @@ md = solve(md, 'Transient') +nlayers=np.size(md.results.TransientSolution[-1].SmbT,1) + #Fields and tolerances to track changes -field_names = ['SmbDz1', 'SmbT1', 'SmbD1', 'SmbRe1', 'SmbGdn1', 'SmbGsp1', 'SmbA1', 'SmbEC1', 'SmbMassBalance1', 'SmbMAdd1', 'SmbDzAdd1', 'SmbFAC1', +field_names = ['Layers','SmbDz1', 'SmbT1', 'SmbD1', 'SmbRe1', 'SmbGdn1', 'SmbGsp1', 'SmbA1', 'SmbEC1', 'SmbMassBalance1', 'SmbMAdd1', 'SmbDzAdd1', 'SmbFAC1', 'SmbDz2', 'SmbT2', 'SmbD2', 'SmbRe2', 'SmbGdn2', 'SmbGsp2', 'SmbA2', 'SmbEC2', 'SmbMassBalance2', 'SmbMAdd2', 'SmbDzAdd2', 'SmbFAC2', 'SmbDz3', 'SmbT3', 'SmbD3', 'SmbRe3', 'SmbGdn3', 'SmbGsp3', 'SmbA3', 'SmbEC3', 'SmbMassBalance3', 'SmbMAdd3', 'SmbDzAdd3', 'SmbFAC3', 'SmbDz4', 'SmbT4', 'SmbD4', 'SmbRe4', 'SmbGdn4', 'SmbGsp4', 'SmbA4', 'SmbEC4', 'SmbMassBalance4', 'SmbMAdd4', 'SmbDzAdd4', 'SmbFAC4'] -field_tolerances = [1e-11, 1e-12, 1e-11, 2e-11, 1e-11, 1e-11, 1e-12, 1e-12, 1e-12, 1e-12, 1e-12, 1e-11, - 1e-11, 1e-12, 1e-11, 8e-10, 1e-11, 4e-11, 1e-12, 4e-12, 1e-12, 1e-12, 1e-12, 2e-11, - 1e-11, 1e-12, 1e-11, 8e-10, 1e-11, 4e-11, 1e-12, 6e-12, 6e-12, 1e-12, 1e-12, 2e-11, - 2e-9, 1e-11, 3e-9, 1e-9, 1e-11, 2e-7, 7e-11, 2e-12, 3e-9, 1e-12, 1e-12, 9e-11] +field_tolerances = [1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12, + 1e-12,1e-12,1e-11,1e-10,1e-11,1e-11,1e-12,4e-12,1e-12,1e-12,1e-12,1e-11, + 1e-12,1e-12,2e-12,2e-11,1e-11,1e-11,1e-12,1e-11,1e-11,1e-12,1e-12,1e-11, + 1e-11,1e-11,4e-11,1e-11,1e-12,3e-11,1e-12,1e-12,1e-10,1e-12,1e-12,2e-11] + #shape is different in python solution (fixed using reshape) which can cause test failure: -field_values = [md.results.TransientSolution[0].SmbDz[0, 0:230].reshape(1, -1), - md.results.TransientSolution[0].SmbT[0, 0:230].reshape(1, -1), - md.results.TransientSolution[0].SmbD[0, 0:230].reshape(1, -1), - md.results.TransientSolution[0].SmbRe[0, 0:230].reshape(1, -1), - md.results.TransientSolution[0].SmbGdn[0, 0:230].reshape(1, -1), - md.results.TransientSolution[0].SmbGsp[0, 0:230].reshape(1, -1), - md.results.TransientSolution[0].SmbA[0, 0:230].reshape(1, -1), +field_values = [nlayers,md.results.TransientSolution[0].SmbDz[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbT[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbD[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbRe[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbGdn[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbGsp[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbA[0, 0:nlayers].reshape(1, -1), md.results.TransientSolution[0].SmbEC[0], md.results.TransientSolution[0].SmbMassBalance[0], @@ -93,11 +96,11 @@ md.results.TransientSolution[0].SmbDzAdd[0], md.results.TransientSolution[0].SmbFAC[0], - md.results.TransientSolution[145].SmbDz[0, 0:230].reshape(1, -1), - md.results.TransientSolution[145].SmbT[0, 0:230].reshape(1, -1), - md.results.TransientSolution[145].SmbD[0, 0:230].reshape(1, -1), - md.results.TransientSolution[145].SmbRe[0, 0:230].reshape(1, -1), - md.results.TransientSolution[145].SmbGdn[0, 0:230].reshape(1, -1), - md.results.TransientSolution[145].SmbGsp[0, 0:230].reshape(1, -1), - md.results.TransientSolution[145].SmbA[0, 0:230].reshape(1, -1), + md.results.TransientSolution[145].SmbDz[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbT[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbD[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbRe[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbGdn[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbGsp[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbA[0, 0:nlayers].reshape(1, -1), md.results.TransientSolution[145].SmbEC[0], md.results.TransientSolution[145].SmbMassBalance[0], @@ -105,11 +108,11 @@ md.results.TransientSolution[145].SmbDzAdd[0], md.results.TransientSolution[145].SmbFAC[0], - md.results.TransientSolution[146].SmbDz[0, 0:230].reshape(1, -1), - md.results.TransientSolution[146].SmbT[0, 0:230].reshape(1, -1), - md.results.TransientSolution[146].SmbD[0, 0:230].reshape(1, -1), - md.results.TransientSolution[146].SmbRe[0, 0:230].reshape(1, -1), - md.results.TransientSolution[146].SmbGdn[0, 0:230].reshape(1, -1), - md.results.TransientSolution[146].SmbGsp[0, 0:230].reshape(1, -1), - md.results.TransientSolution[146].SmbA[0, 0:230].reshape(1, -1), + md.results.TransientSolution[146].SmbDz[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbT[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbD[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbRe[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbGdn[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbGsp[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbA[0, 0:nlayers].reshape(1, -1), md.results.TransientSolution[146].SmbEC[0], md.results.TransientSolution[146].SmbMassBalance[0], @@ -117,11 +120,11 @@ md.results.TransientSolution[146].SmbDzAdd[0], md.results.TransientSolution[146].SmbFAC[0], - md.results.TransientSolution[-1].SmbDz[0, 0:230].reshape(1, -1), - md.results.TransientSolution[-1].SmbT[0, 0:230].reshape(1, -1), - md.results.TransientSolution[-1].SmbD[0, 0:230].reshape(1, -1), - md.results.TransientSolution[-1].SmbRe[0, 0:230].reshape(1, -1), - md.results.TransientSolution[-1].SmbGdn[0, 0:230].reshape(1, -1), - md.results.TransientSolution[-1].SmbGsp[0, 0:230].reshape(1, -1), - md.results.TransientSolution[-1].SmbA[0, 0:230].reshape(1, -1), + md.results.TransientSolution[-1].SmbDz[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbT[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbD[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbRe[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbGdn[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbGsp[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbA[0, 0:nlayers].reshape(1, -1), md.results.TransientSolution[-1].SmbEC[0], md.results.TransientSolution[-1].SmbMassBalance[0], Index: /issm/trunk-jpl/test/NightlyRun/test253.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test253.m (revision 24735) +++ /issm/trunk-jpl/test/NightlyRun/test253.m (revision 24735) @@ -0,0 +1,121 @@ +%Test Name: SquareShelfSMBGembClimConstrainT +md=triangle(model(),'../Exp/Square.exp',350000.); +md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareShelf.par'); +md=setflowequation(md,'SSA','all'); +md.materials.rho_ice=910; +md.cluster=generic('name',oshostname(),'np',3); + +% Use of Gemb method for SMB computation +md.smb = SMBgemb(md.mesh,md.geometry); +md.smb.dsnowIdx = 3; +md.smb.aIdx = 2; +md.smb.denIdx = 1; + +%load hourly surface forcing date from 1979 to 2009: +inputs=load('../Data/gemb_input.mat'); + +%setup the inputs: +md.smb.Ta=[repmat(inputs.Ta0',md.mesh.numberofelements,1);inputs.dateN']; +md.smb.V=[repmat(inputs.V0',md.mesh.numberofelements,1);inputs.dateN']; +md.smb.dswrf=[repmat(inputs.dsw0',md.mesh.numberofelements,1);inputs.dateN']; +md.smb.dlwrf=[repmat(inputs.dlw0',md.mesh.numberofelements,1);inputs.dateN']; +md.smb.P=[repmat(inputs.P0',md.mesh.numberofelements,1);inputs.dateN']; +md.smb.eAir=[repmat(inputs.eAir0',md.mesh.numberofelements,1);inputs.dateN']; +md.smb.pAir=[repmat(inputs.pAir0',md.mesh.numberofelements,1);inputs.dateN']; +md.smb.Vz=repmat(inputs.LP.Vz,md.mesh.numberofelements,1); +md.smb.Tz=repmat(inputs.LP.Tz,md.mesh.numberofelements,1); +md.smb.Tmean=repmat(inputs.LP.Tmean,md.mesh.numberofelements,1); +md.smb.C=repmat(inputs.LP.C,md.mesh.numberofelements,1); + +md.smb.Ta=md.smb.Ta(:,1:365*8); +md.smb.V=md.smb.V(:,1:365*8); +md.smb.dswrf=md.smb.dswrf(:,1:365*8); +md.smb.dlwrf=md.smb.dlwrf(:,1:365*8); +md.smb.P=md.smb.P(:,1:365*8); +md.smb.eAir=md.smb.eAir(:,1:365*8); +md.smb.pAir=md.smb.pAir(:,1:365*8); + +md.smb.isclimatology=1; +md.smb.isconstrainsurfaceT=1; + +%smb settings +md.smb.requested_outputs={'SmbDz','SmbT','SmbD','SmbRe','SmbGdn','SmbGsp','SmbEC','SmbA','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC'}; + +%only run smb core: +md.transient.isstressbalance=0; +md.transient.ismasstransport=0; +md.transient.isthermal=0; + +%time stepping: +md.timestepping.start_time=1965.6; +md.timestepping.final_time=1966.6; +md.timestepping.time_step=1/365.0; +md.timestepping.interp_forcings=0; + +%Run transient +md=solve(md,'Transient'); + +nlayers=size(md.results.TransientSolution(1).SmbT,2); + +%Fields and tolerances to track changes +field_names ={'Layers','SmbDz1','SmbT1' ,'SmbD1' ,'SmbRe1','SmbGdn1','SmbGsp1','SmbA1' ,'SmbEC1','SmbMassBalance1','SmbMAdd1','SmbDzAdd1','SmbFAC1',... + 'SmbDz2','SmbT2' ,'SmbD2' ,'SmbRe2','SmbGdn2','SmbGsp2','SmbA2' ,'SmbEC2','SmbMassBalance2','SmbMAdd2','SmbDzAdd2','SmbFAC2',... + 'SmbDz3','SmbT3' ,'SmbD3' ,'SmbRe3','SmbGdn3','SmbGsp3','SmbA3' ,'SmbEC3','SmbMassBalance3','SmbMAdd3','SmbDzAdd3','SmbFAC3',... + 'SmbDz4','SmbT4' ,'SmbD4' ,'SmbRe4','SmbGdn4','SmbGsp4','SmbA4' ,'SmbEC4','SmbMassBalance4','SmbMAdd4','SmbDzAdd4','SmbFAC4'}; +field_tolerances ={1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,... + 1e-12,1e-12,1e-11,1e-10,4e-11,1e-11,1e-12,1e-11,1e-12,1e-12,1e-12,1e-11,... + 1e-12,1e-12,2e-12,2e-11,1e-10,1e-11,1e-12,1e-11,1e-11,1e-12,1e-12,1e-11,... + 1e-11,1e-11,4e-11,1e-11,1e-12,3e-11,1e-12,4e-12,1e-10,1e-12,1e-12,2e-11}; + +field_values={... + (nlayers)... + (md.results.TransientSolution(1).SmbDz(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbT(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbD(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbRe(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbGdn(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbGsp(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbA(1,1:nlayers)),... + (md.results.TransientSolution(1).SmbEC(1)),... + (md.results.TransientSolution(1).SmbMassBalance(1)),... + (md.results.TransientSolution(1).SmbMAdd(1)),... + (md.results.TransientSolution(1).SmbDzAdd(1)),... + (md.results.TransientSolution(1).SmbFAC(1)),... + (md.results.TransientSolution(146).SmbDz(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbT(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbD(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbRe(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbGdn(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbGsp(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbA(1,1:nlayers)),... + (md.results.TransientSolution(146).SmbEC(1)),... + (md.results.TransientSolution(146).SmbMassBalance(1)),... + (md.results.TransientSolution(146).SmbMAdd(1)),... + (md.results.TransientSolution(146).SmbDzAdd(1)),... + (md.results.TransientSolution(146).SmbFAC(1)),... + (md.results.TransientSolution(147).SmbDz(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbT(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbD(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbRe(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbGdn(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbGsp(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbA(1,1:nlayers)),... + (md.results.TransientSolution(147).SmbEC(1)),... + (md.results.TransientSolution(147).SmbMassBalance(1)),... + (md.results.TransientSolution(147).SmbMAdd(1)),... + (md.results.TransientSolution(147).SmbDzAdd(1)),... + (md.results.TransientSolution(147).SmbFAC(1)),... + (md.results.TransientSolution(end).SmbDz(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbT(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbD(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbRe(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbGdn(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbGsp(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbA(1,1:nlayers)),... + (md.results.TransientSolution(end).SmbEC(1)),... + (md.results.TransientSolution(end).SmbMassBalance(1)),... + (md.results.TransientSolution(end).SmbMAdd(1)),... + (md.results.TransientSolution(end).SmbDzAdd(1)),... + (md.results.TransientSolution(end).SmbFAC(1)),... + }; Index: /issm/trunk-jpl/test/NightlyRun/test253.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test253.py (revision 24735) +++ /issm/trunk-jpl/test/NightlyRun/test253.py (revision 24735) @@ -0,0 +1,135 @@ +#Test Name: SquareShelfSMBGembClimConstrainT +import numpy as np +import sys +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from SMBgemb import * + +md = triangle(model(), '../Exp/Square.exp', 350000.) +md = setmask(md, 'all', '') +md = parameterize(md, '../Par/SquareShelf.py') +md = setflowequation(md, 'SSA', 'all') +md.materials.rho_ice = 910 +md.cluster = generic('name', gethostname(), 'np', 3) + +#Use of Gemb method for SMB computation +md.smb = SMBgemb() +md.smb.setdefaultparameters(md.mesh, md.geometry) +md.smb.dsnowIdx = 3 +md.smb.aIdx = 2 +md.smb.denIdx = 1 + +#load hourly surface forcing date from 1979 to 2009: +if sys.version_info.major == 2: + inputs = np.load('../Data/gemb_input.npy', allow_pickle=True).item() +else: + inputs = np.load('../Data/gemb_input.npy', allow_pickle=True, encoding='bytes').item() + +#setup the inputs: +md.smb.Ta = np.append(np.tile(np.conjugate(inputs[b'Ta0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0) +md.smb.V = np.append(np.tile(np.conjugate(inputs[b'V0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0) +md.smb.dswrf = np.append(np.tile(np.conjugate(inputs[b'dsw0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0) +md.smb.dlwrf = np.append(np.tile(np.conjugate(inputs[b'dlw0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0) +md.smb.P = np.append(np.tile(np.conjugate(inputs[b'P0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0) +md.smb.eAir = np.append(np.tile(np.conjugate(inputs[b'eAir0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0) +md.smb.pAir = np.append(np.tile(np.conjugate(inputs[b'pAir0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0) +md.smb.Vz = np.tile(np.conjugate(inputs[b'LP']['Vz']), (md.mesh.numberofelements, 1)).flatten() +md.smb.Tz = np.tile(np.conjugate(inputs[b'LP']['Tz']), (md.mesh.numberofelements, 1)).flatten() +md.smb.Tmean = np.tile(np.conjugate(inputs[b'LP']['Tmean']), (md.mesh.numberofelements, 1)).flatten() +md.smb.C = np.tile(np.conjugate(inputs[b'LP']['C']), (md.mesh.numberofelements, 1)).flatten() + +md.smb.Ta = md.smb.Ta[:, 0:365 * 8] +md.smb.V = md.smb.V[:, 0:365 * 8] +md.smb.dswrf = md.smb.dswrf[:, 0:365 * 8] +md.smb.dlwrf = md.smb.dlwrf[:, 0:365 * 8] +md.smb.P = md.smb.P[:, 0:365 * 8] +md.smb.eAir = md.smb.eAir[:, 0:365 * 8] +md.smb.pAir = md.smb.pAir[:, 0:365 * 8] + +md.smb.isclimatology = 1 +md.smb.isconstrainsurfaceT = 1 + +#smb settings +md.smb.requested_outputs = ['SmbDz', 'SmbT', 'SmbD', 'SmbRe', 'SmbGdn', 'SmbGsp', 'SmbEC', 'SmbA', 'SmbMassBalance', 'SmbMAdd', 'SmbDzAdd', 'SmbFAC'] + +#only run smb core: +md.transient.isstressbalance = 0 +md.transient.ismasstransport = 0 +md.transient.isthermal = 0 + +#time stepping: +md.timestepping.start_time = 1965.6 +md.timestepping.final_time = 1966.6 +md.timestepping.time_step = 1. / 365.0 +md.timestepping.interp_forcings = 0. + +#Run transient +md = solve(md, 'Transient') + +nlayers=np.size(md.results.TransientSolution[0].SmbT,1) + +#Fields and tolerances to track changes + +field_names = ['Layers','SmbDz1', 'SmbT1', 'SmbD1', 'SmbRe1', 'SmbGdn1', 'SmbGsp1', 'SmbA1', 'SmbEC1', 'SmbMassBalance1', 'SmbMAdd1', 'SmbDzAdd1', 'SmbFAC1', + 'SmbDz2', 'SmbT2', 'SmbD2', 'SmbRe2', 'SmbGdn2', 'SmbGsp2', 'SmbA2', 'SmbEC2', 'SmbMassBalance2', 'SmbMAdd2', 'SmbDzAdd2', 'SmbFAC2', + 'SmbDz3', 'SmbT3', 'SmbD3', 'SmbRe3', 'SmbGdn3', 'SmbGsp3', 'SmbA3', 'SmbEC3', 'SmbMassBalance3', 'SmbMAdd3', 'SmbDzAdd3', 'SmbFAC3', + 'SmbDz4', 'SmbT4', 'SmbD4', 'SmbRe4', 'SmbGdn4', 'SmbGsp4', 'SmbA4', 'SmbEC4', 'SmbMassBalance4', 'SmbMAdd4', 'SmbDzAdd4', 'SmbFAC4'] +field_tolerances = [1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12, + 1e-12,1e-12,1e-11,1e-10,4e-11,1e-11,1e-12,1e-11,1e-12,1e-12,1e-12,1e-11, + 1e-12,1e-12,2e-12,2e-11,1e-10,1e-11,1e-12,1e-11,1e-11,1e-12,1e-12,1e-11, + 1e-11,1e-11,4e-11,1e-11,1e-12,3e-11,1e-12,4e-12,1e-10,1e-12,1e-12,2e-11] + +#shape is different in python solution (fixed using reshape) which can cause test failure: +field_values = [nlayers,md.results.TransientSolution[0].SmbDz[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbT[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbD[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbRe[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbGdn[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbGsp[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbA[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[0].SmbEC[0], + md.results.TransientSolution[0].SmbMassBalance[0], + md.results.TransientSolution[0].SmbMAdd[0], + md.results.TransientSolution[0].SmbDzAdd[0], + md.results.TransientSolution[0].SmbFAC[0], + md.results.TransientSolution[145].SmbDz[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbT[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbD[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbRe[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbGdn[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbGsp[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbA[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[145].SmbEC[0], + md.results.TransientSolution[145].SmbMassBalance[0], + md.results.TransientSolution[145].SmbMAdd[0], + md.results.TransientSolution[145].SmbDzAdd[0], + md.results.TransientSolution[145].SmbFAC[0], + md.results.TransientSolution[146].SmbDz[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbT[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbD[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbRe[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbGdn[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbGsp[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbA[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[146].SmbEC[0], + md.results.TransientSolution[146].SmbMassBalance[0], + md.results.TransientSolution[146].SmbMAdd[0], + md.results.TransientSolution[146].SmbDzAdd[0], + md.results.TransientSolution[146].SmbFAC[0], + md.results.TransientSolution[-1].SmbDz[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbT[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbD[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbRe[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbGdn[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbGsp[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbA[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbEC[0], + md.results.TransientSolution[-1].SmbMassBalance[0], + md.results.TransientSolution[-1].SmbMAdd[0], + md.results.TransientSolution[-1].SmbDzAdd[0], + md.results.TransientSolution[-1].SmbFAC[0]]