Changeset 17087

Timestamp:

01/09/14 15:06:03 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added SMBhenning fonction, to be completed by Henning (dummy for now)

Location:

issm/trunk-jpl/src

Files:

• c/analyses/MasstransportAnalysis.cpp (modified) (1 diff)
• c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (2 diffs)
• m/classes/SMBhenning.m (modified) (4 diffs)

issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

@@ -110 +110 @@
                         break;
                 case SMBhenningEnum:
-                        iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum,0.);
+                        iomodel->FetchDataToInput(elements,SurfaceforcingsSmbrefEnum,0.);
                         break;
                 default:

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

@@ -37 +37 @@
                         if(VerboseSolution())_printf_("  call smb Henning module\n");
                         SmbHenningx(femmodel);
+                        break;
                 default:
                         _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
@@ -160 +161 @@
 void SmbHenningx(FemModel* femmodel){/*{{{*/
 
-        // void SurfaceMassBalancex(hd,agd,ni){
-        //    INPUT parameters: ni: working size of arrays
-        //    INPUT: surface elevation (m): hd(NA)
-        //    OUTPUT: mass-balance (m/yr ice): agd(NA)
-
+        /*Intermediaries*/
+        IssmDouble  z_critical = 200.;
+        IssmDouble  a = 2.;
+        IssmDouble  b = 3.;
+        IssmDouble  smb,smbref,z;
+
+        /*Loop over all the elements of this partition*/
         for(int i=0;i<femmodel->elements->Size();i++){
                 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-        }
-
-}/*}}}*/
+
+                /*Get reference SMB (uncorrected) and allocate all arrays*/
+                int         numvertices = element->GetNumberOfVertices();
+                IssmDouble* surfacelist = xNew<IssmDouble>(numvertices);
+                IssmDouble* smblistref  = xNew<IssmDouble>(numvertices);
+                IssmDouble* smblist     = xNew<IssmDouble>(numvertices);
+                element->GetInputListOnVertices(surfacelist,SurfaceEnum);
+                element->GetInputListOnVertices(smblistref,SurfaceforcingsSmbrefEnum);
+
+                /*Loop over all vertices of element and correct SMB as a function of altitude z*/
+                for(int v=0;v<numvertices;v++){
+
+                        /*Get vertex elevation, reference smb*/
+                        z      = surfacelist[v];
+                        smbref = smblistref[v];
+
+                        /*Compute corrected smb*/
+                        if(z<z_critical){
+                                smb = smbref;
+                        }
+                        else{
+                                smb = a*smbref*smbref+b+z;
+                        }
+
+                        /*Update array accordingly*/
+                        smblist[v] = smb;
+
+                }
+
+                /*Add input to element and Free memory*/
+                element->AddInput(SurfaceforcingsMassBalanceEnum,smblist,P1Enum);
+                xDelete<IssmDouble>(surfacelist);
+                xDelete<IssmDouble>(smblistref);
+                xDelete<IssmDouble>(smblist);
+        }
+
+}/*}}}*/

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

@@ -6 +6 @@
 classdef SMBhenning
         properties (SetAccess=public) 
-                mass_balance = NaN;
+                smbref = NaN;
         end
         methods
@@ -18 +18 @@
                 function self = extrude(self,md) % {{{
 
-                        self.mass_balance=project3d(md,'vector',self.mass_balance,'type','node');
+                        self.smbref=project3d(md,'vector',self.smbref,'type','node');
 
                 end % }}}
                 function self = initialize(self,md) % {{{
 
-                        if isnan(self.mass_balance)
-                                self.mass_balance=zeros(md.mesh.numberofvertices,1);
-                                disp('      no surfaceforcings.mass_balance specified: values set as zero');
+                        if isnan(self.smbref)
+                                self.smbref=zeros(md.mesh.numberofvertices,1);
+                                disp('      no surfaceforcings.smbref specified: values set as zero');
                         end
 
@@ -32 +32 @@
 
                         if ismember(MasstransportAnalysisEnum(),analyses),
-                                md = checkfield(md,'fieldname','surfaceforcings.mass_balance','forcing',1,'NaN',1);
+                                md = checkfield(md,'fieldname','surfaceforcings.smbref','forcing',1,'NaN',1);
                         end
                         if ismember(BalancethicknessAnalysisEnum(),analyses),
-                                md = checkfield(md,'fieldname','surfaceforcings.mass_balance','size',[md.mesh.numberofvertices 1],'NaN',1);
+                                md = checkfield(md,'fieldname','surfaceforcings.smbref','size',[md.mesh.numberofvertices 1],'NaN',1);
                         end
                 end % }}}
                 function disp(obj) % {{{
                         disp(sprintf('   surface forcings parameters:'));
-                        fielddisplay(obj,'mass_balance','surface mass balance [m/yr ice eq]');
+                        fielddisplay(obj,'smbref','reference smb from which deviation is calculated [m/yr ice eq]');
                 end % }}}
                 function marshall(obj,md,fid) % {{{
@@ -47 +47 @@
 
                         WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBhenningEnum(),'format','Integer');
-                        WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','mass_balance','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
+                        WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','smbref','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
                 end % }}}
         end