Changeset 27203

Timestamp:

08/16/22 11:38:38 (3 years ago)

Author:

Cheng Gong

Message:

CHG: add a scaling of calving rate based on ice velocity

Location:

issm/trunk-jpl/src

Files:

• c/analyses/LevelsetAnalysis.cpp (modified) (1 diff)
• c/classes/Elements/Tria.cpp (modified) (5 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• m/classes/calvingparameterization.m (modified) (4 diffs)

issm/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_threshold",CalvingVelThresholdEnum));
                         break;
                 default:

issm/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  gamma, theta, alpha, xoffset, yoffset, vel_threshold, vrate;
 
         /* Get node coordinates and dof list: */
@@ -870 +870 @@
         this->FindParam(&xoffset, CalvingXoffsetEnum);
         this->FindParam(&yoffset, CalvingYoffsetEnum);
+        this->FindParam(&vel_threshold, CalvingVelThresholdEnum);
 
         /* Start looping on the number of vertices: */
@@ -883 +884 @@
                 gr_input->GetInputValue(&groundedice,gauss);
                 bs_input->GetInputValue(&bed,gauss);
+                H_input->GetInputValue(&thickness,gauss);
                 vel=sqrt(vx*vx+vy*vy)+1.e-14;
                 sl_input->GetInputValue(&sealevel,gauss);
                 arate_input->GetInputValue(&arate,gauss);
+                vrate = 1.0;
+                if (vel < vel_threshold) vrate = vel / vel_threshold;
 
                 /*Compute strain rate and viscosity: */
@@ -912 +916 @@
                                 gamma = yoffset -  0.5*theta*tanh(alpha*(bed+xoffset));
                                 break;
+                        case 2:
+                                /* 2 tanh(thicknes): f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
+                                gamma = yoffset -  0.5*theta*tanh(alpha*(-thickness+xoffset));
+                                break;
                         case -1:
                                 /* nothing, just the arate*/
@@ -926 +934 @@
 
                 /*-------------------------------------------*/
-                calvingrate[iv] = arate*gamma;
+                calvingrate[iv] = arate*gamma*vrate;
         }
 

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

@@ -116 +116 @@
         CalvingXoffsetEnum,
         CalvingYoffsetEnum,
+        CalvingVelThresholdEnum,
         ConfigurationTypeEnum,
         ConstantsGEnum,

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

@@ -12 +12 @@
                 xoffset = 0;
                 yoffset = 0;
+                vel_threshold = 6000;
         end
         methods
@@ -51 +52 @@
                         % offset in y-axis 
                         self.yoffset = 0;
+                        % velocity threshold to reduce calving rate
+                        vel_threshold = 6000; % m/a
                 end % }}}
                 function md = checkconsistency(self,md,solution,analyses) % {{{
@@ -57 +60 @@
 
                         md = checkfield(md,'fieldname','calving.min_thickness','>=',0,'NaN',1,'Inf',1,'numel',1);
-                        md = checkfield(md,'fieldname','calving.use_param','values',[-1, 0, 1]);
+                        md = checkfield(md,'fieldname','calving.use_param','values',[-1, 0, 1, 2]);
                         md = checkfield(md,'fieldname','calving.theta','NaN',1,'Inf',1,'numel',1);
                         md = checkfield(md,'fieldname','calving.alpha','NaN',1,'Inf',1,'numel',1);
                         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_threshold','NaN',1,'Inf',1,'numel',1);
                 end % }}}
                 function disp(self) % {{{
                         disp(sprintf('   Calving test parameters:'));
                         fielddisplay(self,'min_thickness','minimum thickness below which no ice is allowed [m]');
-                        fielddisplay(self,'use_param','-1 - just use frontal ablation rate, 0 - f(x) = y_{o} + \alpha (x+x_{o}), 1 - f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o}))');
+                        fielddisplay(self,'use_param','-1 - just use frontal ablation rate, 0 - f(x) = y_{o} + \alpha (x+x_{o}), 1 - f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})), 2 - tanh(thickness)');
                         fielddisplay(self,'theta','the amplifier');
                         fielddisplay(self,'alpha','the slope');
                         fielddisplay(self,'xoffset','offset in x-axis');
                         fielddisplay(self,'yoffset','offset in y-axis');
+                        fielddisplay(self,'vel_threshold','threshold of ice velocity to reduce the calving rate');
                 end % }}}
                 function marshall(self,prefix,md,fid) % {{{
@@ -81 +86 @@
                         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_threshold','format','Double','scale', 1./yts);
                 end % }}}
         end