Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26809) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26810) @@ -2256,4 +2256,41 @@ xDelete(depths); +}/*}}}*/ +void Element::LapseRateBasinSMB(IssmDouble* lapseratepos, IssmDouble* lapserateneg,IssmDouble* refelevation){/*{{{*/ + + /*Variable declaration*/ + const int numvertices = this->GetNumberOfVertices(); + bool isadjustsmb = true; + int basinid; + IssmDouble ela,refelevation_b,lapseratepos_b,lapserateneg_b; + IssmDouble* surfacelist = xNew(numvertices); + IssmDouble* smblist = xNew(numvertices); + + /*Retrieve SMB values non-adjusted for SMB lapse rate*/ + this->GetInputListOnVertices(smblist,SmbMassBalanceEnum); + /*Get basin-specific parameters of the element*/ + this->GetInputValue(&basinid,SmbBasinsIdEnum); + refelevation_b = refelevation[basinid]; + lapseratepos_b = lapseratepos[basinid]; + lapserateneg_b = lapserateneg[basinid]; + /*Get surface elevation on vertices*/ + this->GetInputListOnVertices(surfacelist,SurfaceEnum); + + for(int v=0;v=0 && lapseratepos_b!=0) ela = refelevation_b-smblist[v]/lapseratepos_b; + else if(smblist[v]<0 && lapserateneg_b!=0) ela = refelevation_b-smblist[v]/lapserateneg_b; + /*Lapse rate is zero: SMB not adjusted*/ + else isadjustsmb = false; + /*Adjust SMB value*/ + if(isadjustsmb) smblist[v] = lapseratepos_b*max(surfacelist[v]-ela,0.0)+lapserateneg_b*min(surfacelist[v]-ela,0.0); + } + + /*Add input to element*/ + this->AddInput(SmbMassBalanceEnum,smblist,P1Enum); + + /*Cleanup*/ + xDelete(surfacelist); + xDelete(smblist); }/*}}}*/ void Element::LinearFloatingiceMeltingRate(){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26809) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26810) @@ -153,4 +153,5 @@ bool IsLandInElement(); void Ismip6FloatingiceMeltingRate(); + void LapseRateBasinSMB(IssmDouble* lapseratepos, IssmDouble* lapserateneg,IssmDouble* refelevation); void LinearFloatingiceMeltingRate(); void SpatialLinearFloatingiceMeltingRate(); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26809) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26810) @@ -412,5 +412,14 @@ parameters->AddObject(new DoubleMatParam(SmbPhiEnum,transparam,M,N)); xDelete(transparam); - break; + iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserate_pos"); + parameters->AddObject(new DoubleVecParam(SmbLapseRatePosEnum,transparam,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserate_neg"); + parameters->AddObject(new DoubleVecParam(SmbLapseRateNegEnum,transparam,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.smb.refelevation"); + parameters->AddObject(new DoubleVecParam(SmbRefElevationEnum,transparam,N)); + xDelete(transparam); + break; case SMBgembEnum: parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum)); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26809) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26810) @@ -202,12 +202,18 @@ femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum); IssmDouble tinit_ar; - IssmDouble* beta0 = NULL; - IssmDouble* beta1 = NULL; - IssmDouble* phi = NULL; + IssmDouble* beta0 = NULL; + IssmDouble* beta1 = NULL; + IssmDouble* phi = NULL; + IssmDouble* lapseratepos = NULL; + IssmDouble* lapserateneg = NULL; + IssmDouble* refelevation = NULL; femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum); - femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); + femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); + femmodel->parameters->FindParam(&lapseratepos,&M,SmbLapseRatePosEnum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&lapserateneg,&M,SmbLapseRateNegEnum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum); _assert_(M==numbasins); femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); @@ -228,6 +234,9 @@ for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum); - } + /*Compute autoregression*/ + element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum); + /*Compute lapse rate adjustment*/ + element->LapseRateBasinSMB(lapseratepos,lapserateneg,refelevation); + } /*Cleanup*/ @@ -235,4 +244,7 @@ xDelete(beta1); xDelete(phi); + xDelete(lapseratepos); + xDelete(lapserateneg); + xDelete(refelevation); }/*}}}*/ void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26809) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26810) @@ -421,5 +421,4 @@ syn keyword cConstant StochasticForcingNumFieldsEnum syn keyword cConstant StochasticForcingRandomflagEnum -syn keyword cConstant RotationalPolarMoiEnum syn keyword cConstant SolidearthSettingsReltolEnum syn keyword cConstant SolidearthSettingsSelfAttractionEnum @@ -477,4 +476,6 @@ syn keyword cConstant SmbIsturbulentfluxEnum syn keyword cConstant SmbKEnum +syn keyword cConstant SmbLapseRateNegEnum +syn keyword cConstant SmbLapseRatePosEnum syn keyword cConstant SmbNumBasinsEnum syn keyword cConstant SmbNumRequestedOutputsEnum @@ -482,4 +483,5 @@ syn keyword cConstant SmbPhiEnum syn keyword cConstant SmbRdlEnum +syn keyword cConstant SmbRefElevationEnum syn keyword cConstant SmbRequestedOutputsEnum syn keyword cConstant SmbRlapsEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26809) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26810) @@ -470,4 +470,6 @@ SmbIsturbulentfluxEnum, SmbKEnum, + SmbLapseRateNegEnum, + SmbLapseRatePosEnum, SmbNumBasinsEnum, SmbNumRequestedOutputsEnum, @@ -475,4 +477,5 @@ SmbPhiEnum, SmbRdlEnum, + SmbRefElevationEnum, SmbRequestedOutputsEnum, SmbRlapsEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26809) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26810) @@ -478,4 +478,6 @@ case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux"; case SmbKEnum : return "SmbK"; + case SmbLapseRateNegEnum : return "SmbLapseRateNeg"; + case SmbLapseRatePosEnum : return "SmbLapseRatePos"; case SmbNumBasinsEnum : return "SmbNumBasins"; case SmbNumRequestedOutputsEnum : return "SmbNumRequestedOutputs"; @@ -483,4 +485,5 @@ case SmbPhiEnum : return "SmbPhi"; case SmbRdlEnum : return "SmbRdl"; + case SmbRefElevationEnum : return "SmbRefElevation"; case SmbRequestedOutputsEnum : return "SmbRequestedOutputs"; case SmbRlapsEnum : return "SmbRlaps"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26809) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26810) @@ -469,4 +469,6 @@ syn keyword juliaConstC SmbIsturbulentfluxEnum syn keyword juliaConstC SmbKEnum +syn keyword juliaConstC SmbLapseRateNegEnum +syn keyword juliaConstC SmbLapseRatePosEnum syn keyword juliaConstC SmbNumBasinsEnum syn keyword juliaConstC SmbNumRequestedOutputsEnum @@ -474,4 +476,5 @@ syn keyword juliaConstC SmbPhiEnum syn keyword juliaConstC SmbRdlEnum +syn keyword juliaConstC SmbRefElevationEnum syn keyword juliaConstC SmbRequestedOutputsEnum syn keyword juliaConstC SmbRlapsEnum Index: /issm/trunk-jpl/src/m/classes/SMBautoregression.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBautoregression.m (revision 26809) +++ /issm/trunk-jpl/src/m/classes/SMBautoregression.m (revision 26810) @@ -14,4 +14,7 @@ phi = NaN; basin_id = NaN; + lapserate_pos = NaN; + lapserate_neg = NaN; + refelevation = NaN; steps_per_step = 1; averaging = 0; @@ -70,4 +73,15 @@ md = checkfield(md,'fieldname','smb.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep md = checkfield(md,'fieldname','smb.phi','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]); + + if (any(isnan(md.smb.lapserate_pos)==0) || numel(md.smb.lapserate_pos)>1) + md = checkfield(md,'fieldname','smb.lapserate_pos','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); + end + if (any(isnan(md.smb.lapserate_neg)==0) || numel(md.smb.lapserate_neg)>1) + md = checkfield(md,'fieldname','smb.lapserate_neg','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); + end + if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1) + md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); + end + end md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]); @@ -85,4 +99,7 @@ fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]'); fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]'); + fielddisplay(self,'lapserate_pos','basin-specific SMB lapse rates applied in range of SMB>=0 [m ice eq yr^-1 m^-1] (default: no lapse rate)'); + fielddisplay(self,'lapserate_neg','basin-specific SMB lapse rates applied in range of SMB<0 [m ice eq yr^-1 m^-1] (default: no lapse rate'); + fielddisplay(self,'refelevation','basin-specific reference elevations on which SMB lapse rates are applied (default: basin mean elevation) [m]'); fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'); fielddisplay(self, 'averaging', 'averaging methods from short to long steps'); @@ -97,4 +114,31 @@ yts=md.constants.yts; + templapserate_pos = md.smb.lapserate_pos; + templapserate_neg = md.smb.lapserate_neg; + temprefelevation = md.smb.refelevation; + if(any(isnan(md.smb.lapserate_pos))) + templapserate_pos = zeros(1,md.smb.num_basins); + disp(' smb.lapserate_pos not specified: set to 0'); + end + if(any(isnan(md.smb.lapserate_neg))) + templapserate_neg = zeros(1,md.smb.num_basins); + disp(' smb.lapserate_neg not specified: set to 0'); + end + if(any(isnan(md.smb.refelevation))) + temprefelevation = zeros(1,md.smb.num_basins); + areas = GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y); + for ii=1:md.smb.num_basins + indices = find(md.smb.basin_id==ii); + elemsh = zeros(numel(indices),1); + for jj=1:numel(indices) + elemsh(jj) = mean(md.geometry.surface(md.mesh.elements(indices(jj),:))); + end + temprefelevation(ii) = sum(areas(indices).*elemsh)/sum(areas(indices)); + end + if(any([templapserate_pos,templapserate_neg])~=0) + disp(' smb.refelevation not specified: Reference elevations set to mean surface elevation of basins'); + end + end + WriteData(fid,prefix,'name','md.smb.model','data',55,'format','Integer'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','num_basins','format','Integer'); @@ -106,4 +150,7 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta1','format','DoubleMat','name','md.smb.beta1','scale',1./(yts^2),'yts',yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',yts); + WriteData(fid,prefix,'data',templapserate_pos,'format','DoubleMat','name','md.smb.lapserate_pos','scale',1./yts,'yts',yts); + WriteData(fid,prefix,'data',templapserate_neg,'format','DoubleMat','name','md.smb.lapserate_neg','scale',1./yts,'yts',yts); + WriteData(fid,prefix,'data',temprefelevation,'format','DoubleMat','name','md.smb.refelevation'); WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer'); Index: /issm/trunk-jpl/src/m/classes/SMBautoregression.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBautoregression.py (revision 26809) +++ /issm/trunk-jpl/src/m/classes/SMBautoregression.py (revision 26810) @@ -5,5 +5,5 @@ from project3d import * from WriteData import * - +from GetAreas import * class SMBautoregression(object): @@ -23,4 +23,7 @@ self.phi = np.nan self.basin_id = np.nan + self.lapserate_pos = np.nan + self.lapserate_neg = np.nan + self.refelevation = np.nan self.steps_per_step = 1 self.averaging = 0 @@ -44,6 +47,7 @@ s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]')) s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'covmat', 'inter-basin covariance matrix for multivariate normal noise at each time step [m^2 ice eq. yr^(-2)]')) - s += '{}\n'.format(fielddisplay(self, 'randomflag', 'whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)')) + s += '{}\n'.format(fielddisplay(self, 'lapserate_pos', 'basin-specific SMB lapse rates applied in range of SMB>=0 [m ice eq yr^-1 m^-1] (default: no lapse rate)')) + s += '{}\n'.format(fielddisplay(self, 'lapserate_neg', 'basin-specific SMB lapse rates applied in range of SMB<0 [m ice eq yr^-1 m^-1] (default: no lapse rate)')) + s += '{}\n'.format(fielddisplay(self, 'refelevation', 'basin-specific reference elevations on which SMB lapse rates are applied (default: basin mean elevation) [m]')) s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step')) s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) @@ -97,4 +101,11 @@ md = checkfield(md, 'fieldname', 'smb.ar_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep md = checkfield(md, 'fieldname', 'smb.phi', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ar_order]) + if(np.any(np.isnan(md.smb.lapserate_pos)==False) or np.size(md.smb.lapserate_pos)>1): + md = checkfield(md, 'fieldname', 'smb.lapserate_pos', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) + if(np.any(np.isnan(md.smb.lapserate_neg)==False) or np.size(md.smb.lapserate_neg)>1): + md = checkfield(md, 'fieldname', 'smb.lapserate_neg', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) + if(np.any(np.isnan(md.smb.refelevation)==False) or np.size(md.smb.refelevation)>1): + md = checkfield(md, 'fieldname', 'smb.refelevation', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) + md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1]) md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2]) @@ -105,4 +116,25 @@ def marshall(self, prefix, md, fid): # {{{ yts = md.constants.yts + + templapserate_pos = np.copy(md.smb.lapserate_pos) + templapserate_neg = np.copy(md.smb.lapserate_neg) + temprefelevation = np.copy(md.smb.lapserate_neg) + if(np.any(np.isnan(md.smb.lapserate_pos))): + templapserate_pos = np.zeros((md.smb.num_basins)) + print(' smb.lapserate_pos not specified: set to 0') + if(np.any(np.isnan(md.smb.lapserate_neg))): + templapserate_neg = np.zeros((md.smb.num_basins)) + print(' smb.lapserate_neg not specified: set to 0') + if(np.any(np.isnan(md.smb.refelevation))): + temprefelevation = np.zeros((md.smb.num_basins)) + areas = GetAreas(md.mesh.elements, md.mesh.x, md.mesh.y) + for ii in range(int(md.smb.num_basins)): + indices = np.where(md.smb.basin_id==ii)[0] + elemsh = np.zeros((len(indices))) + for jj in range(len(indices)): + elemsh[jj] = np.mean(md.geometry.surface[md.mesh.elements[indices[jj,:]]]) + temprefelevation[ii] = np.sum(areas[indices]*elemsh)/np.sum(areas[indices]) + if(np.any(templapserate_pos!=0) or np.any(templapserate_neg!=0)): + print(' smb.refelevation not specified: Reference elevations set to mean surface elevation of basins') WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 55, 'format', 'Integer') @@ -115,4 +147,7 @@ WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.smb.beta1', 'scale', 1 / (yts ** 2), 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.smb.phi', 'yts', yts) + WriteData(fid, prefix, 'data', templapserate_pos, 'name', 'md.smb.lapserate_pos', 'format', 'DoubleMat','scale',1/yts,'yts',yts) + WriteData(fid, prefix, 'data', templapserate_neg, 'name', 'md.smb.lapserate_neg', 'format', 'DoubleMat','scale',1/yts,'yts',yts) + WriteData(fid, prefix, 'data', temprefelevation, 'name', 'md.smb.refelevation', 'format', 'DoubleMat') WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer') Index: /issm/trunk-jpl/test/NightlyRun/test257.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26809) +++ /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26810) @@ -45,4 +45,6 @@ md.smb.ar_timestep = 2.0; %timestep of the autoregressive model [yr] md.smb.phi = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1];[0.4,-0.4,0.1,-0.1]]; +md.smb.lapserate_pos = [0.01,0.0,0.0]; +md.smb.lapserate_neg = [0.01,0.0,0.0]; %Stochastic forcing Index: /issm/trunk-jpl/test/NightlyRun/test257.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.py (revision 26809) +++ /issm/trunk-jpl/test/NightlyRun/test257.py (revision 26810) @@ -50,4 +50,6 @@ md.smb.ar_timestep = 2.0 #timestep of the autoregressive model [yr] md.smb.phi = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1], [0.4, -0.4, 0.1, -0.1]]) +md.smb.lapserate_pos = np.array([[0.01,0,0]]) +md.smb.lapserate_neg = np.array([[0.01,0,0]]) # Stochastic forcing