source: issm/oecreview/Archive/23390-24306/ISSM-24063-24064.diff@ 24307

../trunk-jpl/src/c/classes/Loads/Channel.cpp

@@ -54 +54 @@
         int e1 = iomodel->faces[4*index+2];
         int e2 = iomodel->faces[4*index+3];
 
+        if(e2==-1){
+                this->boundary = true;
+        }
+        else{
+                this->boundary = false;
+        }
+
         /*Set Element hook (4th column may be -1 for boundary edges)*/
         this->helement  = new Hook(&e1,1);
 
@@ -544 +551 @@
                 Ngrad   = fabs(dphids);
                 if(Ngrad<1.e-12) Ngrad = 1.e-12;
 
-                /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
-                IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.);
-                IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(Ngrad,BETA_S-2.);
+                /*Compute the effective conductivity Ks = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
+                IssmDouble Ks = ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.);
 
                 /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
                 qc = - Ks * dphids;
@@ -584 +590 @@
 
         /*Initialize Element matrix and return if necessary*/
         Tria*  tria=(Tria*)element;
-        if(!tria->IsIceInElement()){
+        if(!tria->IsIceInElement() || this->boundary){
                 this->S = 0.;
                 return;
         }
@@ -673 +679 @@
         IssmDouble beta = 1./(rho_ice*L)*( fabs(lc*qc*dphids) + C*fFactor*dPw );
 
         /*Solve ODE*/
-        this->S = ODE1(alpha,beta,this->S,dt,0);
+        this->S = ODE1(alpha,beta,this->S,dt,2);
 
+        /*Make sure Area > 0*/
+        if(this->S<0.) this->S = 0.;
+
         /*Clean up and return*/
         delete gauss;
 }

../trunk-jpl/src/c/classes/Loads/Channel.h

@@ -19 +19 @@
 
         private: 
                 IssmDouble S;
+                bool       boundary;
 
         public: 
                 int sid;