Changeset 27289

Timestamp:

09/24/22 10:01:51 (3 years ago)

Author:

sehrenfe

Message:

partially floating and part no ice elements have zero sheet thickness and zero N

File:

• issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (modified) (8 diffs)

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

@@ -518 +518 @@
         IssmDouble  A,B,n,phi,phi_0;
         IssmDouble  alpha,beta;
+        IssmDouble  oceanLS,iceLS;
 
         /*Fetch number vertices for this element*/
@@ -548 +549 @@
         Input* n_input  = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
         Input* phi_input = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
+        Input* oceanLS_input = element->GetInput(MaskOceanLevelsetEnum); _assert_(oceanLS_input);
+        Input* iceLS_input = element->GetInput(MaskIceLevelsetEnum); _assert_(iceLS_input);
+
 
         /* Start  looping on the number of gaussian points: */
@@ -564 +568 @@
                 b_input->GetInputValue(&b,gauss);
                 H_input->GetInputValue(&H,gauss);
+                oceanLS_input->GetInputValue(&oceanLS,gauss);
+                iceLS_input->GetInputValue(&iceLS,gauss);
+
+                /*Set sheet thickness to zero if floating or no ice*/
+                if(oceanLS<0. || iceLS>0.){
+                        h_new[iv] = 0.;
+                }
+                else{
 
                 /*Get values for a few potentials*/
@@ -591 +603 @@
                 if(h_new[iv]<AEPS) h_new[iv] = AEPS;
         }
-        
+        }
+
         /*Force floating ice to have zero sheet thickness*/
-        if(element->IsAllFloating() || !element->IsIceInElement()){
+        /*if(!element->IsAllGrounded() || !element->IsIceInElement()){
                 for(int iv=0;iv<numvertices;iv++) h_new[iv] = 0.;
-        }
+        }*/
 
         element->AddInput(HydrologySheetThicknessEnum,h_new,P1Enum);
@@ -616 +629 @@
         IssmDouble phi_0, phi_m, p_i;
         IssmDouble H,b,phi;
+        IssmDouble oceanLS,iceLS;
 
         int numnodes = element->GetNumberOfNodes();
@@ -627 +641 @@
         Input* b_input   = element->GetInput(BedEnum); _assert_(b_input);
         Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input);
+        Input* oceanLS_input = element->GetInput(MaskOceanLevelsetEnum); _assert_(oceanLS_input);
+        Input* iceLS_input = element->GetInput(MaskIceLevelsetEnum); _assert_(iceLS_input);
 
         /*Set to 0 if inactive element*/
@@ -645 +661 @@
                 b_input->GetInputValue(&b,gauss);
                 phi_input->GetInputValue(&phi,gauss);
+                oceanLS_input->GetInputValue(&oceanLS,gauss);
+                iceLS_input->GetInputValue(&iceLS,gauss);
 
                 /*Elevation potential*/
@@ -657 +675 @@
                 /*Calculate effective pressure*/
                 N[iv] = phi_0 - phi;
+
+                /*Make sure that all floating ice and ice free areas have zero effective pressure*/
+                if(oceanLS<0.0) N[iv] = 0.0;
+                if(iceLS>0.0) N[iv] = 0.0;
 
                 if(xIsNan<IssmDouble>(N[iv])) _error_("NaN found in solution vector");