Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 23316) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 23317) @@ -175,4 +175,5 @@ ./shared/Elements/PrintArrays.cpp\ ./shared/Elements/PddSurfaceMassBalance.cpp\ + ./shared/Elements/PddSurfaceMassBalanceSicopolis.cpp\ ./shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp\ ./shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp\ Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 23316) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 23317) @@ -43,6 +43,7 @@ int smb_model; iomodel->FindConstant(&smb_model,"md.smb.model"); - if(smb_model==SMBpddEnum) isdynamic=true; - if(smb_model==SMBd18opddEnum) isdynamic=true; + if(smb_model==SMBpddEnum) isdynamic=true; + if(smb_model==SMBd18opddEnum) isdynamic=true; + if(smb_model==SMBpddSicopolisEnum) isdynamic=true; } @@ -1200,4 +1201,6 @@ _assert_((Hc+Ht)>0.); lambda = Hc/(Hc+Ht); + _assert_(lambda>=0.); + _assert_(lambda<=1.); kappa = kappa_c*kappa_t/(lambda*kappa_t+(1.-lambda)*kappa_c); // ==(lambda/kappa_c + (1.-lambda)/kappa_t)^-1 } Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23316) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23317) @@ -5,4 +5,7 @@ #include "../modules/modules.h" +// FIX +#include "./shared/io/Print/Print.h" + /*Model processing*/ void SmbAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ @@ -21,5 +24,5 @@ int smb_model; - bool isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled; + bool isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled,isfirnwarming; /*Update elements: */ @@ -88,4 +91,14 @@ } break; + case SMBpddSicopolisEnum: + iomodel->FetchDataToInput(elements,"md.smb.s0p",SmbS0pEnum); + iomodel->FetchDataToInput(elements,"md.smb.s0t",SmbS0tEnum); + iomodel->FindConstant(&isfirnwarming,"md.smb.isfirnwarming"); + iomodel->FetchDataToInput(elements,"md.smb.smb_corr",SmbSmbCorrEnum); + iomodel->FetchDataToInput(elements,"md.smb.precipitation",SmbPrecipitationEnum); + iomodel->FetchDataToInput(elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum); + iomodel->FetchDataToInput(elements,"md.smb.precipitation_anomaly",SmbPrecipitationsAnomalyEnum); + iomodel->FetchDataToInput(elements,"md.smb.temperature_anomaly",SmbTemperaturesAnomalyEnum); + break; case SMBd18opddEnum: iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled"); @@ -148,5 +161,5 @@ int numoutputs; char** requestedoutputs = NULL; - bool isdelta18o,ismungsm,isd18opd,issetpddfac,interp; + bool isdelta18o,ismungsm,isd18opd,issetpddfac,interp,isfirnwarming; int smb_model; IssmDouble *temp = NULL; @@ -215,4 +228,7 @@ } break; + case SMBpddSicopolisEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.smb.isfirnwarming",SmbIssetpddfacEnum)); + break; case SMBd18opddEnum: parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum)); @@ -288,4 +304,8 @@ PositiveDegreeDayx(femmodel); break; + case SMBpddSicopolisEnum: + if(VerboseSolution()) _printf0_(" call SICOPOLIS positive degree day module\n"); + PositiveDegreeDaySicopolisx(femmodel); + break; case SMBd18opddEnum: bool isd18opd; @@ -296,5 +316,5 @@ if(VerboseSolution()) _printf0_(" call positive degree day module\n"); PositiveDegreeDayx(femmodel); - } + } break; case SMBgradientsEnum: Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 23316) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 23317) @@ -33,6 +33,7 @@ int smb_model; iomodel->FindConstant(&smb_model,"md.smb.model"); - if(smb_model==SMBpddEnum) isdynamic=true; - if(smb_model==SMBd18opddEnum) isdynamic=true; + if(smb_model==SMBpddEnum) isdynamic=true; + if(smb_model==SMBd18opddEnum) isdynamic=true; + if(smb_model==SMBpddSicopolisEnum) isdynamic=true; } else{ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23316) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23317) @@ -2559,4 +2559,169 @@ xDelete(tmp); +} +/*}}}*/ +void Element::PositiveDegreeDaySicopolis(bool isfirnwarming){/*{{{*/ + + /* General FIXMEs: get Tmelting point, pddicefactor, pddsnowfactor, sigma from parameters/user input */ + + int numvertices = this->GetNumberOfVertices(); + + int i; + IssmDouble* smb=xNew(numvertices); // surface mass balance + IssmDouble* melt=xNew(numvertices); // melting comp. of surface mass balance + IssmDouble* accu=xNew(numvertices); // accuumulation comp. of surface mass balance + IssmDouble* melt_star=xNew(numvertices); + IssmDouble* monthlytemperatures=xNew(12*numvertices); + IssmDouble* monthlyprec=xNew(12*numvertices); + IssmDouble* yearlytemperatures=xNew(numvertices); memset(yearlytemperatures, 0., numvertices*sizeof(IssmDouble)); + IssmDouble* s=xNew(numvertices); // actual surface height + IssmDouble* s0p=xNew(numvertices); // reference elevation for precip. + IssmDouble* s0t=xNew(numvertices); // reference elevation for temperature + IssmDouble* smbcorr=xNew(numvertices); // surface mass balance correction; will be added after pdd call + IssmDouble* p_ampl=xNew(numvertices); // precip anomaly + IssmDouble* t_ampl=xNew(numvertices); // remperature anomaly + IssmDouble rho_water,rho_ice,desfac,rlaps; + IssmDouble inv_twelve=1./12.; //factor for monthly average + IssmDouble time,yts,time_yr; + + /*Get material parameters :*/ + rho_water=this->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + rho_ice=this->matpar->GetMaterialParameter(MaterialsRhoIceEnum); + + /*Get parameters for height corrections*/ + desfac=this->matpar->GetMaterialParameter(SmbDesfacEnum); + rlaps=this->matpar->GetMaterialParameter(SmbRlapsEnum); + + /*Recover monthly temperatures and precipitation*/ + Input* input=this->inputs->GetInput(SmbMonthlytemperaturesEnum); _assert_(input); + Input* input2=this->inputs->GetInput(SmbPrecipitationEnum); _assert_(input2); + /*Recover smb correction term */ + Input* input3=this->inputs->GetInput(SmbSmbCorrEnum); _assert_(input3); + + /* Get time */ + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + time_yr=floor(time/yts)*yts; + + /* Set parameters for finrnwarming */ + IssmDouble MU_0 = 9.7155; //Firn-warming correction, in (d*deg C)/(mm WE) + IssmDouble mu = MU_0*(1000.0*86400.0)*(rho_ice/rho_water); // (d*deg C)/(mm WE) --> (s*deg C)/(m IE) + + /*loop over vertices: */ + Gauss* gauss=this->NewGauss(); + for(int month=0;month<12;month++){ + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&monthlytemperatures[iv*12+month],gauss,(month+1)/12.*yts); + monthlytemperatures[iv*12+month]=monthlytemperatures[iv*12+month]-273.15; // conversion from Kelvin to celcius for PDD module + input2->GetInputValue(&monthlyprec[iv*12+month],gauss,(month+1)/12.*yts); + monthlyprec[iv*12+month]=monthlyprec[iv*12+month]*yts; + } + } + + /*Recover info at the vertices: */ + GetInputListOnVertices(&s[0],SurfaceEnum); + GetInputListOnVertices(&s0p[0],SmbS0pEnum); + GetInputListOnVertices(&s0t[0],SmbS0tEnum); + GetInputListOnVertices(&smbcorr[0],SmbSmbCorrEnum); + GetInputListOnVertices(&t_ampl[0],SmbTemperaturesAnomalyEnum); + GetInputListOnVertices(&p_ampl[0],SmbPrecipitationsAnomalyEnum); + + /*measure the surface mass balance*/ + for (int iv = 0; iv=melt[iv]){ + yearlytemperatures[iv]= yearlytemperatures[iv]+mu*(melt_star[iv]-melt[iv]); + } + else{ + yearlytemperatures[iv]= yearlytemperatures[iv]; + } + } + if (yearlytemperatures[iv]>273.15) yearlytemperatures[iv]=273.15; + } + + switch(this->ObjectEnum()){ + case TriaEnum: + // this->inputs->AddInput(new TriaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); + this->inputs->AddInput(new TriaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum)); + this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&smb[0],P1Enum)); + this->inputs->AddInput(new TriaInput(SmbAccumulationEnum,&accu[0],P1Enum)); + this->inputs->AddInput(new TriaInput(SmbMeltEnum,&melt[0],P1Enum)); + break; + case PentaEnum: + bool isthermal; + this->parameters->FindParam(&isthermal,TransientIsthermalEnum); + if(isthermal){ + bool isenthalpy; + this->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); + if(IsOnSurface()){ + /*Here, we want to change the BC of the thermal model, keep + * the temperatures as they are for the base of the penta and + * use yearlytemperatures for the top*/ + GetInputListOnVertices(&s[0],TemperatureEnum); + yearlytemperatures[0] = s[0]; + yearlytemperatures[1] = s[1]; + yearlytemperatures[2] = s[2]; + this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); + if(isenthalpy){ + /*Convert that to enthalpy for the enthalpy model*/ + IssmDouble enthalpy[6]; + GetInputListOnVertices(&enthalpy[0],EnthalpyEnum); + ThermalToEnthalpy(&enthalpy[3],yearlytemperatures[3],0.,0.); + ThermalToEnthalpy(&enthalpy[4],yearlytemperatures[4],0.,0.); + ThermalToEnthalpy(&enthalpy[5],yearlytemperatures[5],0.,0.); + this->inputs->AddInput(new PentaInput(EnthalpyEnum,&enthalpy[0],P1Enum)); + } + } + } + this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&smb[0],P1Enum)); + this->inputs->AddInput(new PentaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum)); + this->inputs->AddInput(new PentaInput(SmbAccumulationEnum,&accu[0],P1Enum)); + this->inputs->AddInput(new PentaInput(SmbMeltEnum,&melt[0],P1Enum)); + this->InputExtrude(TemperaturePDDEnum,-1); + this->InputExtrude(SmbMassBalanceEnum,-1); + this->InputExtrude(SmbAccumulationEnum,-1); + this->InputExtrude(SmbMeltEnum,-1); + break; + case TetraEnum: + if(IsOnSurface()){ + GetInputListOnVertices(&s[0],TemperatureEnum); + yearlytemperatures[0] = s[0]; + yearlytemperatures[1] = s[1]; + yearlytemperatures[2] = s[2]; + this->inputs->AddInput(new TetraInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); + } + this->inputs->AddInput(new TetraInput(SmbMassBalanceEnum,&smb[0],P1Enum)); + this->inputs->AddInput(new TetraInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum)); + this->InputExtrude(TemperaturePDDEnum,-1); + this->InputExtrude(SmbMassBalanceEnum,-1); + break; + default: _error_("Not implemented yet"); + } + + /*clean-up*/ + delete gauss; + xDelete(monthlytemperatures); + xDelete(monthlyprec); + xDelete(smb); + xDelete(melt); + xDelete(accu); + xDelete(yearlytemperatures); + xDelete(s); + xDelete(s0t); + xDelete(s0p); + xDelete(t_ampl); + xDelete(p_ampl); + xDelete(smbcorr); + xDelete(melt_star); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 23316) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 23317) @@ -145,4 +145,5 @@ ElementVector* NewElementVector(int approximation_enum=NoneApproximationEnum); void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac); + void PositiveDegreeDaySicopolis(bool isfirnwarming); IssmDouble PureIceEnthalpy(IssmDouble pressure); void ResultInterpolation(int* pinterpolation,int*nodesperelement,int* parray_size, int output_enum); Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 23316) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 23317) @@ -106,4 +106,8 @@ iomodel->FindConstant(&this->rlaps,"md.smb.rlaps"); iomodel->FindConstant(&this->rlapslgm,"md.smb.rlapslgm"); + break; + case SMBpddSicopolisEnum: + iomodel->FindConstant(&this->desfac,"md.smb.desfac"); + iomodel->FindConstant(&this->rlaps,"md.smb.rlaps"); break; case SMBd18opddEnum: @@ -210,4 +214,8 @@ iomodel->FindConstant(&this->rlaps,"md.smb.rlaps"); iomodel->FindConstant(&this->rlapslgm,"md.smb.rlapslgm"); + break; + case SMBpddSicopolisEnum: + iomodel->FindConstant(&this->desfac,"md.smb.desfac"); + iomodel->FindConstant(&this->rlaps,"md.smb.rlaps"); break; case SMBd18opddEnum: Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 23316) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 23317) @@ -371,10 +371,12 @@ femmodel->parameters->FindParam(&smb_model,SmbEnum); if(isenthalpy){ - if(smb_model==SMBpddEnum) ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum); - if(smb_model==SMBd18opddEnum) ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum); + if(smb_model==SMBpddEnum) ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum); + if(smb_model==SMBd18opddEnum) ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum); + if(smb_model==SMBpddSicopolisEnum) ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum); } else{ - if(smb_model==SMBpddEnum) ResetBoundaryConditions(femmodel,ThermalAnalysisEnum); - if(smb_model==SMBd18opddEnum) ResetBoundaryConditions(femmodel,ThermalAnalysisEnum); + if(smb_model==SMBpddEnum) ResetBoundaryConditions(femmodel,ThermalAnalysisEnum); + if(smb_model==SMBd18opddEnum) ResetBoundaryConditions(femmodel,ThermalAnalysisEnum); + if(smb_model==SMBpddSicopolisEnum) ResetBoundaryConditions(femmodel,ThermalAnalysisEnum); } } Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 23316) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 23317) @@ -262,4 +262,15 @@ xDelete(pds); }/*}}}*/ +void PositiveDegreeDaySicopolisx(FemModel* femmodel){/*{{{*/ + + bool isfirnwarming; + femmodel->parameters->FindParam(&isfirnwarming,SmbIssetpddfacEnum); + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->PositiveDegreeDaySicopolis(isfirnwarming); + } + +}/*}}}*/ void SmbHenningx(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 23316) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 23317) @@ -16,4 +16,5 @@ void Delta18opdParameterizationx(FemModel* femmodel); void PositiveDegreeDayx(FemModel* femmodel); +void PositiveDegreeDaySicopolisx(FemModel* femmodel); void SmbHenningx(FemModel* femmodel); void SmbComponentsx(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalanceSicopolis.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalanceSicopolis.cpp (revision 23317) +++ /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalanceSicopolis.cpp (revision 23317) @@ -0,0 +1,131 @@ +/* file: PddSurfaceMassBlanceSicopolis.cpp + Calculating the surface mass balance using the adapted PDD routine from SICOPOLIS. + */ + +#include "./elements.h" +#include "../Numerics/numerics.h" + +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){ + + int imonth; // month counter + IssmDouble B; // output: surface mass balance (m/a IE), melt+accumulation + IssmDouble frac_solid, snowfall, rainfall, runoff; + IssmDouble saccu; // yearly surface accumulation (m/a IE) + IssmDouble smelt; // yearly melt (m/a IE) + IssmDouble smelt_star; // yearly ... + IssmDouble precip; // total precipitation during 1 year + IssmDouble sconv; //rhow_rain/rhoi / 12 months + IssmDouble st; // elevation between altitude of the temp record and current altitude + IssmDouble sp; // elevation between altitude of the prec record and current altitude + IssmDouble q; // q is desert/elev. fact + 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) + + /* initalize fields */ + precip=0.0; + tstar=0.0; + snowfall=0.0; + pdd=0.0; + /* seasonal loop */ + for(imonth=0;imonth<12;imonth++){ + + /********* 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; + monthlytemperatures[imonth]=monthlytemperatures[imonth]-rlaps*st;//*max(st,1e-3); + tstar=monthlytemperatures[imonth]+t_ampl[0]; + + /********* Precipitation correction *************/ + /* Ref: Vizcaino et al 2010; DOI 10.1007/s00382-009-0591-y */ + if(s0p<2000.0) + q=exp(desfac*(max(s,2000.0)-2000.0)); + else + q=exp(desfac*(max(s,2000.0)-s0p)); + + precip_star=q*monthlyprec[imonth]*sconv*p_ampl[0]*yts; // convert precip from m/s -> m/a + precip=precip+precip_star*inv_twelve; + + /********* compute PDD **************************/ + /* Ref: Calov & Greve 2005 Journal of Glaciology, Vol. 51, No. 172, 2005, Correspondence */ + IssmDouble s_stat=5.0; + IssmDouble inv_sqrt2pi =1.0/sqrt(2.0*PI); + IssmDouble inv_s_stat =1.0/s_stat; + IssmDouble inv_sqrt2 =1.0/sqrt(2.0); + pdd=pdd+(s_stat*inv_sqrt2pi*exp(-0.5*pow(tstar*inv_s_stat,2)) + +0.5*tstar*erfc(-tstar*inv_s_stat*inv_sqrt2))*inv_twelve; + + /*Partition of precip in solid and liquid parts, Bales et al. (2009) */ + IssmDouble temp_rain=7.2; // Threshold monthly mean temperature for + // precipitation = 101% rain, in deg C + IssmDouble temp_snow=-11.6; // Threshold monthly mean temperature for + // precipitation = 100% snow, in deg C + + IssmDouble coeff1=5.4714e-01; // Coefficients + IssmDouble coeff2=-9.1603e-02; // of + IssmDouble coeff3=-3.314e-03; // the + IssmDouble coeff4= 4.66e-04; // fifth-order + IssmDouble coeff5=3.8e-05; // polynomial + IssmDouble coeff6=6.0e-07; // fit + + if(tstar>=temp_rain) + frac_solid = 0.0; + else if(tstar<=temp_snow) + frac_solid = 1.0; + else{ + frac_solid=coeff1+tstar*(coeff2 + +tstar*(coeff3+tstar*(coeff4+tstar*(coeff5+tstar*coeff6)))); + } + + snowfall=snowfall+precip_star*frac_solid*inv_twelve; + } + /* end of seasonal loop */ + + rainfall=precip-snowfall; + if(snowfall<0.0) snowfall=0.0; // correction of + if(rainfall<0.0) rainfall=0.0; // negative values + + if(rainfall<=(Pmax*snowfall)){ + if((rainfall+beta1*pdd)<=(Pmax*snowfall)) { + smelt_star = rainfall+beta1*pdd; + smelt = 0.0; + runoff = smelt; + } + else{ + smelt_star = Pmax*snowfall; + smelt = beta2*(pdd-(smelt_star-rainfall)/beta1); + runoff = smelt; + } + } + else{ + smelt_star = Pmax*snowfall; + smelt = beta2*pdd; + runoff = smelt+rainfall-Pmax*snowfall; + } + + saccu = precip; + + /* asign output*/ + melt[0]=runoff/yts; + accu[0]=saccu/yts; + melt_star[0]=smelt_star/yts; + B=(saccu-runoff)/yts; + + return B; +} Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 23316) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 23317) @@ -24,4 +24,8 @@ IssmDouble TdiffTime,IssmDouble sealevTime,IssmDouble pddsnowfac,IssmDouble pddicefac, IssmDouble rho_water, IssmDouble rho_ice); +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); void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime, IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23316) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23317) @@ -1157,4 +1157,9 @@ XYEnum, XYZEnum, + SMBpddSicopolisEnum, + SmbIsfirnwarmingEnum, + SmbSmbCorrEnum, + SmbTemperaturesAnomalyEnum, + SmbPrecipitationsAnomalyEnum, /*}}}*/ /*Unused?{{{*/ Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 23316) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 23317) @@ -173,4 +173,5 @@ case 8: return SMBgembEnum; case 9: return SMBgradientselaEnum; + case 10: return SMBpddSicopolisEnum; default: _error_("Marshalled SMB code \""<=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','smb.s0t','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','smb.rlaps','>=',0,'numel',1); + md = checkfield(md,'fieldname','smb.monthlytemperatures','timeseries',1,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices+1 12]); + md = checkfield(md,'fieldname','smb.precipitation','timeseries',1,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices+1 12]); + end + md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1); + + end % }}} + function disp(self) % {{{ + disp(sprintf(' surface forcings parameters:')); + + disp(sprintf('\n SICOPOLIS PDD scheme (Calov & Greve, 2005) :')); + fielddisplay(self,'monthlytemperatures','monthly surface temperatures [K]'); + fielddisplay(self,'precipitation','monthly surface precipitation [m/yr water eq]'); + fielddisplay(self,'temperature_anomaly','anomaly to monthly reference temperature (additive; [K])'); + fielddisplay(self,'precipitation_anomaly','anomaly to monthly precipitation (multiplicative, e.g. q=q0*exp(0.070458*DeltaT) after Huybrechts (2002)); [no unit])'); + fielddisplay(self,'smb_corr','correction of smb after PDD call [m/a]'); + fielddisplay(self,'s0p','should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]'); + fielddisplay(self,'s0t','should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]'); + fielddisplay(self,'rlaps','present day lapse rate (default is 7.4 degree/km)'); + fielddisplay(self,'desfac','desertification elevation factor (default is -log(2.0)/1000)'); + fielddisplay(self,'isfirnwarming','is firnwarming (Reeh 1991) activated (0 or 1, default is 1)'); + fielddisplay(self,'requested_outputs','additional outputs requested (TemperaturePDD, SmbAccumulation, SmbMelt)'); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + yts=md.constants.yts; + + WriteData(fid,prefix,'name','md.smb.model','data',10,'format','Integer'); + + WriteData(fid,prefix,'object',self,'class','smb','fieldname','isfirnwarming','format','Boolean'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfac','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0p','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0t','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','rlaps','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','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'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','smb_corr','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + + %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/test/NightlyRun/test245.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test245.m (revision 23317) +++ /issm/trunk-jpl/test/NightlyRun/test245.m (revision 23317) @@ -0,0 +1,47 @@ +%Test Name: SquareShelfTranIspddSicopolisSSA2d +md=triangle(model(),'../Exp/Square.exp',150000.); +md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareShelf.par'); + +% Use of SMBpddSicopolis +md.smb = SMBpddSicopolis(); +% initalize pdd fields +md.smb=initialize(md.smb,md); +md.smb.s0p=md.geometry.surface; +md.smb.s0t=md.geometry.surface; + +% +md.smb.monthlytemperatures=[]; md.smb.precipitation=[]; md.smb.precipitation=[]; +temp_ma_present=-10*ones(md.mesh.numberofvertices,1)-md.smb.rlaps*md.geometry.surface/1000; +temp_mj_present=10*ones(md.mesh.numberofvertices,1)-md.smb.rlaps*md.geometry.surface/1000; +precipitation=5*ones(md.mesh.numberofvertices,1); +for imonth=0:11 + md.smb.monthlytemperatures(1:md.mesh.numberofvertices,imonth+1)=md.materials.meltingpoint+temp_ma_present+(temp_mj_present-temp_ma_present)*sin(double(imonth+1-4)*pi/6.0); + md.smb.monthlytemperatures(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); + md.smb.precipitation(1:md.mesh.numberofvertices,imonth+1)=precipitation; + md.smb.precipitation(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); +end + +% time steps and resolution +md.timestepping.time_step=1; +md.settings.output_frequency=1; +md.timestepping.final_time=2; + +md.transient.issmb=1; +md.transient.ismasstransport=1; +md.transient.isstressbalance=0; +md.transient.isthermal=0; + +md.transient.requested_outputs={'default','TemperaturePDD'}; +md.cluster=generic('name',oshostname(),'np',1); % 3 for the cluster +md=solve(md,'Transient'); + +%Fields and tolerances to track changes +field_names ={'TemperaturePDD1','SmbMassBalance1','TemperaturePDD2','SmbMassBalance2'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13}; +field_values={... + (md.results.TransientSolution(1).TemperaturePDD),... + (md.results.TransientSolution(1).SmbMassBalance),... + (md.results.TransientSolution(2).TemperaturePDD),... + (md.results.TransientSolution(2).SmbMassBalance),... + };