Sync upstream -- 2026.02.26

src/c/analyses/SmbAnalysis.cpp

@@ -338,6 +338,8 @@ void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
 			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isfirnwarming",SmbIsfirnwarmingEnum));
 			parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
 			parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
+			parameters->AddObject(iomodel->CopyConstantObject("md.smb.pdd_fac_ice",PddfacIceEnum));
+			parameters->AddObject(iomodel->CopyConstantObject("md.smb.pdd_fac_snow",PddfacSnowEnum));
 			break;
 		case SMBdebrisEvattEnum:
 			parameters->AddObject(iomodel->CopyConstantObject("md.smb.qlaps",SmbDesfacEnum));

src/c/classes/Elements/Element.cpp

@@ -3695,6 +3695,7 @@ void       Element::PositiveDegreeDaySicopolis(bool isfirnwarming){/*{{{*/
 	IssmDouble* p_ampl=xNew<IssmDouble>(NUM_VERTICES);	// precip anomaly
 	IssmDouble* t_ampl=xNew<IssmDouble>(NUM_VERTICES);	// remperature anomaly
 	IssmDouble rho_water,rho_ice,desfac,rlaps;
+	IssmDouble pdd_fac_ice,pdd_fac_snow;
 	IssmDouble inv_twelve=1./12.;								//factor for monthly average
 	IssmDouble time,yts,time_yr;
 
@@ -3709,6 +3710,10 @@ void       Element::PositiveDegreeDaySicopolis(bool isfirnwarming){/*{{{*/
 	desfac=this->FindParam(SmbDesfacEnum);
 	rlaps=this->FindParam(SmbRlapsEnum);
 
+	/*Get pdd melt factors*/
+	pdd_fac_ice=this->FindParam(PddfacIceEnum);
+	pdd_fac_snow=this->FindParam(PddfacSnowEnum);
+
 	/* Get time */
 	this->parameters->FindParam(&time,TimeEnum);
 	this->parameters->FindParam(&yts,ConstantsYtsEnum);
@@ -3747,7 +3752,7 @@ void       Element::PositiveDegreeDaySicopolis(bool isfirnwarming){/*{{{*/
 	for (int iv = 0; iv<NUM_VERTICES; iv++){
 		smb[iv]=PddSurfaceMassBalanceSicopolis(&monthlytemperatures[iv*12], &monthlyprec[iv*12],
 					&melt[iv], &accu[iv], &melt_star[iv], &t_ampl[iv], &p_ampl[iv], yts, s[iv],
-					desfac, s0t[iv], s0p[iv],rlaps,rho_water,rho_ice);
+					desfac, s0t[iv], s0p[iv],rlaps,rho_water,rho_ice,pdd_fac_ice,pdd_fac_snow);
 
 		/* make correction */
 		smb[iv] = smb[iv]+smbcorr[iv];
@@ -5659,17 +5664,24 @@ void       Element::SmbGemb(IssmDouble timeinputs, int count, int steps){/*{{{*/
 		}
 		else pAir=pparam;
 
-		//Hold relative humidity constant and calculte new saturation vapor pressure
+		//Hold relative humidity constant, calculte new saturation vapor pressure,
+		// and new saturation specific humidity to scale precipitation
 		//https://cran.r-project.org/web/packages/humidity/vignettes/humidity-measures.html
 		//es over ice calculation
 		//Ding et al., 2019 after Bolton, 1980
 		//ea37 = rh37*100*6.112.*exp((17.67*(t237-273.15))./(t237-29.65));
-		rhparam=eaparam/6.112/exp((17.67*(taparam-273.15))/(taparam-29.65));
-		eAir=fmax(rhparam*6.112*exp((17.67*(Ta-273.15))/(Ta-29.65)),0.0);
-
-		if (isprecipmap && (eaparam>0) && (eAir>0)){
-			P=prparam*eAir/eaparam;
-			C=C*eAir/eaparam;
+		//ea37s (hPa) = 6.1121*exp(22.587*(t-273.15)/(t-273.15+273.86)); (with respect to ice)
+		IssmDouble esparam, es, qsparam, qs;
+		esparam=6.112*exp((17.67*(taparam-273.15))/(taparam-29.65));
+		es=6.112*exp((17.67*(Ta-273.15))/(Ta-29.65));
+		rhparam=eaparam/esparam;
+		eAir=fmax(rhparam*es,0.0);
+		qs=fmax(0.622*es/(pAir/100 - 0.378*es),0);
+		qsparam=fmax(0.622*esparam/(pparam/100 - 0.378*esparam),0);
+
+		if ((isprecipmap) && (qsparam>0)){ 
+			P=fmax(prparam*qs/qsparam,0.0);
+			C=fmax(C*qs/qsparam,0.0);
 		}
 		else P=prparam;
 

src/c/classes/Loads/Friction.cpp

@@ -641,7 +641,7 @@ void Friction::GetAlpha2Josh(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
 
 	/*Intermediaries: */
 	IssmDouble  T,Tpmp,deltaT,deltaTref,pressure,diff,drag_coefficient;
-	IssmDouble  alpha2,time,gamma,ref,alp_new,alphascaled;
+	IssmDouble  alpha2,time,gamma,ref,alp_new,alphascaled,max_coefficient;
 	const IssmDouble yts = 365*24*3600.;
 
 	/*Get viscous part*/
@@ -655,6 +655,8 @@ void Friction::GetAlpha2Josh(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
 	 *    element->GetInputValue(&Tpdd,gauss,TemperaturePDDEnum);
 	 *       */
 
+	element->parameters->FindParam(&max_coefficient,FrictionMaxCoefficientEnum);
+
 	/*Compute delta T*/
 	element->GetInputValue(&T,gauss,TemperatureEnum);
 	element->GetInputValue(&pressure,gauss,PressureEnum);
@@ -669,7 +671,7 @@ void Friction::GetAlpha2Josh(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
 	alp_new = ref/exp(deltaT/gamma);
 
 	alphascaled = sqrt(alp_new)*drag_coefficient;
-	if (alphascaled > 300) alp_new = (300/drag_coefficient)*(300/drag_coefficient);
+	if (alphascaled > max_coefficient) alp_new = (max_coefficient/drag_coefficient)*(max_coefficient/drag_coefficient);
 
 	alp_new=alp_new*alpha2;
 
@@ -1419,6 +1421,7 @@ void FrictionUpdateParameters(Parameters* parameters,IoModel* iomodel){/*{{{*/
 		case 9:
 			parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
 			parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
+			parameters->AddObject(iomodel->CopyConstantObject("md.friction.coefficient_max",FrictionMaxCoefficientEnum));
 			parameters->AddObject(new IntParam(FrictionCouplingEnum,2));/*comment this line to use effective pressure from Beuler and Pelt (2015)*/
 			break;
 		case 10:

src/c/shared/Elements/PddSurfaceMassBalanceSicopolis.cpp

@@ -11,7 +11,7 @@ IssmDouble PddSurfaceMassBalanceSicopolis(IssmDouble* monthlytemperatures, IssmD
 				 IssmDouble* melt, IssmDouble* accu, IssmDouble* melt_star, IssmDouble* t_ampl, IssmDouble* p_ampl,
 				 IssmDouble yts, IssmDouble s, IssmDouble desfac,
 				 IssmDouble s0t, IssmDouble s0p, IssmDouble rlaps,
-				 IssmDouble rho_water,IssmDouble rho_ice){
+				 IssmDouble rho_water,IssmDouble rho_ice,IssmDouble pdd_fac_ice,IssmDouble pdd_fac_snow){
 
   int			imonth;				// month counter
   IssmDouble B;					// output: surface mass balance (m/a IE), melt+accumulation
@@ -27,16 +27,13 @@ IssmDouble PddSurfaceMassBalanceSicopolis(IssmDouble* monthlytemperatures, IssmD
   IssmDouble pdd;					// pdd factor (a * degC)
   IssmDouble tstar;				// monthly temp. after lapse rate correction (degC)
   IssmDouble precip_star;		// monthly precip after correction (m/a IE)
-  IssmDouble beta1 = 2.73;		// 3 mm IE/(d*deg C),  ablation factor for snow per positive degree day.
-  IssmDouble beta2 = 7.28;		// 8 mm IE/(d*deg C),  ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002).
   IssmDouble Pmax = 0.6;
   IssmDouble inv_twelve=1./12.; 
 
   sconv=(rho_water/rho_ice);		//rhow_rain/rhoi
 
-  /* FIXME betas shoud be user input */
-  beta1=beta1*(0.001*365)*sconv; // (mm WE)/(d*deg C) --> (m IE)/(a*deg C)
-  beta2=beta2*(0.001*365)*sconv; // (mm WE)/(d*deg C) --> (m IE)/(a*deg C)
+  pdd_fac_snow=pdd_fac_snow*(0.001*365)*sconv; // (mm WE)/(d*deg C) --> (m IE)/(a*deg C)
+  pdd_fac_ice=pdd_fac_ice*(0.001*365)*sconv; // (mm WE)/(d*deg C) --> (m IE)/(a*deg C)
 
   /* initalize fields */
   precip=0.0;
@@ -49,8 +46,7 @@ IssmDouble PddSurfaceMassBalanceSicopolis(IssmDouble* monthlytemperatures, IssmD
     /********* Surface temperature correction *******/    
     st=(s-s0t)/1000.;
 
-    // FIXME rlaps ??
-	 rlaps=-6.309e-03+(-5.426e-03-(-6.309e-03))*sin((imonth+1-4)*PI/6.0)*1000.0;
+    /******** Monhtly temperature correction *******/
     monthlytemperatures[imonth]=monthlytemperatures[imonth]-rlaps*st;//*max(st,1e-3);
     tstar=monthlytemperatures[imonth]+t_ampl[0];
 
@@ -109,20 +105,20 @@ IssmDouble PddSurfaceMassBalanceSicopolis(IssmDouble* monthlytemperatures, IssmD
   if(rainfall<0.0) rainfall=0.0;   // negative values
 
   if(rainfall<=(Pmax*snowfall)){
-	  if((rainfall+beta1*pdd)<=(Pmax*snowfall)) {
-		  smelt_star = rainfall+beta1*pdd;
+	  if((rainfall+pdd_fac_snow*pdd)<=(Pmax*snowfall)) {
+		  smelt_star = rainfall+pdd_fac_snow*pdd;
 		  smelt      = 0.0;
 		  runoff     = smelt;
 	  }
 	  else{
 		  smelt_star = Pmax*snowfall;
-		  smelt      = beta2*(pdd-(smelt_star-rainfall)/beta1);
+		  smelt      = pdd_fac_ice*(pdd-(smelt_star-rainfall)/pdd_fac_snow);
 		  runoff     = smelt;
 	  }
   } 
   else{
 	  smelt_star = Pmax*snowfall;
-	  smelt      = beta2*pdd;
+	  smelt      = pdd_fac_ice*pdd;
 	  runoff     = smelt+rainfall-Pmax*snowfall;
   }
 

src/c/shared/Elements/elements.h

@@ -28,7 +28,7 @@ IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures,  IssmDouble* m
 IssmDouble PddSurfaceMassBalanceSicopolis(IssmDouble* monthlytemperatures,  IssmDouble* monthlyprec,
 				 IssmDouble* melt, IssmDouble* accu, IssmDouble* melt_star, IssmDouble* t_ampl, IssmDouble* p_ampl,
 				 IssmDouble yts, IssmDouble s, IssmDouble desfac,IssmDouble s0t,
-				 IssmDouble s0p, IssmDouble rlaps, IssmDouble rho_water, IssmDouble rho_ice);
+				 IssmDouble s0p, IssmDouble rlaps, IssmDouble rho_water, IssmDouble rho_ice, IssmDouble pdd_fac_ice, IssmDouble pdd_fac_snow);
 void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime,
 					     IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime,
 					     IssmDouble* PrecipitationsPresentday,

src/c/shared/Enum/Enum.vim

@@ -221,6 +221,7 @@ syn keyword cConstant FrictionFEnum
 syn keyword cConstant FrictionGammaEnum
 syn keyword cConstant FrictionLawEnum
 syn keyword cConstant FrictionLinearizeEnum
+syn keyword cConstant FrictionMaxCoefficientEnum
 syn keyword cConstant FrictionPseudoplasticityExponentEnum
 syn keyword cConstant FrictionU0Enum
 syn keyword cConstant FrictionThresholdSpeedEnum
@@ -433,6 +434,8 @@ syn keyword cConstant OutputBufferSizePointerEnum
 syn keyword cConstant OutputFileNameEnum
 syn keyword cConstant OutputFilePointerEnum
 syn keyword cConstant OutputdefinitionEnum
+syn keyword cConstant PddfacIceEnum
+syn keyword cConstant PddfacSnowEnum
 syn keyword cConstant QmuErrNameEnum
 syn keyword cConstant QmuInNameEnum
 syn keyword cConstant QmuIsdakotaEnum

src/c/shared/Enum/EnumDefinitions.h

@@ -215,6 +215,7 @@ enum definitions{
 	FrictionGammaEnum,
 	FrictionLawEnum,
 	FrictionLinearizeEnum,
+	FrictionMaxCoefficientEnum,
 	FrictionPseudoplasticityExponentEnum,
 	FrictionU0Enum,
 	FrictionThresholdSpeedEnum,
@@ -427,6 +428,8 @@ enum definitions{
 	OutputFileNameEnum,
 	OutputFilePointerEnum,
 	OutputdefinitionEnum,
+	PddfacIceEnum,
+	PddfacSnowEnum,
 	QmuErrNameEnum,
 	QmuInNameEnum,
 	QmuIsdakotaEnum,

src/c/shared/Enum/EnumToStringx.cpp

@@ -223,6 +223,7 @@ const char* EnumToStringx(int en){
 		case FrictionGammaEnum : return "FrictionGamma";
 		case FrictionLawEnum : return "FrictionLaw";
 		case FrictionLinearizeEnum : return "FrictionLinearize";
+		case FrictionMaxCoefficientEnum : return "FrictionMaxCoefficient";
 		case FrictionPseudoplasticityExponentEnum : return "FrictionPseudoplasticityExponent";
 		case FrictionU0Enum : return "FrictionU0";
 		case FrictionThresholdSpeedEnum : return "FrictionThresholdSpeed";
@@ -435,6 +436,8 @@ const char* EnumToStringx(int en){
 		case OutputFileNameEnum : return "OutputFileName";
 		case OutputFilePointerEnum : return "OutputFilePointer";
 		case OutputdefinitionEnum : return "Outputdefinition";
+		case PddfacIceEnum : return "PddfacIce";
+		case PddfacSnowEnum : return "PddfacSnow";
 		case QmuErrNameEnum : return "QmuErrName";
 		case QmuInNameEnum : return "QmuInName";
 		case QmuIsdakotaEnum : return "QmuIsdakota";

src/c/shared/Enum/Enumjl.vim

@@ -214,6 +214,7 @@ syn keyword juliaConstC FrictionFEnum
 syn keyword juliaConstC FrictionGammaEnum
 syn keyword juliaConstC FrictionLawEnum
 syn keyword juliaConstC FrictionLinearizeEnum
+syn keyword juliaConstC FrictionMaxCoefficientEnum
 syn keyword juliaConstC FrictionPseudoplasticityExponentEnum
 syn keyword juliaConstC FrictionU0Enum
 syn keyword juliaConstC FrictionThresholdSpeedEnum
@@ -426,6 +427,8 @@ syn keyword juliaConstC OutputBufferSizePointerEnum
 syn keyword juliaConstC OutputFileNameEnum
 syn keyword juliaConstC OutputFilePointerEnum
 syn keyword juliaConstC OutputdefinitionEnum
+syn keyword juliaConstC PddfacIceEnum
+syn keyword juliaConstC PddfacSnowEnum
 syn keyword juliaConstC QmuErrNameEnum
 syn keyword juliaConstC QmuInNameEnum
 syn keyword juliaConstC QmuIsdakotaEnum