Changeset 26550


Ignore:
Timestamp:
11/05/21 15:15:23 (3 years ago)
Author:
schlegel
Message:

CHG: update how we deal with emissivity, giving user more options to change value on wet snow/ice, and when grain size exceeds a threshold value. Add option to run with a bias correction to outgoing long wave radiation.

Location:
issm/trunk-jpl/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

    r26477 r26550  
    7575                        iomodel->FetchDataToInput(inputs,elements,"md.smb.Sizeini",SmbSizeiniEnum);
    7676                        iomodel->FetchDataToInput(inputs,elements,"md.smb.aValue",SmbAValueEnum);
     77                        iomodel->FetchDataToInput(inputs,elements,"md.smb.dulwrfValue",SmbDulwrfValueEnum);
    7778                        iomodel->FetchDataToInput(inputs,elements,"md.smb.teValue",SmbTeValueEnum);
    7879                        iomodel->FetchDataToInput(inputs,elements,"md.smb.szaValue",SmbSzaValueEnum);
    … …  
    242243                case SMBgembEnum:
    243244                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIdx",SmbAIdxEnum));
     245                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.eIdx",SmbEIdxEnum));
    244246                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.swIdx",SmbSwIdxEnum));
    245247                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.denIdx",SmbDenIdxEnum));
    … …  
    261263                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.isturbulentflux",SmbIsturbulentfluxEnum));
    262264                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.isconstrainsurfaceT",SmbIsconstrainsurfaceTEnum));
     265                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdeltaLWup",SmbIsdeltaLWupEnum));
    263266                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
    264267                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum));
    265268                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum));
     269                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.teThresh",SmbTeThreshEnum));
    266270                        break;
    267271                case SMBpddEnum:

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r26526 r26550  
    37783778        IssmDouble teValue=1.0;
    37793779        IssmDouble aValue=0.0;
     3780        IssmDouble dulwrfValue=0.0;
    37803781        IssmDouble szaValue=0.0;
    37813782        IssmDouble cotValue=0.0;
    … …  
    37843785        IssmDouble dt,time,smb_dt;
    37853786        int        aIdx=0;
     3787        int        eIdx=0;
    37863788        int        denIdx=0;
    37873789        int        dsnowIdx=0;
    … …  
    38133815        bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux;
    38143816        bool isconstrainsurfaceT=false;
     3817        bool isdeltaLWup=false;
    38153818        IssmDouble init_scaling=0.0;
    38163819        IssmDouble thermo_scaling=1.0;
    38173820        IssmDouble adThresh=1023.0;
     3821        IssmDouble teThresh=10;
    38183822        /*}}}*/
    38193823        /*Output variables:{{{ */
    … …  
    38653869        parameters->FindParam(&smb_dt,SmbDtEnum);                     /*time period for the smb solution,  usually smaller than the glaciological dt*/
    38663870        parameters->FindParam(&aIdx,SmbAIdxEnum);
     3871        parameters->FindParam(&eIdx,SmbEIdxEnum);
    38673872        parameters->FindParam(&denIdx,SmbDenIdxEnum);
    38683873        parameters->FindParam(&swIdx,SmbSwIdxEnum);
    … …  
    38813886        parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum);
    38823887        parameters->FindParam(&isconstrainsurfaceT,SmbIsconstrainsurfaceTEnum);
     3888        parameters->FindParam(&isdeltaLWup,SmbIsdeltaLWupEnum);
    38833889        parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum);
    38843890        parameters->FindParam(&thermo_scaling,SmbThermoDeltaTScalingEnum);
    38853891        parameters->FindParam(&adThresh,SmbAdThreshEnum);
     3892        parameters->FindParam(&teThresh,SmbTeThreshEnum);
    38863893        /*}}}*/
    38873894        /*Retrieve inputs: {{{*/
    … …  
    40734080        Input *teValue_input= this->GetInput(SmbTeValueEnum,timeinputs); _assert_(teValue_input);
    40744081        Input *aValue_input= this->GetInput(SmbAValueEnum,timeinputs); _assert_(aValue_input);
     4082        Input *dulwrfValue_input= this->GetInput(SmbDulwrfValueEnum,timeinputs); _assert_(dulwrfValue_input);
    40754083        Input *szaValue_input= this->GetInput(SmbSzaValueEnum,timeinputs); _assert_(szaValue_input);
    40764084        Input *cotValue_input= this->GetInput(SmbCotValueEnum,timeinputs); _assert_(cotValue_input);
    … …  
    40884096        pAir_input->GetInputValue(&pAir,gauss);  // screen level air pressure [Pa]
    40894097        teValue_input->GetInputValue(&teValue,gauss);  // Emissivity [0-1]
     4098        dulwrfValue_input->GetInputValue(&dulwrfValue,gauss);  // LWup perturbation [W m-2]
    40904099        aValue_input->GetInputValue(&aValue,gauss);  // Albedo [0 1]
    40914100        szaValue_input->GetInputValue(&szaValue,gauss);  // Solar Zenith Angle [degree]
    … …  
    41134122        }
    41144123        /*Thermal profile computation:*/
    4115         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);
     4124        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);
    41164125
    41174126        /*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell.

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

    r26468 r26550  
    694694
    695695}  /*}}}*/
    696 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) { /*{{{*/
     696void 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) { /*{{{*/
    697697
    698698        /* ENGLACIAL THERMODYNAMICS*/
    … …  
    10101010
    10111011                // upward longwave contribution
    1012                 //ulw = - (SB * pow(Ts,4.0)* teValue) * dt; // - deltatest here
    1013                 IssmDouble deltatest=0;
    1014                 ulw = - (SB * pow(Ts,4.0)* teValue - deltatest) * dt; // - deltatest here
     1012                IssmDouble deltaULW=0.0;
     1013                IssmDouble emissivity=1.0;
     1014                //If user wants to set a upward long wave bias
     1015                if(isdeltaLWup) deltaULW = dulwrfValue;
     1016                //If user wants to directly set emissivity, or grain radius is larger than the
     1017                // threshold, or eIdx is 2 and we have wet snow or ice, use prescribed emissivity
     1018                if(eIdx==0 || (teThresh - re[0])<Gdntol || (eIdx==2 && z0>0.001)) emissivity = teValue;
     1019                ulw = - (SB * pow(Ts,4.0)* emissivity + deltaULW) * dt;
    10151020                ulwrf = ulwrf - ulw/dt0;
    10161021
    … …  
    22022207                                        if (fabs(adThresh - 820) < Dtol){
    22032208                                                // ERA5 new aIdx=1, swIdx=0, MODIS 820
    2204                                                 M0 = max(1.8230 - (0.1645 * log(C)),0.25);
    2205                                                 M1 = max(2.5134 - (0.3244 * log(C)),0.25);
     2209                                                //M0 = max(1.8230 - (0.1645 * log(C)),0.25);
     2210                                                //M1 = max(2.5134 - (0.3244 * log(C)),0.25);
     2211                                                // ERA5 new aIdx=1, swIdx=0, MODIS 820, p90 new (ERA5, 40 and e97)
     2212                                                M0 = max(1.3045 - (0.0988 * log(C)),0.25);
     2213                                                M1 = max(1.3694 - (0.1354 * log(C)),0.25);
    22062214                                        }
    22072215                                        else{
    … …  
    22462254                                        if (fabs(adThresh - 820) < Dtol){
    22472255                                                // ERA5 new aIdx=1, swIdx=0, MODIS 820
    2248                                                 M0 = max(1.4174 - (0.1037 * log(C)),0.25);
    2249                                                 M1 = max(2.2010 - (0.2460 * log(C)),0.25);
     2256                                                //M0 = max(1.4174 - (0.1037 * log(C)),0.25);
     2257                                                //M1 = max(2.2010 - (0.2460 * log(C)),0.25);
     2258                                                // ERA5 new aIdx=1, swIdx=0, MODIS 820, p90 new (ERA5 lwt, 40 and e97)
     2259                                                M0 = max(1.2138 - (0.1057 * log(C)),0.25);
     2260                                                M1 = max(1.4946 - (0.1607 * log(C)),0.25);
    22502261                                        }
    22512262                                        else{

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h

    r26477 r26550  
    3434void 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);
    3535void 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);
    36 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);
     36void 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);
    3737void 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);
    3838void 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);

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r26526 r26550  
    443443syn keyword cConstant SmbDtEnum
    444444syn keyword cConstant SmbEnum
     445syn keyword cConstant SmbEIdxEnum
    445446syn keyword cConstant SmbFEnum
    446447syn keyword cConstant SmbInitDensityScalingEnum
    … …  
    451452syn keyword cConstant SmbIsdelta18oEnum
    452453syn keyword cConstant SmbIsdensificationEnum
     454syn keyword cConstant SmbIsdeltaLWupEnum
    453455syn keyword cConstant SmbIsfirnwarmingEnum
    454456syn keyword cConstant SmbIsgraingrowthEnum
    … …  
    478480syn keyword cConstant SmbT0dryEnum
    479481syn keyword cConstant SmbT0wetEnum
     482syn keyword cConstant SmbTeThreshEnum
    480483syn keyword cConstant SmbTdiffEnum
    481484syn keyword cConstant SmbThermoDeltaTScalingEnum
    … …  
    913916syn keyword cConstant SmbDiniEnum
    914917syn keyword cConstant SmbDlwrfEnum
     918syn keyword cConstant SmbDulwrfValueEnum
    915919syn keyword cConstant SmbDswrfEnum
    916920syn keyword cConstant SmbDswdiffrfEnum

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r26526 r26550  
    437437        SmbDtEnum,
    438438        SmbEnum,
     439        SmbEIdxEnum,
    439440        SmbFEnum,
    440441        SmbInitDensityScalingEnum,
    … …  
    445446        SmbIsdelta18oEnum,
    446447        SmbIsdensificationEnum,
     448        SmbIsdeltaLWupEnum,
    447449        SmbIsfirnwarmingEnum,
    448450        SmbIsgraingrowthEnum,
    … …  
    472474        SmbT0dryEnum,
    473475        SmbT0wetEnum,
     476        SmbTeThreshEnum,
    474477        SmbTdiffEnum,
    475478        SmbThermoDeltaTScalingEnum,
    … …  
    909912        SmbDiniEnum,
    910913        SmbDlwrfEnum,
     914        SmbDulwrfValueEnum,
    911915        SmbDswrfEnum,
    912916        SmbDswdiffrfEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r26526 r26550  
    445445                case SmbDtEnum : return "SmbDt";
    446446                case SmbEnum : return "Smb";
     447                case SmbEIdxEnum : return "SmbEIdx";
    447448                case SmbFEnum : return "SmbF";
    448449                case SmbInitDensityScalingEnum : return "SmbInitDensityScaling";
    … …  
    453454                case SmbIsdelta18oEnum : return "SmbIsdelta18o";
    454455                case SmbIsdensificationEnum : return "SmbIsdensification";
     456                case SmbIsdeltaLWupEnum : return "SmbIsdeltaLWup";
    455457                case SmbIsfirnwarmingEnum : return "SmbIsfirnwarming";
    456458                case SmbIsgraingrowthEnum : return "SmbIsgraingrowth";
    … …  
    480482                case SmbT0dryEnum : return "SmbT0dry";
    481483                case SmbT0wetEnum : return "SmbT0wet";
     484                case SmbTeThreshEnum : return "SmbTeThresh";
    482485                case SmbTdiffEnum : return "SmbTdiff";
    483486                case SmbThermoDeltaTScalingEnum : return "SmbThermoDeltaTScaling";
    … …  
    915918                case SmbDiniEnum : return "SmbDini";
    916919                case SmbDlwrfEnum : return "SmbDlwrf";
     920                case SmbDulwrfValueEnum : return "SmbDulwrfValue";
    917921                case SmbDswrfEnum : return "SmbDswrf";
    918922                case SmbDswdiffrfEnum : return "SmbDswdiffrf";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r26526 r26550  
    454454              else if (strcmp(name,"SmbDt")==0) return SmbDtEnum;
    455455              else if (strcmp(name,"Smb")==0) return SmbEnum;
     456              else if (strcmp(name,"SmbEIdx")==0) return SmbEIdxEnum;
    456457              else if (strcmp(name,"SmbF")==0) return SmbFEnum;
    457458              else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum;
    … …  
    462463              else if (strcmp(name,"SmbIsdelta18o")==0) return SmbIsdelta18oEnum;
    463464              else if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum;
     465              else if (strcmp(name,"SmbIsdeltaLWup")==0) return SmbIsdeltaLWupEnum;
    464466              else if (strcmp(name,"SmbIsfirnwarming")==0) return SmbIsfirnwarmingEnum;
    465467              else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum;
    … …  
    489491              else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum;
    490492              else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum;
     493              else if (strcmp(name,"SmbTeThresh")==0) return SmbTeThreshEnum;
    491494              else if (strcmp(name,"SmbTdiff")==0) return SmbTdiffEnum;
    492495              else if (strcmp(name,"SmbThermoDeltaTScaling")==0) return SmbThermoDeltaTScalingEnum;
    … …  
    503506              else if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum;
    504507              else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum;
    505               else if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum;
    506               else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum;
    507               else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum;
     511              if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum;
     512              else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum;
     513              else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum;
     514              else if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum;
    512515              else if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum;
    513516              else if (strcmp(name,"StressbalanceRequestedOutputs")==0) return StressbalanceRequestedOutputsEnum;
    … …  
    626629              else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
    627630              else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
    628               else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
    629               else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum;
    630               else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
     634              if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
     635              else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum;
     636              else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
     637              else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
    635638              else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
    636639              else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
    … …  
    749752              else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
    750753              else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
    751               else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum;
    752               else if (strcmp(name,"HydrologyTws")==0) return HydrologyTwsEnum;
    753               else if (strcmp(name,"HydrologyTwsSpc")==0) return HydrologyTwsSpcEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"HydrologyTwsAnalysis")==0) return HydrologyTwsAnalysisEnum;
     757              if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum;
     758              else if (strcmp(name,"HydrologyTws")==0) return HydrologyTwsEnum;
     759              else if (strcmp(name,"HydrologyTwsSpc")==0) return HydrologyTwsSpcEnum;
     760              else if (strcmp(name,"HydrologyTwsAnalysis")==0) return HydrologyTwsAnalysisEnum;
    758761              else if (strcmp(name,"HydrologyWatercolumnMax")==0) return HydrologyWatercolumnMaxEnum;
    759762              else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
    … …  
    872875              else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
    873876              else if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum;
    874               else if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum;
    875               else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum;
    876               else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum;
     880              if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum;
     881              else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum;
     882              else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum;
     883              else if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum;
    881884              else if (strcmp(name,"SealevelchangeGUsubelOcean")==0) return SealevelchangeGUsubelOceanEnum;
    882885              else if (strcmp(name,"SealevelchangeGEsubelOcean")==0) return SealevelchangeGEsubelOceanEnum;
    … …  
    936939              else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum;
    937940              else if (strcmp(name,"SmbDlwrf")==0) return SmbDlwrfEnum;
     941              else if (strcmp(name,"SmbDulwrfValue")==0) return SmbDulwrfValueEnum;
    938942              else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum;
    939943              else if (strcmp(name,"SmbDswdiffrf")==0) return SmbDswdiffrfEnum;
    … …  
    994998              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    995999              else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
    996               else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
    9971004              else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
    9981005              else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
    9991006              else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
     1007              else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
    10041008              else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum;
    10051009              else if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
    … …  
    11171121              else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
    11181122              else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum;
    1119               else if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
    11201127              else if (strcmp(name,"Outputdefinition30")==0) return Outputdefinition30Enum;
    11211128              else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum;
    11221129              else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
     1130              else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
    11271131              else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
    11281132              else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
    … …  
    12401244              else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum;
    12411245              else if (strcmp(name,"Channel")==0) return ChannelEnum;
    1242               else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
    12431250              else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum;
    12441251              else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum;
    12451252              else if (strcmp(name,"Closed")==0) return ClosedEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"Colinear")==0) return ColinearEnum;
     1253              else if (strcmp(name,"Colinear")==0) return ColinearEnum;
    12501254              else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
    12511255              else if (strcmp(name,"Contact")==0) return ContactEnum;
    … …  
    13631367              else if (strcmp(name,"Intersect")==0) return IntersectEnum;
    13641368              else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
    1365               else if (strcmp(name,"J")==0) return JEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"J")==0) return JEnum;
    13661373              else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum;
    13671374              else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum;
    13681375              else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
     1376              else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
    13731377              else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
    13741378              else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
    … …  
    14861490              else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
    14871491              else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
    1488               else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
    14891496              else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
    14901497              else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;
    14911498              else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;
    1492          else stage=13;
    1493    }
    1494    if(stage==13){
    1495               if (strcmp(name,"SMBpdd")==0) return SMBpddEnum;
     1499              else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum;
    14961500              else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum;
    14971501              else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum;

  • issm/trunk-jpl/src/m/classes/SMBgemb.m

    r26358 r26550  
    2424                isconstrainsurfaceT;
    2525                isclimatology;
     26                isdeltaLWup;
    2627
    2728                %inputs:
    … …  
    4243                %optional inputs:
    4344                aValue = NaN; %Albedo forcing at every element.  Used only if aIdx == 0, or density exceeds adThresh
    44                 teValue = NaN; %Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)
     45                teValue = NaN; %Outward longwave radiation thermal emissivity forcing at every element (default in code is 1).
     46                               %Used only if eIdx== 0, or effective grain radius exceeds teThresh
     47                dulwrfValue = NaN; %Delta with which to perturn the long wave radiation upwards. Use if isdeltaLWup is true; 
    4548
    4649                % Initialization of snow properties
    … …  
    5962                %settings:
    6063                aIdx   = NaN; %method for calculating albedo and subsurface absorption (default is 1)
    61                 % 0: direct input from aValue parameter
     64                % 0: direct input from aValue parameter, no use of adThresh
    6265                % 1: effective grain radius [Gardner & Sharp, 2009]
    6366                % 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)
    6467                % 3: density and cloud amount [Greuell & Konzelmann, 1994]
    6568                % 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
     69
     70                eIdx   = NaN; %method for calculating emissivity (default is 1)
     71                % 0: direct input from teValue parameter, no use of teThresh
     72                % 1: default value of 1, in areas with grain radius below teThresh
     73                % 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
    6674
    6775                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))
    … …  
    111119                %or else apply direct input value from aValue, allowing albedo to be altered.
    112120                %Default value is rho water (1023 kg m-3).
     121                teThresh = NaN; %Apply eIdx method to all areas with grain radii below this value,
     122                %or else apply direct input value from teValue, allowing emissivity to be altered.
     123                %Default value is a effective grain radius of 10 mm.
    113124
    114125                %densities:
    … …  
    152163                        fielddisplay(self,'isdensification','run densification module (default true)');
    153164                        fielddisplay(self,'isturbulentflux','run turbulant heat fluxes module (default true)');
    154                         fielddisplay(self,'isconstrainsurfaceT','constrain surface temperatures to air temperature, turn off EC and surface flux contribution to surface temperature change');
     165                        fielddisplay(self,'isconstrainsurfaceT','constrain surface temperatures to air temperature, turn off EC and surface flux contribution to surface temperature change (default false)');
     166                        fielddisplay(self,'isdeltaLWup','set to true to invoke a bias in the long wave upward spatially, specified by dulwrfValue (default false)');
    155167                        fielddisplay(self,'Ta','2 m air temperature, in Kelvin');
    156168                        fielddisplay(self,'V','wind speed (m s-1)');
    … …  
    183195                                '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'})
    184196
     197                        fielddisplay(self,'dulwrfValue','Specified bias to be applied to the outward long wave radiation every element (W/m-2, +upward)');
    185198                        fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)');
     199                        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.'});
     200                        fielddisplay(self,'eIdx',{'method for calculating emissivity (default is 1)',...
     201                                '0: direct input from teValue parameter, no use of teThresh',...
     202                                '1: default value of 1, in areas with grain radius below teThresh',...
     203                                '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'});
    186204
    187205                        %snow properties init
    … …  
    299317                                self.cciceValue=project3d(md,'vector',self.cciceValue,'type','element');
    300318                        end
    301                         if (aIdx == 0) & ~isnan(self.aValue)
     319                        if ~isnan(self.aValue)
    302320                                self.aValue=project3d(md,'vector',self.aValue,'type','element');
    303321                        end
    … …  
    305323                                self.teValue=project3d(md,'vector',self.teValue,'type','element');
    306324                        end
    307 
    308325
    309326                end % }}}
    … …  
    322339                        self.isturbulentflux=1;
    323340                        self.isconstrainsurfaceT=0;
     341                        self.isdeltaLWup=0;
    324342
    325343                        self.aIdx = 1;
     344                        self.eIdx = 1;
    326345                        self.swIdx = 1;
    327346                        self.denIdx = 2;
    … …  
    348367                        self.K = 7;
    349368                        self.adThresh = 1023;
     369                        self.teThresh = 10;
    350370
    351371                        self.teValue = ones(mesh.numberofelements,1);
    352372                        self.aValue = self.aSnow*ones(mesh.numberofelements,1);
     373                        self.dulwrfValue = zeros(mesh.numberofelements,1);
    353374
    354375                        self.dswdiffrf=0.0*ones(mesh.numberofelements,1);
    … …  
    384405                        md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]);
    385406                        md = checkfield(md,'fieldname','smb.isconstrainsurfaceT','values',[0 1]);
     407                        md = checkfield(md,'fieldname','smb.isdeltaLWup','values',[0 1]);
    386408
    387409                        md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value
    … …  
    400422
    401423                        md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
     424                        md = checkfield(md,'fieldname','smb.dulwrfValue','timeseries',1,'NaN',1,'Inf',1);
    402425
    403426                        md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4]);
     427                        md = checkfield(md,'fieldname','smb.eIdx','NaN',1,'Inf',1,'values',[0,1,2]);
    404428                        md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]);
    405429                        md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5,6,7]);
    … …  
    414438                        md = checkfield(md,'fieldname','smb.ThermoDeltaTScaling','NaN',1,'Inf',1,'>=',0,'<=',1);
    415439                        md = checkfield(md,'fieldname','smb.adThresh','NaN',1,'Inf',1,'>=',0);
    416 
     440                        md = checkfield(md,'fieldname','smb.teThresh','NaN',1,'Inf',1,'>=',0);
     441
     442                        md = checkfield(md,'fieldname','smb.aValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
    417443                        switch self.aIdx,
    418                                 case 0
    419                                         md = checkfield(md,'fieldname','smb.aValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1);
    420444                                case {1 2}
    421445                                        md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'>=',.64,'<=',.89);
    … …  
    460484                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','isturbulentflux','format','Boolean');
    461485                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','isconstrainsurfaceT','format','Boolean');
     486                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','isdeltaLWup','format','Boolean');
    462487
    463488                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    … …  
    483508
    484509                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','aIdx','format','Integer');
     510                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','eIdx','format','Integer');
    485511                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','swIdx','format','Integer');
    486512                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','denIdx','format','Integer');
    … …  
    496522                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double');
    497523                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','adThresh','format','Double');
     524                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','teThresh','format','Double');
    498525
    499526                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','aValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    500527                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','teValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
     528                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','dulwrfValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    501529                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','szaValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    502530                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','cotValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    … …  
    533561                                error('If GEMB forcing dswdiffrf is transient, it must have the same time steps as input Ta in the final row!');
    534562                        end
     563                        if size(md.smb.dulwrfValue,2)>1 & any(md.smb.dulwrfValue(end,:) - md.smb.Ta(end,:) ~= 0)
     564                                error('If GEMB forcing dulwrfValue is transient, it must have the same time steps as input Ta in the final row!');
     565                        end
    535566                        if size(md.smb.aValue,2)>1 & any(md.smb.aValue(end,:) - md.smb.Ta(end,:) ~= 0)
    536567                                error('If GEMB forcing aValue is transient, it must have the same time steps as input Ta in the final row!');

  • issm/trunk-jpl/src/m/classes/SMBgemb.py

    r26359 r26550  
    3131        self.isturbulentflux = 0
    3232        self.isconstrainsurfaceT = 0
     33        self.isdeltaLWup = 0
    3334        self.isclimatology = np.nan
    3435
    … …  
    5051        # Optional inputs
    5152        self.aValue                 = np.nan    # Albedo forcing at every element. Used only if aIdx == 0, or density exceeds adThresh.
    52         self.teValue                = np.nan    # Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)
     53        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
     54        dulwrfValue                 = np.nan    #Delta with which to perturn the long wave radiation upwards. Use if isdeltaLWup is true
    5355
    5456        # Initialization of snow properties
    … …  
    6769        # Settings
    6870        self.aIdx                   = np.nan    # method for calculating albedo and subsurface absorption (default is 1)
    69         # 0: direct input from aValue parameter
     71        # 0: direct input from aValue parameter, no use of adThresh
    7072        # 1: effective grain radius [Gardner & Sharp, 2009]
    7173        # 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)
    7274        # 3: density and cloud amount [Greuell & Konzelmann, 1994]
    7375        # 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
     76       
     77        self.eIdx                   = np.nan    #method for calculating emissivity (default is 1)
     78        # 0: direct input from teValue parameter, no use of teThresh
     79        # 1: default value of 1, in areas with grain radius below teThresh
     80        # 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
    7481
    7582        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))
    … …  
    119126        # or else apply direct input value from aValue, allowing albedo to be altered.
    120127        # Default value is rho water (1023 kg m-3).
     128        teThresh                    = np.nan    #Apply eIdx method to all areas with grain radii below this value,
     129        #or else apply direct input value from teValue, allowing emissivity to be altered.
     130        #Default value is a effective grain radius of 10 mm.
    121131
    122132        # Densities
    … …  
    158168        s += '{}\n'.format(fielddisplay(self, 'isdensification', 'run densification module (default true)'))
    159169        s += '{}\n'.format(fielddisplay(self, 'isturbulentflux', 'run turbulant heat fluxes module (default true)'))
    160         s += '{}\n'.format(fielddisplay(self, 'isconstrainsurfaceT', 'constrain surface temperatures to air temperature, turn off EC and surface flux contribution to surface temperature change'))
     170        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)'))
     171        s += '{}\n'.format(fielddisplay(self, 'isdeltaLWup', 'set to true to invoke a bias in the long wave upward spatially, specified by dulwrfValue (default false)'))
    161172        s += '{}\n'.format(fielddisplay(self, 'Ta', '2 m air temperature, in Kelvin'))
    162173        s += '{}\n'.format(fielddisplay(self, 'V', 'wind speed (m s-1)'))
    … …  
    188199            '3: density and cloud amount [Greuell & Konzelmann, 1994]',
    189200            '4: exponential time decay & wetness [Bougamont & Bamber, 2005]']))
     201
     202        s += '{}\n'.format(fielddisplay(self, 'dulwrfValue', 'Specified bias to be applied to the outward long wave radiation every element (W/m-2, +upward)'))
    190203        s += '{}\n'.format(fielddisplay(self, 'teValue', 'Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)'))
     204        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.']))
     205        s += '{}\n'.format(fielddisplay(self, 'eIdx', ['method for calculating emissivity (default is 1)',
     206            '0: direct input from teValue parameter, no use of teThresh',
     207            '1: default value of 1, in areas with grain radius below teThresh',
     208            '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']))
    191209
    192210        # Snow properties init
    … …  
    286304            self.cciceValue=project3d(md,'vector',self.cciceValue,'type','element');
    287305
    288         if (self.aIdx == 0) and (not np.isnan(self.aValue)):
     306        if not np.isnan(self.aValue):
    289307            self.aValue = project3d(md, 'vector', self.aValue, 'type', 'element')
    290308        if not np.isnan(self.teValue):
    … …  
    308326        self.isturbulentflux = 1
    309327        self.isconstrainsurfaceT = 0
     328        self.isdeltaLWup = 0
    310329
    311330        self.aIdx = 1
     331        self.eIdx = 1
    312332        self.swIdx = 1
    313333        self.denIdx = 2
    … …  
    334354        self.K = 7
    335355        self.adThresh = 1023
     356        self.teThresh = 10
    336357
    337358        self.teValue = np.ones((mesh.numberofelements,))
    338359        self.aValue = self.aSnow * np.ones(mesh.numberofelements,)
     360        self.dulwrfValue = np.zeros((mesh.numberofelements,))
    339361
    340362        self.dswdiffrf = 0.0 * np.ones(mesh.numberofelements,)
    … …  
    373395        md = checkfield(md, 'fieldname', 'smb.isdensification', 'values', [0, 1])
    374396        md = checkfield(md, 'fieldname', 'smb.isturbulentflux', 'values', [0, 1])
     397        md = checkfield(md, 'fieldname', 'smb.isdeltaLWup', 'values',[0, 1])
    375398        md = checkfield(md, 'fieldname', 'smb.isconstrainsurfaceT', 'values', [0, 1])
    376399
    … …  
    390413
    391414        md = checkfield(md, 'fieldname', 'smb.teValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
     415        md = checkfield(md, 'fieldname', 'smb.dulwrfValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    392416
    393417        md = checkfield(md, 'fieldname', 'smb.aIdx', 'NaN', 1, 'Inf', 1, 'values', [0, 1, 2, 3, 4])
     418        md = checkfield(md, 'fieldname', 'smb.eIdx', 'NaN', 1, 'Inf', 1, 'values', [0, 1, 2])
    394419        md = checkfield(md, 'fieldname', 'smb.swIdx', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
    395420        md = checkfield(md, 'fieldname', 'smb.denIdx', 'NaN', 1, 'Inf', 1, 'values', [1, 2, 3, 4, 5, 6, 7])
    … …  
    404429        md = checkfield(md, 'fieldname', 'smb.ThermoDeltaTScaling', 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1)
    405430        md = checkfield(md, 'fieldname', 'smb.adThresh', 'NaN', 1, 'Inf', 1, '>=', 0)
    406 
    407         if self.aIdx == 0:
    408             md = checkfield(md, 'fieldname', 'smb.aValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
    409         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)):
     431        md = checkfield(md, 'fieldname', 'smb.teThresh', 'NaN', 1, 'Inf',1,'>=',0)
     432       
     433        md = checkfield(md, 'fieldname', 'smb.aValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
     434        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)):
    410435            md = checkfield(md, 'fieldname', 'smb.aSnow', 'NaN', 1, 'Inf', 1, '> = ', .64, '< = ', .89)
    411436            md = checkfield(md, 'fieldname', 'smb.aIce', 'NaN', 1, 'Inf', 1, '> = ', .27, '< = ', .58)
    … …  
    446471        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isturbulentflux', 'format', 'Boolean')
    447472        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isconstrainsurfaceT', 'format', 'Boolean')
     473        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isdeltaLWup','format','Boolean')
    448474
    449475        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Ta', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
    … …  
    469495
    470496        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'aIdx', 'format', 'Integer')
     497        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'eIdx', 'format', 'Integer')
    471498        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'swIdx', 'format', 'Integer')
    472499        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'denIdx', 'format', 'Integer')
    … …  
    482509        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'K', 'format', 'Double')
    483510        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'adThresh', 'format', 'Double')
     511        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'teThresh', 'format', 'Double')
    484512
    485513        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'aValue', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
    486514        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'teValue', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
     515        WriteData(fid,prefix,'object',self,'class','smb','fieldname','dulwrfValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
    487516        WriteData(fid,prefix,'object',self,'class','smb','fieldname','szaValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
    488517        WriteData(fid,prefix,'object',self,'class','smb','fieldname','cotValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
    … …  
    514543        if ((np.ndim(self.aValue)>1) & np.any(self.aValue[-1] - self.Ta[-1] != 0)):
    515544            raise IOError('If GEMB forcing aValue is transient, it must have the same time steps as input Ta in the final row!')
     545        if ((np.ndim(self.dulwrfValue)>1) & np.any(self.dulwrfValue[-1] - self.Ta[-1] != 0)):
     546            raise IOError('If GEMB forcing dulwrfValue is transient, it must have the same time steps as input Ta in the final row!')
    516547        if ((np.ndim(self.szaValue)>1) & np.any(self.szaValue[-1] - self.Ta[-1] != 0)):
    517548            raise IOError('If GEMB forcing szaValue is transient, it must have the same time steps as input Ta in the final row!')
Note: See TracChangeset for help on using the changeset viewer.