source: issm/oecreview/Archive/27230-27343/ISSM-27251-27252.diff@ 27344

../trunk-jpl/src/m/classes/calvingparameterization.m

@@ -11 +11 @@
                 alpha = 0;
                 xoffset = 0;
                 yoffset = 0;
-                vel_upperbound = 6000;
+                vel_upperbound = 8000;
+                vel_threshold = 6000;
                 vel_lowerbound = 0;
         end
         methods
@@ -67 +68 @@
                         md = checkfield(md,'fieldname','calving.xoffset','NaN',1,'Inf',1,'numel',1);
                         md = checkfield(md,'fieldname','calving.yoffset','NaN',1,'Inf',1,'numel',1);
                         md = checkfield(md,'fieldname','calving.vel_lowerbound','NaN',1,'Inf',1,'numel',1);
+                        md = checkfield(md,'fieldname','calving.vel_threshold','NaN',1,'Inf',1,'numel',1);
                         md = checkfield(md,'fieldname','calving.vel_upperbound','NaN',1,'Inf',1,'numel',1);
                 end % }}}
                 function disp(self) % {{{
@@ -78 +80 @@
                         fielddisplay(self,'xoffset','offset in x-axis');
                         fielddisplay(self,'yoffset','offset in y-axis');
                         fielddisplay(self,'vel_lowerbound','lowerbound of ice velocity to reduce the calving rate [m/a]');
+                        fielddisplay(self,'vel_threshold','threshold of ice velocity to reduce the calving rate [m/a]');
                         fielddisplay(self,'vel_upperbound','upperbound of ice velocity to reduce the calving rate [m/a]');
                 end % }}}
                 function marshall(self,prefix,md,fid) % {{{
@@ -90 +93 @@
                         WriteData(fid,prefix,'object',self,'fieldname','xoffset','format','Double');
                         WriteData(fid,prefix,'object',self,'fieldname','yoffset','format','Double');
                         WriteData(fid,prefix,'object',self,'fieldname','vel_lowerbound','format','Double','scale', 1./yts);
+                        WriteData(fid,prefix,'object',self,'fieldname','vel_threshold','format','Double','scale', 1./yts);
                         WriteData(fid,prefix,'object',self,'fieldname','vel_upperbound','format','Double','scale', 1./yts);
                 end % }}}
         end

../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

@@ -225 +225 @@
                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.xoffset",CalvingXoffsetEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.yoffset",CalvingYoffsetEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.vel_lowerbound",CalvingVelLowerboundEnum));
+                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.vel_threshold",CalvingVelThresholdEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.vel_upperbound",CalvingVelUpperboundEnum));
                         break;
                 default:

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

@@ -844 +844 @@
         IssmDouble  arate, rho_ice, rho_water, thickness;
         int                     use_parameter=-1;
         IssmDouble  gamma, theta, alpha, xoffset, yoffset;
-        IssmDouble  vel_lower, vel_upper, vrate, truncateVrate;
+        IssmDouble  vel_lower, vel_upper, vel_threshold, vrate, truncateVrate;
 
         /* Get node coordinates and dof list: */
         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
@@ -872 +872 @@
         this->FindParam(&yoffset, CalvingYoffsetEnum);
         this->FindParam(&vel_lower, CalvingVelLowerboundEnum);
         this->FindParam(&vel_upper, CalvingVelUpperboundEnum);
+        this->FindParam(&vel_threshold, CalvingVelThresholdEnum);
 
         /* Start looping on the number of vertices: */
         GaussTria* gauss=new GaussTria();
@@ -890 +891 @@
                 sl_input->GetInputValue(&sealevel,gauss);
                 arate_input->GetInputValue(&arate,gauss);
                 vrate = 1.0;
-                if (vel < vel_upper) vrate = vel / vel_upper;
+                if (vel < vel_threshold) vrate = vel / vel_threshold;
 
                 /*Compute strain rate and viscosity: */
                 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
@@ -923 +924 @@
                                 break;
                         case 3:
                                 /* 3 tanh(normal vel): f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
-                                _error_("The normalized velocity is not supported yet!");
-                                gamma = yoffset -  0.5*theta*tanh(alpha*(vel+xoffset));
+                                truncateVrate = (min(vel_upper, max(vel_lower, vel))-vel_lower) / (vel_upper - vel_lower);
+                                gamma = yoffset -  0.5*theta*tanh(alpha*(truncateVrate+xoffset));
                                 break;
                         case 4:
-                                /* 4 tanh(truncated vel): f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
-                                truncateVrate = (min(vel_upper, max(vel_lower, vel))-vel_lower) / vel_upper;
+                                /* 4 tanh(normal vel), cross (0,0): f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
+                                truncateVrate = (min(vel_upper, max(vel_lower, vel))-vel_lower) / (vel_upper - vel_lower);
                                 gamma = 0.5*theta*(tanh(alpha*xoffset) - tanh(alpha*(truncateVrate+xoffset)));
                                 break;
                         case -1:
@@ -945 +946 @@
                 if (bed >= sealevel) gamma = 0.0;
 
                 /*-------------------------------------------*/
-                if (use_parameter < 3) {
-                        calvingrate[iv] = arate*gamma*vrate;
-                } 
-                else {
-                        calvingrate[iv] = arate*gamma;
-                }
+                calvingrate[iv] = arate*gamma*vrate;
         }
         /*Add input*/
         this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);

../trunk-jpl/src/c/cores/movingfront_core.cpp

@@ -19 +19 @@
         int  domaintype, num_extrapol_vars, index,reinit_frequency,step;
         int* extrapol_vars=NULL;
         Analysis  *analysis=NULL;
+        IssmDouble maxVel;
 
         /* recover parameters */
         femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
@@ -50 +51 @@
         /* smoothen slope of lsf for computation of normal on ice domain*/
         levelsetfunctionslope_core(femmodel);
 
+        /* compute the maximal velocity over the whole domain */
+        femmodel->MaxVelx(&maxVel);
+        femmodel->parameters->SetParam(maxVel, CalvingVelThresholdEnum);
+
         /* start the work from here */
         if(VerboseSolution()) _printf0_("   computing calving and undercutting\n");
         Calvingx(femmodel);

../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -118 +118 @@
         CalvingXoffsetEnum,
         CalvingYoffsetEnum,
         CalvingVelLowerboundEnum,
+        CalvingVelThresholdEnum,
         CalvingVelUpperboundEnum,
         ConfigurationTypeEnum,
         ConstantsGEnum,