Changeset 26464

Timestamp:

09/30/21 05:38:15 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixing crevasse depth law: missing factor 2

Location:

issm/trunk-jpl

Files:

• src/c/classes/Elements/Tria.cpp (modified) (7 diffs)
• test/Archives/Archive809.arch (modified) ( previous)

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

@@ -29 +29 @@
 #define NUMVERTICES   3
 #define NUMVERTICES1D 2
+//#define ISMICI        1
 
 /*Constructors/destructor/copy*/
@@ -406 +407 @@
         IssmDouble  calvingrate[NUMVERTICES];
         IssmDouble  vx,vy;
-        IssmDouble  water_height, bed,Ho,thickness,surface;
+        IssmDouble  water_height, bed,Hab,thickness,surface;
         IssmDouble  surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal;
         IssmDouble  strainparallel, straineffective,B,n;
         IssmDouble  s_xx,s_xy,s_yy,s1,s2,stmp;
         int         crevasse_opening_stress;
+
+
+        /*reset if no ice in element*/
+        if(!this->IsIceInElement()){
+                for(int i=0;i<NUMVERTICES;i++){
+                        surface_crevasse[i] = 0.;
+                        basal_crevasse[i] = 0.;
+                        crevasse_depth[i] = 0.;
+                }
+                this->AddInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1DGEnum);
+                this->AddInput(BasalCrevasseEnum,&basal_crevasse[0],P1DGEnum);
+                this->AddInput(CrevasseDepthEnum,&crevasse_depth[0],P1DGEnum);
+                return;
+        }
+
 
         /*retrieve the type of crevasse_opening_stress*/
@@ -462 +478 @@
          /*Benn2017,Todd2018: maximum principal stress */
                         Matrix2x2Eigen(&s1,&s2,NULL,NULL,s_xx,s_xy,s_yy);
-                        if(fabs(s2)>fabs(s1)){
-            stmp=s2; s2=s1; s1=stmp;
-         }
+                        s1 = max(0.,max(s1,s2));
                 }
       else{
@@ -471 +485 @@
 
       /*Surface crevasse: sigma'_xx - rho_i g d + rho_fw g d_w = 0*/
-      surface_crevasse[iv] = s1 / (rho_ice*constant_g) + (rho_freshwater/rho_ice)*water_height;
+      surface_crevasse[iv] = 2*s1 / (rho_ice*constant_g) + (rho_freshwater/rho_ice)*water_height;
                 if(surface_crevasse[iv]<0.){
                         surface_crevasse[iv]=0.;
@@ -482 +496 @@
       }
       else{
-         Ho = thickness - (rho_seawater/rho_ice) * (-bed);
-         if(Ho<0.)  Ho=0.;
-         basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
+         Hab = thickness - (rho_seawater/rho_ice) * (-bed);
+         if(Hab<0.)  Hab=0.;
+         basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (2*s1/ (rho_ice*constant_g)-Hab);
          if(basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
       }
@@ -4323 +4337 @@
         }
 
+        #ifdef MICI
         /**************************************  MICI  START ************************************/
         /*MICI from Crawford et al. 2021:
@@ -4331 +4346 @@
          * Tice = 5C;  Bn  ->  I = 1.9e-16; α = 7.3
          */
-        //Input* surface_input = this->GetInput(SurfaceEnum); _assert_(surface_input);
-        //Input* bed_input     = this->GetInput(BedEnum);     _assert_(bed_input);
-        //Input* ls_input      = this->GetInput(MaskIceLevelsetEnum);   _assert_(ls_input);
-        //IssmDouble Hc,bed,ls;
-        //for(int iv=0;iv<NUMVERTICES;iv++){
-        //      gauss.GaussVertex(iv);
-
-        //      surface_input->GetInputValue(&Hc,&gauss);
-        //      bed_input->GetInputValue(&bed,&gauss);
-        //      ls_input->GetInputValue(&ls,&gauss);
-
-        //      if(Hc>135. && bed<0. && fabs(ls)<100.e3){
-
-        //              lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
-        //              norm_dlsf=0.;
-        //              for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
-        //              norm_dlsf=sqrt(norm_dlsf);
-
-        //              /*5C Bn (worst case scenario)*/
-        //              IssmDouble I     = 1.9e-16;
-        //              IssmDouble alpha = 7.3;
-        //              IssmDouble C = I*pow(Hc,alpha)/(24*3600.); /*convert from m/day to m/s*/
-        //              movingfrontvx[iv] = -C*dlsf[0]/norm_dlsf;
-        //              movingfrontvy[iv] = -C*dlsf[1]/norm_dlsf;
-        //      }
-        //}
+        Input* surface_input = this->GetInput(SurfaceEnum); _assert_(surface_input);
+        Input* bed_input     = this->GetInput(BedEnum);     _assert_(bed_input);
+        Input* ls_input      = this->GetInput(MaskIceLevelsetEnum);   _assert_(ls_input);
+        IssmDouble Hc,bed,ls;
+        for(int iv=0;iv<NUMVERTICES;iv++){
+                gauss.GaussVertex(iv);
+
+                surface_input->GetInputValue(&Hc,&gauss);
+                bed_input->GetInputValue(&bed,&gauss);
+                ls_input->GetInputValue(&ls,&gauss);
+
+                /*Do we assume that the calving front does not move?*/
+                //movingfrontvx[iv] = 0.;
+                //movingfrontvy[iv] = 0.;
+
+                //if(Hc>80. && bed<0. && fabs(ls)<100.e3){ //Pollard & De Conto
+                if(Hc>135. && bed<0. && fabs(ls)<100.e3){ // Crawford et all
+
+                        lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
+                        norm_dlsf=0.;
+                        for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
+                        norm_dlsf=sqrt(norm_dlsf);
+
+                        /*use vel direction instead of LSF*/
+                        vx_input->GetInputValue(&v[0],&gauss);
+                        vy_input->GetInputValue(&v[1],&gauss);
+                        vel=sqrt(v[0]*v[0] + v[1]*v[1]);
+                        norm_dlsf = max(vel,1.e-10);
+                        dlsf[0] = v[0];
+                        dlsf[1] = v[1];
+
+                        /*5C Bn (worst case scenario)*/
+                        IssmDouble I     = 1.9e-16;
+                        IssmDouble alpha = 7.3;
+                        IssmDouble C = min(2000.,I*pow(Hc,alpha))/(24*3600.); /*convert from m/day to m/s*/
+                        //IssmDouble C = (min(max(Hc,80.),100.) - 80.)/20. * 10./(24*3600.); /*Original MICI! convert from m/day to m/s*/
+                        movingfrontvx[iv] = -C*dlsf[0]/norm_dlsf;
+                        movingfrontvy[iv] = -C*dlsf[1]/norm_dlsf;
+                }
+        }
+        #endif
         /**************************************  END MICI  *************************************/