Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22538) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22539) @@ -2816,4 +2816,5 @@ IssmDouble sumMassAdd=0.0; IssmDouble sumdz_add=0.0; + IssmDouble fac=0.0; IssmDouble sumMass=0.0; IssmDouble dMass=0.0; @@ -3117,8 +3118,12 @@ sumP = P + sumP; sumEC = sumEC + EC; // evap (-)/cond(+) - sumdz_add=dz_add+sumdz_add; /*Calculate total system mass:*/ - sumMass=0; for(int i=0;iAddInput(new DoubleInput(SmbRunoffEnum,sumR/dt/rho_ice)); this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/dt/rho_ice)); - this->AddInput(new DoubleInput(SmbDz_addEnum,sumdz_add/yts)); - this->AddInput(new DoubleInput(SmbM_addEnum,sumMassAdd/dt)); + this->AddInput(new DoubleInput(SmbDzAddEnum,sumdz_add)); + this->AddInput(new DoubleInput(SmbMAddEnum,sumMassAdd/dt)); + this->AddInput(new DoubleInput(SmbFACEnum,fac/1000)); // output in meters /*Free allocations:{{{*/ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 22538) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 22539) @@ -583,10 +583,7 @@ /*outputs:*/ IssmDouble EC=0.0; - IssmDouble* T=NULL; + IssmDouble* T=*pT; if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" thermal module\n"); - - /*Recover pointers: */ - T=*pT; ds = d[0]; // density of top grid cell @@ -686,5 +683,5 @@ } dt=max_fdt; - + // determine mean (harmonic mean) of K/dz for u, d, & p Au = xNew(m); @@ -783,5 +780,5 @@ // calculate heat/wind 'coef_H' stability factor - if (Ri < -0.03-Ttol) coefH = 1.3 * coefM; + if (Ri <= -0.03+Ttol) coefH = coefM/1.3; else coefH = coefM; @@ -790,5 +787,5 @@ shf = C * CA * (Ta - Ts); - // adjust using Monin–Obukhov stability theory + // adjust using Monin-Obukhov stability theory shf = shf / (coefM * coefH); @@ -796,18 +793,20 @@ // determine if snow pack is melting & calcualte surface vapour pressure over ice or liquid water if (Ts >= CtoK-Ttol){ - L = LV; - + L = LV; //for liquid water at 273.15 k to vapor + + //for liquid surface (assume liquid on surface when Ts == 0 deg C) + // Wright (1997), US Meteorological Handbook from Murphy and Koop, 2005 Appendix A + eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK)); + } + else{ + L = LS; // latent heat of sublimation + // for an ice surface Murphy and Koop, 2005 [Equation 7] eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts); } - else{ - L = LS; // latent heat of sublimation for liquid surface (assume liquid on surface when Ts == 0 deg C) - // Wright (1997), US Meteorological Handbook from Murphy and Koop, 2005 Appendix A - eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK)); - } - + // Latent heat flux [W m-2] lhf = C * L * (eAir - eS) * 0.622 / pAir; - + // adjust using Monin-Obukhov stability theory (if lhf '+' then there is energy and mass gained at the surface, // if '-' then there is mass and energy loss at the surface. @@ -832,5 +831,5 @@ // temperature diffusion - for(int j=0;j(K); @@ -1680,5 +1679,5 @@ // WORKS FINE BUT HAS BEEN DISABLED FOR CONVIENCE OF MODEL OUTPUT // INTERPRETATION - // // calculate total model depth + // calculate total model depth Ztot = cellsum(dz,n); @@ -1993,5 +1992,5 @@ // if V = 0 goes to infinity therfore if V = 0 change - if(V < 0.01) V=0.01; + if(V<0.01-Dtol)V=0.01; // calculate the Bulk Richardson Number (Ri) @@ -1999,16 +1998,21 @@ // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H' - + // do not allow Ri to exceed 0.19 - if(Ri>0.19-Ttol)Ri= 0.19; - - // calculate momentum 'coef_M' stability factor - if (Ri > 0.0) coef_M = pow(1.0-5.2*Ri,-1.0); // if stable - else coef_M = pow(1.0-18*Ri,-0.25); + Ri = fmin(Ri, 0.19); + + // calculate momentum 'coefM' stability factor + if (Ri > 0.0+Ttol){ + // if stable + coef_M = 1.0/(1.0-5.2*Ri); + } + else { + coef_M =pow (1.0-18*Ri,-0.25); + } // calculate heat/wind 'coef_H' stability factor - if (Ri < -0.03) coef_H = 1.3 * coef_M; + if (Ri <= -0.03+Ttol) coef_H = coef_M/1.3; else coef_H = coef_M; - + //// Bulk-transfer coefficient An = pow(0.4,2) / pow(log(Tz/z0),2); // Bulk-transfer coefficient @@ -2022,20 +2026,20 @@ shf = shf / (coef_M * coef_H); - //// Latent Heat - // determine if snow pack is melting & calcualte surface vapour pressure - // over ice or liquid water + // Latent Heat + // determine if snow pack is melting & calcualte surface vapour pressure over ice or liquid water if (Ts >= CtoK-Ttol){ - L = LV; - - // for an ice surface Murphy and Koop, 2005 [Equation 7] - eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts)- 0.00728332 * Ts); + L = LV; //for liquid water at 273.15 k to vapor + + //for liquid surface (assume liquid on surface when Ts == 0 deg C) + // Wright (1997), US Meteorological Handbook from Murphy and Koop, 2005 Appendix A + eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK)); } else{ L = LS; // latent heat of sublimation - // for liquid surface (assume liquid on surface when Ts == 0 deg C) - // Wright (1997), US Meteorological Handbook from Murphy and Koop, - // 2005 Apendix A - eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK)); - } + + // for an ice surface Murphy and Koop, 2005 [Equation 7] + eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts); + } + // Latent heat flux [W m-2] Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22538) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22539) @@ -484,6 +484,7 @@ SmbIsdensificationEnum, SmbIsturbulentfluxEnum, - SmbDz_addEnum, - SmbM_addEnum, + SmbDzAddEnum, + SmbMAddEnum, + SmbFACEnum, /*SMBpdd*/ SMBpddEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22538) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22539) @@ -486,6 +486,7 @@ case SmbIsdensificationEnum : return "SmbIsdensification"; case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux"; - case SmbDz_addEnum : return "SmbDz_add"; - case SmbM_addEnum : return "SmbM_add"; + case SmbDzAddEnum : return "SmbDzAdd"; + case SmbMAddEnum : return "SmbMAdd"; + case SmbFACEnum : return "SmbFAC"; case SMBpddEnum : return "SMBpdd"; case SmbDelta18oEnum : return "SmbDelta18o"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22538) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22539) @@ -495,6 +495,7 @@ else if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum; else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum; - else if (strcmp(name,"SmbDz_add")==0) return SmbDz_addEnum; - else if (strcmp(name,"SmbM_add")==0) return SmbM_addEnum; + else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum; + else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; + else if (strcmp(name,"SmbFAC")==0) return SmbFACEnum; else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"SmbIsd18opd")==0) return SmbIsd18opdEnum; else if (strcmp(name,"SmbIstemperaturescaled")==0) return SmbIstemperaturescaledEnum; - else if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum; + if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum; + else if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum; else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum; else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; - else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum; + if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; + else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum; else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; else if (strcmp(name,"LambdaS")==0) return LambdaSEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum; else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum; - else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum; else stage=7; } if(stage==7){ - if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum; + if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum; + else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum; else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum; else if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum; @@ -874,9 +875,9 @@ else if (strcmp(name,"SealevelriseRequestedOutputs")==0) return SealevelriseRequestedOutputsEnum; else if (strcmp(name,"SealevelriseNumRequestedOutputs")==0) return SealevelriseNumRequestedOutputsEnum; - else if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"LoveFrequencies")==0) return LoveFrequenciesEnum; + if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum; + else if (strcmp(name,"LoveFrequencies")==0) return LoveFrequenciesEnum; else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; else if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum; @@ -997,9 +998,9 @@ else if (strcmp(name,"Tria")==0) return TriaEnum; else if (strcmp(name,"TriaInput")==0) return TriaInputEnum; - else if (strcmp(name,"Tetra")==0) return TetraEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"TetraInput")==0) return TetraInputEnum; + if (strcmp(name,"Tetra")==0) return TetraEnum; + else if (strcmp(name,"TetraInput")==0) return TetraInputEnum; else if (strcmp(name,"Penta")==0) return PentaEnum; else if (strcmp(name,"PentaInput")==0) return PentaInputEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"P1xP3")==0) return P1xP3Enum; else if (strcmp(name,"P1xP4")==0) return P1xP4Enum; - else if (strcmp(name,"P2xP4")==0) return P2xP4Enum; else stage=10; } if(stage==10){ - if (strcmp(name,"P1P1")==0) return P1P1Enum; + if (strcmp(name,"P2xP4")==0) return P2xP4Enum; + else if (strcmp(name,"P1P1")==0) return P1P1Enum; else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum; else if (strcmp(name,"MINI")==0) return MINIEnum; Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 22538) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 22539) @@ -199,4 +199,8 @@ elseif strcmp(fieldname,'SmbRunoff'), field = field*yts; + elseif strcmp(fieldname,'SmbEvaporation'), + field = field*yts; + elseif strcmp(fieldname,'SmbRefreeze'), + field = field*yts; elseif strcmp(fieldname,'SmbEC'), field = field*yts; @@ -205,8 +209,6 @@ elseif strcmp(fieldname,'SmbMelt'), field = field*yts; - elseif strcmp(fieldname,'SmbDz_add'), - field = field*yts; - elseif strcmp(fieldname,'SmbM_add'), - field = field*yts; + elseif strcmp(fieldname,'SmbMAdd'), + field = field*yts; elseif strcmp(fieldname,'CalvingCalvingrate'), field = field*yts; Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 22538) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 22539) @@ -209,4 +209,8 @@ elif fieldname=='SmbRunoff': field = field*yts + elif fieldname=='SmbEvaporation': + field = field*yts; + elif fieldname=='SmbRefreeze': + field = field*yts; elif fieldname=='SmbEC': field = field*yts @@ -215,7 +219,5 @@ elif fieldname=='SmbMelt': field = field*yts - elif fieldname=='SmbDz_add': - field = field*yts - elif fieldname=='SmbM_add': + elif fieldname=='SmbMAdd': field = field*yts elif fieldname=='CalvingCalvingrate': Index: /issm/trunk-jpl/test/NightlyRun/test243.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test243.m (revision 22538) +++ /issm/trunk-jpl/test/NightlyRun/test243.m (revision 22539) @@ -27,5 +27,5 @@ %smb settings -md.smb.requested_outputs={'SmbDz','SmbT','SmbD','SmbRe','SmbGdn','SmbGsp','SmbEC','SmbA','SmbMassBalance'}; +md.smb.requested_outputs={'SmbDz','SmbT','SmbD','SmbRe','SmbGdn','SmbGsp','SmbEC','SmbA','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC'}; %only run smb core: @@ -44,6 +44,6 @@ %Fields and tolerances to track changes -field_names ={'SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance'}; -field_tolerances ={1e-11,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-12,1e-12}; +field_names ={'SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC'}; +field_tolerances ={1e-11,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12}; field_values={... @@ -57,3 +57,6 @@ (md.results.TransientSolution(end).SmbEC(1)),... (md.results.TransientSolution(end).SmbMassBalance(1)),... + (md.results.TransientSolution(end).SmbMAdd(1)),... + (md.results.TransientSolution(end).SmbDzAdd(1)),... + (md.results.TransientSolution(end).SmbFAC(1)) }; Index: /issm/trunk-jpl/test/NightlyRun/test243.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test243.py (revision 22538) +++ /issm/trunk-jpl/test/NightlyRun/test243.py (revision 22539) @@ -40,5 +40,5 @@ #smb settings -md.smb.requested_outputs = ['SmbDz','SmbT','SmbD','SmbRe','SmbGdn','SmbGsp','SmbEC','SmbA','SmbMassBalance'] +md.smb.requested_outputs = ['SmbDz','SmbT','SmbD','SmbRe','SmbGdn','SmbGsp','SmbEC','SmbA','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC'] #only run smb core: @@ -57,6 +57,6 @@ #Fields and tolerances to track changes -field_names = ['SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance'] -field_tolerances = [1e-11,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-12,1e-12] +field_names = ['SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC'] +field_tolerances = [1e-11,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12] #shape is different in python solution (fixed using reshape) which can cause test failure: field_values = [ @@ -70,3 +70,6 @@ md.results.TransientSolution[-1].SmbEC[0], md.results.TransientSolution[-1].SmbMassBalance[0], + md.results.TransientSolution[-1].SmbMAdd[0], + md.results.TransientSolution[-1].SmbDzAdd[0], + md.results.TransientSolution[-1].SmbFAC[0], ]