diff --git a/src/c/analyses/SmbAnalysis.cpp b/src/c/analyses/SmbAnalysis.cpp index 3c17be2be..02052f30d 100644 --- a/src/c/analyses/SmbAnalysis.cpp +++ b/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)); diff --git a/src/c/classes/Elements/Element.cpp b/src/c/classes/Elements/Element.cpp index 1f902fc7d..2b7fcccd4 100644 --- a/src/c/classes/Elements/Element.cpp +++ b/src/c/classes/Elements/Element.cpp @@ -3695,6 +3695,7 @@ void Element::PositiveDegreeDaySicopolis(bool isfirnwarming){/*{{{*/ IssmDouble* p_ampl=xNew(NUM_VERTICES); // precip anomaly IssmDouble* t_ampl=xNew(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; iv0) && (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; diff --git a/src/c/classes/Loads/Friction.cpp b/src/c/classes/Loads/Friction.cpp index 6a0a3cc72..76f8f8587 100644 --- a/src/c/classes/Loads/Friction.cpp +++ b/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: diff --git a/src/c/shared/Elements/PddSurfaceMassBalanceSicopolis.cpp b/src/c/shared/Elements/PddSurfaceMassBalanceSicopolis.cpp index a94f72f68..380f479ef 100644 --- a/src/c/shared/Elements/PddSurfaceMassBalanceSicopolis.cpp +++ b/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; } diff --git a/src/c/shared/Elements/elements.h b/src/c/shared/Elements/elements.h index 5f66ef92a..f11a8677b 100644 --- a/src/c/shared/Elements/elements.h +++ b/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, diff --git a/src/c/shared/Enum/Enum.vim b/src/c/shared/Enum/Enum.vim index 0c6266e60..bb95a37d1 100644 --- a/src/c/shared/Enum/Enum.vim +++ b/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 diff --git a/src/c/shared/Enum/EnumDefinitions.h b/src/c/shared/Enum/EnumDefinitions.h index 70557519f..f941e36b5 100644 --- a/src/c/shared/Enum/EnumDefinitions.h +++ b/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, diff --git a/src/c/shared/Enum/EnumToStringx.cpp b/src/c/shared/Enum/EnumToStringx.cpp index 2d4d382b0..14e78bb32 100644 --- a/src/c/shared/Enum/EnumToStringx.cpp +++ b/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"; diff --git a/src/c/shared/Enum/Enumjl.vim b/src/c/shared/Enum/Enumjl.vim index 14c41d125..ce3c81591 100644 --- a/src/c/shared/Enum/Enumjl.vim +++ b/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 diff --git a/src/c/shared/Enum/StringToEnumx.cpp b/src/c/shared/Enum/StringToEnumx.cpp index f31d143ae..b75f68d5c 100644 --- a/src/c/shared/Enum/StringToEnumx.cpp +++ b/src/c/shared/Enum/StringToEnumx.cpp @@ -226,6 +226,7 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"FrictionGamma")==0) return FrictionGammaEnum; else if (strcmp(name,"FrictionLaw")==0) return FrictionLawEnum; else if (strcmp(name,"FrictionLinearize")==0) return FrictionLinearizeEnum; + else if (strcmp(name,"FrictionMaxCoefficient")==0) return FrictionMaxCoefficientEnum; else if (strcmp(name,"FrictionPseudoplasticityExponent")==0) return FrictionPseudoplasticityExponentEnum; else if (strcmp(name,"FrictionU0")==0) return FrictionU0Enum; else if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum; @@ -258,11 +259,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"FrontalForcingsSdNumberofParams")==0) return FrontalForcingsSdNumberofParamsEnum; else if (strcmp(name,"FrontalForcingsSdpolyparams")==0) return FrontalForcingsSdpolyparamsEnum; else if (strcmp(name,"GrdModel")==0) return GrdModelEnum; - else if (strcmp(name,"GroundinglineFrictionInterpolation")==0) return GroundinglineFrictionInterpolationEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"GroundinglineMeltInterpolation")==0) return GroundinglineMeltInterpolationEnum; + if (strcmp(name,"GroundinglineFrictionInterpolation")==0) return GroundinglineFrictionInterpolationEnum; + else if (strcmp(name,"GroundinglineMeltInterpolation")==0) return GroundinglineMeltInterpolationEnum; else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum; else if (strcmp(name,"GroundinglineNumRequestedOutputs")==0) return GroundinglineNumRequestedOutputsEnum; else if (strcmp(name,"GroundinglineRequestedOutputs")==0) return GroundinglineRequestedOutputsEnum; @@ -381,11 +382,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"LoveIntegrationScheme")==0) return LoveIntegrationSchemeEnum; else if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum; else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum; - else if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"LoveNTemporalIterations")==0) return LoveNTemporalIterationsEnum; + if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum; + else if (strcmp(name,"LoveNTemporalIterations")==0) return LoveNTemporalIterationsEnum; else if (strcmp(name,"LoveNYiEquations")==0) return LoveNYiEquationsEnum; else if (strcmp(name,"LoveR0")==0) return LoveR0Enum; else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; @@ -444,6 +445,8 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum; else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum; + else if (strcmp(name,"PddfacIce")==0) return PddfacIceEnum; + else if (strcmp(name,"PddfacSnow")==0) return PddfacSnowEnum; else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum; else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum; else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum; @@ -502,13 +505,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"LovePolarMotionTransferFunctionColinear")==0) return LovePolarMotionTransferFunctionColinearEnum; else if (strcmp(name,"LovePolarMotionTransferFunctionOrthogonal")==0) return LovePolarMotionTransferFunctionOrthogonalEnum; else if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum; - else if (strcmp(name,"TidalLoveK")==0) return TidalLoveKEnum; - else if (strcmp(name,"TidalLoveL")==0) return TidalLoveLEnum; - else if (strcmp(name,"TidalLoveK2Secular")==0) return TidalLoveK2SecularEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"LoadLoveH")==0) return LoadLoveHEnum; + if (strcmp(name,"TidalLoveK")==0) return TidalLoveKEnum; + else if (strcmp(name,"TidalLoveL")==0) return TidalLoveLEnum; + else if (strcmp(name,"TidalLoveK2Secular")==0) return TidalLoveK2SecularEnum; + else if (strcmp(name,"LoadLoveH")==0) return LoadLoveHEnum; else if (strcmp(name,"LoadLoveK")==0) return LoadLoveKEnum; else if (strcmp(name,"LoadLoveL")==0) return LoadLoveLEnum; else if (strcmp(name,"LoveTimeFreq")==0) return LoveTimeFreqEnum; @@ -625,13 +628,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"SmbIsmappedforcing")==0) return SmbIsmappedforcingEnum; else if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum; else if (strcmp(name,"SmbIsmungsm")==0) return SmbIsmungsmEnum; - else if (strcmp(name,"SmbIsprecipforcingremapped")==0) return SmbIsprecipforcingremappedEnum; - else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum; - else if (strcmp(name,"SmbIssetpddfac")==0) return SmbIssetpddfacEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum; + if (strcmp(name,"SmbIsprecipforcingremapped")==0) return SmbIsprecipforcingremappedEnum; + else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum; + else if (strcmp(name,"SmbIssetpddfac")==0) return SmbIssetpddfacEnum; + else if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum; else if (strcmp(name,"SmbIstemperaturescaled")==0) return SmbIstemperaturescaledEnum; else if (strcmp(name,"SmbIsthermal")==0) return SmbIsthermalEnum; else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum; @@ -748,13 +751,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; - else if (strcmp(name,"TransientIsmmemasstransport")==0) return TransientIsmmemasstransportEnum; - else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; - else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; + if (strcmp(name,"TransientIsmmemasstransport")==0) return TransientIsmmemasstransportEnum; + else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; + else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; + else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; else if (strcmp(name,"TransientIssampling")==0) return TransientIssamplingEnum; else if (strcmp(name,"TransientIsslc")==0) return TransientIsslcEnum; else if (strcmp(name,"TransientIssmb")==0) return TransientIssmbEnum; @@ -871,13 +874,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Dsl")==0) return DslEnum; else if (strcmp(name,"DeltaStr")==0) return DeltaStrEnum; else if (strcmp(name,"StrOld")==0) return StrOldEnum; - else if (strcmp(name,"Str")==0) return StrEnum; - else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; - else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; + if (strcmp(name,"Str")==0) return StrEnum; + else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; + else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; + else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; @@ -994,13 +997,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"DebrisMaskNodeActivation")==0) return DebrisMaskNodeActivationEnum; else if (strcmp(name,"Ice")==0) return IceEnum; else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; - else if (strcmp(name,"Input")==0) return InputEnum; - else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; - else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum; + if (strcmp(name,"Input")==0) return InputEnum; + else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; + else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; + else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum; else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; @@ -1117,13 +1120,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"SealevelchangeAzimuthIndexOcean")==0) return SealevelchangeAzimuthIndexOceanEnum; else if (strcmp(name,"SealevelchangeAzimuthIndexIce")==0) return SealevelchangeAzimuthIndexIceEnum; else if (strcmp(name,"SealevelchangeAzimuthIndexHydro")==0) return SealevelchangeAzimuthIndexHydroEnum; - else if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum; - else if (strcmp(name,"SealevelchangeViscousSG")==0) return SealevelchangeViscousSGEnum; - else if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum; + if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum; + else if (strcmp(name,"SealevelchangeViscousSG")==0) return SealevelchangeViscousSGEnum; + else if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum; + else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum; else if (strcmp(name,"SealevelchangeViscousE")==0) return SealevelchangeViscousEEnum; else if (strcmp(name,"CouplingTransferCount")==0) return CouplingTransferCountEnum; else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; @@ -1240,13 +1243,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"SmbPddfacIce")==0) return SmbPddfacIceEnum; else if (strcmp(name,"SmbPddfacSnow")==0) return SmbPddfacSnowEnum; else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum; - else if (strcmp(name,"SmbPrecipitationSubstep")==0) return SmbPrecipitationSubstepEnum; - else if (strcmp(name,"SmbPrecipitationsAnomaly")==0) return SmbPrecipitationsAnomalyEnum; - else if (strcmp(name,"SmbDsradiationAnomaly")==0) return SmbDsradiationAnomalyEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"SmbDlradiationAnomaly")==0) return SmbDlradiationAnomalyEnum; + if (strcmp(name,"SmbPrecipitationSubstep")==0) return SmbPrecipitationSubstepEnum; + 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; @@ -1363,13 +1366,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; else if (strcmp(name,"VxBase")==0) return VxBaseEnum; else if (strcmp(name,"VxDebris")==0) return VxDebrisEnum; - else if (strcmp(name,"Vx")==0) return VxEnum; - else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; - else if (strcmp(name,"VxObs")==0) return VxObsEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"VxShear")==0) return VxShearEnum; + if (strcmp(name,"Vx")==0) return VxEnum; + else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; + else if (strcmp(name,"VxObs")==0) return VxObsEnum; + else if (strcmp(name,"VxShear")==0) return VxShearEnum; else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum; else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; else if (strcmp(name,"VyBase")==0) return VyBaseEnum; @@ -1486,13 +1489,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum; else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum; else if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum; - else if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum; - else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum; - else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; else stage=13; } if(stage==13){ - if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum; + if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum; + else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum; + else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; + else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum; else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum; else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum; else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum; @@ -1609,13 +1612,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition211")==0) return Outputdefinition211Enum; else if (strcmp(name,"Outputdefinition212")==0) return Outputdefinition212Enum; else if (strcmp(name,"Outputdefinition213")==0) return Outputdefinition213Enum; - else if (strcmp(name,"Outputdefinition214")==0) return Outputdefinition214Enum; - else if (strcmp(name,"Outputdefinition215")==0) return Outputdefinition215Enum; - else if (strcmp(name,"Outputdefinition216")==0) return Outputdefinition216Enum; else stage=14; } if(stage==14){ - if (strcmp(name,"Outputdefinition217")==0) return Outputdefinition217Enum; + if (strcmp(name,"Outputdefinition214")==0) return Outputdefinition214Enum; + else if (strcmp(name,"Outputdefinition215")==0) return Outputdefinition215Enum; + else if (strcmp(name,"Outputdefinition216")==0) return Outputdefinition216Enum; + else if (strcmp(name,"Outputdefinition217")==0) return Outputdefinition217Enum; else if (strcmp(name,"Outputdefinition218")==0) return Outputdefinition218Enum; else if (strcmp(name,"Outputdefinition219")==0) return Outputdefinition219Enum; else if (strcmp(name,"Outputdefinition220")==0) return Outputdefinition220Enum; @@ -1732,13 +1735,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition330")==0) return Outputdefinition330Enum; else if (strcmp(name,"Outputdefinition331")==0) return Outputdefinition331Enum; else if (strcmp(name,"Outputdefinition332")==0) return Outputdefinition332Enum; - else if (strcmp(name,"Outputdefinition333")==0) return Outputdefinition333Enum; - else if (strcmp(name,"Outputdefinition334")==0) return Outputdefinition334Enum; - else if (strcmp(name,"Outputdefinition335")==0) return Outputdefinition335Enum; else stage=15; } if(stage==15){ - if (strcmp(name,"Outputdefinition336")==0) return Outputdefinition336Enum; + if (strcmp(name,"Outputdefinition333")==0) return Outputdefinition333Enum; + else if (strcmp(name,"Outputdefinition334")==0) return Outputdefinition334Enum; + else if (strcmp(name,"Outputdefinition335")==0) return Outputdefinition335Enum; + else if (strcmp(name,"Outputdefinition336")==0) return Outputdefinition336Enum; else if (strcmp(name,"Outputdefinition337")==0) return Outputdefinition337Enum; else if (strcmp(name,"Outputdefinition338")==0) return Outputdefinition338Enum; else if (strcmp(name,"Outputdefinition339")==0) return Outputdefinition339Enum; @@ -1855,13 +1858,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition449")==0) return Outputdefinition449Enum; else if (strcmp(name,"Outputdefinition404")==0) return Outputdefinition404Enum; else if (strcmp(name,"Outputdefinition450")==0) return Outputdefinition450Enum; - else if (strcmp(name,"Outputdefinition451")==0) return Outputdefinition451Enum; - else if (strcmp(name,"Outputdefinition452")==0) return Outputdefinition452Enum; - else if (strcmp(name,"Outputdefinition453")==0) return Outputdefinition453Enum; else stage=16; } if(stage==16){ - if (strcmp(name,"Outputdefinition454")==0) return Outputdefinition454Enum; + if (strcmp(name,"Outputdefinition451")==0) return Outputdefinition451Enum; + else if (strcmp(name,"Outputdefinition452")==0) return Outputdefinition452Enum; + else if (strcmp(name,"Outputdefinition453")==0) return Outputdefinition453Enum; + else if (strcmp(name,"Outputdefinition454")==0) return Outputdefinition454Enum; else if (strcmp(name,"Outputdefinition455")==0) return Outputdefinition455Enum; else if (strcmp(name,"Outputdefinition456")==0) return Outputdefinition456Enum; else if (strcmp(name,"Outputdefinition457")==0) return Outputdefinition457Enum; @@ -1978,13 +1981,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition567")==0) return Outputdefinition567Enum; else if (strcmp(name,"Outputdefinition568")==0) return Outputdefinition568Enum; else if (strcmp(name,"Outputdefinition569")==0) return Outputdefinition569Enum; - else if (strcmp(name,"Outputdefinition506")==0) return Outputdefinition506Enum; - else if (strcmp(name,"Outputdefinition570")==0) return Outputdefinition570Enum; - else if (strcmp(name,"Outputdefinition571")==0) return Outputdefinition571Enum; else stage=17; } if(stage==17){ - if (strcmp(name,"Outputdefinition572")==0) return Outputdefinition572Enum; + if (strcmp(name,"Outputdefinition506")==0) return Outputdefinition506Enum; + else if (strcmp(name,"Outputdefinition570")==0) return Outputdefinition570Enum; + else if (strcmp(name,"Outputdefinition571")==0) return Outputdefinition571Enum; + else if (strcmp(name,"Outputdefinition572")==0) return Outputdefinition572Enum; else if (strcmp(name,"Outputdefinition573")==0) return Outputdefinition573Enum; else if (strcmp(name,"Outputdefinition574")==0) return Outputdefinition574Enum; else if (strcmp(name,"Outputdefinition575")==0) return Outputdefinition575Enum; @@ -2101,13 +2104,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition685")==0) return Outputdefinition685Enum; else if (strcmp(name,"Outputdefinition686")==0) return Outputdefinition686Enum; else if (strcmp(name,"Outputdefinition687")==0) return Outputdefinition687Enum; - else if (strcmp(name,"Outputdefinition688")==0) return Outputdefinition688Enum; - else if (strcmp(name,"Outputdefinition689")==0) return Outputdefinition689Enum; - else if (strcmp(name,"Outputdefinition608")==0) return Outputdefinition608Enum; else stage=18; } if(stage==18){ - if (strcmp(name,"Outputdefinition690")==0) return Outputdefinition690Enum; + if (strcmp(name,"Outputdefinition688")==0) return Outputdefinition688Enum; + else if (strcmp(name,"Outputdefinition689")==0) return Outputdefinition689Enum; + else if (strcmp(name,"Outputdefinition608")==0) return Outputdefinition608Enum; + else if (strcmp(name,"Outputdefinition690")==0) return Outputdefinition690Enum; else if (strcmp(name,"Outputdefinition691")==0) return Outputdefinition691Enum; else if (strcmp(name,"Outputdefinition692")==0) return Outputdefinition692Enum; else if (strcmp(name,"Outputdefinition693")==0) return Outputdefinition693Enum; @@ -2224,13 +2227,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition811")==0) return Outputdefinition811Enum; else if (strcmp(name,"Outputdefinition812")==0) return Outputdefinition812Enum; else if (strcmp(name,"Outputdefinition813")==0) return Outputdefinition813Enum; - else if (strcmp(name,"Outputdefinition814")==0) return Outputdefinition814Enum; - else if (strcmp(name,"Outputdefinition815")==0) return Outputdefinition815Enum; - else if (strcmp(name,"Outputdefinition816")==0) return Outputdefinition816Enum; else stage=19; } if(stage==19){ - if (strcmp(name,"Outputdefinition817")==0) return Outputdefinition817Enum; + if (strcmp(name,"Outputdefinition814")==0) return Outputdefinition814Enum; + else if (strcmp(name,"Outputdefinition815")==0) return Outputdefinition815Enum; + else if (strcmp(name,"Outputdefinition816")==0) return Outputdefinition816Enum; + else if (strcmp(name,"Outputdefinition817")==0) return Outputdefinition817Enum; else if (strcmp(name,"Outputdefinition818")==0) return Outputdefinition818Enum; else if (strcmp(name,"Outputdefinition819")==0) return Outputdefinition819Enum; else if (strcmp(name,"Outputdefinition820")==0) return Outputdefinition820Enum; @@ -2347,13 +2350,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition930")==0) return Outputdefinition930Enum; else if (strcmp(name,"Outputdefinition931")==0) return Outputdefinition931Enum; else if (strcmp(name,"Outputdefinition932")==0) return Outputdefinition932Enum; - else if (strcmp(name,"Outputdefinition933")==0) return Outputdefinition933Enum; - else if (strcmp(name,"Outputdefinition934")==0) return Outputdefinition934Enum; - else if (strcmp(name,"Outputdefinition935")==0) return Outputdefinition935Enum; else stage=20; } if(stage==20){ - if (strcmp(name,"Outputdefinition936")==0) return Outputdefinition936Enum; + if (strcmp(name,"Outputdefinition933")==0) return Outputdefinition933Enum; + else if (strcmp(name,"Outputdefinition934")==0) return Outputdefinition934Enum; + else if (strcmp(name,"Outputdefinition935")==0) return Outputdefinition935Enum; + else if (strcmp(name,"Outputdefinition936")==0) return Outputdefinition936Enum; else if (strcmp(name,"Outputdefinition937")==0) return Outputdefinition937Enum; else if (strcmp(name,"Outputdefinition938")==0) return Outputdefinition938Enum; else if (strcmp(name,"Outputdefinition939")==0) return Outputdefinition939Enum; @@ -2470,13 +2473,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1049")==0) return Outputdefinition1049Enum; else if (strcmp(name,"Outputdefinition1004")==0) return Outputdefinition1004Enum; else if (strcmp(name,"Outputdefinition1050")==0) return Outputdefinition1050Enum; - else if (strcmp(name,"Outputdefinition1051")==0) return Outputdefinition1051Enum; - else if (strcmp(name,"Outputdefinition1052")==0) return Outputdefinition1052Enum; - else if (strcmp(name,"Outputdefinition1053")==0) return Outputdefinition1053Enum; else stage=21; } if(stage==21){ - if (strcmp(name,"Outputdefinition1054")==0) return Outputdefinition1054Enum; + if (strcmp(name,"Outputdefinition1051")==0) return Outputdefinition1051Enum; + else if (strcmp(name,"Outputdefinition1052")==0) return Outputdefinition1052Enum; + else if (strcmp(name,"Outputdefinition1053")==0) return Outputdefinition1053Enum; + else if (strcmp(name,"Outputdefinition1054")==0) return Outputdefinition1054Enum; else if (strcmp(name,"Outputdefinition1055")==0) return Outputdefinition1055Enum; else if (strcmp(name,"Outputdefinition1056")==0) return Outputdefinition1056Enum; else if (strcmp(name,"Outputdefinition1057")==0) return Outputdefinition1057Enum; @@ -2593,13 +2596,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1167")==0) return Outputdefinition1167Enum; else if (strcmp(name,"Outputdefinition1168")==0) return Outputdefinition1168Enum; else if (strcmp(name,"Outputdefinition1169")==0) return Outputdefinition1169Enum; - else if (strcmp(name,"Outputdefinition1106")==0) return Outputdefinition1106Enum; - else if (strcmp(name,"Outputdefinition1170")==0) return Outputdefinition1170Enum; - else if (strcmp(name,"Outputdefinition1171")==0) return Outputdefinition1171Enum; else stage=22; } if(stage==22){ - if (strcmp(name,"Outputdefinition1172")==0) return Outputdefinition1172Enum; + if (strcmp(name,"Outputdefinition1106")==0) return Outputdefinition1106Enum; + else if (strcmp(name,"Outputdefinition1170")==0) return Outputdefinition1170Enum; + else if (strcmp(name,"Outputdefinition1171")==0) return Outputdefinition1171Enum; + else if (strcmp(name,"Outputdefinition1172")==0) return Outputdefinition1172Enum; else if (strcmp(name,"Outputdefinition1173")==0) return Outputdefinition1173Enum; else if (strcmp(name,"Outputdefinition1174")==0) return Outputdefinition1174Enum; else if (strcmp(name,"Outputdefinition1175")==0) return Outputdefinition1175Enum; @@ -2716,13 +2719,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1285")==0) return Outputdefinition1285Enum; else if (strcmp(name,"Outputdefinition1286")==0) return Outputdefinition1286Enum; else if (strcmp(name,"Outputdefinition1287")==0) return Outputdefinition1287Enum; - else if (strcmp(name,"Outputdefinition1288")==0) return Outputdefinition1288Enum; - else if (strcmp(name,"Outputdefinition1289")==0) return Outputdefinition1289Enum; - else if (strcmp(name,"Outputdefinition1208")==0) return Outputdefinition1208Enum; else stage=23; } if(stage==23){ - if (strcmp(name,"Outputdefinition1290")==0) return Outputdefinition1290Enum; + if (strcmp(name,"Outputdefinition1288")==0) return Outputdefinition1288Enum; + else if (strcmp(name,"Outputdefinition1289")==0) return Outputdefinition1289Enum; + else if (strcmp(name,"Outputdefinition1208")==0) return Outputdefinition1208Enum; + else if (strcmp(name,"Outputdefinition1290")==0) return Outputdefinition1290Enum; else if (strcmp(name,"Outputdefinition1291")==0) return Outputdefinition1291Enum; else if (strcmp(name,"Outputdefinition1292")==0) return Outputdefinition1292Enum; else if (strcmp(name,"Outputdefinition1293")==0) return Outputdefinition1293Enum; @@ -2839,13 +2842,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1411")==0) return Outputdefinition1411Enum; else if (strcmp(name,"Outputdefinition1412")==0) return Outputdefinition1412Enum; else if (strcmp(name,"Outputdefinition1413")==0) return Outputdefinition1413Enum; - else if (strcmp(name,"Outputdefinition1414")==0) return Outputdefinition1414Enum; - else if (strcmp(name,"Outputdefinition1415")==0) return Outputdefinition1415Enum; - else if (strcmp(name,"Outputdefinition1416")==0) return Outputdefinition1416Enum; else stage=24; } if(stage==24){ - if (strcmp(name,"Outputdefinition1417")==0) return Outputdefinition1417Enum; + if (strcmp(name,"Outputdefinition1414")==0) return Outputdefinition1414Enum; + else if (strcmp(name,"Outputdefinition1415")==0) return Outputdefinition1415Enum; + else if (strcmp(name,"Outputdefinition1416")==0) return Outputdefinition1416Enum; + else if (strcmp(name,"Outputdefinition1417")==0) return Outputdefinition1417Enum; else if (strcmp(name,"Outputdefinition1418")==0) return Outputdefinition1418Enum; else if (strcmp(name,"Outputdefinition1419")==0) return Outputdefinition1419Enum; else if (strcmp(name,"Outputdefinition1420")==0) return Outputdefinition1420Enum; @@ -2962,13 +2965,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1530")==0) return Outputdefinition1530Enum; else if (strcmp(name,"Outputdefinition1531")==0) return Outputdefinition1531Enum; else if (strcmp(name,"Outputdefinition1532")==0) return Outputdefinition1532Enum; - else if (strcmp(name,"Outputdefinition1533")==0) return Outputdefinition1533Enum; - else if (strcmp(name,"Outputdefinition1534")==0) return Outputdefinition1534Enum; - else if (strcmp(name,"Outputdefinition1535")==0) return Outputdefinition1535Enum; else stage=25; } if(stage==25){ - if (strcmp(name,"Outputdefinition1536")==0) return Outputdefinition1536Enum; + if (strcmp(name,"Outputdefinition1533")==0) return Outputdefinition1533Enum; + else if (strcmp(name,"Outputdefinition1534")==0) return Outputdefinition1534Enum; + else if (strcmp(name,"Outputdefinition1535")==0) return Outputdefinition1535Enum; + else if (strcmp(name,"Outputdefinition1536")==0) return Outputdefinition1536Enum; else if (strcmp(name,"Outputdefinition1537")==0) return Outputdefinition1537Enum; else if (strcmp(name,"Outputdefinition1538")==0) return Outputdefinition1538Enum; else if (strcmp(name,"Outputdefinition1539")==0) return Outputdefinition1539Enum; @@ -3085,13 +3088,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1649")==0) return Outputdefinition1649Enum; else if (strcmp(name,"Outputdefinition1604")==0) return Outputdefinition1604Enum; else if (strcmp(name,"Outputdefinition1650")==0) return Outputdefinition1650Enum; - else if (strcmp(name,"Outputdefinition1651")==0) return Outputdefinition1651Enum; - else if (strcmp(name,"Outputdefinition1652")==0) return Outputdefinition1652Enum; - else if (strcmp(name,"Outputdefinition1653")==0) return Outputdefinition1653Enum; else stage=26; } if(stage==26){ - if (strcmp(name,"Outputdefinition1654")==0) return Outputdefinition1654Enum; + if (strcmp(name,"Outputdefinition1651")==0) return Outputdefinition1651Enum; + else if (strcmp(name,"Outputdefinition1652")==0) return Outputdefinition1652Enum; + else if (strcmp(name,"Outputdefinition1653")==0) return Outputdefinition1653Enum; + else if (strcmp(name,"Outputdefinition1654")==0) return Outputdefinition1654Enum; else if (strcmp(name,"Outputdefinition1655")==0) return Outputdefinition1655Enum; else if (strcmp(name,"Outputdefinition1656")==0) return Outputdefinition1656Enum; else if (strcmp(name,"Outputdefinition1657")==0) return Outputdefinition1657Enum; @@ -3208,13 +3211,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1767")==0) return Outputdefinition1767Enum; else if (strcmp(name,"Outputdefinition1768")==0) return Outputdefinition1768Enum; else if (strcmp(name,"Outputdefinition1769")==0) return Outputdefinition1769Enum; - else if (strcmp(name,"Outputdefinition1706")==0) return Outputdefinition1706Enum; - else if (strcmp(name,"Outputdefinition1770")==0) return Outputdefinition1770Enum; - else if (strcmp(name,"Outputdefinition1771")==0) return Outputdefinition1771Enum; else stage=27; } if(stage==27){ - if (strcmp(name,"Outputdefinition1772")==0) return Outputdefinition1772Enum; + if (strcmp(name,"Outputdefinition1706")==0) return Outputdefinition1706Enum; + else if (strcmp(name,"Outputdefinition1770")==0) return Outputdefinition1770Enum; + else if (strcmp(name,"Outputdefinition1771")==0) return Outputdefinition1771Enum; + else if (strcmp(name,"Outputdefinition1772")==0) return Outputdefinition1772Enum; else if (strcmp(name,"Outputdefinition1773")==0) return Outputdefinition1773Enum; else if (strcmp(name,"Outputdefinition1774")==0) return Outputdefinition1774Enum; else if (strcmp(name,"Outputdefinition1775")==0) return Outputdefinition1775Enum; @@ -3331,13 +3334,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1885")==0) return Outputdefinition1885Enum; else if (strcmp(name,"Outputdefinition1886")==0) return Outputdefinition1886Enum; else if (strcmp(name,"Outputdefinition1887")==0) return Outputdefinition1887Enum; - else if (strcmp(name,"Outputdefinition1888")==0) return Outputdefinition1888Enum; - else if (strcmp(name,"Outputdefinition1889")==0) return Outputdefinition1889Enum; - else if (strcmp(name,"Outputdefinition1808")==0) return Outputdefinition1808Enum; else stage=28; } if(stage==28){ - if (strcmp(name,"Outputdefinition1890")==0) return Outputdefinition1890Enum; + if (strcmp(name,"Outputdefinition1888")==0) return Outputdefinition1888Enum; + else if (strcmp(name,"Outputdefinition1889")==0) return Outputdefinition1889Enum; + else if (strcmp(name,"Outputdefinition1808")==0) return Outputdefinition1808Enum; + else if (strcmp(name,"Outputdefinition1890")==0) return Outputdefinition1890Enum; else if (strcmp(name,"Outputdefinition1891")==0) return Outputdefinition1891Enum; else if (strcmp(name,"Outputdefinition1892")==0) return Outputdefinition1892Enum; else if (strcmp(name,"Outputdefinition1893")==0) return Outputdefinition1893Enum; @@ -3454,13 +3457,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"AdaptiveTimestepping")==0) return AdaptiveTimesteppingEnum; else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum; else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; - else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; - else if (strcmp(name,"AgeAnalysis")==0) return AgeAnalysisEnum; - else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum; else stage=29; } if(stage==29){ - if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum; + if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; + else if (strcmp(name,"AgeAnalysis")==0) return AgeAnalysisEnum; + else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum; + else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum; else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum; else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum; else if (strcmp(name,"AutodiffJacobian")==0) return AutodiffJacobianEnum; @@ -3577,13 +3580,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"GLheightadvectionAnalysis")==0) return GLheightadvectionAnalysisEnum; else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum; else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; - else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; - else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; - else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; else stage=30; } if(stage==30){ - if (strcmp(name,"GenericParam")==0) return GenericParamEnum; + if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; + else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; + else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; + else if (strcmp(name,"GenericParam")==0) return GenericParamEnum; else if (strcmp(name,"GenericExternalResult")==0) return GenericExternalResultEnum; else if (strcmp(name,"Gradient1")==0) return Gradient1Enum; else if (strcmp(name,"Gradient2")==0) return Gradient2Enum; @@ -3700,13 +3703,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum; else if (strcmp(name,"MeshX")==0) return MeshXEnum; else if (strcmp(name,"MeshY")==0) return MeshYEnum; - else if (strcmp(name,"MinVel")==0) return MinVelEnum; - else if (strcmp(name,"MinVx")==0) return MinVxEnum; - else if (strcmp(name,"MinVy")==0) return MinVyEnum; else stage=31; } if(stage==31){ - if (strcmp(name,"MinVz")==0) return MinVzEnum; + if (strcmp(name,"MinVel")==0) return MinVelEnum; + else if (strcmp(name,"MinVx")==0) return MinVxEnum; + else if (strcmp(name,"MinVy")==0) return MinVyEnum; + else if (strcmp(name,"MinVz")==0) return MinVzEnum; else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum; else if (strcmp(name,"Moulin")==0) return MoulinEnum; else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; @@ -3823,13 +3826,13 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"IntrusionMelt")==0) return IntrusionMeltEnum; else if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum; else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum; - else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; - else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; - else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; else stage=32; } if(stage==32){ - if (strcmp(name,"Tetra")==0) return TetraEnum; + if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; + else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; + else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; + else if (strcmp(name,"Tetra")==0) return TetraEnum; else if (strcmp(name,"TetraInput")==0) return TetraInputEnum; else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum; else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum; diff --git a/src/c/shared/Enum/issmenums.jl b/src/c/shared/Enum/issmenums.jl index 666e47e26..98ba8ce2d 100644 --- a/src/c/shared/Enum/issmenums.jl +++ b/src/c/shared/Enum/issmenums.jl @@ -210,6 +210,7 @@ FrictionGammaEnum FrictionLawEnum FrictionLinearizeEnum + FrictionMaxCoefficientEnum FrictionPseudoplasticityExponentEnum FrictionU0Enum FrictionThresholdSpeedEnum @@ -422,6 +423,8 @@ OutputFileNameEnum OutputFilePointerEnum OutputdefinitionEnum + PddfacIceEnum + PddfacSnowEnum QmuErrNameEnum QmuInNameEnum QmuIsdakotaEnum @@ -3986,6 +3989,7 @@ function EnumToString(enum::IssmEnum) if(enum==FrictionGammaEnum) return "FrictionGamma" end if(enum==FrictionLawEnum) return "FrictionLaw" end if(enum==FrictionLinearizeEnum) return "FrictionLinearize" end + if(enum==FrictionMaxCoefficientEnum) return "FrictionMaxCoefficient" end if(enum==FrictionPseudoplasticityExponentEnum) return "FrictionPseudoplasticityExponent" end if(enum==FrictionU0Enum) return "FrictionU0" end if(enum==FrictionThresholdSpeedEnum) return "FrictionThresholdSpeed" end @@ -4198,6 +4202,8 @@ function EnumToString(enum::IssmEnum) if(enum==OutputFileNameEnum) return "OutputFileName" end if(enum==OutputFilePointerEnum) return "OutputFilePointer" end if(enum==OutputdefinitionEnum) return "Outputdefinition" end + if(enum==PddfacIceEnum) return "PddfacIce" end + if(enum==PddfacSnowEnum) return "PddfacSnow" end if(enum==QmuErrNameEnum) return "QmuErrName" end if(enum==QmuInNameEnum) return "QmuInName" end if(enum==QmuIsdakotaEnum) return "QmuIsdakota" end @@ -7762,6 +7768,7 @@ function StringToEnum(name::String) if(name=="FrictionGamma") return FrictionGammaEnum end if(name=="FrictionLaw") return FrictionLawEnum end if(name=="FrictionLinearize") return FrictionLinearizeEnum end + if(name=="FrictionMaxCoefficient") return FrictionMaxCoefficientEnum end if(name=="FrictionPseudoplasticityExponent") return FrictionPseudoplasticityExponentEnum end if(name=="FrictionU0") return FrictionU0Enum end if(name=="FrictionThresholdSpeed") return FrictionThresholdSpeedEnum end @@ -7974,6 +7981,8 @@ function StringToEnum(name::String) if(name=="OutputFileName") return OutputFileNameEnum end if(name=="OutputFilePointer") return OutputFilePointerEnum end if(name=="Outputdefinition") return OutputdefinitionEnum end + if(name=="PddfacIce") return PddfacIceEnum end + if(name=="PddfacSnow") return PddfacSnowEnum end if(name=="QmuErrName") return QmuErrNameEnum end if(name=="QmuInName") return QmuInNameEnum end if(name=="QmuIsdakota") return QmuIsdakotaEnum end diff --git a/src/m/classes/SMBpddSicopolis.m b/src/m/classes/SMBpddSicopolis.m index 153264197..a94cee8a4 100644 --- a/src/m/classes/SMBpddSicopolis.m +++ b/src/m/classes/SMBpddSicopolis.m @@ -16,6 +16,8 @@ s0t = NaN; rlaps = 0; isfirnwarming = 0; + pdd_fac_ice = 0; + pdd_fac_snow = 0; steps_per_step = 1 averaging = 0 requested_outputs = {}; @@ -71,6 +73,8 @@ self.isfirnwarming = 1; self.desfac = -log(2.0)/1000; self.rlaps = 7.4; + self.pdd_fac_ice = 7.28; + self.pdd_fac_snow = 2.73; self.requested_outputs={'default'}; end % }}} @@ -85,6 +89,8 @@ md = checkfield(md,'fieldname','smb.rlaps','>=',0,'numel',1); md = checkfield(md,'fieldname','smb.monthlytemperatures','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 12]); md = checkfield(md,'fieldname','smb.precipitation','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 12]); + md = checkfield(md,'fieldname','smb.pdd_fac_ice','>',0,'numel',1); + md = checkfield(md,'fieldname','smb.pdd_fac_snow','>',0,'numel',1); end md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]); @@ -108,6 +114,8 @@ function disp(self) % {{{ fielddisplay(self,'isfirnwarming','is firnwarming (Reeh 1991) activated (0 or 1, default is 1)'); fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'); fielddisplay(self,'averaging','averaging methods from short to long steps'); + fielddisplay(self,'pdd_fac_ice','Pdd factor for ice for all the domain [mm ice equiv/day/degree C]'); + fielddisplay(self,'pdd_fac_snow','Pdd factor for snow for all the domain [mm ice equiv/day/degree C]'); disp(sprintf('%51s 0: Arithmetic (default)',' ')); disp(sprintf('%51s 1: Geometric',' ')); disp(sprintf('%51s 2: Harmonic',' ')); @@ -124,6 +132,8 @@ function marshall(self,prefix,md,fid) % {{{ 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','pdd_fac_ice','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','pdd_fac_snow','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); diff --git a/src/m/classes/SMBpddSicopolis.py b/src/m/classes/SMBpddSicopolis.py index 9d627a5e6..502f1a785 100644 --- a/src/m/classes/SMBpddSicopolis.py +++ b/src/m/classes/SMBpddSicopolis.py @@ -26,6 +26,8 @@ def __init__(self, *args): # {{{ self.s0t = np.nan self.rlaps = 0 self.isfirnwarming = 0 + self.pdd_fac_ice = 0 + self.pdd_fac_snow = 0 self.steps_per_step = 1 self.averaging = 0 self.requested_outputs = [] @@ -50,7 +52,9 @@ def __repr__(self): # {{{ s += '{}\n'.format(fielddisplay(self, 'desfac', 'desertification elevation factor (default is -log(2.0)/1000)')) s += '{}\n'.format(fielddisplay(self, 'isfirnwarming', 'is firnwarming (Reeh 1991) activated (0 or 1, default is 1)')) s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step')) - s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) + s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) + s += '{}\n'.format(fielddisplay(self, 'pdd_fac_ice', 'Pdd factor for ice for all the domain [mm ice equiv/day/degree C]')) + s += '{}\n'.format(fielddisplay(self, 'pdd_fac_snow', 'Pdd factor for snow for all the domain [mm ice equiv/day/degree C]')) s += '\t\t{}\n'.format('0: Arithmetic (default)') s += '\t\t{}\n'.format('1: Geometric') s += '\t\t{}\n'.format('2: Harmonic') @@ -100,6 +104,8 @@ def setdefaultparameters(self): # {{{ self.isfirnwarming = 1 self.desfac = -np.log(2.0) / 1000 self.rlaps = 7.4 + self.pdd_fac_ice = 7.28 + self.pdd_fac_snow = 2.73 self.requested_outputs = ['default'] return self # }}} @@ -114,6 +120,8 @@ def checkconsistency(self, md, solution, analyses): # {{{ md = checkfield(md, 'fieldname', 'smb.rlaps', '>=', 0, 'numel', 1) md = checkfield(md, 'fieldname', 'smb.monthlytemperatures', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 12]) md = checkfield(md, 'fieldname', 'smb.precipitation', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 12]) + md = checkfield(md, 'fieldname', 'smb.pdd_fac_ice', '>', 0, 'numel', 1) + md = checkfield(md, 'fieldname', 'smb.pdd_fac_snow', '>', 0, 'numel', 1) 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) @@ -128,6 +136,8 @@ def marshall(self, prefix, md, fid): # {{{ 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', 'pdd_fac_ice', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'pdd_fac_snow', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'monthlytemperatures', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'precipitation', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'temperature_anomaly', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) diff --git a/src/m/classes/frictionjosh.m b/src/m/classes/frictionjosh.m index 011e468b3..57686ff68 100644 --- a/src/m/classes/frictionjosh.m +++ b/src/m/classes/frictionjosh.m @@ -5,10 +5,11 @@ classdef frictionjosh properties (SetAccess=public) - coefficient = NaN; + coefficient = NaN; pressure_adjusted_temperature = NaN; - gamma = 0.; - effective_pressure_limit = 0; + gamma = 0.0; + effective_pressure_limit = 0.0; + coefficient_max = 0.0; end methods function self = extrude(self,md) % {{{ @@ -30,8 +31,11 @@ %Default gamma: 1 self.gamma = 1.; + % Default max friction coefficient: 300 + self.coefficient_max = 300.0; + %Default 0 - self.effective_pressure_limit = 0; + self.effective_pressure_limit = 0.0; end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -43,6 +47,7 @@ md = checkfield(md,'fieldname','friction.pressure_adjusted_temperature','NaN',1,'Inf',1); md = checkfield(md,'fieldname','friction.gamma','numel',1,'NaN',1,'Inf',1,'>',0.); md = checkfield(md,'fieldname','friction.effective_pressure_limit','numel',[1],'>=',0); + md = checkfield(md,'fieldname','friction.coefficient_max','numel',1,'NaN',1,'Inf',1,'>',0.); %Check that temperature is provided md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'Inf',1,'size','universal'); @@ -53,14 +58,16 @@ function disp(self) % {{{ fielddisplay(self,'pressure_adjusted_temperature','friction pressure_adjusted_temperature (T - Tpmp) [K]'); fielddisplay(self,'gamma','(T - Tpmp)/gamma [K]'); fielddisplay(self,'effective_pressure_limit','Neff do not allow to fall below a certain limit: effective_pressure_limit*rho_ice*g*thickness (default 0)'); + fielddisplay(self,'coefficient_max', 'effective friction C = min(coefficient_max, sqrt(exp(T_b(modern) - T_b(t))/gamma) * coefficient)'); end % }}} function marshall(self,prefix,md,fid) % {{{ - WriteData(fid,prefix,'name','md.friction.law','data',9,'format','Integer'); - WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'class','friction','object',self,'fieldname','pressure_adjusted_temperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'class','friction','object',self,'fieldname','gamma','format','Double'); - WriteData(fid,prefix,'object',self,'class','friction','fieldname','effective_pressure_limit','format','Double'); + WriteData(fid,prefix, 'name', 'md.friction.law', 'data',9, 'format', 'Integer'); + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype',1, 'timeserieslength',md.mesh.numberofvertices+1, 'yts',md.constants.yts); + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'pressure_adjusted_temperature', 'format', 'DoubleMat', 'mattype',1, 'timeserieslength',md.mesh.numberofvertices+1, 'yts',md.constants.yts); + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'gamma', 'format', 'Double'); + WriteData(fid,prefix, 'object',self, 'class', 'friction', 'fieldname', 'effective_pressure_limit', 'format', 'Double'); + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'coefficient_max', 'format', 'Double'); end % }}} end end diff --git a/src/m/classes/frictionjosh.py b/src/m/classes/frictionjosh.py index a333bc361..45a56f9cc 100644 --- a/src/m/classes/frictionjosh.py +++ b/src/m/classes/frictionjosh.py @@ -14,13 +14,13 @@ class frictionjosh(object): """ def __init__(self): # {{{ - self.coefficient = np.nan + self.coefficient = np.nan self.pressure_adjusted_temperature = np.nan - self.gamma = 0 - self.effective_pressure_limit = 0 + self.gamma = 0. + self.effective_pressure_limit = 0. + self.coefficient_max = 0. self.setdefaultparameters() - #self.requested_outputs = [] # }}} def __repr__(self): # {{{ @@ -30,7 +30,7 @@ def __repr__(self): # {{{ s += '{}\n'.format(fielddisplay(self, 'pressure_adjusted_temperature', 'friction pressure_adjusted_temperature (T - Tpmp) [K]')) s += '{}\n'.format(fielddisplay(self, 'gamma', '(T - Tpmp)/gamma [K]')) s += '{}\n'.format(fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)')) - #s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested')) + s += '{}\n'.format(fielddisplay(self, 'coefficient_max', 'effective friction C = min(coefficient_max, sqrt(exp(T_b(modern) - T_b(t))/gamma) * coefficient)')) return s # }}} @@ -41,34 +41,41 @@ def extrude(self, md): # {{{ # }}} def setdefaultparameters(self): # {{{ + + #Default gamma: 1 self.gamma = 1. - #self.requested_outputs = ['default'] - self.effective_pressure_limit = 0 - return self - # }}} - def defaultoutputs(self, md): # {{{ - list = [] - return list + #Default 0. + self.effective_pressure_limit = 0.0 + + #Default max friction coefficient: 300 + self.coefficient_max = 300. + + return self # }}} def checkconsistency(self, md, solution, analyses): # {{{ + # Early return if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses: return md + md = checkfield(md, 'fieldname', 'friction.coefficient', 'timeseries', 1, 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'friction.pressure_adjusted_temperature','NaN',1,'Inf',1) md = checkfield(md, 'fieldname', 'friction.gamma','numel',1,'NaN',1,'Inf',1,'>',0.) md = checkfield(md, 'fieldname', 'friction.effective_pressure_limit', 'numel', [1], '>=', 0) + md = checkfield(md,'fieldname', 'friction.coefficient_max', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0.) + # Check that temperature is provided md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'Inf',1,'size','universal') return md # }}} def marshall(self, prefix, md, fid): # {{{ - WriteData(fid,prefix,'name','md.friction.law','data',9,'format','Integer') - WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - WriteData(fid,prefix,'class','friction','object',self,'fieldname','pressure_adjusted_temperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - WriteData(fid,prefix,'class','friction','object',self,'fieldname','gamma','format','Double') - WriteData(fid,prefix,'object',self,'class','friction','fieldname','effective_pressure_limit','format','Double') + WriteData(fid,prefix, 'name', 'md.friction.law', 'data',9, 'format', 'Integer') + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype',1, 'timeserieslength',md.mesh.numberofvertices+1, 'yts',md.constants.yts) + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'pressure_adjusted_temperature', 'format', 'DoubleMat', 'mattype',1, 'timeserieslength',md.mesh.numberofvertices+1, 'yts',md.constants.yts) + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'gamma', 'format', 'Double') + WriteData(fid,prefix, 'object',self, 'class', 'friction', 'fieldname', 'effective_pressure_limit', 'format', 'Double') + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'coefficient_max', 'format', 'Double') # }}} diff --git a/src/m/classes/model.m b/src/m/classes/model.m index 1586acabc..8623f0eaf 100644 --- a/src/m/classes/model.m +++ b/src/m/classes/model.m @@ -198,7 +198,12 @@ %2022 Oct 28 if ~isa(md.debris,'debris'); md.debris=debris(); end %Mmetransport: Jun 2022: - if ~isa(md.mmemasstransport,'mmemasstransport'); md.mmemasstransport=mmemasstransport(); end; + if ~isa(md.mmemasstransport,'mmemasstransport'); md.mmemasstransport=mmemasstransport(); end + % 2026 February 18 + if isa(md.friction, 'frictionjosh') && md.friction.coefficient_max==0; md.friction.coefficient_max=300; end + % 2026 February 20 + if isa(md.smb, 'SMBpddSicopolis') && md.smb.pdd_fac_ice==0; md.smb.pdd_fac_ice=7.28; end + if isa(md.smb, 'SMBpddSicopolis') && md.smb.pdd_fac_snow==0; md.smb.pdd_fac_snow=2.73; end end% }}} end methods @@ -598,7 +603,7 @@ function disp(self) % {{{ %copy model md1=md; - %recover options: + %recover optoins: options=pairoptions(varargin{:}); %some checks @@ -662,13 +667,13 @@ function disp(self) % {{{ %loop over model fields model_fields=fields(md1); - for i=1:length(model_fields) + for i=1:length(model_fields), %get field field=md1.(model_fields{i}); fieldsize=size(field); if isobject(field), %recursive call object_fields=fields(md1.(model_fields{i})); - for j=1:length(object_fields) + for j=1:length(object_fields), %get field field=md1.(model_fields{i}).(object_fields{j}); fieldsize=size(field); @@ -677,7 +682,7 @@ function disp(self) % {{{ md2.(model_fields{i}).(object_fields{j})=field(pos_node,:); elseif (fieldsize(1)==numberofvertices1+1) md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)]; - %size = number of elements * n + %size = number of elements * n elseif fieldsize(1)==numberofelements1 md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); elseif (fieldsize(1)==numberofelements1+1) @@ -690,7 +695,7 @@ function disp(self) % {{{ md2.(model_fields{i})=field(pos_node,:); elseif (fieldsize(1)==numberofvertices1+1) md2.(model_fields{i})=[field(pos_node,:); field(end,:)]; - %size = number of elements * n + %size = number of elements * n elseif fieldsize(1)==numberofelements1 md2.(model_fields{i})=field(pos_elem,:); elseif (fieldsize(1)==numberofelements1+1) @@ -706,13 +711,6 @@ function disp(self) % {{{ md2.mesh.numberofvertices=numberofvertices2; md2.mesh.elements=elements_2; - % Extract ISMIP6 basal tf field - if isa(md1.basalforcings, 'basalforcingsismip6') - for i=1:numel(md.basalforcings.tf) - md2.basalforcings.tf{i} = [md1.basalforcings.tf{i}(pos_node); md1.basalforcings.tf{i}(end)]; - end - end - %mesh.uppervertex mesh.lowervertex if isa(md1.mesh,'mesh3dprisms'), md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node); diff --git a/src/m/interp/ProfileValues.py b/src/m/interp/ProfileValues.py new file mode 100644 index 000000000..5eee099fa --- /dev/null +++ b/src/m/interp/ProfileValues.py @@ -0,0 +1,35 @@ +#!/usr/bin/env python3 +from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d +from InterpFromMeshToMesh3d import InterpFromMeshToMesh3d +from project2d import project2d +import numpy as np +from numpy import ceil, mean, linspace, arange, transpose, reshape, ones, shape + +def ProfileValues(md,data,xprof,yprof,resolution): + """ + PROFILEVALUES - compute the value of a field on a vertical profile + + This routine gets the value of a given field of the model on + a point given by its coordinates + + Usage: + z,data=ProfileValues(md,data,xcoord,ycoord,resolution) + """ + + #Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system + offset=10e-3 + bed=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.base,1),xprof,yprof)+offset + surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),xprof,yprof)-offset + + #Some useful parameters + layers=ceil(mean(md.geometry.thickness)/resolution) + Z=arange(bed,surface+resolution,resolution) + Z=transpose(Z) + X=xprof*ones(shape(Z)) + Y=yprof*ones(shape(Z)) + data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X,Y,Z,np.nan) + + #Get np.ndarray in list + data_interp=data_interp[0] + + return Z, data_interp diff --git a/src/m/mech/basalstress.m b/src/m/mech/basalstress.m index 4350c7dbe..68fdc38b9 100644 --- a/src/m/mech/basalstress.m +++ b/src/m/mech/basalstress.m @@ -16,9 +16,10 @@ sealevel = 0; p_ice = g*rho_ice*md.geometry.thickness; -if isfield(md.friction, 'coupling') +if isprop(md.friction, 'coupling') coupling = md.friction.coupling; else + warning(sprintf('md.friction.coupling is not found. Default coupling is set to 0.')); coupling = 0; end diff --git a/src/m/mech/basalstress.py b/src/m/mech/basalstress.py index a99f2eeeb..6b14d1c7a 100644 --- a/src/m/mech/basalstress.py +++ b/src/m/mech/basalstress.py @@ -1,4 +1,5 @@ from friction import friction +from frictionweertman import frictionweertman from frictionschoof import frictionschoof from averaging import averaging import numpy as np @@ -16,7 +17,7 @@ def basalstress(md): ''' #Check md.friction class - if not isinstance(md.friction,friction): + if not (isinstance(md.friction,friction) | isinstance(md.friction,frictionschoof) | isinstance(md.friction,frictionweertman)): raise Exception('Error: md.friction only supports "friction.m" class.') #Compute effective pressure @@ -74,6 +75,12 @@ def basalstress(md): Cmax=averaging(md,md.friction.Cmax,0) alpha2 = (C**2 * ub**(m-1))/(1 + (C**2/(Cmax*N))**(1/m)*ub)**m + + elif isinstance(md.friction,frictionweertman): + m = averaging(md,md.friction.m,0.0) + C = md.friction.C + alpha2 = C**2 * ub**(1/m-1) + else: raise Exception('not supported yet') diff --git a/src/m/modules/InterpFromMeshToMesh3d.py b/src/m/modules/InterpFromMeshToMesh3d.py index 743ded6cf..b5805b06c 100644 --- a/src/m/modules/InterpFromMeshToMesh3d.py +++ b/src/m/modules/InterpFromMeshToMesh3d.py @@ -18,10 +18,10 @@ def InterpFromMeshToMesh3d(index, x, y, z, data, x_prime, y_prime, z_prime, defa """ # Check usage - nargs = len(args) - if nargs != 9: - print(InterpFromMeshToMesh3d.__doc__) - raise Exception('Wrong usage (see above)') + #nargs = len(args) + #if nargs != 9: + # print(InterpFromMeshToMesh3d.__doc__) + # raise Exception('Wrong usage (see above)') # Call Python module data_prime = InterpFromMeshToMesh3d_python(index, x, y, z, data, x_prime, y_prime, z_prime, default_value) diff --git a/src/m/plot/applyoptions.py b/src/m/plot/applyoptions.py index 80e6b0862..e4074a73d 100644 --- a/src/m/plot/applyoptions.py +++ b/src/m/plot/applyoptions.py @@ -231,12 +231,12 @@ def applyoptions(md, data, options, fig, axgrid, gridindex): if options.exist('log'): formatter = mpl.ticker.LogFormatterSciNotation(base=options.getfieldvalue('log')) - try: - #NOTE: axis made by ImageGrid contains "cax". If axis generated by "plt.subplots()", we have to manually generate cax... + #NOTE: axis made by ImageGrid contains "cax". If axis generated by "plt.subplots()", we have to manually generate cax... + if hasattr(ax, 'cax'): # type(ax) == mpl_toolkits.axes_grid1.mpl_axes.Axes cb = mpl.colorbar.ColorbarBase(ax.cax, cmap=cmap, norm=norm, extend=extend, format=formatter) - except: + else: from mpl_toolkits.axes_grid1 import make_axes_locatable driver = make_axes_locatable(ax) cax = driver.append_axes(options.getfieldvalue('colorbarpos','right'), diff --git a/src/m/plot/plot_landsat.m b/src/m/plot/plot_landsat.m index 56dba6ce7..147f6fabd 100644 --- a/src/m/plot/plot_landsat.m +++ b/src/m/plot/plot_landsat.m @@ -1,152 +1,170 @@ -% Explain -% This funtion loads Landsat Image Mosaic Antarctica (LIMA) for background image. -% -% Usage -% plot_landsat(md,data,options,plotlines,plotcols,i), -% function plot_landsat(md,data,options,plotlines,plotcols,i), + % Explain + % This funtion loads Landsat Image Mosaic Antarctica (LIMA) for background image. + % + % Usage + % plot_landsat(md,data,options,plotlines,plotcols,i), + % + + %process mesh and data + [x2d y2d z2d elements2d is2d isplanet]=processmesh(md,[],options); + [data datatype]=processdata(md,data,options); + + %check is2d + if ~is2d, + error('buildgridded error message: gridded not supported for 3d meshes, project on a layer'); + end -%process mesh and data -[x2d y2d z2d elements2d is2d isplanet]=processmesh(md,[],options); -[data datatype]=processdata(md,data,options); - -%check is2d -if ~is2d, - error('buildgridded error message: gridded not supported for 3d meshes, project on a layer'); -end - -%Get some options -transparency = getfieldvalue(options,'transparency',.3); -highres = getfieldvalue(options,'highres',0); - -% get xlim, and ylim -xlim=getfieldvalue(options,'xlim',[min(x2d) max(x2d)])/getfieldvalue(options,'unit',1); -ylim=getfieldvalue(options,'ylim',[min(y2d) max(y2d)])/getfieldvalue(options,'unit',1); - -if md.mesh.epsg == 3031 & (isempty(md.radaroverlay.pwr) | isempty(md.radaroverlay.x) | isempty(md.radaroverlay.y) | length(size(md.radaroverlay.pwr)) < 3), % Antarctica region {{{ - if highres, - disp(' LIMA with geotiff'), % {{{ - disp('WARNING : this image shoud be collected with geocoded tif file'); - % find merged mosaic landsat image {{{ - limapath = {'/drive/project_inwoo/issm/Data/LIMA/AntarcticaLandsat.tif'}; - pos = zeros(length(limapath),1); - for ii = 1:length(limapath) - if exist(limapath{ii}), pos(ii) = 1; end + %Get some options + transparency = getfieldvalue(options,'transparency',.3); + highres = getfieldvalue(options,'highres',0); + + %Get xlim, and ylim + xlim=getfieldvalue(options,'xlim',[min(x2d) max(x2d)])/getfieldvalue(options,'unit',1); + ylim=getfieldvalue(options,'ylim',[min(y2d) max(y2d)])/getfieldvalue(options,'unit',1); + + pwr = md.radaroverlay.pwr; + xm = md.radaroverlay.x; + ym = md.radaroverlay.y; + nx = numel(xm); + ny = numel(ym); + + if md.mesh.epsg == 3031 & (isempty(md.radaroverlay.pwr) | isempty(md.radaroverlay.x) | isempty(md.radaroverlay.y) | length(size(md.radaroverlay.pwr)) < 3), % Antarctica region {{{ + if highres, + disp(' LIMA with geotiff'), % {{{ + disp('WARNING : this image shoud be collected with geocoded tif file'); + % find merged mosaic landsat image {{{ + limapath = {'/drive/project_inwoo/issm/Data/LIMA/AntarcticaLandsat.tif'}; + pos = zeros(length(limapath),1); + for ii = 1:length(limapath) + if exist(limapath{ii}), pos(ii) = 1; end + end + limapath = limapath{find(pos)}; + fprintf(' LIMA path is %s\n', limapath); + % }}} + + % read image + im = imread(limapath); + + % Region of LIMA data set + info = gdalinfo(limapath); % get geotiff info + xm = info.xmin + info.dx*[0:info.nx-1]; + ym = info.ymax - info.dy*[0:info.ny-1]; + + % find region of model at LIMA + offset = 1e+4; + posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset)); + posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset)); + % }}} + else + disp(' LIMA with reduced tiff'), + % find merged mosaic landsat image {{{ + limapath = {'/drive/project_inwoo/issm/Data/LIMA/tiff_90pct/00000-20080319-092059124.tif'}; + pos = zeros(length(limapath),1); + for ii = 1:length(limapath) + if exist(limapath{ii}), pos(ii) = 1; end + end + + if sum(pos) == 0, + fprintf('download website : https://lima.usgs.gov/fullcontinent.php\n'); + error('Landsat image at Antarctic region should be downloaded at above website'); + end + limapath = limapath{find(pos)}; + fprintf(' LIMA path is %s\n', limapath); + % }}} + + % read image + im = imread(limapath); + + % Region of LIMA data set + info = gdalinfo(limapath); % get geotiff info + xm = info.xmin + info.dx*[0:info.nx-1]; + ym = info.ymax - info.dy*[0:info.ny-1]; + + % find region of model at LIMA + offset = 1e+4; + posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset)); + posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset)); end - limapath = limapath{find(pos)}; - fprintf(' LIMA path is %s\n', limapath); + + % update region of radaroverlay + md.radaroverlay.x = xm(posx); + md.radaroverlay.y = ym(posy); + md.radaroverlay.pwr = im(posy, posx,:); % }}} + elseif length(size(md.radaroverlay.pwr)) == 3 + % it already contains LIMA image. + elseif md.mesh.epsg == 3431 & (isempty(md.radaroverlay.pwr) | isempty(md.radaroverlay.x) | isempty(md.radaroverlay.y) | length(size(md.radaroverlay.pwr)) < 3), % Greenladn region + error('Greenland region is not yet available.'); + else + error('Check md.mesh.epsg, available Landsat regeion is at Antarctica (EPSG:3031)'); + end + + %Check dataset + assert(ndims(pwr) ~= 3, 'Error: Check ndims(md.radaroverlay.pwr) should be equal to 3.'); + assert(size(pwr) == [nx, ny, 3], 'Error: Given md.radaroverlay.pwr shoule be equal to (nx, ny, 3)'); - % read image - im = imread(limapath); + if any(diff(xm)) < 0 + disp('WARNING: md.radaroverlay.x should be increasing order.') + xm = flip(xm); + pwr= flip(pwr, 1); + end + if any(diff(ym)) < 0 + disp('WARNING: md.radaroverlay.y should be increasing order.') + ym = flip(ym); + pwr= flip(pwr, 2); + end - % Region of LIMA data set - info = gdalinfo(limapath); % get geotiff info - xm = info.xmin + info.dx*[0:info.nx-1]; - ym = info.ymax - info.dy*[0:info.ny-1]; + %Process image from model + final = double(pwr)/double(max(md.radaroverlay.pwr(:))); %rescale between 0 and 1 - % find region of model at LIMA - offset = 1e+4; - posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset)); - posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset)); - % }}} + %Prepare grid + if size(md.radaroverlay.x,1)==1 | size(md.radaroverlay.x,2)==1, + data_grid=InterpFromMeshToGrid(elements2d,x2d/getfieldvalue(options,'unit',1),y2d/getfieldvalue(options,'unit',1),data,xm,ym,NaN); + %data_grid=InterpFromMeshToGrid(md.mesh.elements,md.mesh.x/getfieldvalue(options,'unit',1),md.mesh.y/getfieldvalue(options,'unit',1),data,x_m,y_m,NaN); else - disp(' LIMA with reduced tiff'), - % find merged mosaic landsat image {{{ - limapath = {'/drive/project_inwoo/issm/Data/LIMA/tiff_90pct/00000-20080319-092059124.tif'}; - pos = zeros(length(limapath),1); - for ii = 1:length(limapath) - if exist(limapath{ii}), pos(ii) = 1; end - end - - if sum(pos) == 0, - fprintf('download website : https://lima.usgs.gov/fullcontinent.php\n'); - error('Landsat image at Antarctic region should be downloaded at above website'); - end - limapath = limapath{find(pos)}; - fprintf(' LIMA path is %s\n', limapath); - % }}} + X = md.radaroverlay.x; + Y = md.radaroverlay.y; + data_grid=InterpFromMeshToMesh2d(elements2d,x2d,y2d,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X)); + %data_grid=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X)); + xm=X(1,:); ym=Y(:,1); + end - % read image - im = imread(limapath); + data_nan=isnan(data_grid); + if exist(options,'caxis'), + caxis_opt=getfieldvalue(options,'caxis'); + data_grid(find(data_gridcaxis_opt(2)))=caxis_opt(2); + data_min=caxis_opt(1); + data_max=caxis_opt(2); + else + data_min=min(data_grid(:)); + data_max=max(data_grid(:)); + end + colorm = getcolormap(options); + image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(colorm)/(data_max-data_min))),colorm); - % Region of LIMA data set - info = gdalinfo(limapath); % get geotiff info - xm = info.xmin + info.dx*[0:info.nx-1]; - ym = info.ymax - info.dy*[0:info.ny-1]; + alpha=ones(size(data_grid)); + alpha(find(~data_nan))=transparency; + alpha=repmat(alpha,[1 1 3]); - % find region of model at LIMA - offset = 1e+4; - posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset)); - posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset)); - % }}} + final=alpha.*final+(1-alpha).*image_rgb; + + %Select plot area + subplotmodel(plotlines,plotcols,i,options); + + h=imagesc(xm*getfieldvalue(options,'unit',1),ym*getfieldvalue(options,'unit',1),final); + + %last step: mesh gridded? + if exist(options,'edgecolor'), + A=elements(:,1); B=elements(:,2); C=elements(:,3); + patch('Faces',[A B C],'Vertices', [x y z],'FaceVertexCData',data_grid(1)*ones(size(x)),'FaceColor','none','EdgeColor',getfieldvalue(options,'edgecolor')); end - % update region of radaroverlay - md.radaroverlay.x = xm(posx); - md.radaroverlay.y = ym(posy); - md.radaroverlay.pwr = im(posy, posx,:); - % }}} -elseif length(size(md.radaroverlay.pwr)) == 3 - % it already contains LIMA image. -elseif md.mesh.epsg == 3431 & (isempty(md.radaroverlay.pwr) | isempty(md.radaroverlay.x) | isempty(md.radaroverlay.y) | length(size(md.radaroverlay.pwr)) < 3), % Greenladn region - error('Greenland region is not yet available.'); -else - error('Check md.mesh.epsg, available Landsat regeion is at Antarctica (EPSG:3031)'); -end - -%Process image from model -final = double(md.radaroverlay.pwr)/double(max(md.radaroverlay.pwr(:))); %rescale between 0 and 1 - -%Prepare grid -if size(md.radaroverlay.x,1)==1 | size(md.radaroverlay.x,2)==1, - x_m = md.radaroverlay.x; - y_m = md.radaroverlay.y; - data_grid=InterpFromMeshToGrid(elements2d,x2d/getfieldvalue(options,'unit',1),y2d/getfieldvalue(options,'unit',1),data,x_m,y_m,NaN); - %data_grid=InterpFromMeshToGrid(md.mesh.elements,md.mesh.x/getfieldvalue(options,'unit',1),md.mesh.y/getfieldvalue(options,'unit',1),data,x_m,y_m,NaN); -else - X = md.radaroverlay.x; - Y = md.radaroverlay.y; - data_grid=InterpFromMeshToMesh2d(elements2d,x2d,y2d,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X)); - %data_grid=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X)); - x_m=X(1,:); y_m=Y(:,1); -end - -data_nan=isnan(data_grid); -if exist(options,'caxis'), - caxis_opt=getfieldvalue(options,'caxis'); - data_grid(find(data_gridcaxis_opt(2)))=caxis_opt(2); - data_min=caxis_opt(1); - data_max=caxis_opt(2); -else - data_min=min(data_grid(:)); - data_max=max(data_grid(:)); -end -colorm = getcolormap(options); -image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(colorm)/(data_max-data_min))),colorm); - -alpha=ones(size(data_grid)); -alpha(find(~data_nan))=transparency; -alpha=repmat(alpha,[1 1 3]); - -final=alpha.*final+(1-alpha).*image_rgb; - -%Select plot area -subplotmodel(plotlines,plotcols,i,options); - -h=imagesc(x_m*getfieldvalue(options,'unit',1),y_m*getfieldvalue(options,'unit',1),final); - -%last step: mesh gridded? -if exist(options,'edgecolor'), - A=elements(:,1); B=elements(:,2); C=elements(:,3); - patch('Faces',[A B C],'Vertices', [x y z],'FaceVertexCData',data_grid(1)*ones(size(x)),'FaceColor','none','EdgeColor',getfieldvalue(options,'edgecolor')); -end - -%Apply options -if ~isnan(data_min), - options=changefieldvalue(options,'caxis',[data_min data_max]); % force caxis so that the colorbar is ready -end -options=addfielddefault(options,'axis','xy equal off'); % default axis -applyoptions(md,data,options); -end + %Apply options + if ~isnan(data_min), + options=changefieldvalue(options,'caxis',[data_min data_max]); % force caxis so that the colorbar is ready + end + options=addfielddefault(options,'axis','xy equal off'); % default axis + applyoptions(md,data,options); + end diff --git a/src/m/plot/plot_landsat.py b/src/m/plot/plot_landsat.py index c54fe4959..5a0ed494b 100644 --- a/src/m/plot/plot_landsat.py +++ b/src/m/plot/plot_landsat.py @@ -45,8 +45,8 @@ def plot_landsat(md,data,options,fig,axgrid,gridindex): isunit = options.getfieldvalue('unit',1) #Get xlim, and ylim - xlim=options.getfieldvalue('xlim',np.array([min(x2d),max(x2d)]))/isunit - ylim=options.getfieldvalue('ylim',np.array([min(y2d),max(y2d)]))/isunit + xlim=np.array(options.getfieldvalue('xlim',[min(x2d),max(x2d)]))/isunit + ylim=np.array(options.getfieldvalue('ylim',[min(y2d),max(y2d)]))/isunit pwr = md.radaroverlay.pwr xm = md.radaroverlay.x @@ -113,8 +113,17 @@ def plot_landsat(md,data,options,fig,axgrid,gridindex): if np.ndim(pwr) != 3: raise Exception('Error: Check np.ndim(md.radaroverlay.pwr) should be equal to 3.') + if np.any(np.diff(xm) < 0): + print('WARNING: md.radaroverlay.x should be increasing order.') + xm = np.flip(xm) + pwr= np.flip(pwr,axis=0) + if np.any(np.diff(md.radaroverlay.y) < 0): + print('WARNING: md.radaroverlay.y should be increasing order.') + ym = np.flip(ym) + pwr= np.flip(pwr,axis=1) + #Check image size - #shape of image should be (ny, nx, band) + #shape of image should be (nx, ny, band) nx = len(xm) ny = len(ym) if (np.shape(pwr)[0]==nx) & (np.shape(pwr)[1]==ny): @@ -146,11 +155,19 @@ def plot_landsat(md,data,options,fig,axgrid,gridindex): else: data_min=np.nanmin(data_grid) data_max=np.nanmax(data_grid) + caxis_opt=[data_min, data_max] # Back-up caxis for 'applyoptions' options.addfielddefault('colormap',plt.cm.viridis) cmap = getcolormap(copy.deepcopy(options)) #TODO: Matlab version #image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(colorm)/(data_max-data_min))),colorm); + #Set log scale + if options.exist('log'): + #NOTE: Tricy part for log scale dataset. "log" scale option does not rely on "processdata.py" function. + data_grid=np.log(data_grid)/np.log(options.getfieldvalue('log')) + data_min =np.log(data_min)/np.log(options.getfieldvalue('log')) + data_max =np.log(data_max)/np.log(options.getfieldvalue('log')) + #NOTE: Python version for ind2rgb image_rgb = cmap((data_grid-data_min)/(data_max-data_min)) @@ -167,6 +184,7 @@ def plot_landsat(md,data,options,fig,axgrid,gridindex): xmax = max(xm)/isunit ymin = min(ym)/isunit ymax = max(ym)/isunit + #Draw RGB image h=ax.imshow(final, extent=[xmin, xmax, ymin, ymax], origin='lower') @@ -179,6 +197,6 @@ def plot_landsat(md,data,options,fig,axgrid,gridindex): #Apply options if ~np.isnan(data_min): - options.changefieldvalue('caxis',[data_min, data_max]) # force caxis so that the colorbar is ready + options.changefieldvalue('caxis',caxis_opt) # force caxis so that the colorbar is ready options.addfielddefault('axis','xy equal off') # default axis applyoptions(md,data,options,fig,axgrid,gridindex) diff --git a/src/m/plot/plot_unit.py b/src/m/plot/plot_unit.py index 85c11b73b..0706d18a9 100644 --- a/src/m/plot/plot_unit.py +++ b/src/m/plot/plot_unit.py @@ -174,10 +174,10 @@ def plot_unit(x, y, z, elements, data, is2d, isplanet, datatype, options, fig, a if is2d: if np.ma.is_masked(data): if hasattr(np, 'isin'): #Numpy 2017+ - tmp = np.isin(index, np.where(data.mask)) + tmp = np.isin(range(len(data)), np.where(data.mask)) else: #For backward compatibility - tmp = np.in1d(index, np.where(data.mask)) - EltMask = np.asarray([np.any(tmp) for index in elements]) + tmp = np.in1d(range(len(data)), np.where(data.mask)) + EltMask = np.asarray([np.any(tmp[index]) for index in elements]) triangles = mpl.tri.Triangulation(x, y, elements, EltMask) else: triangles = mpl.tri.Triangulation(x, y, elements) diff --git a/src/m/plot/processdata.py b/src/m/plot/processdata.py index e770adc23..dfb455890 100644 --- a/src/m/plot/processdata.py +++ b/src/m/plot/processdata.py @@ -158,9 +158,9 @@ def processdata(md, data, options): # }}} - # layer procesiong? {{{ + # layer processing? {{{ if options.getfieldvalue('layer',0)>=1: - procdata=project2d(md,data,options.getfieldvalue('layer')) + procdata=project2d(md,procdata,options.getfieldvalue('layer')) if options.getfieldvalue('depthaverage',0): raise Exception('Error: "depthaverage" option is not supported in Python.') diff --git a/test/Archives/Archive245.arch b/test/Archives/Archive245.arch index afe512437..4dafb9726 100644 Binary files a/test/Archives/Archive245.arch and b/test/Archives/Archive245.arch differ diff --git a/test/Archives/Archive258.arch b/test/Archives/Archive258.arch index 48c6811b5..92461413d 100644 Binary files a/test/Archives/Archive258.arch and b/test/Archives/Archive258.arch differ