Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 19524) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 19524) @@ -0,0 +1,199 @@ +%SMBgemb Class definition. +% This is the class that hosts all the inputs for the Alberta Glacier Surface Mass Balance Model +% Alex Gardner, University of Alberta. +% +% Usage: +% SMBgemb=SMBgemb(); + +classdef SMBgemb + properties (SetAccess=public) + % {{{ + %each one of these properties is a transient forcing to the GEMB model, loaded from meteorological data derived + %from an automatic weather stations (AWS). Each property is therefore a matrix, of size (numberofvertices x number + %of time steps. ) + + %inputs: + Ta = NaN; %2 m air temperature, in Celsius + V = NaN; %wind speed (m/s-1) + dswrf = NaN; %downward shortwave radiation flux [W/m^2] + dlwrf = NaN; %downward longwave radiation flux [W/m^2] + RH = NaN; %relative humidity [%] + P = NaN; %precipitation [mm w.e. / m^2] + eAir = NaN; %screen level vapor pressure [Pa] + pAir = NaN; %surface pressure [Pa] + + Tmean = NaN; %mean annual temperature [K] + C = NaN; %mean annual snow accumulation [kg m-2 yr-1] + Tz = NaN; %height above ground at which temperature (T) was sampled [m] + Vz = NaN; %height above ground at which wind (V) eas sampled [m] + + %settings: + spinUp = NaN; %number of cycles of met data run before output is calculated (default is 0) + aIdx = NaN; %method for calculating albedo and subsurface absorption (default is 1) + % 1: effective grain radius [Gardner & Sharp, 2009] + % 2: effective grain radius [Brun et al., 2009] + % 3: density and cloud amount [Greuell & Konzelmann, 1994] + % 4: exponential time decay & wetness [Bougamont & Bamber, 2005] + swIdx = NaN; % apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1) + + denIdx = NaN; %densification model to use (default is 2): + % 1 = emperical model of Herron and Langway (1980) + % 2 = semi-emerical model of Anthern et al. (2010) + % 3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010) + % 4 = DO NOT USE: emperical model of Li and Zwally (2004) + % 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008) + + 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] + dzMin = NaN; % initial min vertical allowable grid spacing (default dzMin/2) [m] + zY = NaN; % strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)] + zMax = NaN; %initial max model depth (default is min(thickness,500)) [m] + zMin = NaN; %initial min model depth (default is min(thickness,30)) [m] + outputFreq = NaN; %output frequency in days (default is monthly, 30) + + %specific albedo parameters: + %Method 1 and 2: + aSnow = NaN; % new snow albedo (0.64 - 0.89) + aIce = NaN; % range 0.27-0.58 for old snow + %Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo + cldFrac = NaN; % average cloud amount + %Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005) + t0wet = NaN; % time scale for wet snow (15-21.9) + t0dry = NaN; % warm snow timescale (30) + K = NaN; % time scale temperature coef. (7) + + %Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM. + %dateN: that's the last row of the above fields. + %dt: included in dateN. Not an input. + %elev: this is taken from the ISSM surface itself. + + end % }}} + methods + function self = SMBgemb(varargin) % {{{ + switch nargin + case 2 + mesh=varargin{1}; + geometry=varargin{2}; + self=setdefaultparameters(self,mesh,geometry); + otherwise + error('constructor not supported'); + end + end % }}} + function self = extrude(self,md) % {{{ + + self.Ta=project3d(md,'vector',self.Ta,'type','node'); + self.V=project3d(md,'vector',self.V,'type','node'); + self.dswrf=project3d(md,'vector',self.dswrf,'type','node'); + self.dswrf=project3d(md,'vector',self.dswrf,'type','node'); + self.RH=project3d(md,'vector',self.RH,'type','node'); + self.P=project3d(md,'vector',self.P,'type','node'); + self.eAir=project3d(md,'vector',self.eAir,'type','node'); + self.pAir=project3d(md,'vector',self.pAir,'type','node'); + + end % }}} + function self = setdefaultparameters(self,mesh,geometry) % {{{ + + + self.spinUp=0; + self.aIdx = 1; + self.swIdx = 1; + self.denIdx = 2; + self.zTop=10*ones(mesh.numberofelements,1); + self.dzTop = .05* ones (mesh.numberofelements,1); + self.dzMin = self.dzTop/2; + + he=sum(geometry.thickness(mesh.elements),2)/size(mesh.elements,2); + self.zMax=min(500,he); + self.zMin=min(30,he); + self.zY = 1.10*ones(mesh.numberofelements,1); + self.outputFreq = 30; + + %additional albedo parameters + self.aSnow = 0.85; + self.aIce = 0.48; + self.cldFrac = 0.1; + self.t0wet = 15; + self.t0dry = 30; + self.K = 7; + + end % }}} +function md = checkconsistency(self,md,solution,analyses) % {{{ + + md = checkfield(md,'fieldname','surfaceforcings.Ta','timeseries',1,'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.V','timeseries',1,'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.dswrf','timeseries',1,'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.dlwrf','timeseries',1,'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.RH','timeseries',1,'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.P','timeseries',1,'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.eAir','timeseries',1,'NaN',1); + + %check zTop is < local thickness: + he=sum(md.geometry.thickness(md.mesh.elements),2)/size(md.mesh.elements,2); + if any(he