Changeset 26780

Timestamp:

01/13/22 13:17:50 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: allow for both MICI as in Crawford et al. and DeConto

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (3 diffs)

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

@@ -29 +29 @@
 #define NUMVERTICES   3
 #define NUMVERTICES1D 2
-//#define ISMICI        1
+//#define MICI          1 //1 = DeConto & Pollard, 2 = DOMINOS
 
 /*Constructors/destructor/copy*/
@@ -4381 +4381 @@
                 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];
+                /*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];
+
+                /*Do we assume that the calving front does not move if MICI is not engaged?*/
+                movingfrontvx[iv] = 0.;
+                movingfrontvy[iv] = 0.;
+                //movingfrontvx[iv] = -2000./(365*24*3600.)*dlsf[0]/norm_dlsf;
+                //movingfrontvy[iv] = -2000./(365*24*3600.)*dlsf[1]/norm_dlsf;
+
+                if(MICI==1 && Hc>80. && bed<0. && fabs(ls)<100.e3){ //Pollard & De Conto
+                        IssmDouble C = (min(max(Hc,80.),100.) - 80.)/20. * 10./(24*3600.); /*Original MICI! convert from m/day to m/s*/
+
+                        /*Front motion = ice speed (v) - calving rate*/
+                        movingfrontvx[iv] = v[0]-C*dlsf[0]/norm_dlsf;
+                        movingfrontvy[iv] = v[1]-C*dlsf[1]/norm_dlsf;
+                }
+                else if (MICI==2 && Hc>135. && bed<0. && fabs(ls)<100.e3){ // Crawford et all
 
                         /*5C Bn (worst case scenario)*/
@@ -4405 +4408 @@
                         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;
+
+                        /*Front motion = ice speed (v) - calving rate*/
+                        movingfrontvx[iv] = v[0] -C*dlsf[0]/norm_dlsf;
+                        movingfrontvy[iv] = v[1] -C*dlsf[1]/norm_dlsf;
                 }
         }