Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 27846) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 27847) @@ -253,4 +253,5 @@ ./modules/StochasticForcingx/StochasticForcingx.cpp \ ./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp \ + ./modules/ComputeMeanElax/ComputeMeanElax.cpp \ ./cores/ProcessArguments.cpp \ ./cores/ResetBoundaryConditions.cpp \ Index: /issm/trunk-jpl/src/c/analyses/DebrisAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DebrisAnalysis.cpp (revision 27846) +++ /issm/trunk-jpl/src/c/analyses/DebrisAnalysis.cpp (revision 27847) @@ -10,5 +10,7 @@ #define FINITEELEMENT P1Enum -#define EPS 1e-14 +#define EPS 1e-8 +#define RHO 2 +//#define KAPPA 5 /*Model processing*/ @@ -115,4 +117,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.debris.min_thickness",DebrisMinThicknessEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.debris.packingfraction",DebrisPackingFractionEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.debris.max_displacementvelocity",DebrisMaxdisplacementvelocityEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.debris.removal_slope_threshold",DebrisRemovalSlopeThresholdEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.debris.removal_stress_threshold",DebrisRemovalStressThresholdEnum)); @@ -128,11 +131,10 @@ void DebrisAnalysis::Core(FemModel* femmodel){/*{{{*/ - //PreProcessing(femmodel); - //femmodel->parameters->SetParam(VxDebrisEnum,InputToExtrudeEnum); - //extrudefromtop_core(femmodel); + PreProcessing(femmodel); femmodel->SetCurrentConfiguration(DebrisAnalysisEnum); solutionsequence_linear(femmodel); + femmodel->inputs->DuplicateInput(DebrisThicknessEnum,DebrisThicknessOldEnum); PostProcessing(femmodel); - + }/*}}}*/ void DebrisAnalysis::PreCore(FemModel* femmodel){/*{{{*/ @@ -197,5 +199,11 @@ Input* vy_input=NULL; if(dim>1){vy_input = topelement->GetInput(VyDebrisEnum); _assert_(vy_input);} + Input* slopex_input = topelement->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input); + Input* slopey_input = topelement->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input); + Input* thickness_input = topelement->GetInput(DebrisThicknessEnum); _assert_(thickness_input); h=topelement->CharacteristicLength(); + //IssmDouble hx,hy,hz; + //topelement->ElementSizes(&hx,&hy,&hz); + //h=sqrt(hx*hx+hy*hy); /* Start looping on the number of gaussian points: */ @@ -209,5 +217,5 @@ vx_input->GetInputValue(&vx,gauss); if(dim==2) vy_input->GetInputValue(&vy,gauss); - + D_scalar=gauss->weight*Jdet; for(int i=0;ivalues[i*numnodes+j] += D_scalar*basis[i]*basis[j]; @@ -215,8 +223,8 @@ /*Advection terms*/ vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + dvxdx=dvx[0]; D_scalar=dt*gauss->weight*Jdet; if(dim==2){ vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); - dvxdx=dvx[0]; dvydy=dvy[1]; for(int i=0;iFindParam(&step,StepEnum); + IssmDouble slope_threshold; + topelement->FindParam(&slope_threshold,DebrisRemovalSlopeThresholdEnum); + IssmDouble rad2deg=180./M_PI; + IssmDouble kappa,f,smb,debristhickness;//,slopex; + IssmDouble M=1,C=-0.2,m,alpha_t_deg=25.0; + IssmDouble alpha_t_rad=tan(alpha_t_deg/rad2deg); + IssmDouble Diff,fraction; + Diff=3.2/3e7; + //Input* slopex_input=topelement->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input); + Input* smb_input=topelement->GetInput(SmbMassBalanceEnum); _assert_(smb_input); + Input* debristhickness_input=topelement->GetInput(DebrisThicknessEnum); _assert_(debristhickness_input); + + /*if(isdisplacement){ + + slopex_input->GetInputValue(&slopex, gauss); + smb_input->GetInputValue(&smb, gauss); + debristhickness_input->GetInputValue(&debristhickness, gauss); + if((atan(fabs(slopex))*rad2deg)=alpha_t_deg & fabs(slopex)<=((1-C)/M)){ + m=1/(((1-C)/M)-alpha_t_rad); + f=m*fabs(slopex)-m*alpha_t_rad; + }else{ + f=1; + } + //f=1; + //kappa=-5.6e16*smb*debristhickness*f; + //kappa=debristhickness/h*4e9*f; + //kappa=14.2809e8*f; // 25° + kappa=247*f; // 35° + if(dim==2){ + for(int i=0;ivalues[i*numnodes+j] += D_scalar*kappa*( + dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i] + ); + } + } + }else{ + for(int i=0;ivalues[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i]); + } + } + } + } + /*******************************************************************/ + + IssmDouble rho,KAPPA,slopex,slopey, thickness; if(FINITEELEMENT==P1Enum){ - rho=2.; + rho=RHO; }else if(FINITEELEMENT==P2Enum){ - rho=4.; + rho=RHO*2.; } @@ -249,5 +309,6 @@ /*SSA*/ if(dim==1){ - vx_input->GetInputValue(&vx,gauss); + //vx_input->GetInputValue(&vx,gauss); + vx_input->GetInputAverage(&vx); D[0]=h/rho*fabs(vx); }else{ @@ -262,9 +323,12 @@ if(dim==1){ vx_input->GetInputValue(&vx,gauss); + //vx_input->GetInputAverage(&vx); vel=fabs(vx)+EPS; D[0] = h/(rho*vel)*vx*vx; }else{ vx_input->GetInputValue(&vx,gauss); + //vx_input->GetInputAverage(&vx); vy_input->GetInputValue(&vy,gauss); + //vy_input->GetInputAverage(&vy); vel=sqrt(vx*vx+vy*vy)+EPS; D[0*dim+0]=h/(rho*vel)*vx*vx; @@ -277,9 +341,28 @@ if(dim==1){ vx_input->GetInputValue(&vx,gauss); - tau=h/(rho*max(fabs(vx),EPS)); + //vx_input->GetInputAverage(&vx); + //tau=h/(rho*max(fabs(vx),EPS)); + tau=h/(rho*(fabs(vx)+EPS)); }else{ vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); - tau=h/(rho*sqrt(vx*vx+vy*vy)+EPS); + vx_input->GetInputAverage(&vx); + vy_input->GetInputAverage(&vy); + vel=sqrt(vx*vx+vy*vy)+EPS; + + KAPPA=vel; + // adapt diffusion coefficient from here: http://websrv.cs.umt.edu/isis/index.php/Solving_the_equation_for_thickness_evolution + slopex_input->GetInputValue(&slopex,gauss); + slopey_input->GetInputValue(&slopey,gauss); + thickness_input->GetInputValue(&thickness,gauss); + //KAPPA=vel*thickness/(pow(pow(slopex,2)+pow(slopey,2),0.5)+EPS); + tau=h/(rho*vel)*(cosh(vel*h/2/KAPPA)/sinh(vel*h/2/KAPPA)-1./(vel*h/2/KAPPA)); + + /*kappa=KAPPA/yts; + if(vel*h/(3*2*kappa)<1){ + tau=pow(h,2)/(3*2*2*kappa); + } + else tau=h/(2*vel);*/ + //tau=dt*3; } } @@ -327,5 +410,4 @@ } - /*Advection matrix - part 2, B*/ for(int i=0;ivalues[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx)*(basis[i]*dvxdx);; + Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx)*(basis[i]*dvxdx); } } @@ -442,4 +524,6 @@ Input* smb_input = topelement->GetInput(SmbMassBalanceEnum); _assert_(smb_input); Input* thickness_input = topelement->GetInput(DebrisThicknessEnum); _assert_(thickness_input); + Input* slopex_input = topelement->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input); + Input* slopey_input = topelement->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input); Input* vx_input = topelement->GetInput(VxDebrisEnum); _assert_(vx_input); Input* vy_input=NULL; @@ -448,10 +532,13 @@ } h=topelement->CharacteristicLength(); - - IssmDouble rho; + //IssmDouble hx,hy,hz; + //topelement->ElementSizes(&hx,&hy,&hz); + //h=sqrt(hx*hx+hy*hy); + + IssmDouble rho,KAPPA,slopex,slopey; if(FINITEELEMENT==P1Enum){ - rho=2.; + rho=RHO; }else if(FINITEELEMENT==P2Enum){ - rho=4.; + rho=RHO*2.; } @@ -466,9 +553,9 @@ thickness_input->GetInputValue(&thickness,gauss); if(smb>0.){ - psi=thickness-0.*dt*smb*pf; + psi=thickness-0*dt*smb*pf; }else{ psi=thickness-dt*smb*pf; } - + //psi=thickness-dt*smb*pf; for(int i=0;ivalues[i]+=Jdet*gauss->weight*psi*basis[i]; @@ -480,5 +567,7 @@ if(dim==1){ vx_input->GetInputValue(&vx,gauss); - tau=h/(rho*max(fabs(vx),EPS)); + //vx_input->GetInputAverage(&vx); + //tau=h/(rho*max(fabs(vx),EPS)); + tau=h/(rho*(fabs(vx)+EPS)); /*Force vector - part 2*/ @@ -494,8 +583,24 @@ vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); + vx_input->GetInputAverage(&vx); + vy_input->GetInputAverage(&vy); vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); - vel=sqrt(vx*vx+vy*vy); - dvydy=dvy[1]; - tau=h/(rho*vel+EPS); + dvydy=dvy[1]; + + vel=sqrt(vx*vx+vy*vy)+EPS; + // tuned value + KAPPA=vel; + // adapt diffusion coefficient from here: http://websrv.cs.umt.edu/isis/index.php/Solving_the_equation_for_thickness_evolution + slopex_input->GetInputValue(&slopex,gauss); + slopey_input->GetInputValue(&slopey,gauss); + //KAPPA=vel*thickness/(pow(pow(slopex,2)+pow(slopey,2),0.5)+EPS); + tau=h/(rho*vel)*(cosh(vel*h/2/KAPPA)/sinh(vel*h/2/KAPPA)-1./(vel*h/2/KAPPA)); + + /*kappa=KAPPA/yts; + if(vel*h/(3*2*kappa)<1){ + tau=pow(h,2)/(3*2*2*kappa); + } + else tau=h/(2*vel);*/ + //tau=dt*3; /*Force vector - part 2*/ @@ -531,8 +636,11 @@ int numnodes = element->GetNumberOfNodes(); IssmDouble* newthickness = xNew(numnodes); + IssmDouble* icethickness = xNew(numnodes); + /*Use the dof list to index into the solution vector: */ IssmDouble minthickness = element->FindParam(DebrisMinThicknessEnum); element->GetDofListLocal(&ddoflist,NoneApproximationEnum,GsetEnum); + element->GetInputListOnVertices(&icethickness[0],ThicknessEnum); for(int i=0;i(newthickness); + xDelete(icethickness); xDelete(ddoflist); //*/ @@ -599,5 +708,4 @@ femmodel->parameters->FindParam(&removalmodel,DebrisRemovalmodelEnum); Element* element= NULL; - Element* topelement = NULL; if(removalmodel==0){ @@ -616,9 +724,9 @@ IssmDouble* slopey = xNew(numnodes); IssmDouble* onsurface = xNew(numnodes); - IssmDouble* ls_active = xNew(numnodes); - element->GetInputListOnNodes(debristhickness,DebrisThicknessEnum); + //IssmDouble* ls_active = xNew(numnodes); + element->GetInputListOnNodes(debristhickness,DebrisThicknessOldEnum); element->GetInputListOnNodes(icethickness,ThicknessEnum); element->GetInputListOnNodes(onsurface,MeshVertexonsurfaceEnum); - element->GetInputListOnNodes(ls_active,DebrisMaskNodeActivationEnum); + //element->GetInputListOnNodes(ls_active,DebrisMaskNodeActivationEnum); dim=1; @@ -643,6 +751,7 @@ if(icethickness[k]<=(iceminthickness+0.00001)) kk++; } - isminthicknessinelement=true; - if(kkAddInput(DebrisThicknessEnum,debristhickness,P1Enum); + element->AddInput(DebrisThicknessEnum,debristhickness,P1DGEnum); + + xDelete(debristhickness); + xDelete(onsurface); + xDelete(icethickness); + xDelete(slopex); + xDelete(slopey); + //xDelete(ls_active); + break; + } + case 2:{ + IssmDouble stress_threshold=element->FindParam(DebrisRemovalStressThresholdEnum); + IssmDouble gravity=element->FindParam(ConstantsGEnum); + IssmDouble stress,rhod=1900.; + int kk=0; + for(k=0; kstress_threshold) debristhickness[k]=0; + if(stress>stress_threshold) remove_debris=true; + } + if(remove_debris){ + for(k=0; kAddInput(DebrisThicknessEnum,debristhickness,P1DGEnum); xDelete(debristhickness); @@ -667,41 +811,4 @@ break; } - case 2:{ - IssmDouble stress_threshold = element->FindParam(DebrisRemovalStressThresholdEnum); - IssmDouble gravity = element->FindParam(ConstantsGEnum); - IssmDouble stress,rhod=1900.; - int kk=0; - for(k=0; kstress_threshold) remove_debris=true; - //if(stress>stress_threshold) debristhickness[k]=0.; - } - //if((stress_sum/double(kk))>stress_threshold) remove_debris=true; - if(remove_debris){ - for(k=0; kAddInput(DebrisThicknessEnum,debristhickness,P1Enum); - - xDelete(debristhickness); - xDelete(icethickness); - xDelete(slopex); - xDelete(slopey); - xDelete(ls_active); - break; - } default: _error_("removalmodel "<parameters->FindParam(&removalmodel,DebrisRemovalmodelEnum); + femmodel->parameters->FindParam(&displacementmodel,DebrisDisplacementmodelEnum); Element* element= NULL; @@ -718,66 +829,121 @@ element=xDynamicCast(object); - int numvertices = element->GetNumberOfVertices(); + int numvertices=element->GetNumberOfVertices(); + int numnodes=element->GetNumberOfNodes(); + numvertices=numnodes; IssmDouble* vx = xNew(numvertices); + IssmDouble* vy = xNew(numvertices); + IssmDouble* vx_new = xNew(numvertices); + IssmDouble* vy_new = xNew(numvertices); IssmDouble* debristhickness = xNew(numvertices); IssmDouble* slopex = xNew(numvertices); - IssmDouble* onsurface = xNew(numvertices); - IssmDouble* icethickness = xNew(numvertices); - + IssmDouble* slopey = xNew(numvertices); + IssmDouble* icethickness = xNew(numvertices); + int dim=1,domaintype; + femmodel->parameters->FindParam(&domaintype,DomainTypeEnum); + if(domaintype!=Domain2DverticalEnum){ + dim=2; + } element->GetInputListOnVertices(&debristhickness[0],DebrisThicknessEnum); - element->GetInputListOnVertices(&vx[0],VxDebrisEnum); + element->GetInputListOnVertices(&vx[0],VxEnum); element->GetInputListOnVertices(&slopex[0],SurfaceSlopeXEnum); - element->GetInputListOnVertices(&onsurface[0],MeshVertexonsurfaceEnum); element->GetInputListOnVertices(&icethickness[0],ThicknessEnum); - - IssmDouble slope,rad2deg=180./M_PI; //=57.2958 - IssmDouble vslipx,rhod=1900.; + if(dim==2){ + element->GetInputListOnVertices(&vy[0],VyEnum); + element->GetInputListOnVertices(&slopey[0],SurfaceSlopeYEnum); + } + + IssmDouble rad2deg=180./M_PI; //=57.2958 + IssmDouble vslipx,vslipy,stress,rhod=1900.; + IssmDouble taud_plus, taud_minus, taudx, taudy, taud; IssmDouble gravity=element->FindParam(ConstantsGEnum); IssmDouble slope_threshold=element->FindParam(DebrisRemovalSlopeThresholdEnum); + IssmDouble stress_threshold = element->FindParam(DebrisRemovalStressThresholdEnum); + IssmDouble maxv = element->FindParam(DebrisMaxdisplacementvelocityEnum); IssmDouble iceminthickness=element->FindParam(MasstransportMinThicknessEnum); + IssmDouble hx,hy,hz; + element->ElementSizes(&hx,&hy,&hz); int step; - IssmDouble dt, maxv; + IssmDouble dt; IssmDouble yts=31536000.; femmodel->parameters->FindParam(&step,StepEnum); femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - bool isminthicknessinelement=false; - for(int i=0;i.5){ - slope=fabs(slopex[i]); - //if((atan(slope)*rad2deg)>25.){ - //if(debristhickness[i]>0.01){ - vslipx=1.0/yts; - //maxv=10.0/2./dt; - //vslipx=-slope_threshold*rhod*gravity*debristhickness[i]*slopex[i]/yts; - vx[i]=vx[i]+vslipx; - //debristhickness[i]=debristhickness[i]; - //if(vx[i]>maxv) vx[i]=maxv; - //} - //} - //} - } - } - //if(step%100==0) - element->AddInput(VxDebrisEnum,vx,P1Enum); - + //isminthicknessinelement=true; + //if(icethickness[i]<=(iceminthickness+0.01)){ + //if(kk=taud_plus){ + vx[i]=vx[i]+vslipx; + if(dim==2) vy[i]=vy[i]+vslipy; + }else if(taud>=taud_minus & taudmaxv/yts) vx[i]=maxv/yts*fabs(vx[i])/(vx[i]+1e-10); + if(fabs(vy[i])>maxv/yts) vy[i]=maxv/yts*fabs(vy[i])/(vy[i]+1e-10); + } + } + element->AddInput(VxDebrisEnum,vx,P1Enum); + if(dim==2) element->AddInput(VyDebrisEnum,vy,P1Enum); + } /* Free resources */ xDelete(debristhickness); xDelete(icethickness); xDelete(vx); + xDelete(vx_new); xDelete(slopex); - xDelete(onsurface); - } - -}/*}}}*/ + xDelete(vy); + xDelete(vy_new); + xDelete(slopey); + } +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 27846) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 27847) @@ -242,7 +242,21 @@ } break; - case SMBdebrisMLEnum: - iomodel->FetchDataToInput(inputs,elements,"md.initialization.debris",DebrisThicknessEnum); - break; + case SMBdebrisEvattEnum: + iomodel->FetchDataToInput(inputs,elements,"md.initialization.debris",DebrisThicknessEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.snowheight",SmbSnowheightEnum); + iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum); + iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlyprecipitation",SmbPrecipitationEnum); + iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlydsradiation",SmbMonthlydsradiationEnum); + iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlydlradiation",SmbMonthlydlradiationEnum); + iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlywindspeed",SmbMonthlywindspeedEnum); + iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlyairhumidity",SmbMonthlyairhumidityEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.precipitation_anomaly",SmbPrecipitationsAnomalyEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.temperature_anomaly",SmbTemperaturesAnomalyEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.dsradiation_anomaly",SmbDsradiationAnomalyEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.dlradiation_anomaly",SmbDlradiationAnomalyEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.windspeed_anomaly",SmbWindspeedAnomalyEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.airhumidity_anomaly",SmbAirhumidityAnomalyEnum); + break; default: _error_("Surface mass balance model "< /*}}}*/ #define MAXVERTICES 6 /*Maximum number of vertices per element, currently Penta, to avoid dynamic mem allocation*/ @@ -2350,7 +2351,7 @@ IssmDouble minthickness = this->FindParam(MasstransportMinThicknessEnum); - + Input* input=this->GetInput(ThicknessEnum); _assert_(input); - if(input->GetInputMax()<=(minthickness+0.00000001)){ + if(input->GetInputMin()<=(minthickness+0.00000001)){ return true; } @@ -3400,5 +3401,5 @@ monthlytemperatures[iv*12+month]=monthlytemperatures[iv*12+month]-273.15; // conversion from Kelvin to celcius for PDD module dinput2->GetInputValue(&monthlyprec[iv*12+month],gauss,month); - monthlyprec[iv*12+month]=monthlyprec[iv*12+month]*yts; + monthlyprec[iv*12+month]=monthlyprec[iv*12+month];//*yts; } } @@ -3594,5 +3595,4 @@ Gauss* gauss=this->NewGauss(); for(int month=0;month<12;month++){ - for(int iv=0;ivGaussVertex(iv); @@ -3727,4 +3727,334 @@ xDelete(smbcorr); xDelete(melt_star); +} +/*}}}*/ +void Element::SmbDebrisEvatt(){/*{{{*/ + + const int NUM_VERTICES = this->GetNumberOfVertices(); + const int NUM_VERTICES_MONTHS_PER_YEAR = NUM_VERTICES * 365; + + int i,vertexlids[MAXVERTICES];; + IssmDouble* smb=xNew(NUM_VERTICES); // surface mass balance + IssmDouble* melt=xNew(NUM_VERTICES); // melting comp. of surface mass balance + IssmDouble* summermelt=xNew(NUM_VERTICES); + IssmDouble* albedo=xNew(NUM_VERTICES); // melting comp. of surface mass balance + IssmDouble* summeralbedo=xNew(NUM_VERTICES); + IssmDouble* accu=xNew(NUM_VERTICES); // accuumulation comp. of surface mass balance + IssmDouble* monthlytemperatures=xNew(NUM_VERTICES_MONTHS_PER_YEAR); + IssmDouble* monthlyprecip=xNew(NUM_VERTICES_MONTHS_PER_YEAR); + IssmDouble* monthlylw=xNew(NUM_VERTICES_MONTHS_PER_YEAR); + IssmDouble* monthlysw=xNew(NUM_VERTICES_MONTHS_PER_YEAR); + IssmDouble* monthlywind=xNew(NUM_VERTICES_MONTHS_PER_YEAR); + IssmDouble* monthlyhumidity=xNew(NUM_VERTICES_MONTHS_PER_YEAR); + IssmDouble* yearlytemperatures=xNew(NUM_VERTICES); memset(yearlytemperatures, 0., NUM_VERTICES*sizeof(IssmDouble)); + IssmDouble* p_ampl=xNew(NUM_VERTICES); // precip anomaly + IssmDouble* t_ampl=xNew(NUM_VERTICES); // temperature anomaly + IssmDouble* lw_ampl=xNew(NUM_VERTICES); + IssmDouble* sw_ampl=xNew(NUM_VERTICES); + IssmDouble* wind_ampl=xNew(NUM_VERTICES); + IssmDouble* humidity_ampl=xNew(NUM_VERTICES); + + IssmDouble* surface=xNew(NUM_VERTICES); // actual surface height + IssmDouble* s0t=xNew(NUM_VERTICES); // reference elevation + IssmDouble* snowheight=xNew(NUM_VERTICES); + IssmDouble* debriscover=xNew(NUM_VERTICES); + IssmDouble rho_water,rho_ice,Tf,debris,debris_here; + IssmDouble qlaps,rlaps,dslaps,dllaps,windspeedlaps,humiditylaps,Tm; + IssmDouble inv_twelve=1./365.; //factor for monthly average + IssmDouble time,yts,time_yr,lambda; + IssmDouble MonthlyMelt,CleanIceMonthlyMelt, CumMonthlyMelt=0,CleanIceMelt,CumMonthlySummerMelt=0; + IssmDouble MeanAlbedo=0, MeanSummerAlbedo=0; + bool isdebris,ismonthly,isAnderson,iscryokarst; + this->parameters->FindParam(&isdebris,TransientIsdebrisEnum); + this->parameters->FindParam(&ismonthly,SmbDebrisIsMonthlyEnum); + this->parameters->FindParam(&isAnderson,SmbDebrisIsAndersonEnum); + this->parameters->FindParam(&iscryokarst,SmbDebrisIsCryokarstEnum); + IssmDouble PhiD=0.,p; + IssmDouble icealbedo=this->FindParam(SmbIcealbedoEnum); + IssmDouble snowalbedo=this->FindParam(SmbSnowalbedoEnum); + IssmDouble debrisalbedo=this->FindParam(SmbDebrisalbedoEnum); + IssmDouble Lm=this->FindParam(MaterialsLatentheatEnum); + IssmDouble D0=this->FindParam(SmbDebrisAndersonD0Enum); + int step; + this->FindParam(&step,StepEnum); + + // cryokarst + int dim=1,domaintype; + this->parameters->FindParam(&domaintype,DomainTypeEnum); + if(domaintype!=Domain2DverticalEnum) dim=2; + IssmDouble taud_plus=110e3, taud_minus=60e3; + IssmDouble taud, slope, gravity, taudx, taudy; + this->parameters->FindParam(&gravity,ConstantsGEnum); + IssmDouble* slopex = xNew(NUM_VERTICES); + IssmDouble* slopey = xNew(NUM_VERTICES); + IssmDouble* icethickness = xNew(NUM_VERTICES); + + /*Get material parameters :*/ + rho_water=this->FindParam(MaterialsRhoSeawaterEnum); + rho_ice=this->FindParam(MaterialsRhoIceEnum); + IssmDouble sconv=(rho_water/rho_ice); + Tf=this->FindParam(MaterialsMeltingpointEnum); + + /*Get parameters for height corrections*/ + qlaps=this->FindParam(SmbDesfacEnum); // comment MR; on alpine galciers we dont have the dertification effect + rlaps=this->FindParam(SmbRlapsEnum); + dslaps=this->FindParam(SmbSWlapsEnum); + dllaps=this->FindParam(SmbLWlapsEnum); + windspeedlaps=this->FindParam(SmbWindspeedlapsEnum); + humiditylaps=this->FindParam(SmbHumiditylapsEnum); + + /* Get time */ + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + time_yr=floor(time/yts)*yts; + + /*Get inputs*/ + DatasetInput* tempmonthly =this->GetDatasetInput(SmbMonthlytemperaturesEnum); _assert_(tempmonthly); + DatasetInput* precipmonthly =this->GetDatasetInput(SmbPrecipitationEnum); _assert_(precipmonthly); + DatasetInput* lwmonthly =this->GetDatasetInput(SmbMonthlydlradiationEnum); _assert_(lwmonthly); + DatasetInput* swmonthly =this->GetDatasetInput(SmbMonthlydsradiationEnum); _assert_(swmonthly); + DatasetInput* windmonthly =this->GetDatasetInput(SmbMonthlywindspeedEnum); _assert_(windmonthly); + DatasetInput* humiditymonthly =this->GetDatasetInput(SmbMonthlyairhumidityEnum); _assert_(humiditymonthly); + + /*loop over vertices: */ + Gauss* gauss=this->NewGauss(); + for(int month=0;month<365;month++){ + for(int iv=0;ivGaussVertex(iv); + tempmonthly->GetInputValue(&monthlytemperatures[iv*365+month],gauss,month); + monthlytemperatures[iv*365+month]=monthlytemperatures[iv*365+month]-Tf; // conversion from Kelvin to celcius for PDD module + precipmonthly->GetInputValue(&monthlyprecip[iv*365+month],gauss,month); + monthlyprecip[iv*365+month]=monthlyprecip[iv*365+month]*yts; // from m/s to m/a + lwmonthly->GetInputValue(&monthlylw[iv*365+month],gauss,month); + swmonthly->GetInputValue(&monthlysw[iv*365+month],gauss,month); + windmonthly->GetInputValue(&monthlywind[iv*365+month],gauss,month); + humiditymonthly->GetInputValue(&monthlyhumidity[iv*365+month],gauss,month); + } + } + + /*Recover info at the vertices: */ + GetInputListOnVertices(&surface[0],SurfaceEnum); + GetInputListOnVertices(&s0t[0],SmbS0tEnum); + GetInputListOnVertices(&snowheight[0],SmbSnowheightEnum); + GetInputListOnVertices(&debriscover[0],DebrisThicknessEnum); + GetInputListOnVertices(&t_ampl[0],SmbTemperaturesAnomalyEnum); + GetInputListOnVertices(&p_ampl[0],SmbPrecipitationsAnomalyEnum); + GetInputListOnVertices(&lw_ampl[0],SmbDsradiationAnomalyEnum); + GetInputListOnVertices(&sw_ampl[0],SmbDlradiationAnomalyEnum); + GetInputListOnVertices(&wind_ampl[0],SmbWindspeedAnomalyEnum); + GetInputListOnVertices(&humidity_ampl[0],SmbAirhumidityAnomalyEnum); + if(iscryokarst){ + GetInputListOnVertices(&slopex[0],SurfaceSlopeXEnum); + GetInputListOnVertices(&icethickness[0],ThicknessEnum); + if(dim==2){ + GetInputListOnVertices(&slopey[0],SurfaceSlopeYEnum); + } + taudx=rho_ice*gravity*icethickness[i]*slopex[i]; + if(dim==2) taudy=rho_ice*gravity*icethickness[i]*slopey[i]; + taud=sqrt(taudx*taudx+taudy*taudy); + } + IssmDouble Alphaeff,Alphaeff_cleanice; + + /*measure the surface mass balance*/ + for (int iv = 0; iv i.e + /*Partition of precip in solid and liquid parts */ + IssmDouble temp_plus=1; + IssmDouble temp_minus=-1.; + IssmDouble frac_solid; + if(Tm>=temp_plus){ + frac_solid=0; + }else if(Tm<=temp_minus){ + frac_solid=1; + }else{ + frac_solid=1*(1-cos(PI*(temp_plus-Tm)/(temp_plus-temp_minus)))/2; + } + + /*Get yearly temperatures and accumulation */ + yearlytemperatures[iv]=yearlytemperatures[iv]+((monthlytemperatures[iv*365+imonth]-rlaps*st+Tf+t_ampl[iv]))/365; // Has to be in Kelvin + accu[iv]=accu[iv]+P*frac_solid; + if(yearlytemperatures[iv]>Tf) yearlytemperatures[iv]=Tf; + + CleanIceMonthlyMelt=((In-(Eps*Sigma*(Tf*Tf*Tf*Tf))+ + Q*(1.-Alphaeff)+ + (Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))*Tm)/((1-PhiD)*rho_ice*Lm)/(1.+ + ((Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/ + K*debris)-(Lv*Ustar*Ustar*((Humidity))*(exp(-Gamma*Xr)))/((1.-PhiD)* + rho_ice*Lm*Ustar)/(((Um + -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris))); + if(CleanIceMonthlyMelt<0) CleanIceMonthlyMelt=0.; + MonthlyMelt=CleanIceMonthlyMelt; + + if(ismb==1){ + snowheight[iv]=snowheight[iv]+(P-CleanIceMonthlyMelt*yts/365); + if(snowheight[iv]<=0) snowheight[iv]=0.; + if(snowheight[iv]<=0.0001){ + p=debris_here*PhiD/(2*0.2*0.01); //Eq. 51 from Evatt et al 2015 without source term g*t + if(p>1.) p=1.; + if(p>=0.999){ + Alphaeff=debrisalbedo; + } else { + Alphaeff=Alphaeff_cleanice+p*(debrisalbedo-Alphaeff_cleanice); + } + debris=debris_here; + MonthlyMelt=((In-(Eps*Sigma*(Tf*Tf*Tf*Tf))+ + Q*(1.-Alphaeff)+ + (Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))*Tm)/((1-PhiD)*rho_ice*Lm)/(1.+ + ((Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/ + K*debris)-(Lv*Ustar*Ustar*((Humidity))*(exp(-Gamma*Xr)))/((1.-PhiD)* + rho_ice*Lm*Ustar)/(((Um-2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris))); + if(MonthlyMelt<0) MonthlyMelt=0.; + MeanSummerAlbedo=MeanSummerAlbedo+Alphaeff; + CumMonthlySummerMelt=CumMonthlySummerMelt+MonthlyMelt/365; + } + } + CumMonthlyMelt=CumMonthlyMelt+MonthlyMelt/365; + } + MeanAlbedo=MeanAlbedo+Alphaeff; + if(ismb==0) CleanIceMelt=CumMonthlyMelt; + } + + if(iscryokarst){ + if(taud>=taud_plus){ + lambda=0; + }else if(taud>=taud_minus & taudAddInput(SmbMassBalanceEnum,smb,P1Enum); + this->AddInput(SmbAccumulationEnum,accu,P1Enum); + this->AddInput(SmbMeltEnum,melt,P1Enum); + this->AddInput(SmbSummerMeltEnum,summermelt,P1Enum); + this->AddInput(SmbSnowheightEnum,snowheight,P1Enum); + this->AddInput(SmbAlbedoEnum,albedo,P1Enum); + this->AddInput(SmbSummerAlbedoEnum,summeralbedo,P1Enum); + this->AddInput(TemperaturePDDEnum,yearlytemperatures,P1Enum); // TemperaturePDD is wrong here, but don't want to create new Enum ... + + /*clean-up*/ + xDelete(monthlytemperatures); + xDelete(monthlyprecip); + xDelete(monthlylw); + xDelete(monthlysw); + xDelete(monthlywind); + xDelete(monthlyhumidity); + xDelete(smb); + xDelete(surface); + xDelete(melt); + xDelete(summermelt); + xDelete(albedo); + xDelete(summeralbedo); + xDelete(accu); + xDelete(yearlytemperatures); + xDelete(s0t); + xDelete(snowheight); + xDelete(debriscover); + xDelete(t_ampl); + xDelete(p_ampl); + xDelete(lw_ampl); + xDelete(sw_ampl); + xDelete(humidity_ampl); + xDelete(wind_ampl); + xDelete(slopex); + xDelete(slopey); + xDelete(icethickness); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27846) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27847) @@ -177,4 +177,5 @@ void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac); void PositiveDegreeDaySicopolis(bool isfirnwarming); + void SmbDebrisEvatt(); void RignotMeltParameterization(); void ResultInterpolation(int* pinterpolation,int*nodesperelement,int* parray_size, int output_enum); @@ -256,4 +257,5 @@ virtual void ComputeStressTensor(void)=0; virtual void ComputeEsaStrainAndVorticity(void)=0; + virtual void ComputeMeanEla(IssmDouble* paltitude,int* pcounter)=0; virtual void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters,Inputs* inputsin)=0; virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int M,int N,int interp)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 27846) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 27847) @@ -5283,2 +5283,26 @@ /*}}}*/ #endif +void Penta::ComputeMeanEla(IssmDouble* paltitude, int* pcounter){/*{{{*/ + + if(!this->IsOnSurface() & !this->IsIceInElement()) return; + int ipos=0,ineg=0, counter=0; + IssmDouble mean_element_surface=0; + + /* input */ + IssmDouble smb[NUMVERTICES],surface[NUMVERTICES]; + Element::GetInputListOnVertices(&smb[0],SmbMassBalanceEnum); + Element::GetInputListOnVertices(&surface[0],SurfaceEnum); + + for(int iv=0;iv0) ipos=ipos+1; + if(smb[iv]<0) ineg=ineg+1; + } + + /* we define ELA if an element has pos and neg smb */ + if(ipos>0 & ineg>0){ + for(int iv=0;ivparameters->FindParam(&domaintype,DomainTypeEnum); femmodel->parameters->FindParam(&save_results,SaveResultsEnum); femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); @@ -43,12 +44,4 @@ femmodel->inputs->DuplicateInput(VyEnum,VyDebrisEnum); } - femmodel->parameters->SetParam(VxEnum,InputToDepthaverageInEnum); - femmodel->parameters->SetParam(VxAverageEnum,InputToDepthaverageOutEnum); - depthaverage_core(femmodel); - if(domaintype!=Domain2DverticalEnum){ - femmodel->parameters->SetParam(VyEnum,InputToDepthaverageInEnum); - femmodel->parameters->SetParam(VyAverageEnum,InputToDepthaverageOutEnum); - depthaverage_core(femmodel); - } debris_analysis = new DebrisAnalysis(); @@ -59,4 +52,9 @@ extrudefromtop_core(femmodel); + IssmDouble mean_ela; + ComputeMeanElax(&mean_ela,femmodel); + femmodel->parameters->SetParam(mean_ela,SmbMeanElaEnum); + _printf_("core ELA: " << mean_ela << "\n"); + if(save_results) femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); if(solution_type==DebrisSolutionEnum)femmodel->RequestedDependentsx(); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 27846) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 27847) @@ -489,5 +489,18 @@ parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismethod",SmbSemicMethodEnum)); break; - case SMBdebrisMLEnum: + case SMBdebrisEvattEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.smb.qlaps",SmbDesfacEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.dslaps",SmbSWlapsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.dllaps",SmbLWlapsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.windspeedlaps",SmbWindspeedlapsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.humiditylaps",SmbHumiditylapsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.icealbedo",SmbIcealbedoEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.snowalbedo",SmbSnowalbedoEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.debrisalbedo",SmbDebrisalbedoEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismonthly",SmbDebrisIsMonthlyEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.isAnderson",SmbDebrisIsAndersonEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.iscryokarst",SmbDebrisIsCryokarstEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.AndersonD0",SmbDebrisAndersonD0Enum)); break; default: Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 27846) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 27847) @@ -504,145 +504,9 @@ }/*}}}*/ -void SmbDebrisMLx(FemModel* femmodel){/*{{{*/ - - // The function is based on: - // Evatt GW, Abrahams ID, Heil M, Mayer C, Kingslake J, Mitchell SL, et al. Glacial melt under a porous debris layer. Journal of Glaciology 61 (2015) 825–836, doi:10.3189/2 - // Constants/Values are taken from Mayer, Licciulli (2021): https://www.frontiersin.org/articles/10.3389/feart.2021.710276/full#B7 - // function taken from https://github.com/carlolic/DebrisExp/blob/main/USFs/USF_DebrisCoverage.f90 - - /*Intermediaries*/ - // altitude gradients of the crucial parameters (radiation from Marty et al., TaAClimat; 2002) - IssmDouble LW=2.9; // W/m^2 /100m 2.9 - IssmDouble SW=1.3; // W/m^2 /100m 1.3 - IssmDouble HumidityG=0; // % /100m rough estimate - IssmDouble AirTemp=0.7; // C /100m - IssmDouble WindSpeed=0.02; // m/s /100m rough estimate 0.2 - - // accumulation follows a linear increase above the ELA up to a plateau - IssmDouble AccG=0.1; // m w.e. /100m - IssmDouble AccMax=1.; // m w.e. - IssmDouble ReferenceElevation; - IssmDouble AblationDays=120.; // - - IssmDouble In=100.; // Wm^-2 incoming long wave - IssmDouble Q=500.; // Wm^-2 incoming short wave - IssmDouble K=0.585; // Wm^-1K^-1 thermal conductivity 0.585 - IssmDouble Qm=0.0012; // kg m^-3 measured humiditiy level - IssmDouble Qh=0.006 ; // kg m^-3 saturated humidity level - IssmDouble Tm=2.; // C air temperature - IssmDouble Rhoaa=1.22; // kgm^-3 air densitiy - IssmDouble Um=1.5; // ms^-1 measured wind speed - IssmDouble Xm=1.5; // ms^-1 measurement height - IssmDouble Xr=0.01; // ms^-1 surface roughness 0.01 - IssmDouble Alphad=0.07; // debris albedo 0.07 - IssmDouble Alphai=0.4; // ice ablbedo - IssmDouble Alphaeff; - IssmDouble Ustar=0.16; // ms^-1 friction velocity 0.16 - IssmDouble Ca=1000.; // jkg^-1K^-1 specific heat capacity of air - IssmDouble Lm;//=3.34E+05; // jkg^-1K^-1 latent heat of ice melt - IssmDouble Lv=2.50E+06; // jkg^-1K^-1 latent heat of evaporation - IssmDouble Tf=273.; // K water freeezing temperature - IssmDouble Eps=0.95; // thermal emissivity - IssmDouble Rhoi=900.; // kgm^-3 ice density - IssmDouble Sigma=5.67E-08; // Wm^-2K^-4 Stefan Boltzmann constant - IssmDouble Kstar=0.4; // von kármán constant - IssmDouble Gamma=180.; // m^-1 wind speed attenuation 234 - IssmDouble PhiD;//=0.005; // debris packing fraction 0.01 - IssmDouble Humidity=0.2; // relative humidity - - IssmDouble smb,yts,z,debris; - IssmDouble MassBalanceCmDayDebris,MassBalanceMYearDebris; - bool isdebris; - int domaintype; - femmodel->parameters->FindParam(&isdebris,TransientIsdebrisEnum); - - /*Get material parameters and constants */ - //femmodel->parameters->FindParam(&Rhoi,MaterialsRhoIceEnum); // Note Carlo's model used as benchmark was run with different densities for debris and FS - femmodel->parameters->FindParam(&Lm,MaterialsLatentheatEnum); - femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); - PhiD=0.; - if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum); - - /* Loop over all the elements of this partition */ - for(Object* & object : femmodel->elements->objects){ - Element* element=xDynamicCast(object); - - /* Allocate all arrays */ - int numvertices=element->GetNumberOfVertices(); - IssmDouble* surfacelist=xNew(numvertices); - IssmDouble* smb=xNew(numvertices); - IssmDouble* debriscover=xNew(numvertices); - element->GetInputListOnVertices(surfacelist,SurfaceEnum); - - /* Get inputs */ - element->GetInputListOnVertices(debriscover,DebrisThicknessEnum); - element->FindParam(&domaintype,DomainTypeEnum); - - /*Loop over all vertices of element and calculate SMB as function of Debris Cover and z */ - for(int v=0;vnmax){ - Alphaeff=Alphad; - } else { - Alphaeff=Alphai+n*(Alphad-Alphai)/nmax; - } - ReferenceElevation=3200.; // m HEF - - - Alphaeff=Alphad; - ReferenceElevation=2200.; // m M&L - - MassBalanceCmDayDebris=(((In-(z-ReferenceElevation)*LW/100.)-(Eps*Sigma*(Tf*Tf*Tf*Tf))+ - (Q+(z-ReferenceElevation)*SW/100.)*(1.-Alphaeff)+ - (Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* - WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))*(Tm-(z- - ReferenceElevation)*AirTemp/100.))/((1-PhiD)*Rhoi*Lm)/(1.+ - ((Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* - WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/ - K*debris)-(Lv*Ustar*Ustar*(Qh-(Qh*(Humidity-(z- - ReferenceElevation)*HumidityG/100.)))*(exp(-Gamma*Xr)))/((1.-PhiD)* - Rhoi*Lm*Ustar)/((((Um-(z-ReferenceElevation)*WindSpeed/100.) - -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)))*100.*24.*60.*60.; - } - - /* account form ablation days, and convert to m/s */ - MassBalanceMYearDebris=-MassBalanceCmDayDebris/100.*AblationDays/yts; - - /*Update array accordingly*/ - smb[v]=MassBalanceMYearDebris; - } - - /*Add input to element and Free memory*/ - element->AddInput(SmbMassBalanceEnum,smb,P1Enum); - xDelete(surfacelist); - xDelete(smb); - xDelete(debriscover); - } +void SmbDebrisEvattx(FemModel* femmodel){/*{{{*/ + for(Object* & object : femmodel->elements->objects){ + Element* element=xDynamicCast(object); + element->SmbDebrisEvatt(); + } }/*}}}*/ void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 27846) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 27847) @@ -23,5 +23,5 @@ void SmbMeltComponentsx(FemModel* femmodel); void SmbGradientsComponentsx(FemModel* femmodel); -void SmbDebrisMLx(FemModel* femmodel); +void SmbDebrisEvattx(FemModel* femmodel); /* SEMIC: */ void SmbSemicx(FemModel* femmodel, int ismethod); Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 27846) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 27847) @@ -86,4 +86,5 @@ #include "./ResetConstraintsx/ResetConstraintsx.h" #include "./ResetFSBasalBoundaryConditionx/ResetFSBasalBoundaryConditionx.h" +#include "./ComputeMeanElax/ComputeMeanElax.h" #include "./RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h" #include "./RheologyBAbsGradientx/RheologyBAbsGradientx.h" Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27846) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27847) @@ -176,4 +176,5 @@ syn keyword cConstant DebrisRemovalStressThresholdEnum syn keyword cConstant DebrisPackingFractionEnum +syn keyword cConstant DebrisMaxdisplacementvelocityEnum syn keyword cConstant DebugProfilingEnum syn keyword cConstant DomainDimensionEnum @@ -539,4 +540,9 @@ syn keyword cConstant SmbARMAmaOrderEnum syn keyword cConstant SmbAveragingEnum +syn keyword cConstant SmbDebrisalbedoEnum +syn keyword cConstant SmbDebrisIsMonthlyEnum +syn keyword cConstant SmbDebrisIsAndersonEnum +syn keyword cConstant SmbDebrisIsCryokarstEnum +syn keyword cConstant SmbDebrisAndersonD0Enum syn keyword cConstant SmbDesfacEnum syn keyword cConstant SmbDesfacElevEnum @@ -552,5 +558,7 @@ syn keyword cConstant SmbEIdxEnum syn keyword cConstant SmbFEnum +syn keyword cConstant SmbHumiditylapsEnum syn keyword cConstant SmbInitDensityScalingEnum +syn keyword cConstant SmbIcealbedoEnum syn keyword cConstant SmbIsaccumulationEnum syn keyword cConstant SmbIsalbedoEnum @@ -572,4 +580,8 @@ syn keyword cConstant SmbKEnum syn keyword cConstant SmbLapseRatesEnum +syn keyword cConstant SmbLWlapsEnum +syn keyword cConstant SmbSWlapsEnum +syn keyword cConstant SmbSnowalbedoEnum +syn keyword cConstant SmbMeanElaEnum syn keyword cConstant SmbNumBasinsEnum syn keyword cConstant SmbNumBreaksEnum @@ -612,4 +624,5 @@ syn keyword cConstant SmbTemperaturesReconstructedYearsEnum syn keyword cConstant SmbPrecipitationsReconstructedYearsEnum +syn keyword cConstant SmbWindspeedlapsEnum syn keyword cConstant SmoothThicknessMultiplierEnum syn keyword cConstant SolutionTypeEnum @@ -779,4 +792,5 @@ syn keyword cConstant DamageFEnum syn keyword cConstant DebrisThicknessEnum +syn keyword cConstant DebrisThicknessOldEnum syn keyword cConstant DegreeOfChannelizationEnum syn keyword cConstant DepthBelowSurfaceEnum @@ -1107,4 +1121,8 @@ syn keyword cConstant SmbMeltEnum syn keyword cConstant SmbMonthlytemperaturesEnum +syn keyword cConstant SmbMonthlydsradiationEnum +syn keyword cConstant SmbMonthlydlradiationEnum +syn keyword cConstant SmbMonthlywindspeedEnum +syn keyword cConstant SmbMonthlyairhumidityEnum syn keyword cConstant SmbMSurfEnum syn keyword cConstant SmbNetLWEnum @@ -1116,4 +1134,8 @@ syn keyword cConstant SmbPrecipitationEnum syn keyword cConstant SmbPrecipitationsAnomalyEnum +syn keyword cConstant SmbDsradiationAnomalyEnum +syn keyword cConstant SmbDlradiationAnomalyEnum +syn keyword cConstant SmbWindspeedAnomalyEnum +syn keyword cConstant SmbAirhumidityAnomalyEnum syn keyword cConstant SmbPrecipitationsLgmEnum syn keyword cConstant SmbPrecipitationsPresentdayEnum @@ -1135,4 +1157,7 @@ syn keyword cConstant SmbSmbrefEnum syn keyword cConstant SmbSzaValueEnum +syn keyword cConstant SmbSummerMeltEnum +syn keyword cConstant SmbSummerAlbedoEnum +syn keyword cConstant SmbSnowheightEnum syn keyword cConstant SmbTEnum syn keyword cConstant SmbTaEnum @@ -1654,5 +1679,5 @@ syn keyword cConstant SMBarmaEnum syn keyword cConstant SMBcomponentsEnum -syn keyword cConstant SMBdebrisMLEnum +syn keyword cConstant SMBdebrisEvattEnum syn keyword cConstant SMBd18opddEnum syn keyword cConstant SMBforcingEnum @@ -1780,4 +1805,5 @@ syn keyword cType Cfsurfacesquaretransient syn keyword cType Channel +syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1787,6 +1813,6 @@ syn keyword cType ControlParam syn keyword cType Covertree +syn keyword cType DatasetInput syn keyword cType DataSetParam -syn keyword cType DatasetInput syn keyword cType Definition syn keyword cType DependentObject @@ -1801,6 +1827,6 @@ syn keyword cType ElementInput syn keyword cType ElementMatrix +syn keyword cType Elements syn keyword cType ElementVector -syn keyword cType Elements syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1809,9 +1835,10 @@ syn keyword cType Friction syn keyword cType Gauss +syn keyword cType GaussianVariogram +syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria -syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1830,4 +1857,5 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface +syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1840,4 +1868,5 @@ syn keyword cType Matice syn keyword cType Matlitho +syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1852,6 +1881,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType Options syn keyword cType OptionUtilities -syn keyword cType Options syn keyword cType Param syn keyword cType Parameters @@ -1867,11 +1896,11 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType Riftfront syn keyword cType RiftStruct -syn keyword cType Riftfront syn keyword cType SealevelGeometry syn keyword cType Seg syn keyword cType SegInput +syn keyword cType Segment syn keyword cType SegRef -syn keyword cType Segment syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1892,8 +1921,4 @@ syn keyword cType Vertex syn keyword cType Vertices -syn keyword cType classes -syn keyword cType gaussobjects -syn keyword cType krigingobjects -syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27846) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27847) @@ -170,4 +170,5 @@ DebrisRemovalStressThresholdEnum, DebrisPackingFractionEnum, + DebrisMaxdisplacementvelocityEnum, DebugProfilingEnum, DomainDimensionEnum, @@ -530,12 +531,17 @@ SmbAlbedoLandEnum, SmbARMATimestepEnum, - SmbARMAarOrderEnum, - SmbARMAmaOrderEnum, + SmbARMAarOrderEnum, + SmbARMAmaOrderEnum, SmbAveragingEnum, + SmbDebrisalbedoEnum, + SmbDebrisIsMonthlyEnum, + SmbDebrisIsAndersonEnum, + SmbDebrisIsCryokarstEnum, + SmbDebrisAndersonD0Enum, SmbDesfacEnum, SmbDesfacElevEnum, SmbDpermilEnum, SmbDsnowIdxEnum, - SmbElevationBinsEnum, + SmbElevationBinsEnum, SmbCldFracEnum, SmbDelta18oEnum, @@ -546,5 +552,7 @@ SmbEIdxEnum, SmbFEnum, + SmbHumiditylapsEnum, SmbInitDensityScalingEnum, + SmbIcealbedoEnum, SmbIsaccumulationEnum, SmbIsalbedoEnum, @@ -565,5 +573,9 @@ SmbIsturbulentfluxEnum, SmbKEnum, - SmbLapseRatesEnum, + SmbLapseRatesEnum, + SmbLWlapsEnum, + SmbSWlapsEnum, + SmbSnowalbedoEnum, + SmbMeanElaEnum, SmbNumBasinsEnum, SmbNumBreaksEnum, @@ -606,4 +618,5 @@ SmbTemperaturesReconstructedYearsEnum, SmbPrecipitationsReconstructedYearsEnum, + SmbWindspeedlapsEnum, SmoothThicknessMultiplierEnum, SolutionTypeEnum, @@ -775,4 +788,5 @@ DamageFEnum, DebrisThicknessEnum, + DebrisThicknessOldEnum, DegreeOfChannelizationEnum, DepthBelowSurfaceEnum, @@ -1043,5 +1057,5 @@ SmbAdiffiniEnum, SmbAiniEnum, - SmbARMANoiseEnum, + SmbARMANoiseEnum, SmbBasinsIdEnum, SmbBMaxEnum, @@ -1104,4 +1118,8 @@ SmbMInitnum, SmbMonthlytemperaturesEnum, + SmbMonthlydsradiationEnum, + SmbMonthlydlradiationEnum, + SmbMonthlywindspeedEnum, + SmbMonthlyairhumidityEnum, SmbMSurfEnum, SmbNetLWEnum, @@ -1113,4 +1131,8 @@ SmbPrecipitationEnum, SmbPrecipitationsAnomalyEnum, + SmbDsradiationAnomalyEnum, + SmbDlradiationAnomalyEnum, + SmbWindspeedAnomalyEnum, + SmbAirhumidityAnomalyEnum, SmbPrecipitationsLgmEnum, SmbPrecipitationsPresentdayEnum, @@ -1132,4 +1154,7 @@ SmbSmbrefEnum, SmbSzaValueEnum, + SmbSummerMeltEnum, + SmbSummerAlbedoEnum, + SmbSnowheightEnum, SmbTEnum, SmbTaEnum, @@ -1653,5 +1678,5 @@ SMBarmaEnum, SMBcomponentsEnum, - SMBdebrisMLEnum, + SMBdebrisEvattEnum, SMBd18opddEnum, SMBforcingEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27846) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27847) @@ -178,4 +178,5 @@ case DebrisRemovalStressThresholdEnum : return "DebrisRemovalStressThreshold"; case DebrisPackingFractionEnum : return "DebrisPackingFraction"; + case DebrisMaxdisplacementvelocityEnum : return "DebrisMaxdisplacementvelocity"; case DebugProfilingEnum : return "DebugProfiling"; case DomainDimensionEnum : return "DomainDimension"; @@ -541,4 +542,9 @@ case SmbARMAmaOrderEnum : return "SmbARMAmaOrder"; case SmbAveragingEnum : return "SmbAveraging"; + case SmbDebrisalbedoEnum : return "SmbDebrisalbedo"; + case SmbDebrisIsMonthlyEnum : return "SmbDebrisIsMonthly"; + case SmbDebrisIsAndersonEnum : return "SmbDebrisIsAnderson"; + case SmbDebrisIsCryokarstEnum : return "SmbDebrisIsCryokarst"; + case SmbDebrisAndersonD0Enum : return "SmbDebrisAndersonD0"; case SmbDesfacEnum : return "SmbDesfac"; case SmbDesfacElevEnum : return "SmbDesfacElev"; @@ -554,5 +560,7 @@ case SmbEIdxEnum : return "SmbEIdx"; case SmbFEnum : return "SmbF"; + case SmbHumiditylapsEnum : return "SmbHumiditylaps"; case SmbInitDensityScalingEnum : return "SmbInitDensityScaling"; + case SmbIcealbedoEnum : return "SmbIcealbedo"; case SmbIsaccumulationEnum : return "SmbIsaccumulation"; case SmbIsalbedoEnum : return "SmbIsalbedo"; @@ -574,4 +582,8 @@ case SmbKEnum : return "SmbK"; case SmbLapseRatesEnum : return "SmbLapseRates"; + case SmbLWlapsEnum : return "SmbLWlaps"; + case SmbSWlapsEnum : return "SmbSWlaps"; + case SmbSnowalbedoEnum : return "SmbSnowalbedo"; + case SmbMeanElaEnum : return "SmbMeanEla"; case SmbNumBasinsEnum : return "SmbNumBasins"; case SmbNumBreaksEnum : return "SmbNumBreaks"; @@ -614,4 +626,5 @@ case SmbTemperaturesReconstructedYearsEnum : return "SmbTemperaturesReconstructedYears"; case SmbPrecipitationsReconstructedYearsEnum : return "SmbPrecipitationsReconstructedYears"; + case SmbWindspeedlapsEnum : return "SmbWindspeedlaps"; case SmoothThicknessMultiplierEnum : return "SmoothThicknessMultiplier"; case SolutionTypeEnum : return "SolutionType"; @@ -781,4 +794,5 @@ case DamageFEnum : return "DamageF"; case DebrisThicknessEnum : return "DebrisThickness"; + case DebrisThicknessOldEnum : return "DebrisThicknessOld"; case DegreeOfChannelizationEnum : return "DegreeOfChannelization"; case DepthBelowSurfaceEnum : return "DepthBelowSurface"; @@ -1109,4 +1123,8 @@ case SmbMeltEnum : return "SmbMelt"; case SmbMonthlytemperaturesEnum : return "SmbMonthlytemperatures"; + case SmbMonthlydsradiationEnum : return "SmbMonthlydsradiation"; + case SmbMonthlydlradiationEnum : return "SmbMonthlydlradiation"; + case SmbMonthlywindspeedEnum : return "SmbMonthlywindspeed"; + case SmbMonthlyairhumidityEnum : return "SmbMonthlyairhumidity"; case SmbMSurfEnum : return "SmbMSurf"; case SmbNetLWEnum : return "SmbNetLW"; @@ -1118,4 +1136,8 @@ case SmbPrecipitationEnum : return "SmbPrecipitation"; case SmbPrecipitationsAnomalyEnum : return "SmbPrecipitationsAnomaly"; + case SmbDsradiationAnomalyEnum : return "SmbDsradiationAnomaly"; + case SmbDlradiationAnomalyEnum : return "SmbDlradiationAnomaly"; + case SmbWindspeedAnomalyEnum : return "SmbWindspeedAnomaly"; + case SmbAirhumidityAnomalyEnum : return "SmbAirhumidityAnomaly"; case SmbPrecipitationsLgmEnum : return "SmbPrecipitationsLgm"; case SmbPrecipitationsPresentdayEnum : return "SmbPrecipitationsPresentday"; @@ -1137,4 +1159,7 @@ case SmbSmbrefEnum : return "SmbSmbref"; case SmbSzaValueEnum : return "SmbSzaValue"; + case SmbSummerMeltEnum : return "SmbSummerMelt"; + case SmbSummerAlbedoEnum : return "SmbSummerAlbedo"; + case SmbSnowheightEnum : return "SmbSnowheight"; case SmbTEnum : return "SmbT"; case SmbTaEnum : return "SmbTa"; @@ -1656,5 +1681,5 @@ case SMBarmaEnum : return "SMBarma"; case SMBcomponentsEnum : return "SMBcomponents"; - case SMBdebrisMLEnum : return "SMBdebrisML"; + case SMBdebrisEvattEnum : return "SMBdebrisEvatt"; case SMBd18opddEnum : return "SMBd18opdd"; case SMBforcingEnum : return "SMBforcing"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27846) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27847) @@ -169,4 +169,5 @@ syn keyword juliaConstC DebrisRemovalStressThresholdEnum syn keyword juliaConstC DebrisPackingFractionEnum +syn keyword juliaConstC DebrisMaxdisplacementvelocityEnum syn keyword juliaConstC DebugProfilingEnum syn keyword juliaConstC DomainDimensionEnum @@ -532,4 +533,9 @@ syn keyword juliaConstC SmbARMAmaOrderEnum syn keyword juliaConstC SmbAveragingEnum +syn keyword juliaConstC SmbDebrisalbedoEnum +syn keyword juliaConstC SmbDebrisIsMonthlyEnum +syn keyword juliaConstC SmbDebrisIsAndersonEnum +syn keyword juliaConstC SmbDebrisIsCryokarstEnum +syn keyword juliaConstC SmbDebrisAndersonD0Enum syn keyword juliaConstC SmbDesfacEnum syn keyword juliaConstC SmbDesfacElevEnum @@ -545,5 +551,7 @@ syn keyword juliaConstC SmbEIdxEnum syn keyword juliaConstC SmbFEnum +syn keyword juliaConstC SmbHumiditylapsEnum syn keyword juliaConstC SmbInitDensityScalingEnum +syn keyword juliaConstC SmbIcealbedoEnum syn keyword juliaConstC SmbIsaccumulationEnum syn keyword juliaConstC SmbIsalbedoEnum @@ -565,4 +573,8 @@ syn keyword juliaConstC SmbKEnum syn keyword juliaConstC SmbLapseRatesEnum +syn keyword juliaConstC SmbLWlapsEnum +syn keyword juliaConstC SmbSWlapsEnum +syn keyword juliaConstC SmbSnowalbedoEnum +syn keyword juliaConstC SmbMeanElaEnum syn keyword juliaConstC SmbNumBasinsEnum syn keyword juliaConstC SmbNumBreaksEnum @@ -605,4 +617,5 @@ syn keyword juliaConstC SmbTemperaturesReconstructedYearsEnum syn keyword juliaConstC SmbPrecipitationsReconstructedYearsEnum +syn keyword juliaConstC SmbWindspeedlapsEnum syn keyword juliaConstC SmoothThicknessMultiplierEnum syn keyword juliaConstC SolutionTypeEnum @@ -772,4 +785,5 @@ syn keyword juliaConstC DamageFEnum syn keyword juliaConstC DebrisThicknessEnum +syn keyword juliaConstC DebrisThicknessOldEnum syn keyword juliaConstC DegreeOfChannelizationEnum syn keyword juliaConstC DepthBelowSurfaceEnum @@ -1100,4 +1114,8 @@ syn keyword juliaConstC SmbMeltEnum syn keyword juliaConstC SmbMonthlytemperaturesEnum +syn keyword juliaConstC SmbMonthlydsradiationEnum +syn keyword juliaConstC SmbMonthlydlradiationEnum +syn keyword juliaConstC SmbMonthlywindspeedEnum +syn keyword juliaConstC SmbMonthlyairhumidityEnum syn keyword juliaConstC SmbMSurfEnum syn keyword juliaConstC SmbNetLWEnum @@ -1109,4 +1127,8 @@ syn keyword juliaConstC SmbPrecipitationEnum syn keyword juliaConstC SmbPrecipitationsAnomalyEnum +syn keyword juliaConstC SmbDsradiationAnomalyEnum +syn keyword juliaConstC SmbDlradiationAnomalyEnum +syn keyword juliaConstC SmbWindspeedAnomalyEnum +syn keyword juliaConstC SmbAirhumidityAnomalyEnum syn keyword juliaConstC SmbPrecipitationsLgmEnum syn keyword juliaConstC SmbPrecipitationsPresentdayEnum @@ -1128,4 +1150,7 @@ syn keyword juliaConstC SmbSmbrefEnum syn keyword juliaConstC SmbSzaValueEnum +syn keyword juliaConstC SmbSummerMeltEnum +syn keyword juliaConstC SmbSummerAlbedoEnum +syn keyword juliaConstC SmbSnowheightEnum syn keyword juliaConstC SmbTEnum syn keyword juliaConstC SmbTaEnum @@ -1647,5 +1672,5 @@ syn keyword juliaConstC SMBarmaEnum syn keyword juliaConstC SMBcomponentsEnum -syn keyword juliaConstC SMBdebrisMLEnum +syn keyword juliaConstC SMBdebrisEvattEnum syn keyword juliaConstC SMBd18opddEnum syn keyword juliaConstC SMBforcingEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27846) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27847) @@ -181,4 +181,5 @@ else if (strcmp(name,"DebrisRemovalStressThreshold")==0) return DebrisRemovalStressThresholdEnum; else if (strcmp(name,"DebrisPackingFraction")==0) return DebrisPackingFractionEnum; + else if (strcmp(name,"DebrisMaxdisplacementvelocity")==0) return DebrisMaxdisplacementvelocityEnum; else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum; else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"Hydrologyarmadatebreaks")==0) return HydrologyarmadatebreaksEnum; else if (strcmp(name,"Hydrologyarmamalagcoefs")==0) return HydrologyarmamalagcoefsEnum; - else if (strcmp(name,"HydrologyarmamaOrder")==0) return HydrologyarmamaOrderEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"HydrologyarmaMonthlyFactors")==0) return HydrologyarmaMonthlyFactorsEnum; + if (strcmp(name,"HydrologyarmamaOrder")==0) return HydrologyarmamaOrderEnum; + else if (strcmp(name,"HydrologyarmaMonthlyFactors")==0) return HydrologyarmaMonthlyFactorsEnum; else if (strcmp(name,"HydrologyarmaNumBreaks")==0) return HydrologyarmaNumBreaksEnum; else if (strcmp(name,"HydrologyarmaNumParams")==0) return HydrologyarmaNumParamsEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum; else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum; - else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum; + if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum; + else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum; else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum; else if (strcmp(name,"MaterialsEarthDensity")==0) return MaterialsEarthDensityEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum; else if (strcmp(name,"SolidearthSettingsTimeAcc")==0) return SolidearthSettingsTimeAccEnum; - else if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"SolidearthSettingsMaxiter")==0) return SolidearthSettingsMaxiterEnum; + if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum; + else if (strcmp(name,"SolidearthSettingsMaxiter")==0) return SolidearthSettingsMaxiterEnum; else if (strcmp(name,"SolidearthSettingsGrdOcean")==0) return SolidearthSettingsGrdOceanEnum; else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum; @@ -553,4 +554,9 @@ else if (strcmp(name,"SmbARMAmaOrder")==0) return SmbARMAmaOrderEnum; else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum; + else if (strcmp(name,"SmbDebrisalbedo")==0) return SmbDebrisalbedoEnum; + else if (strcmp(name,"SmbDebrisIsMonthly")==0) return SmbDebrisIsMonthlyEnum; + else if (strcmp(name,"SmbDebrisIsAnderson")==0) return SmbDebrisIsAndersonEnum; + else if (strcmp(name,"SmbDebrisIsCryokarst")==0) return SmbDebrisIsCryokarstEnum; + else if (strcmp(name,"SmbDebrisAndersonD0")==0) return SmbDebrisAndersonD0Enum; else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum; else if (strcmp(name,"SmbDesfacElev")==0) return SmbDesfacElevEnum; @@ -566,5 +572,7 @@ else if (strcmp(name,"SmbEIdx")==0) return SmbEIdxEnum; else if (strcmp(name,"SmbF")==0) return SmbFEnum; + else if (strcmp(name,"SmbHumiditylaps")==0) return SmbHumiditylapsEnum; else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum; + else if (strcmp(name,"SmbIcealbedo")==0) return SmbIcealbedoEnum; else if (strcmp(name,"SmbIsaccumulation")==0) return SmbIsaccumulationEnum; else if (strcmp(name,"SmbIsalbedo")==0) return SmbIsalbedoEnum; @@ -586,4 +594,8 @@ else if (strcmp(name,"SmbK")==0) return SmbKEnum; else if (strcmp(name,"SmbLapseRates")==0) return SmbLapseRatesEnum; + else if (strcmp(name,"SmbLWlaps")==0) return SmbLWlapsEnum; + else if (strcmp(name,"SmbSWlaps")==0) return SmbSWlapsEnum; + else if (strcmp(name,"SmbSnowalbedo")==0) return SmbSnowalbedoEnum; + else if (strcmp(name,"SmbMeanEla")==0) return SmbMeanElaEnum; else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum; else if (strcmp(name,"SmbNumBreaks")==0) return SmbNumBreaksEnum; @@ -617,5 +629,8 @@ else if (strcmp(name,"SmbSemicTmax")==0) return SmbSemicTmaxEnum; else if (strcmp(name,"SmbStepsPerStep")==0) return SmbStepsPerStepEnum; - else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum; else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum; else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum; @@ -626,11 +641,9 @@ else if (strcmp(name,"SmbTemperaturesReconstructedYears")==0) return SmbTemperaturesReconstructedYearsEnum; else if (strcmp(name,"SmbPrecipitationsReconstructedYears")==0) return SmbPrecipitationsReconstructedYearsEnum; + else if (strcmp(name,"SmbWindspeedlaps")==0) return SmbWindspeedlapsEnum; else if (strcmp(name,"SmoothThicknessMultiplier")==0) return SmoothThicknessMultiplierEnum; else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum; else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum; + else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum; else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum; else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; @@ -739,5 +752,8 @@ else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum; else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum; - else if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum; else if (strcmp(name,"BasalforcingsSpatialUpperwaterMeltingRate")==0) return BasalforcingsSpatialUpperwaterMeltingRateEnum; else if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum; @@ -752,8 +768,5 @@ else if (strcmp(name,"BasalforcingsPicoOverturningCoeff")==0) return BasalforcingsPicoOverturningCoeffEnum; else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum; + else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum; else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum; else if (strcmp(name,"BasalStressx")==0) return BasalStressxEnum; @@ -799,4 +812,5 @@ else if (strcmp(name,"DamageF")==0) return DamageFEnum; else if (strcmp(name,"DebrisThickness")==0) return DebrisThicknessEnum; + else if (strcmp(name,"DebrisThicknessOld")==0) return DebrisThicknessOldEnum; else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum; else if (strcmp(name,"DepthBelowSurface")==0) return DepthBelowSurfaceEnum; @@ -861,5 +875,8 @@ else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum; else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum; - else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; @@ -875,8 +892,5 @@ else if (strcmp(name,"UGia")==0) return UGiaEnum; else if (strcmp(name,"UGiaRate")==0) return UGiaRateEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"Gradient")==0) return GradientEnum; + else if (strcmp(name,"Gradient")==0) return GradientEnum; else if (strcmp(name,"GroundinglineHeight")==0) return GroundinglineHeightEnum; else if (strcmp(name,"HydraulicPotential")==0) return HydraulicPotentialEnum; @@ -984,5 +998,8 @@ else if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; - else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum; else if (strcmp(name,"SatGraviGRD")==0) return SatGraviGRDEnum; else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum; @@ -998,8 +1015,5 @@ else if (strcmp(name,"SealevelBarystaticHydroLatbar")==0) return SealevelBarystaticHydroLatbarEnum; else if (strcmp(name,"SealevelBarystaticHydroLongbar")==0) return SealevelBarystaticHydroLongbarEnum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum; + else if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum; else if (strcmp(name,"SealevelBarystaticBpMask")==0) return SealevelBarystaticBpMaskEnum; else if (strcmp(name,"SealevelBarystaticBpWeights")==0) return SealevelBarystaticBpWeightsEnum; @@ -1107,5 +1121,8 @@ else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum; else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum; - else if (strcmp(name,"SmbEla")==0) return SmbElaEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"SmbEla")==0) return SmbElaEnum; else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum; else if (strcmp(name,"SmbFAC")==0) return SmbFACEnum; @@ -1121,8 +1138,5 @@ else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum; else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; + else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; else if (strcmp(name,"SmbMassBalanceSnow")==0) return SmbMassBalanceSnowEnum; else if (strcmp(name,"SmbMassBalanceIce")==0) return SmbMassBalanceIceEnum; @@ -1136,4 +1150,8 @@ else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum; else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum; + else if (strcmp(name,"SmbMonthlydsradiation")==0) return SmbMonthlydsradiationEnum; + else if (strcmp(name,"SmbMonthlydlradiation")==0) return SmbMonthlydlradiationEnum; + else if (strcmp(name,"SmbMonthlywindspeed")==0) return SmbMonthlywindspeedEnum; + else if (strcmp(name,"SmbMonthlyairhumidity")==0) return SmbMonthlyairhumidityEnum; else if (strcmp(name,"SmbMSurf")==0) return SmbMSurfEnum; else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum; @@ -1145,4 +1163,8 @@ else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum; else if (strcmp(name,"SmbPrecipitationsAnomaly")==0) return SmbPrecipitationsAnomalyEnum; + else if (strcmp(name,"SmbDsradiationAnomaly")==0) return SmbDsradiationAnomalyEnum; + else if (strcmp(name,"SmbDlradiationAnomaly")==0) return SmbDlradiationAnomalyEnum; + else if (strcmp(name,"SmbWindspeedAnomaly")==0) return SmbWindspeedAnomalyEnum; + else if (strcmp(name,"SmbAirhumidityAnomaly")==0) return SmbAirhumidityAnomalyEnum; else if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum; else if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum; @@ -1164,4 +1186,7 @@ else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum; else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum; + else if (strcmp(name,"SmbSummerMelt")==0) return SmbSummerMeltEnum; + else if (strcmp(name,"SmbSummerAlbedo")==0) return SmbSummerAlbedoEnum; + else if (strcmp(name,"SmbSnowheight")==0) return SmbSnowheightEnum; else if (strcmp(name,"SmbT")==0) return SmbTEnum; else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; @@ -1219,5 +1244,8 @@ else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum; else if (strcmp(name,"Surface")==0) return SurfaceEnum; - else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum; else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; @@ -1244,8 +1272,5 @@ else if (strcmp(name,"ThicknessResidual")==0) return ThicknessResidualEnum; else if (strcmp(name,"TransientAccumulatedDeltaIceThickness")==0) return TransientAccumulatedDeltaIceThicknessEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"Vel")==0) return VelEnum; + else if (strcmp(name,"Vel")==0) return VelEnum; else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; else if (strcmp(name,"VxBase")==0) return VxBaseEnum; @@ -1342,5 +1367,8 @@ else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum; else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum; - else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum; else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum; else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum; @@ -1367,8 +1395,5 @@ else if (strcmp(name,"Outputdefinition82")==0) return Outputdefinition82Enum; else if (strcmp(name,"Outputdefinition83")==0) return Outputdefinition83Enum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum; + else if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum; else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum; else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum; @@ -1465,5 +1490,8 @@ else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum; - else if (strcmp(name,"Dense")==0) return DenseEnum; + else stage=13; + } + if(stage==13){ + if (strcmp(name,"Dense")==0) return DenseEnum; else if (strcmp(name,"DependentObject")==0) return DependentObjectEnum; else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; @@ -1490,8 +1518,5 @@ else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum; else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; - else stage=13; - } - if(stage==13){ - if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; + else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum; else if (strcmp(name,"FSSolver")==0) return FSSolverEnum; @@ -1588,5 +1613,8 @@ else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum; else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum; - else if (strcmp(name,"LoveKf")==0) return LoveKfEnum; + else stage=14; + } + if(stage==14){ + if (strcmp(name,"LoveKf")==0) return LoveKfEnum; else if (strcmp(name,"LoveKt")==0) return LoveKtEnum; else if (strcmp(name,"LoveLf")==0) return LoveLfEnum; @@ -1613,8 +1641,5 @@ else if (strcmp(name,"Materials")==0) return MaterialsEnum; else if (strcmp(name,"Matestar")==0) return MatestarEnum; - else stage=14; - } - if(stage==14){ - if (strcmp(name,"Matice")==0) return MaticeEnum; + else if (strcmp(name,"Matice")==0) return MaticeEnum; else if (strcmp(name,"Matlitho")==0) return MatlithoEnum; else if (strcmp(name,"Mathydro")==0) return MathydroEnum; @@ -1695,5 +1720,5 @@ else if (strcmp(name,"SMBarma")==0) return SMBarmaEnum; else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; - else if (strcmp(name,"SMBdebrisML")==0) return SMBdebrisMLEnum; + else if (strcmp(name,"SMBdebrisEvatt")==0) return SMBdebrisEvattEnum; else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; @@ -1711,5 +1736,8 @@ else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum; - else if (strcmp(name,"Scaled")==0) return ScaledEnum; + else stage=15; + } + if(stage==15){ + if (strcmp(name,"Scaled")==0) return ScaledEnum; else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum; else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum; @@ -1736,8 +1764,5 @@ else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum; else if (strcmp(name,"Sset")==0) return SsetEnum; - else stage=15; - } - if(stage==15){ - if (strcmp(name,"StatisticsSolution")==0) return StatisticsSolutionEnum; + else if (strcmp(name,"StatisticsSolution")==0) return StatisticsSolutionEnum; else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum; else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum; Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 27846) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 27847) @@ -262,5 +262,5 @@ case 12: return SMBsemicEnum; case 13: return SMBarmaEnum; - case 14: return SMBdebrisMLEnum; + case 14: return SMBdebrisEvattEnum; default: _error_("Marshalled SMB code \""<=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + end + md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]); + md = checkfield(md,'fieldname','smb.averaging', 'numel', [1], 'values', [0, 1, 2]); + md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1); + md = checkfield(md,'fieldname','smb.icealbedo','>=',0,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.snowalbedo','>=',0,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.debrisalbedo','>=',0,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.ismonthly','numel', [1], 'values', [0, 1]); + md = checkfield(md,'fieldname','smb.isAnderson','numel', [1], 'values', [0, 1]); + md = checkfield(md,'fieldname','smb.iscryokarst','numel', [1], 'values', [0, 1]); + + end % }}} + function disp(self) % {{{ + + disp(sprintf('Evatt et al. (2015) Debris Model (doi: 10.3189/2015JoG14J235)')); + disp(sprintf('If isAnderson==1 -> Eq. 6 from Ferguson & Vieli (2021) ist taken (https://doi.org/10.5194/tc-15-3377-2021)')); + disp(sprintf('If iscryokarst==1 -> Eqs. 9,10 from Ferguson & Vieli (2021) are taken (https://doi.org/10.5194/tc-15-3377-2021)')); + disp(sprintf('Clean-ice SMB is taken from the Evatt et al. (2015) EBM with debris=0')); + + fielddisplay(self,'ismonthly','do we use monthly forcing (otherwise yearly averaged forcing is applied) (default is 0)'); + fielddisplay(self,'isAnderson','do we use the Anderson parametrization (default is 0)'); + fielddisplay(self,'iscryokarst','do we use a cryokarst parametrization (default is 0)'); + fielddisplay(self,'monthlytemperatures','monthly surface temperatures [K]'); + fielddisplay(self,'monthlyprecipitation','monthly surface precipitation [m/yr water eq]'); + fielddisplay(self,'monthlydsradiation','monthly downwelling shortwave radiation [W m-2]'); + fielddisplay(self,'monthlydlradiation','monthly downwelling longwave radiation [W m-2]'); + fielddisplay(self,'monthlywindspeed','monthly surface wind speed [m s-1]'); + fielddisplay(self,'monthlyairhumidity','monthly near-surface specific humidity [kg kg-1]'); + fielddisplay(self,'temperature_anomaly','anomaly to monthly reference temperature (additive)'); + fielddisplay(self,'precipitation_anomaly','anomaly to monthly precipitation (multiplicative)'); + fielddisplay(self,'dsradiation_anomaly','anomaly to monthly reference downwelling shortwave radiation'); + fielddisplay(self,'dlradiation_anomaly','anomaly to monthly reference downwelling longwave radiation (additive'); + fielddisplay(self,'windspeed_anomaly','anomaly to monthly reference surface wind speed (additive)'); + fielddisplay(self,'airhumidity_anomaly','anomaly to monthly reference near-surface specific humidity (additive)'); + + fielddisplay(self,'s0t','should be set to elevation from RCM/GCM source (between 0 and a few 1000s m, default is 0) [m]'); + fielddisplay(self,'snowheight','guess of snowheight at the end of the summer, will be further evolved'); + fielddisplay(self,'rlaps','present day temperature lapse rate (default is 7.4 degree/km)'); + fielddisplay(self,'dslaps','present day SW lapse rate (default is 1.3 W/m^2/km)'); + fielddisplay(self,'dllaps','present day LW lapse rate (default is 2.9 W/m^2/km)'); + fielddisplay(self,'windspeedlaps','present day wind speed lapse rate (default is 0.02 m/s/km)'); + fielddisplay(self,'humiditylaps','present day humidity lapse rate (default is 0)'); + fielddisplay(self,'qlaps','precip change (default is 0.1/km'); + fielddisplay(self,'icealbedo','albedo for ice (default is 0.3)'); + fielddisplay(self,'snowalbedo','albedo for snow (default is 0.75)'); + fielddisplay(self,'debrisalbedo','albedo for debris (default is 0.07)'); + fielddisplay(self,'AndersonD0','parameter to represent the debris effect (default is 0.5)'); + fielddisplay(self,'steps_per_step', 'number of smb steps per time step'); + fielddisplay(self,'averaging','averaging methods from short to long steps'); + disp(sprintf('%51s 0: Arithmetic (default)',' ')); + disp(sprintf('%51s 1: Geometric',' ')); + disp(sprintf('%51s 2: Harmonic',' ')); + fielddisplay(self,'requested_outputs','additional outputs requested (TemperaturePDD, SmbAccumulation, SmbMelt, SmbSummerMelt, SmbAlbedo, SmbSummerAlbedo, SmbSnowheight)'); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + yts=md.constants.yts; + + WriteData(fid,prefix,'name','md.smb.model','data',14,'format','Integer'); + + WriteData(fid,prefix,'object',self,'class','smb','fieldname','qlaps','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0t','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','snowheight','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','rlaps','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dslaps','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dllaps','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','windspeedlaps','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','humiditylaps','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','icealbedo','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','snowalbedo','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','debrisalbedo','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','ismonthly','format','Boolean'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','isAnderson','format','Boolean'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','iscryokarst','format','Boolean'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','AndersonD0','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','monthlyprecipitation','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','monthlydsradiation','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','monthlydlradiation','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','monthlywindspeed','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','monthlyairhumidity','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','temperature_anomaly','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','precipitation_anomaly','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dsradiation_anomaly','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dlradiation_anomaly','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','windspeed_anomaly','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','airhumidity_anomaly','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer') + + %process requested outputs + outputs = self.requested_outputs; + pos = find(ismember(outputs,'default')); + if ~isempty(pos), + outputs(pos) = []; %remove 'default' from outputs + outputs = [outputs defaultoutputs(self,md)]; %add defaults + end + WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray'); + + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/debris.m =================================================================== --- /issm/trunk-jpl/src/m/classes/debris.m (revision 27846) +++ /issm/trunk-jpl/src/m/classes/debris.m (revision 27847) @@ -12,4 +12,5 @@ removalmodel = 0; displacementmodel = 0; + max_displacementvelocity = 0; removal_slope_threshold = 0; removal_stress_threshold = 0; @@ -66,5 +67,5 @@ function list = defaultoutputs(self,md) % {{{ - list = {'DebrisThickness'}; + list = {'DebrisThickness','DebrisMaskNodeActivation','VxDebris','VyDebris'}; end % }}} @@ -92,4 +93,7 @@ self.removal_stress_threshold=0; + %Max velocity for displacementmodel (1) + self.max_displacementvelocity=0; + %default output self.requested_outputs={'default'}; @@ -101,5 +105,5 @@ md = checkfield(md,'fieldname','debris.spcthickness'); - md = checkfield(md,'fieldname','debris.stabilization','values',[0 1 2 3]); + md = checkfield(md,'fieldname','debris.stabilization','values',[0 1 2 3 4 5]); md = checkfield(md,'fieldname','debris.min_thickness','>=',0); md = checkfield(md,'fieldname','debris.packingfraction','>=',0); @@ -108,4 +112,5 @@ md = checkfield(md,'fieldname','debris.removal_slope_threshold','>=',0); md = checkfield(md,'fieldname','debris.removal_stress_threshold','>=',0); + md = checkfield(md,'fieldname','debris.max_displacementvelocity','>=',0); md = checkfield(md,'fieldname','debris.requested_outputs','stringrow',1); if ~any(isnan(md.stressbalance.vertex_pairing)), @@ -120,5 +125,6 @@ fielddisplay(self,'stabilization','0: no stabilization, 1: artificial diffusion, 2: streamline upwinding, 3: streamline upwind Petrov-Galerkin (SUPG)'); fielddisplay(self,'removalmodel','frontal removal of debris. 0: no removal, 1: Slope-triggered debris removal, 2: driving-stress triggered debris removal'); - fielddisplay(self,'displacementmodel','debris displacement. 0: no displacement, 1: ...'); + fielddisplay(self,'max_displacementvelocity','maximum velocity of debris transport (v_ice + v_displacement) (m/a)'); + fielddisplay(self,'displacementmodel','debris displacement. 0: no displacement, 1: additional debris velocity above the critical slope/stress threshold'); fielddisplay(self,'removal_slope_threshold','critical slope (degrees) for removalmodel (1)'); fielddisplay(self,'removal_stress_threshold','critical stress (Pa) for removalmodel (2)'); @@ -135,4 +141,5 @@ WriteData(fid,prefix,'object',self,'fieldname','removalmodel','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','displacementmodel','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','max_displacementvelocity','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','removal_slope_threshold','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','removal_stress_threshold','format','Double'); @@ -156,4 +163,5 @@ writejsdouble(fid,[modelname '.debris.removalmodel'],self.removalmodel); writejsdouble(fid,[modelname '.debris.displacementmodel'],self.displacementmodel); + writejsdouble(fid,[modelname '.debris.max_displacementvelocity'],self.displacementmodel); writejsdouble(fid,[modelname '.debris.removal_slope_threshold'],self.removal_slope_threshold); writejsdouble(fid,[modelname '.debris.removal_stress_threshold'],self.removal_stress_threshold); Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 27846) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 27847) @@ -1289,4 +1289,5 @@ if isfield(structmd,'cfl_coefficient'), md.timestepping.cfl_coefficient=structmd.cfl_coefficient; end if isfield(structmd,'spcthickness'), md.masstransport.spcthickness=structmd.spcthickness; end + if isfield(structmd,'spcthickness'), md.debris.spcthickness=structmd.spcthickness; end if isfield(structmd,'artificial_diffusivity'), md.masstransport.stabilization=structmd.artificial_diffusivity; end if isfield(structmd,'hydrostatic_adjustment'), md.masstransport.hydrostatic_adjustment=structmd.hydrostatic_adjustment; end Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 27846) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 27847) @@ -240,4 +240,8 @@ field = field*yts; elseif strcmp(fieldname,'VyAverage'), + field = field*yts; + elseif strcmp(fieldname,'VxDebris'), + field = field*yts; + elseif strcmp(fieldname,'VyDebris'), field = field*yts; elseif strcmp(fieldname,'BasalforcingsGroundediceMeltingRate'),