Changeset 22124

Timestamp:

09/27/17 09:15:22 (8 years ago)

Author:

Mathieu Morlighem

Message:

CHG: calving = 0 if grounded above sea level

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Penta.cpp (modified) (4 diffs)
• Tria.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -193 +193 @@
         IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
         IssmDouble  sigma_vm,sigma_max,sigma_max_floating,sigma_max_grounded;
-        IssmDouble  epse_2,groundedice;
+        IssmDouble  epse_2,groundedice,bed;
 
         /* Get node coordinates and dof list: */
@@ -207 +207 @@
         Input* vy_input = inputs->GetInput(VyAverageEnum); _assert_(vy_input);
         Input* gr_input = inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input);
+        Input* bs_input = inputs->GetInput(BaseEnum);                    _assert_(bs_input);
         Input* smax_fl_input = inputs->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
         Input* smax_gr_input = inputs->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
@@ -221 +222 @@
                 vy_input->GetInputValue(&vy,gauss);
                 gr_input->GetInputValue(&groundedice,gauss);
+                bs_input->GetInputValue(&bed,gauss);
                 smax_fl_input->GetInputValue(&sigma_max_floating,gauss);
                 smax_gr_input->GetInputValue(&sigma_max_grounded,gauss);
@@ -248 +250 @@
 
                 /*Assign values*/
-                calvingratex[iv]=vx*sigma_vm/sigma_max;
-                calvingratey[iv]=vy*sigma_vm/sigma_max;
-                calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
+                if(bed>0.){
+                        calvingratex[iv]=0.;
+                        calvingratey[iv]=0.;
+                }
+                else{
+                        calvingratex[iv]=vx*sigma_vm/sigma_max;
+                        calvingratey[iv]=vy*sigma_vm/sigma_max;
+                }
+                calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
         }
 

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -297 +297 @@
 
         IssmDouble  xyz_list[NUMVERTICES][3];
-        GaussTria* gauss=NULL;
         IssmDouble  vx,vy,vel;
         IssmDouble  strainparallel;
-        IssmDouble  propcoeff;
+        IssmDouble  propcoeff,bed;
         IssmDouble  strainperpendicular;
         IssmDouble  calvingratex[NUMVERTICES];
@@ -310 +309 @@
 
         /*Retrieve all inputs and parameters we will need*/
-        Input* vx_input=inputs->GetInput(VxEnum);                                                                                                                                               _assert_(vx_input);
-        Input* vy_input=inputs->GetInput(VyEnum);                                                                                                                                               _assert_(vy_input);
-        Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);                                                           _assert_(strainparallel_input);
-        Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);                                 _assert_(strainperpendicular_input);
-        Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum);                     _assert_(levermanncoeff_input);
+        Input* vx_input=inputs->GetInput(VxEnum);                                                                                                       _assert_(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum);                                                                                                       _assert_(vy_input);
+        Input* bs_input = inputs->GetInput(BaseEnum);                                 _assert_(bs_input);
+        Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);                   _assert_(strainparallel_input);
+        Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);_assert_(strainperpendicular_input);
+        Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum);      _assert_(levermanncoeff_input);
 
         /* Start looping on the number of vertices: */
-        gauss=new GaussTria();
+        GaussTria* gauss=new GaussTria();
         for (int iv=0;iv<NUMVERTICES;iv++){
                 gauss->GaussVertex(iv);
@@ -328 +328 @@
                 strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
                 levermanncoeff_input->GetInputValue(&propcoeff,gauss);
+                bs_input->GetInputValue(&bed,gauss);
 
                 /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
                 calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
-                if(calvingrate[iv]<0){
-                        calvingrate[iv]=0;
+                if(calvingrate[iv]<0. || bed>0.){
+                        calvingrate[iv]=0.;
                 }
                 calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14);