Changeset 26810

Timestamp:

01/25/22 05:41:26 (3 years ago)

Author:

vverjans

Message:

NEW: added SMB lapse rate computation per basin to SMBautoregression

Location:

issm/trunk-jpl

Files:

• src/c/classes/Elements/Element.cpp (modified) (1 diff)
• src/c/classes/Elements/Element.h (modified) (1 diff)
• src/c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (3 diffs)
• src/c/shared/Enum/Enum.vim (modified) (3 diffs)
• src/c/shared/Enum/EnumDefinitions.h (modified) (2 diffs)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (2 diffs)
• src/c/shared/Enum/Enumjl.vim (modified) (2 diffs)
• src/m/classes/SMBautoregression.m (modified) (5 diffs)
• src/m/classes/SMBautoregression.py (modified) (6 diffs)
• test/Archives/Archive257.arch (modified) ( previous)
• test/NightlyRun/test257.m (modified) (1 diff)
• test/NightlyRun/test257.py (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -2256 +2256 @@
         xDelete<IssmDouble>(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<IssmDouble>(numvertices);
+   IssmDouble* smblist     = xNew<IssmDouble>(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<numvertices;v++){
+      /*Find ELA*/
+                if(smblist[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<IssmDouble>(surfacelist);
+   xDelete<IssmDouble>(smblist);
 }/*}}}*/
 void       Element::LinearFloatingiceMeltingRate(){/*{{{*/

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -153 +153 @@
                 bool               IsLandInElement();
                 void               Ismip6FloatingiceMeltingRate();
+                void               LapseRateBasinSMB(IssmDouble* lapseratepos, IssmDouble* lapserateneg,IssmDouble* refelevation);
                 void               LinearFloatingiceMeltingRate();
                 void               SpatialLinearFloatingiceMeltingRate();

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -412 +412 @@
          parameters->AddObject(new DoubleMatParam(SmbPhiEnum,transparam,M,N));
          xDelete<IssmDouble>(transparam);
-         break;
+         iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserate_pos");
+         parameters->AddObject(new DoubleVecParam(SmbLapseRatePosEnum,transparam,N));
+         xDelete<IssmDouble>(transparam);
+         iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserate_neg");
+         parameters->AddObject(new DoubleVecParam(SmbLapseRateNegEnum,transparam,N));
+         xDelete<IssmDouble>(transparam);
+                        iomodel->FetchData(&transparam,&M,&N,"md.smb.refelevation");
+         parameters->AddObject(new DoubleVecParam(SmbRefElevationEnum,transparam,N));
+         xDelete<IssmDouble>(transparam);
+                        break;
                 case SMBgembEnum:
                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

@@ -202 +202 @@
    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 +234 @@
    for(Object* &object:femmodel->elements->objects){
       Element* element = xDynamicCast<Element*>(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 +244 @@
    xDelete<IssmDouble>(beta1);
    xDelete<IssmDouble>(phi);
+   xDelete<IssmDouble>(lapseratepos);
+   xDelete<IssmDouble>(lapserateneg);
+   xDelete<IssmDouble>(refelevation);
 }/*}}}*/
 void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/

issm/trunk-jpl/src/c/shared/Enum/Enum.vim

@@ -421 +421 @@
 syn keyword cConstant StochasticForcingNumFieldsEnum
 syn keyword cConstant StochasticForcingRandomflagEnum
-syn keyword cConstant RotationalPolarMoiEnum
 syn keyword cConstant SolidearthSettingsReltolEnum
 syn keyword cConstant SolidearthSettingsSelfAttractionEnum
@@ -477 +476 @@
 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 +483 @@
 syn keyword cConstant SmbPhiEnum
 syn keyword cConstant SmbRdlEnum
+syn keyword cConstant SmbRefElevationEnum
 syn keyword cConstant SmbRequestedOutputsEnum
 syn keyword cConstant SmbRlapsEnum

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -470 +470 @@
         SmbIsturbulentfluxEnum,
         SmbKEnum,
+        SmbLapseRateNegEnum,
+   SmbLapseRatePosEnum,
         SmbNumBasinsEnum,
         SmbNumRequestedOutputsEnum,
@@ -475 +477 @@
         SmbPhiEnum,
         SmbRdlEnum,
+        SmbRefElevationEnum,
         SmbRequestedOutputsEnum,
         SmbRlapsEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -478 +478 @@
                 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 +485 @@
                 case SmbPhiEnum : return "SmbPhi";
                 case SmbRdlEnum : return "SmbRdl";
+                case SmbRefElevationEnum : return "SmbRefElevation";
                 case SmbRequestedOutputsEnum : return "SmbRequestedOutputs";
                 case SmbRlapsEnum : return "SmbRlaps";

issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim

@@ -469 +469 @@
 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 +476 @@
 syn keyword juliaConstC SmbPhiEnum
 syn keyword juliaConstC SmbRdlEnum
+syn keyword juliaConstC SmbRefElevationEnum
 syn keyword juliaConstC SmbRequestedOutputsEnum
 syn keyword juliaConstC SmbRlapsEnum

issm/trunk-jpl/src/m/classes/SMBautoregression.m

@@ -14 +14 @@
                 phi               = NaN;
                 basin_id          = NaN;
+                lapserate_pos     = NaN;
+                lapserate_neg     = NaN;
+                refelevation      = NaN;
                 steps_per_step    = 1;
                 averaging         = 0;
@@ -70 +73 @@
                                 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 +99 @@
                         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 +114 @@
                         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 +150 @@
                         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');

issm/trunk-jpl/src/m/classes/SMBautoregression.py

@@ -5 +5 @@
 from project3d import *
 from WriteData import *
-
+from GetAreas import *
 
 class SMBautoregression(object):
@@ -23 +23 @@
         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 +47 @@
         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 +101 @@
             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 +116 @@
     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 +147 @@
         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')

issm/trunk-jpl/test/NightlyRun/test257.m

@@ -45 +45 @@
 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

issm/trunk-jpl/test/NightlyRun/test257.py

@@ -50 +50 @@
 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