Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23467) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23468) @@ -57,4 +57,5 @@ iomodel->FetchDataToInput(elements,"md.smb.zMin",SmbZMinEnum); iomodel->FetchDataToInput(elements,"md.smb.Tmean",SmbTmeanEnum); + iomodel->FetchDataToInput(elements,"md.smb.Vmean",SmbVmeanEnum); iomodel->FetchDataToInput(elements,"md.smb.C",SmbCEnum); iomodel->FetchDataToInput(elements,"md.smb.Tz",SmbTzEnum); @@ -184,4 +185,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.smb.swIdx",SmbSwIdxEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.denIdx",SmbDenIdxEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.dsnowIdx",SmbDsnowIdxEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.cldFrac",SmbCldFracEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0wet",SmbT0wetEnum)); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23467) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23468) @@ -3026,4 +3026,5 @@ IssmDouble dzMin=0.0; IssmDouble Tmean=0.0; + IssmDouble Vmean=0.0; IssmDouble C=0.0; IssmDouble Tz,Vz=0.0; @@ -3043,4 +3044,5 @@ int aIdx=0; int denIdx=0; + int dsnowIdx=0; int swIdx=0; IssmDouble cldFrac,t0wet, t0dry, K; @@ -3113,4 +3115,5 @@ parameters->FindParam(&denIdx,SmbDenIdxEnum); parameters->FindParam(&swIdx,SmbSwIdxEnum); + parameters->FindParam(&dsnowIdx,SmbDsnowIdxEnum); parameters->FindParam(&cldFrac,SmbCldFracEnum); parameters->FindParam(&t0wet,SmbT0wetEnum); @@ -3138,4 +3141,5 @@ Input* zY_input=this->GetInput(SmbZYEnum); _assert_(zY_input); Input* Tmean_input=this->GetInput(SmbTmeanEnum); _assert_(Tmean_input); + Input* Vmean_input=this->GetInput(SmbVmeanEnum); _assert_(Vmean_input); Input* C_input=this->GetInput(SmbCEnum); _assert_(C_input); Input* Tz_input=this->GetInput(SmbTzEnum); _assert_(Tz_input); @@ -3161,4 +3165,5 @@ zY_input->GetInputValue(&zY,gauss); Tmean_input->GetInputValue(&Tmean,gauss); + Vmean_input->GetInputValue(&Vmean,gauss); C_input->GetInputValue(&C,gauss); Tz_input->GetInputValue(&Tz,gauss); @@ -3316,5 +3321,5 @@ /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/ - if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, Ta, P, dzMin, aSnow,rho_ice,this->Sid()); + if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, dsnowIdx, Tmean, Ta, P, dzMin, aSnow, C, V, Vmean, rho_ice,this->Sid()); /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K @@ -3327,5 +3332,5 @@ /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every * sub-time step in thermo equations*/ - ulw = 5.67E-8 * pow(T[0],4.0) * teValue; + ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here /*Calculate net longwave [W m-2]*/ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23467) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23468) @@ -820,5 +820,5 @@ // upward longwave contribution - ulw = - (SB * pow(Ts,4.0)* teValue) * dt ; + ulw = - (SB * pow(Ts,4.0)* teValue) * dt ; //+20 dT_ulw = ulw / TCs; @@ -1076,5 +1076,5 @@ } /*}}}*/ -void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid){ /*{{{*/ +void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid){ /*{{{*/ // Adds precipitation and deposition to the model grid @@ -1086,4 +1086,7 @@ adjusts the properties of the top grid cell to account for accumulation T_air & T = Air and top grid cell temperatures [K] + Tmean = average surface temperature [K] + V = wind velocity [m s-1] + C = average accumulation rate [kg m-2 yr-1] dz = topgrid cell length [m] d = density of top grid gell [kg m-3] @@ -1095,5 +1098,5 @@ // MAIN FUNCTION // specify constants - const IssmDouble dSnow = 150; // density of snow [kg m-3] + IssmDouble dSnow = 150; // density of snow [kg m-3] const IssmDouble reNew = 0.1; // new snow grain size [mm] const IssmDouble gdnNew = 1.0; // new snow dendricity @@ -1130,4 +1133,28 @@ gsp=*pgsp; m=*pm; + + //Density of fresh snow [kg m-3] + switch (dsnowIdx){ + case 0: // Default value defined above + break; + + case 1: // Density of Antarctica snow + dSnow = 350; + break; + + case 2: // Density of Greenland snow, Fausto et al., 2008 + dSnow = 315; + break; + + case 3: //Surface snow accumulation density from Kaspers et al., 2004, Antarctica + //dSnow = alpha1 + beta1*T + delta1*C + epsilon1*W + // 7.36x10-2 1.06x10-3 6.69x10-2 4.77x10-3 + dSnow=(7.362e-2 + 1.06e-3*Tmean + 6.69e-2*C/1000 + 4.77e-3*Vmean)*1000; + break; + + case 4: // Kuipers Munneke and others (2015), Greenland + dSnow = 481.0 + 4.834*(Tmean-CtoK); + break; + } // determine initial mass @@ -1915,6 +1942,8 @@ c0arth = 0.07 * H; c1arth = 0.03 * H; - M0 = max(1.435 - (0.151 * log(C)),0.25); - M1 = max(2.366 - (0.293 * log(C)),0.25); + M0 = max(1.6599 - (0.1724 * log(C)),0.25); + M1 = max(2.0102 - (0.2458 * log(C)),0.25); + //M0 = max(1.435 - (0.151 * log(C)),0.25); + //M1 = max(2.366 - (0.293 * log(C)),0.25); c0 = M0*c0arth; c1 = M1*c1arth; @@ -1926,6 +1955,8 @@ c0arth = 0.07 * H; c1arth = 0.03 * H; - M0 = max(1.042 - (0.0916 * log(C)),0.25); - M1 = max(1.734 - (0.2039 * log(C)),0.25); + //M0 = max(1.042 - (0.0916 * log(C)),0.25); + //M1 = max(1.734 - (0.2039 * log(C)),0.25); + M0 = max(1.6201 - (0.1450 * log(C)),0.25); + M1 = max(2.5577 - (0.2899 * log(C)),0.25); c0 = M0*c0arth; c1 = M1*c1arth; Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 23467) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 23468) @@ -29,5 +29,5 @@ void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid); void thermo(IssmDouble* pEC, IssmDouble** T, 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); -void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx,IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid); +void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, 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* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid); void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23467) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23468) @@ -279,4 +279,5 @@ SmbAIceEnum, SmbAIdxEnum, + SmbDsnowIdxEnum, SmbASnowEnum, SmbCldFracEnum, @@ -607,4 +608,5 @@ SmbTiniEnum, SmbTmeanEnum, + SmbVmeanEnum, SmbTzEnum, SmbVEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23467) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23468) @@ -287,4 +287,5 @@ case SmbAIceEnum : return "SmbAIce"; case SmbAIdxEnum : return "SmbAIdx"; + case SmbDsnowIdxEnum : return "SmbDsnowIdx"; case SmbASnowEnum : return "SmbASnow"; case SmbCldFracEnum : return "SmbCldFrac"; @@ -613,4 +614,5 @@ case SmbTiniEnum : return "SmbTini"; case SmbTmeanEnum : return "SmbTmean"; + case SmbVmeanEnum : return "SmbVmean"; case SmbTzEnum : return "SmbTz"; case SmbVEnum : return "SmbV"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23467) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23468) @@ -293,4 +293,5 @@ else if (strcmp(name,"SmbAIce")==0) return SmbAIceEnum; else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum; + else if (strcmp(name,"SmbDsnowIdx")==0) return SmbDsnowIdxEnum; else if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum; else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; - else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"TransientIsslr")==0) return TransientIsslrEnum; + if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; + else if (strcmp(name,"TransientIsslr")==0) return TransientIsslrEnum; else if (strcmp(name,"TransientIssmb")==0) return TransientIssmbEnum; else if (strcmp(name,"TransientIsstressbalance")==0) return TransientIsstressbalanceEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"HydrologyDrainageRate")==0) return HydrologyDrainageRateEnum; else if (strcmp(name,"Ice")==0) return IceEnum; - else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"Input")==0) return InputEnum; + if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; + else if (strcmp(name,"Input")==0) return InputEnum; else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; @@ -625,12 +626,13 @@ else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum; + else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum; else if (strcmp(name,"SmbTz")==0) return SmbTzEnum; else if (strcmp(name,"SmbV")==0) return SmbVEnum; - else if (strcmp(name,"SmbVz")==0) return SmbVzEnum; - else if (strcmp(name,"SmbW")==0) return SmbWEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"SmbWini")==0) return SmbWiniEnum; + if (strcmp(name,"SmbVz")==0) return SmbVzEnum; + else if (strcmp(name,"SmbW")==0) return SmbWEnum; + else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum; else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum; else if (strcmp(name,"SmbZMin")==0) return SmbZMinEnum; @@ -750,10 +752,10 @@ else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum; else if (strcmp(name,"DataSet")==0) return DataSetEnum; - else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; - else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; + if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; + else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; + else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum; else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum; else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum; - else if (strcmp(name,"LambdaS")==0) return LambdaSEnum; - else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; + if (strcmp(name,"LambdaS")==0) return LambdaSEnum; + else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; + else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum; @@ -996,10 +998,10 @@ else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum; - else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum; - else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum; else stage=9; } if(stage==9){ - if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum; + if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum; + else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum; + else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum; else if (strcmp(name,"Outputdefinition3")==0) return Outputdefinition3Enum; else if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"SigmaVM")==0) return SigmaVMEnum; else if (strcmp(name,"SmbAnalysis")==0) return SmbAnalysisEnum; - else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; - else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum; + if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; + else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; + else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum; else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum; Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 23467) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 23468) @@ -33,4 +33,5 @@ Tmean = NaN; %mean annual temperature [K] + Vmean = NaN; %mean annual wind velocity [m s-1] C = NaN; %mean annual snow accumulation [kg m-2 yr-1] Tz = NaN; %height above ground at which temperature (T) was sampled [m] @@ -72,4 +73,11 @@ % 7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015) + dsnowIdx = NaN; %model for fresh snow accumulation density (default is 1): + % 0 = Original GEMB value, 150 kg/m^3 + % 1 = Antarctica value of fresh snow density, 350 kg/m^3 + % 2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008) + % 3 = Antarctica model of Kaspers et al. (2004) + % 4 = Greenland model of Kuipers Munneke et al. (2015) + zTop = NaN; % depth over which grid length is constant at the top of the snopack (default 10) [m] dzTop = NaN; % initial top vertical grid spacing (default .05) [m] @@ -155,4 +163,5 @@ self.swIdx = 1; self.denIdx = 2; + self.dsnowIdx = 1; self.zTop=10*ones(mesh.numberofelements,1); self.dzTop = .05* ones (mesh.numberofelements,1); @@ -160,4 +169,6 @@ self.InitDensityScaling = 1.0; self.ThermoDeltaTScaling = 1/11; + + self.Vmean=10.0*ones(mesh.numberofelements,1); self.zMax=250*ones(mesh.numberofelements,1); @@ -212,4 +223,5 @@ md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0); + md = checkfield(md,'fieldname','smb.Vmean','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000); md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000); @@ -220,4 +232,5 @@ 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]); + md = checkfield(md,'fieldname','smb.dsnowIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4]); md = checkfield(md,'fieldname','smb.zTop','NaN',1,'Inf',1,'>=',0); @@ -265,5 +278,5 @@ fielddisplay(self,'isturbulentflux','run turbulant heat fluxes module (default true)'); fielddisplay(self,'Ta','2 m air temperature, in Kelvin'); - fielddisplay(self,'V','wind speed (m/s-1)'); + fielddisplay(self,'V','wind speed (m s-1)'); fielddisplay(self,'dlwrf','downward shortwave radiation flux [W/m^2]'); fielddisplay(self,'dswrf','downward longwave radiation flux [W/m^2]'); @@ -273,4 +286,5 @@ fielddisplay(self,'Tmean','mean annual temperature [K]'); fielddisplay(self,'C','mean annual snow accumulation [kg m-2 yr-1]'); + fielddisplay(self,'Vmean','mean annual snow accumulation [m s-1] (default 10 m/s)'); fielddisplay(self,'Tz','height above ground at which temperature (T) was sampled [m]'); fielddisplay(self,'Vz','height above ground at which wind (V) eas sampled [m]'); @@ -331,4 +345,11 @@ '6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)',... '7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)'}); + fielddisplay(self,'dsnowIdx',{'model for fresh snow accumulation density (default is 1):',... + '0 = Original GEMB value, 150 kg/m^3',... + '1 = Antarctica value of fresh snow density, 350 kg/m^3',... + '2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)',... + '3 = Antarctica model of Kaspers et al. (2004), Make sure to set Vmean accurately',... + '4 = Greenland model of Kuipers Munneke et al. (2015)'}); + fielddisplay(self,'requested_outputs','additional outputs requested'); @@ -360,4 +381,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tmean','format','DoubleMat','mattype',2); WriteData(fid,prefix,'object',self,'class','smb','fieldname','C','format','DoubleMat','mattype',2); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Vmean','format','DoubleMat','mattype',2); WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tz','format','DoubleMat','mattype',2); WriteData(fid,prefix,'object',self,'class','smb','fieldname','Vz','format','DoubleMat','mattype',2); @@ -372,4 +394,6 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','swIdx','format','Integer'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','denIdx','format','Integer'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dsnowIdx','format','Integer'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dsnowIdx','format','Integer'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','InitDensityScaling','format','Double'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','ThermoDeltaTScaling','format','Double'); Index: /issm/trunk-jpl/src/m/classes/SMBgemb.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 23467) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 23468) @@ -39,4 +39,5 @@ Tmean = float('NaN') #mean annual temperature [K] + Vmean = float('NaN') #mean annual wind velocity [m s-1] C = float('NaN') #mean annual snow accumulation [kg m-2 yr-1] Tz = float('NaN') #height above ground at which temperature (T) was sampled [m] @@ -78,4 +79,11 @@ # 6 = Antarctica semi-emperical model of Ligtenberg et al. (2011) # 7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015) + + dsnowIdx = float('NaN') #model for fresh snow accumulation density (default is 1): + # 0 = Original GEMB value, 150 kg/m^3 + # 1 = Antarctica value of fresh snow density, 350 kg/m^3 + # 2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008) + # 3 = Antarctica model of Kaspers et al. (2004) + # 4 = Greenland model of Kuipers Munneke et al. (2015) zTop = float('NaN') # depth over which grid length is constant at the top of the snopack (default 10) [m] @@ -131,5 +139,5 @@ string = "%s\n%s"%(string,fielddisplay(self,'isturbulentflux','run turbulant heat fluxes module (default true)')) string = "%s\n%s"%(string,fielddisplay(self,'Ta','2 m air temperature, in Kelvin')) - string = "%s\n%s"%(string,fielddisplay(self,'V','wind speed (m/s-1)')) + string = "%s\n%s"%(string,fielddisplay(self,'V','wind speed (m s-1)')) string = "%s\n%s"%(string,fielddisplay(self,'dlwrf','downward shortwave radiation flux [W/m^2]')) string = "%s\n%s"%(string,fielddisplay(self,'dswrf','downward longwave radiation flux [W/m^2]')) @@ -139,4 +147,5 @@ string = "%s\n%s"%(string,fielddisplay(self,'Tmean','mean annual temperature [K]')) string = "%s\n%s"%(string,fielddisplay(self,'C','mean annual snow accumulation [kg m-2 yr-1]')) + string = "%s\n%s"%(string,fielddisplay(self,'Vmean','mean annual temperature [m s-1] (default 10 m/s)')) string = "%s\n%s"%(string,fielddisplay(self,'Tz','height above ground at which temperature (T) was sampled [m]')) string = "%s\n%s"%(string,fielddisplay(self,'Vz','height above ground at which wind (V) eas sampled [m]')) @@ -196,4 +205,12 @@ '6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)', '7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)'])) + + string = "%s\n%s"%(string,fielddisplay(self,'dsnowIdx',['model for fresh snow accumulation density (default is 1):', + '0 = Original GEMB value, 150 kg/m^3', + '1 = Antarctica value of fresh snow density, 350 kg/m^3', + '2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)', + '3 = Antarctica model of Kaspers et al. (2004), Make sure to set Vmean accurately', + '4 = Greenland model of Kuipers Munneke et al. (2015)'])); + string = "%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) return string @@ -234,4 +251,5 @@ self.swIdx = 1 self.denIdx = 2 + self.dsnowIdx = 1 self.zTop = 10*np.ones((mesh.numberofelements,)) self.dzTop = .05* np.ones((mesh.numberofelements,)) @@ -239,4 +257,6 @@ self.InitDensityScaling = 1.0 self.ThermoDeltaTScaling = 1/11.0 + + self.Vmean = 10*np.ones((mesh.numberofelements,)) self.zMax = 250*np.ones((mesh.numberofelements,)) @@ -291,4 +311,5 @@ md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'>',273-100,'<',273+100) #-100/100 celsius min/max value md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0) + md = checkfield(md,'fieldname','smb.Vmean','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0) md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000) md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000) @@ -299,4 +320,5 @@ 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]) + md = checkfield(md,'fieldname','smb.dsnowIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4]) md = checkfield(md,'fieldname','smb.zTop','NaN',1,'Inf',1,'> = ',0) @@ -356,4 +378,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tmean','format','DoubleMat','mattype',2) WriteData(fid,prefix,'object',self,'class','smb','fieldname','C','format','DoubleMat','mattype',2) + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Vmean','format','DoubleMat','mattype',2) WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tz','format','DoubleMat','mattype',2) WriteData(fid,prefix,'object',self,'class','smb','fieldname','Vz','format','DoubleMat','mattype',2) @@ -368,4 +391,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','swIdx','format','Integer') WriteData(fid,prefix,'object',self,'class','smb','fieldname','denIdx','format','Integer') + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dsnowIdx','format','Integer') WriteData(fid,prefix,'object',self,'class','smb','fieldname','InitDensityScaling','format','Double') WriteData(fid,prefix,'object',self,'class','smb','fieldname','ThermoDeltaTScaling','format','Double') Index: /issm/trunk-jpl/test/NightlyRun/test243.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test243.m (revision 23467) +++ /issm/trunk-jpl/test/NightlyRun/test243.m (revision 23468) @@ -9,4 +9,5 @@ % Use of Gemb method for SMB computation md.smb = SMBgemb(md.mesh,md.geometry); +md.smb.dsnowIdx = 0; %load hourly surface forcing date from 1979 to 2009: Index: /issm/trunk-jpl/test/NightlyRun/test243.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test243.py (revision 23467) +++ /issm/trunk-jpl/test/NightlyRun/test243.py (revision 23468) @@ -20,4 +20,5 @@ md.smb = SMBgemb() md.smb.setdefaultparameters(md.mesh,md.geometry) +md.smb.dsnowIdx = 0 #load hourly surface forcing date from 1979 to 2009: Index: /issm/trunk-jpl/test/NightlyRun/test244.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test244.m (revision 23467) +++ /issm/trunk-jpl/test/NightlyRun/test244.m (revision 23468) @@ -10,4 +10,5 @@ % Use of Gemb method for SMB computation md.smb = SMBgemb(md.mesh,md.geometry); +md.smb.dsnowIdx = 0; %load hourly surface forcing date from 1979 to 2009: Index: /issm/trunk-jpl/test/NightlyRun/test244.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test244.py (revision 23467) +++ /issm/trunk-jpl/test/NightlyRun/test244.py (revision 23468) @@ -30,4 +30,5 @@ md.smb = SMBgemb() md.smb.setdefaultparameters(md.mesh,md.geometry) +md.smb.dsnowIdx = 0 #load hourly surface forcing date from 1979 to 2009: