Changeset 17403

Timestamp:

03/10/14 09:14:38 (11 years ago)

Author:

hakesson

Message:

CHG: updating Henning's surface mass balance scheme

File:

• issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (2 diffs)

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

@@ -162 +162 @@
 
         /*Intermediaries*/
-        IssmDouble  z_critical = 200.;
-        IssmDouble  a = 2.;
-        IssmDouble  b = 3.;
-        IssmDouble  smb,smbref,z;
-
+        IssmDouble  z_critical = 1675.;
+        IssmDouble  dz = 0;
+        IssmDouble  a = -15.86;
+        IssmDouble  b = 0.00969;
+        IssmDouble  c = -0.235;
+        IssmDouble  f = 1.;
+        IssmDouble  g = -0.0011;
+        IssmDouble  h = -1.54e-5;
+        IssmDouble  smb,smbref,anomaly,yts,z;
+    
+    /* Get constants */
+    femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
+    /*iomodel->Constant(&yts,ConstantsYtsEnum);*/
+    /*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
+    /*Mathieu original*/
+    /*IssmDouble  smb,smbref,z;*/
+    
         /*Loop over all the elements of this partition*/
         for(int i=0;i<femmodel->elements->Size();i++){
@@ -182 +194 @@
                 for(int v=0;v<numvertices;v++){
 
-                        /*Get vertex elevation, reference smb*/
+                        /*Get vertex elevation, anoma smb*/
                         z      = surfacelist[v];
-                        smbref = smblistref[v];
-
-                        /*Compute corrected smb*/
-                        if(z<z_critical){
-                                smb = smbref;
+                        anomaly = smblistref[v];
+
+            /* Henning edited acc. to Riannes equations*/
+            /* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
+            z_critical = z_critical + dz;
+            
+            /* Calculate smb acc. to the surface elevation z */
+            if(z<z_critical){
+                                smb = a + b*z + c;
                         }
                         else{
-                                smb = a*smbref*smbref+b+z;
+                          smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
                         }
-
+            
+            /* Compute smb including anomaly,
+             correct for number of seconds in a year [s/yr]*/
+            smb = smb + anomaly*yts;
+            
                         /*Update array accordingly*/
                         smblist[v] = smb;