Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 26549) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 26550) @@ -75,4 +75,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.smb.Sizeini",SmbSizeiniEnum); iomodel->FetchDataToInput(inputs,elements,"md.smb.aValue",SmbAValueEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.dulwrfValue",SmbDulwrfValueEnum); iomodel->FetchDataToInput(inputs,elements,"md.smb.teValue",SmbTeValueEnum); iomodel->FetchDataToInput(inputs,elements,"md.smb.szaValue",SmbSzaValueEnum); @@ -242,4 +243,5 @@ case SMBgembEnum: parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIdx",SmbAIdxEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.eIdx",SmbEIdxEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.swIdx",SmbSwIdxEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.denIdx",SmbDenIdxEnum)); @@ -261,7 +263,9 @@ parameters->AddObject(iomodel->CopyConstantObject("md.smb.isturbulentflux",SmbIsturbulentfluxEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.isconstrainsurfaceT",SmbIsconstrainsurfaceTEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdeltaLWup",SmbIsdeltaLWupEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.teThresh",SmbTeThreshEnum)); break; case SMBpddEnum: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26549) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26550) @@ -3778,4 +3778,5 @@ IssmDouble teValue=1.0; IssmDouble aValue=0.0; + IssmDouble dulwrfValue=0.0; IssmDouble szaValue=0.0; IssmDouble cotValue=0.0; @@ -3784,4 +3785,5 @@ IssmDouble dt,time,smb_dt; int aIdx=0; + int eIdx=0; int denIdx=0; int dsnowIdx=0; @@ -3813,7 +3815,9 @@ bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux; bool isconstrainsurfaceT=false; + bool isdeltaLWup=false; IssmDouble init_scaling=0.0; IssmDouble thermo_scaling=1.0; IssmDouble adThresh=1023.0; + IssmDouble teThresh=10; /*}}}*/ /*Output variables:{{{ */ @@ -3865,4 +3869,5 @@ parameters->FindParam(&smb_dt,SmbDtEnum); /*time period for the smb solution, usually smaller than the glaciological dt*/ parameters->FindParam(&aIdx,SmbAIdxEnum); + parameters->FindParam(&eIdx,SmbEIdxEnum); parameters->FindParam(&denIdx,SmbDenIdxEnum); parameters->FindParam(&swIdx,SmbSwIdxEnum); @@ -3881,7 +3886,9 @@ parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum); parameters->FindParam(&isconstrainsurfaceT,SmbIsconstrainsurfaceTEnum); + parameters->FindParam(&isdeltaLWup,SmbIsdeltaLWupEnum); parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum); parameters->FindParam(&thermo_scaling,SmbThermoDeltaTScalingEnum); parameters->FindParam(&adThresh,SmbAdThreshEnum); + parameters->FindParam(&teThresh,SmbTeThreshEnum); /*}}}*/ /*Retrieve inputs: {{{*/ @@ -4073,4 +4080,5 @@ Input *teValue_input= this->GetInput(SmbTeValueEnum,timeinputs); _assert_(teValue_input); Input *aValue_input= this->GetInput(SmbAValueEnum,timeinputs); _assert_(aValue_input); + Input *dulwrfValue_input= this->GetInput(SmbDulwrfValueEnum,timeinputs); _assert_(dulwrfValue_input); Input *szaValue_input= this->GetInput(SmbSzaValueEnum,timeinputs); _assert_(szaValue_input); Input *cotValue_input= this->GetInput(SmbCotValueEnum,timeinputs); _assert_(cotValue_input); @@ -4088,4 +4096,5 @@ pAir_input->GetInputValue(&pAir,gauss); // screen level air pressure [Pa] teValue_input->GetInputValue(&teValue,gauss); // Emissivity [0-1] + dulwrfValue_input->GetInputValue(&dulwrfValue,gauss); // LWup perturbation [W m-2] aValue_input->GetInputValue(&aValue,gauss); // Albedo [0 1] szaValue_input->GetInputValue(&szaValue,gauss); // Solar Zenith Angle [degree] @@ -4113,5 +4122,5 @@ } /*Thermal profile computation:*/ - if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid(),isconstrainsurfaceT); + if(isthermal)thermo(&EC, &T, &ulw, re, dz, d, swf, dlw, Ta, V, eAir, pAir, eIdx, teValue, dulwrfValue, teThresh, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid(),isconstrainsurfaceT,isdeltaLWup); /*Change in thickness of top cell due to evaporation/condensation assuming same density as top cell. Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 26549) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 26550) @@ -694,5 +694,5 @@ } /*}}}*/ -void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* pulwrf, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid, bool isconstrainsurfaceT) { /*{{{*/ +void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* pulwrf, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble eIdx, IssmDouble teValue, IssmDouble dulwrfValue, IssmDouble teThresh, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid, bool isconstrainsurfaceT, bool isdeltaLWup) { /*{{{*/ /* ENGLACIAL THERMODYNAMICS*/ @@ -1010,7 +1010,12 @@ // upward longwave contribution - //ulw = - (SB * pow(Ts,4.0)* teValue) * dt; // - deltatest here - IssmDouble deltatest=0; - ulw = - (SB * pow(Ts,4.0)* teValue - deltatest) * dt; // - deltatest here + IssmDouble deltaULW=0.0; + IssmDouble emissivity=1.0; + //If user wants to set a upward long wave bias + if(isdeltaLWup) deltaULW = dulwrfValue; + //If user wants to directly set emissivity, or grain radius is larger than the + // threshold, or eIdx is 2 and we have wet snow or ice, use prescribed emissivity + if(eIdx==0 || (teThresh - re[0])0.001)) emissivity = teValue; + ulw = - (SB * pow(Ts,4.0)* emissivity + deltaULW) * dt; ulwrf = ulwrf - ulw/dt0; @@ -2202,6 +2207,9 @@ if (fabs(adThresh - 820) < Dtol){ // ERA5 new aIdx=1, swIdx=0, MODIS 820 - M0 = max(1.8230 - (0.1645 * log(C)),0.25); - M1 = max(2.5134 - (0.3244 * log(C)),0.25); + //M0 = max(1.8230 - (0.1645 * log(C)),0.25); + //M1 = max(2.5134 - (0.3244 * log(C)),0.25); + // ERA5 new aIdx=1, swIdx=0, MODIS 820, p90 new (ERA5, 40 and e97) + M0 = max(1.3045 - (0.0988 * log(C)),0.25); + M1 = max(1.3694 - (0.1354 * log(C)),0.25); } else{ @@ -2246,6 +2254,9 @@ if (fabs(adThresh - 820) < Dtol){ // ERA5 new aIdx=1, swIdx=0, MODIS 820 - M0 = max(1.4174 - (0.1037 * log(C)),0.25); - M1 = max(2.2010 - (0.2460 * log(C)),0.25); + //M0 = max(1.4174 - (0.1037 * log(C)),0.25); + //M1 = max(2.2010 - (0.2460 * log(C)),0.25); + // ERA5 new aIdx=1, swIdx=0, MODIS 820, p90 new (ERA5 lwt, 40 and e97) + M0 = max(1.2138 - (0.1057 * log(C)),0.25); + M1 = max(1.4946 - (0.1607 * log(C)),0.25); } else{ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 26549) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 26550) @@ -34,5 +34,5 @@ void albedo(IssmDouble** a, IssmDouble** adiff, int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid); void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble dswdiff, IssmDouble as, IssmDouble asdiff, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid); -void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* pulwrf, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid, bool isconstrainsurfaceT); +void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* pulwrf, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble eIdx, IssmDouble teValue, IssmDouble dulwrfValue, IssmDouble teThresh, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid, bool isconstrainsurfaceT, bool isdeltaLWup); void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* pRa, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid); void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pF, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble Ra, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY, IssmDouble dIce, int sid); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26549) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26550) @@ -443,4 +443,5 @@ syn keyword cConstant SmbDtEnum syn keyword cConstant SmbEnum +syn keyword cConstant SmbEIdxEnum syn keyword cConstant SmbFEnum syn keyword cConstant SmbInitDensityScalingEnum @@ -451,4 +452,5 @@ syn keyword cConstant SmbIsdelta18oEnum syn keyword cConstant SmbIsdensificationEnum +syn keyword cConstant SmbIsdeltaLWupEnum syn keyword cConstant SmbIsfirnwarmingEnum syn keyword cConstant SmbIsgraingrowthEnum @@ -478,4 +480,5 @@ syn keyword cConstant SmbT0dryEnum syn keyword cConstant SmbT0wetEnum +syn keyword cConstant SmbTeThreshEnum syn keyword cConstant SmbTdiffEnum syn keyword cConstant SmbThermoDeltaTScalingEnum @@ -913,4 +916,5 @@ syn keyword cConstant SmbDiniEnum syn keyword cConstant SmbDlwrfEnum +syn keyword cConstant SmbDulwrfValueEnum syn keyword cConstant SmbDswrfEnum syn keyword cConstant SmbDswdiffrfEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26549) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26550) @@ -437,4 +437,5 @@ SmbDtEnum, SmbEnum, + SmbEIdxEnum, SmbFEnum, SmbInitDensityScalingEnum, @@ -445,4 +446,5 @@ SmbIsdelta18oEnum, SmbIsdensificationEnum, + SmbIsdeltaLWupEnum, SmbIsfirnwarmingEnum, SmbIsgraingrowthEnum, @@ -472,4 +474,5 @@ SmbT0dryEnum, SmbT0wetEnum, + SmbTeThreshEnum, SmbTdiffEnum, SmbThermoDeltaTScalingEnum, @@ -909,4 +912,5 @@ SmbDiniEnum, SmbDlwrfEnum, + SmbDulwrfValueEnum, SmbDswrfEnum, SmbDswdiffrfEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26549) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26550) @@ -445,4 +445,5 @@ case SmbDtEnum : return "SmbDt"; case SmbEnum : return "Smb"; + case SmbEIdxEnum : return "SmbEIdx"; case SmbFEnum : return "SmbF"; case SmbInitDensityScalingEnum : return "SmbInitDensityScaling"; @@ -453,4 +454,5 @@ case SmbIsdelta18oEnum : return "SmbIsdelta18o"; case SmbIsdensificationEnum : return "SmbIsdensification"; + case SmbIsdeltaLWupEnum : return "SmbIsdeltaLWup"; case SmbIsfirnwarmingEnum : return "SmbIsfirnwarming"; case SmbIsgraingrowthEnum : return "SmbIsgraingrowth"; @@ -480,4 +482,5 @@ case SmbT0dryEnum : return "SmbT0dry"; case SmbT0wetEnum : return "SmbT0wet"; + case SmbTeThreshEnum : return "SmbTeThresh"; case SmbTdiffEnum : return "SmbTdiff"; case SmbThermoDeltaTScalingEnum : return "SmbThermoDeltaTScaling"; @@ -915,4 +918,5 @@ case SmbDiniEnum : return "SmbDini"; case SmbDlwrfEnum : return "SmbDlwrf"; + case SmbDulwrfValueEnum : return "SmbDulwrfValue"; case SmbDswrfEnum : return "SmbDswrf"; case SmbDswdiffrfEnum : return "SmbDswdiffrf"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26549) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26550) @@ -454,4 +454,5 @@ else if (strcmp(name,"SmbDt")==0) return SmbDtEnum; else if (strcmp(name,"Smb")==0) return SmbEnum; + else if (strcmp(name,"SmbEIdx")==0) return SmbEIdxEnum; else if (strcmp(name,"SmbF")==0) return SmbFEnum; else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum; @@ -462,4 +463,5 @@ else if (strcmp(name,"SmbIsdelta18o")==0) return SmbIsdelta18oEnum; else if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum; + else if (strcmp(name,"SmbIsdeltaLWup")==0) return SmbIsdeltaLWupEnum; else if (strcmp(name,"SmbIsfirnwarming")==0) return SmbIsfirnwarmingEnum; else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum; @@ -489,4 +491,5 @@ else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum; else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum; + else if (strcmp(name,"SmbTeThresh")==0) return SmbTeThreshEnum; else if (strcmp(name,"SmbTdiff")==0) return SmbTdiffEnum; else if (strcmp(name,"SmbThermoDeltaTScaling")==0) return SmbThermoDeltaTScalingEnum; @@ -503,11 +506,11 @@ else if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum; else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum; - else if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum; - else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum; - else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum; + if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum; + else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum; + else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum; + else if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum; else if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum; else if (strcmp(name,"StressbalanceRequestedOutputs")==0) return StressbalanceRequestedOutputsEnum; @@ -626,11 +629,11 @@ else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum; else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum; - else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; - else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum; - else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum; + if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; + else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum; + else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; + else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum; else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum; else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum; @@ -749,11 +752,11 @@ else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum; else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum; - else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum; - else if (strcmp(name,"HydrologyTws")==0) return HydrologyTwsEnum; - else if (strcmp(name,"HydrologyTwsSpc")==0) return HydrologyTwsSpcEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"HydrologyTwsAnalysis")==0) return HydrologyTwsAnalysisEnum; + if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum; + else if (strcmp(name,"HydrologyTws")==0) return HydrologyTwsEnum; + else if (strcmp(name,"HydrologyTwsSpc")==0) return HydrologyTwsSpcEnum; + else if (strcmp(name,"HydrologyTwsAnalysis")==0) return HydrologyTwsAnalysisEnum; else if (strcmp(name,"HydrologyWatercolumnMax")==0) return HydrologyWatercolumnMaxEnum; else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; @@ -872,11 +875,11 @@ else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum; else if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum; - else if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum; - else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum; - else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum; + if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum; + else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum; + else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum; + else if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum; else if (strcmp(name,"SealevelchangeGUsubelOcean")==0) return SealevelchangeGUsubelOceanEnum; else if (strcmp(name,"SealevelchangeGEsubelOcean")==0) return SealevelchangeGEsubelOceanEnum; @@ -936,4 +939,5 @@ else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum; else if (strcmp(name,"SmbDlwrf")==0) return SmbDlwrfEnum; + else if (strcmp(name,"SmbDulwrfValue")==0) return SmbDulwrfValueEnum; else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum; else if (strcmp(name,"SmbDswdiffrf")==0) return SmbDswdiffrfEnum; @@ -994,12 +998,12 @@ else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum; - else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum; else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum; else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum; else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; + else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum; else if (strcmp(name,"SmbTz")==0) return SmbTzEnum; @@ -1117,12 +1121,12 @@ else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum; else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum; - else if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum; else if (strcmp(name,"Outputdefinition30")==0) return Outputdefinition30Enum; else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum; else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum; + else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum; else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum; else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; @@ -1240,12 +1244,12 @@ else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum; else if (strcmp(name,"Channel")==0) return ChannelEnum; - else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum; else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum; else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum; else if (strcmp(name,"Closed")==0) return ClosedEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"Colinear")==0) return ColinearEnum; + else if (strcmp(name,"Colinear")==0) return ColinearEnum; else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; else if (strcmp(name,"Contact")==0) return ContactEnum; @@ -1363,12 +1367,12 @@ else if (strcmp(name,"Intersect")==0) return IntersectEnum; else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum; - else if (strcmp(name,"J")==0) return JEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"J")==0) return JEnum; else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum; else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum; else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum; + else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum; else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum; else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; @@ -1486,12 +1490,12 @@ else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; - else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; + else stage=13; + } + if(stage==13){ + if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; - else stage=13; - } - if(stage==13){ - if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; + else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 26549) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 26550) @@ -24,4 +24,5 @@ isconstrainsurfaceT; isclimatology; + isdeltaLWup; %inputs: @@ -42,5 +43,7 @@ %optional inputs: aValue = NaN; %Albedo forcing at every element. Used only if aIdx == 0, or density exceeds adThresh - teValue = NaN; %Outward longwave radiation thermal emissivity forcing at every element (default in code is 1) + teValue = NaN; %Outward longwave radiation thermal emissivity forcing at every element (default in code is 1). + %Used only if eIdx== 0, or effective grain radius exceeds teThresh + dulwrfValue = NaN; %Delta with which to perturn the long wave radiation upwards. Use if isdeltaLWup is true; % Initialization of snow properties @@ -59,9 +62,14 @@ %settings: aIdx = NaN; %method for calculating albedo and subsurface absorption (default is 1) - % 0: direct input from aValue parameter + % 0: direct input from aValue parameter, no use of adThresh % 1: effective grain radius [Gardner & Sharp, 2009] % 2: effective grain radius [Brun et al., 1992; LeFebre et al., 2003]], with swIdx=1, SW penetration follows grain size in 3 spectral bands (Brun et al., 1992) % 3: density and cloud amount [Greuell & Konzelmann, 1994] % 4: exponential time decay & wetness [Bougamont & Bamber, 2005] + + eIdx = NaN; %method for calculating emissivity (default is 1) + % 0: direct input from teValue parameter, no use of teThresh + % 1: default value of 1, in areas with grain radius below teThresh + % 2: default value of 1, in areas with grain radius below teThresh and areas of dry snow (not bare ice or wet) at the surface swIdx = NaN; % apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1, if aIdx=2 function of effective radius (Brun et al., 1992) or else dependent on snow density (taken from Bassford, 2002)) @@ -111,4 +119,7 @@ %or else apply direct input value from aValue, allowing albedo to be altered. %Default value is rho water (1023 kg m-3). + teThresh = NaN; %Apply eIdx method to all areas with grain radii below this value, + %or else apply direct input value from teValue, allowing emissivity to be altered. + %Default value is a effective grain radius of 10 mm. %densities: @@ -152,5 +163,6 @@ fielddisplay(self,'isdensification','run densification module (default true)'); fielddisplay(self,'isturbulentflux','run turbulant heat fluxes module (default true)'); - fielddisplay(self,'isconstrainsurfaceT','constrain surface temperatures to air temperature, turn off EC and surface flux contribution to surface temperature change'); + fielddisplay(self,'isconstrainsurfaceT','constrain surface temperatures to air temperature, turn off EC and surface flux contribution to surface temperature change (default false)'); + fielddisplay(self,'isdeltaLWup','set to true to invoke a bias in the long wave upward spatially, specified by dulwrfValue (default false)'); fielddisplay(self,'Ta','2 m air temperature, in Kelvin'); fielddisplay(self,'V','wind speed (m s-1)'); @@ -183,5 +195,11 @@ '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'}) + fielddisplay(self,'dulwrfValue','Specified bias to be applied to the outward long wave radiation every element (W/m-2, +upward)'); fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)'); + fielddisplay(self,'teThresh',{'Apply eIdx method to all areas with effective grain radius below this value,','or else apply direct input value from teValue, allowing emissivity to be altered.'}); + fielddisplay(self,'eIdx',{'method for calculating emissivity (default is 1)',... + '0: direct input from teValue parameter, no use of teThresh',... + '1: default value of 1, in areas with grain radius below teThresh',... + '2: default value of 1, in areas with grain radius below teThresh and areas of dry snow (not bare ice or wet) at the surface'}); %snow properties init @@ -299,5 +317,5 @@ self.cciceValue=project3d(md,'vector',self.cciceValue,'type','element'); end - if (aIdx == 0) & ~isnan(self.aValue) + if ~isnan(self.aValue) self.aValue=project3d(md,'vector',self.aValue,'type','element'); end @@ -305,5 +323,4 @@ self.teValue=project3d(md,'vector',self.teValue,'type','element'); end - end % }}} @@ -322,6 +339,8 @@ self.isturbulentflux=1; self.isconstrainsurfaceT=0; + self.isdeltaLWup=0; self.aIdx = 1; + self.eIdx = 1; self.swIdx = 1; self.denIdx = 2; @@ -348,7 +367,9 @@ self.K = 7; self.adThresh = 1023; + self.teThresh = 10; self.teValue = ones(mesh.numberofelements,1); self.aValue = self.aSnow*ones(mesh.numberofelements,1); + self.dulwrfValue = zeros(mesh.numberofelements,1); self.dswdiffrf=0.0*ones(mesh.numberofelements,1); @@ -384,4 +405,5 @@ md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]); md = checkfield(md,'fieldname','smb.isconstrainsurfaceT','values',[0 1]); + md = checkfield(md,'fieldname','smb.isdeltaLWup','values',[0 1]); md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value @@ -400,6 +422,8 @@ md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1); + md = checkfield(md,'fieldname','smb.dulwrfValue','timeseries',1,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4]); + md = checkfield(md,'fieldname','smb.eIdx','NaN',1,'Inf',1,'values',[0,1,2]); md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]); md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5,6,7]); @@ -414,8 +438,8 @@ md = checkfield(md,'fieldname','smb.ThermoDeltaTScaling','NaN',1,'Inf',1,'>=',0,'<=',1); md = checkfield(md,'fieldname','smb.adThresh','NaN',1,'Inf',1,'>=',0); - + md = checkfield(md,'fieldname','smb.teThresh','NaN',1,'Inf',1,'>=',0); + + md = checkfield(md,'fieldname','smb.aValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1); switch self.aIdx, - case 0 - md = checkfield(md,'fieldname','smb.aValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1); case {1 2} md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'>=',.64,'<=',.89); @@ -460,4 +484,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','isturbulentflux','format','Boolean'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','isconstrainsurfaceT','format','Boolean'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','isdeltaLWup','format','Boolean'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); @@ -483,4 +508,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','aIdx','format','Integer'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','eIdx','format','Integer'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','swIdx','format','Integer'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','denIdx','format','Integer'); @@ -496,7 +522,9 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','adThresh','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','teThresh','format','Double'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','aValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','teValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dulwrfValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','szaValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','cotValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); @@ -533,4 +561,7 @@ error('If GEMB forcing dswdiffrf is transient, it must have the same time steps as input Ta in the final row!'); end + if size(md.smb.dulwrfValue,2)>1 & any(md.smb.dulwrfValue(end,:) - md.smb.Ta(end,:) ~= 0) + error('If GEMB forcing dulwrfValue is transient, it must have the same time steps as input Ta in the final row!'); + end if size(md.smb.aValue,2)>1 & any(md.smb.aValue(end,:) - md.smb.Ta(end,:) ~= 0) error('If GEMB forcing aValue is transient, it must have the same time steps as input Ta in the final row!'); Index: /issm/trunk-jpl/src/m/classes/SMBgemb.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 26549) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 26550) @@ -31,4 +31,5 @@ self.isturbulentflux = 0 self.isconstrainsurfaceT = 0 + self.isdeltaLWup = 0 self.isclimatology = np.nan @@ -50,5 +51,6 @@ # Optional inputs self.aValue = np.nan # Albedo forcing at every element. Used only if aIdx == 0, or density exceeds adThresh. - self.teValue = np.nan # Outward longwave radiation thermal emissivity forcing at every element (default in code is 1) + self.teValue = np.nan # Outward longwave radiation thermal emissivity forcing at every element (default in code is 1), Used only if eIdx== 0, or effective grain radius exceeds teThresh + dulwrfValue = np.nan #Delta with which to perturn the long wave radiation upwards. Use if isdeltaLWup is true # Initialization of snow properties @@ -67,9 +69,14 @@ # Settings self.aIdx = np.nan # method for calculating albedo and subsurface absorption (default is 1) - # 0: direct input from aValue parameter + # 0: direct input from aValue parameter, no use of adThresh # 1: effective grain radius [Gardner & Sharp, 2009] # 2: effective grain radius [Brun et al., 1992; LeFebre et al., 2003], with swIdx=1, SW penetration follows grain size in 3 spectral bands (Brun et al., 1992) # 3: density and cloud amount [Greuell & Konzelmann, 1994] # 4: exponential time decay & wetness [Bougamont & Bamber, 2005] + + self.eIdx = np.nan #method for calculating emissivity (default is 1) + # 0: direct input from teValue parameter, no use of teThresh + # 1: default value of 1, in areas with grain radius below teThresh + # 2: default value of 1, in areas with grain radius below teThresh and areas of dry snow (not bare ice or wet) at the surface self.swIdx = np.nan # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1, if aIdx=2 function of effective radius (Brun et al., 1992) or else dependent on snow density (taken from Bassford, 2002)) @@ -119,4 +126,7 @@ # or else apply direct input value from aValue, allowing albedo to be altered. # Default value is rho water (1023 kg m-3). + teThresh = np.nan #Apply eIdx method to all areas with grain radii below this value, + #or else apply direct input value from teValue, allowing emissivity to be altered. + #Default value is a effective grain radius of 10 mm. # Densities @@ -158,5 +168,6 @@ s += '{}\n'.format(fielddisplay(self, 'isdensification', 'run densification module (default true)')) s += '{}\n'.format(fielddisplay(self, 'isturbulentflux', 'run turbulant heat fluxes module (default true)')) - s += '{}\n'.format(fielddisplay(self, 'isconstrainsurfaceT', 'constrain surface temperatures to air temperature, turn off EC and surface flux contribution to surface temperature change')) + s += '{}\n'.format(fielddisplay(self, 'isconstrainsurfaceT', 'constrain surface temperatures to air temperature, turn off EC and surface flux contribution to surface temperature change (default false)')) + s += '{}\n'.format(fielddisplay(self, 'isdeltaLWup', 'set to true to invoke a bias in the long wave upward spatially, specified by dulwrfValue (default false)')) s += '{}\n'.format(fielddisplay(self, 'Ta', '2 m air temperature, in Kelvin')) s += '{}\n'.format(fielddisplay(self, 'V', 'wind speed (m s-1)')) @@ -188,5 +199,12 @@ '3: density and cloud amount [Greuell & Konzelmann, 1994]', '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'])) + + s += '{}\n'.format(fielddisplay(self, 'dulwrfValue', 'Specified bias to be applied to the outward long wave radiation every element (W/m-2, +upward)')) s += '{}\n'.format(fielddisplay(self, 'teValue', 'Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)')) + s += '{}\n'.format(fielddisplay(self, 'teThresh', ['Apply eIdx method to all areas with effective grain radius below this value,', 'or else apply direct input value from teValue, allowing emissivity to be altered.'])) + s += '{}\n'.format(fielddisplay(self, 'eIdx', ['method for calculating emissivity (default is 1)', + '0: direct input from teValue parameter, no use of teThresh', + '1: default value of 1, in areas with grain radius below teThresh', + '2: default value of 1, in areas with grain radius below teThresh and areas of dry snow (not bare ice or wet) at the surface'])) # Snow properties init @@ -286,5 +304,5 @@ self.cciceValue=project3d(md,'vector',self.cciceValue,'type','element'); - if (self.aIdx == 0) and (not np.isnan(self.aValue)): + if not np.isnan(self.aValue): self.aValue = project3d(md, 'vector', self.aValue, 'type', 'element') if not np.isnan(self.teValue): @@ -308,6 +326,8 @@ self.isturbulentflux = 1 self.isconstrainsurfaceT = 0 + self.isdeltaLWup = 0 self.aIdx = 1 + self.eIdx = 1 self.swIdx = 1 self.denIdx = 2 @@ -334,7 +354,9 @@ self.K = 7 self.adThresh = 1023 + self.teThresh = 10 self.teValue = np.ones((mesh.numberofelements,)) self.aValue = self.aSnow * np.ones(mesh.numberofelements,) + self.dulwrfValue = np.zeros((mesh.numberofelements,)) self.dswdiffrf = 0.0 * np.ones(mesh.numberofelements,) @@ -373,4 +395,5 @@ md = checkfield(md, 'fieldname', 'smb.isdensification', 'values', [0, 1]) md = checkfield(md, 'fieldname', 'smb.isturbulentflux', 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'smb.isdeltaLWup', 'values',[0, 1]) md = checkfield(md, 'fieldname', 'smb.isconstrainsurfaceT', 'values', [0, 1]) @@ -390,6 +413,8 @@ md = checkfield(md, 'fieldname', 'smb.teValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1) + md = checkfield(md, 'fieldname', 'smb.dulwrfValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'smb.aIdx', 'NaN', 1, 'Inf', 1, 'values', [0, 1, 2, 3, 4]) + md = checkfield(md, 'fieldname', 'smb.eIdx', 'NaN', 1, 'Inf', 1, 'values', [0, 1, 2]) md = checkfield(md, 'fieldname', 'smb.swIdx', 'NaN', 1, 'Inf', 1, 'values', [0, 1]) md = checkfield(md, 'fieldname', 'smb.denIdx', 'NaN', 1, 'Inf', 1, 'values', [1, 2, 3, 4, 5, 6, 7]) @@ -404,8 +429,8 @@ md = checkfield(md, 'fieldname', 'smb.ThermoDeltaTScaling', 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1) md = checkfield(md, 'fieldname', 'smb.adThresh', 'NaN', 1, 'Inf', 1, '>=', 0) - - if self.aIdx == 0: - md = checkfield(md, 'fieldname', 'smb.aValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1) - elif isinstance(self.aIdx, (list, type(np.array([1, 2])))) and (self.aIdx == [1, 2] or (1 in self.aIdx and 2 in self.aIdx)): + md = checkfield(md, 'fieldname', 'smb.teThresh', 'NaN', 1, 'Inf',1,'>=',0) + + md = checkfield(md, 'fieldname', 'smb.aValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1) + if isinstance(self.aIdx, (list, type(np.array([1, 2])))) and (self.aIdx == [1, 2] or (1 in self.aIdx and 2 in self.aIdx)): md = checkfield(md, 'fieldname', 'smb.aSnow', 'NaN', 1, 'Inf', 1, '> = ', .64, '< = ', .89) md = checkfield(md, 'fieldname', 'smb.aIce', 'NaN', 1, 'Inf', 1, '> = ', .27, '< = ', .58) @@ -446,4 +471,5 @@ WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isturbulentflux', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isconstrainsurfaceT', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isdeltaLWup','format','Boolean') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Ta', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) @@ -469,4 +495,5 @@ WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'aIdx', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'eIdx', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'swIdx', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'denIdx', 'format', 'Integer') @@ -482,7 +509,9 @@ WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'K', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'adThresh', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'teThresh', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'aValue', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'teValue', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dulwrfValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts) WriteData(fid,prefix,'object',self,'class','smb','fieldname','szaValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts) WriteData(fid,prefix,'object',self,'class','smb','fieldname','cotValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts) @@ -514,4 +543,6 @@ if ((np.ndim(self.aValue)>1) & np.any(self.aValue[-1] - self.Ta[-1] != 0)): raise IOError('If GEMB forcing aValue is transient, it must have the same time steps as input Ta in the final row!') + if ((np.ndim(self.dulwrfValue)>1) & np.any(self.dulwrfValue[-1] - self.Ta[-1] != 0)): + raise IOError('If GEMB forcing dulwrfValue is transient, it must have the same time steps as input Ta in the final row!') if ((np.ndim(self.szaValue)>1) & np.any(self.szaValue[-1] - self.Ta[-1] != 0)): raise IOError('If GEMB forcing szaValue is transient, it must have the same time steps as input Ta in the final row!')