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

Index: ../trunk-jpl/src/c/classes/Loads/Channel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Channel.cpp        (revision 24063)
+++ ../trunk-jpl/src/c/classes/Loads/Channel.cpp        (revision 24064)
@@ -54,6 +54,13 @@
        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,9 +551,8 @@
                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,7 +590,7 @@
 
        /*Initialize Element matrix and return if necessary*/
        Tria*  tria=(Tria*)element;
-       if(!tria->IsIceInElement()){
+       if(!tria->IsIceInElement() || this->boundary){
                this->S = 0.;
                return;
        }
@@ -673,8 +679,11 @@
        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;
 }
Index: ../trunk-jpl/src/c/classes/Loads/Channel.h
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Channel.h  (revision 24063)
+++ ../trunk-jpl/src/c/classes/Loads/Channel.h  (revision 24064)
@@ -19,6 +19,7 @@
 
        private: 
                IssmDouble S;
+               bool       boundary;
 
        public: 
                int sid;