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

Index: ../trunk-jpl/src/m/classes/calvingparameterization.m
===================================================================
--- ../trunk-jpl/src/m/classes/calvingparameterization.m        (revision 27251)
+++ ../trunk-jpl/src/m/classes/calvingparameterization.m        (revision 27252)
@@ -11,7 +11,8 @@
                alpha = 0;
                xoffset = 0;
                yoffset = 0;
-               vel_upperbound = 6000;
+               vel_upperbound = 8000;
+               vel_threshold = 6000;
                vel_lowerbound = 0;
        end
        methods
@@ -67,6 +68,7 @@
                        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,6 +80,7 @@
                        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,6 +93,7 @@
                        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
Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp    (revision 27251)
+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp    (revision 27252)
@@ -225,6 +225,7 @@
                        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:
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp        (revision 27251)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp        (revision 27252)
@@ -844,7 +844,7 @@
        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,6 +872,7 @@
        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,7 +891,7 @@
                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,12 +924,12 @@
                                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,12 +946,7 @@
                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);
Index: ../trunk-jpl/src/c/cores/movingfront_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/movingfront_core.cpp       (revision 27251)
+++ ../trunk-jpl/src/c/cores/movingfront_core.cpp       (revision 27252)
@@ -19,6 +19,7 @@
        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,6 +51,10 @@
        /* 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);
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h    (revision 27251)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h    (revision 27252)
@@ -118,6 +118,7 @@
        CalvingXoffsetEnum,
        CalvingYoffsetEnum,
        CalvingVelLowerboundEnum,
+       CalvingVelThresholdEnum,
        CalvingVelUpperboundEnum,
        ConfigurationTypeEnum,
        ConstantsGEnum,