| 1 | %MATICE class definition
 | 
|---|
| 2 | %
 | 
|---|
| 3 | %   Usage:
 | 
|---|
| 4 | %      matice=matice();
 | 
|---|
| 5 | 
 | 
|---|
| 6 | classdef matice
 | 
|---|
| 7 |         properties (SetAccess=public) 
 | 
|---|
| 8 |                 rho_ice                    = 0.;
 | 
|---|
| 9 |                 rho_water                  = 0.;
 | 
|---|
| 10 |                 rho_freshwater             = 0.;
 | 
|---|
| 11 |                 mu_water                   = 0.;
 | 
|---|
| 12 |                 heatcapacity               = 0.;
 | 
|---|
| 13 |                 latentheat                 = 0.;
 | 
|---|
| 14 |                 thermalconductivity        = 0.;
 | 
|---|
| 15 |                 temperateiceconductivity   = 0.;
 | 
|---|
| 16 |                 meltingpoint               = 0.;
 | 
|---|
| 17 |                 beta                       = 0.;
 | 
|---|
| 18 |                 mixed_layer_capacity       = 0.;
 | 
|---|
| 19 |                 thermal_exchange_velocity  = 0.;
 | 
|---|
| 20 |                 rheology_B   = NaN;
 | 
|---|
| 21 |                 rheology_n   = NaN;
 | 
|---|
| 22 |                 rheology_law = '';
 | 
|---|
| 23 | 
 | 
|---|
| 24 |                 %gia: 
 | 
|---|
| 25 |                 lithosphere_shear_modulus  = 0.;
 | 
|---|
| 26 |                 lithosphere_density        = 0.;
 | 
|---|
| 27 |                 mantle_shear_modulus       = 0.;
 | 
|---|
| 28 |                 mantle_density             = 0.;
 | 
|---|
| 29 | 
 | 
|---|
| 30 |                 %slr
 | 
|---|
| 31 |                 earth_density              = 0;
 | 
|---|
| 32 | 
 | 
|---|
| 33 |         end
 | 
|---|
| 34 |         methods
 | 
|---|
| 35 |                 function self = extrude(self,md) % {{{
 | 
|---|
| 36 |                         self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node');
 | 
|---|
| 37 |                         self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element');
 | 
|---|
| 38 |                 end % }}}
 | 
|---|
| 39 |                 function self = matice(varargin) % {{{
 | 
|---|
| 40 |                         switch nargin
 | 
|---|
| 41 |                                 case 0
 | 
|---|
| 42 |                                         self=setdefaultparameters(self);
 | 
|---|
| 43 |                                 case 1
 | 
|---|
| 44 |                                         inputstruct=varargin{1};
 | 
|---|
| 45 |                                         list1 = properties('matice');
 | 
|---|
| 46 |                                         list2 = fieldnames(inputstruct);
 | 
|---|
| 47 |                                         for i=1:length(list1)
 | 
|---|
| 48 |                                                 fieldname = list1{i};
 | 
|---|
| 49 |                                                 if ismember(fieldname,list2),
 | 
|---|
| 50 |                                                         self.(fieldname) = inputstruct.(fieldname);
 | 
|---|
| 51 |                                                 end
 | 
|---|
| 52 |                                         end
 | 
|---|
| 53 |                                 otherwise
 | 
|---|
| 54 |                                         error('constructor not supported');
 | 
|---|
| 55 |                         end
 | 
|---|
| 56 |                 end % }}}
 | 
|---|
| 57 |                 function self = setdefaultparameters(self) % {{{
 | 
|---|
| 58 | 
 | 
|---|
| 59 |                         %ice density (kg/m^3)
 | 
|---|
| 60 |                         self.rho_ice=917.;
 | 
|---|
| 61 | 
 | 
|---|
| 62 |                         %ocean water density (kg/m^3)
 | 
|---|
| 63 |                         self.rho_water=1023.;
 | 
|---|
| 64 | 
 | 
|---|
| 65 |                         %fresh water density (kg/m^3)
 | 
|---|
| 66 |                         self.rho_freshwater=1000.;
 | 
|---|
| 67 | 
 | 
|---|
| 68 |                         %water viscosity (N.s/m^2)
 | 
|---|
| 69 |                         self.mu_water=0.001787;  
 | 
|---|
| 70 | 
 | 
|---|
| 71 |                         %ice heat capacity cp (J/kg/K)
 | 
|---|
| 72 |                         self.heatcapacity=2093.;
 | 
|---|
| 73 | 
 | 
|---|
| 74 |                         %ice latent heat of fusion L (J/kg)
 | 
|---|
| 75 |                         self.latentheat=3.34*10^5;
 | 
|---|
| 76 | 
 | 
|---|
| 77 |                         %ice thermal conductivity (W/m/K)
 | 
|---|
| 78 |                         self.thermalconductivity=2.4;
 | 
|---|
| 79 |                         
 | 
|---|
| 80 |                         %wet ice thermal conductivity (W/m/K)
 | 
|---|
| 81 |                         self.temperateiceconductivity=.24;
 | 
|---|
| 82 | 
 | 
|---|
| 83 |                         %the melting point of ice at 1 atmosphere of pressure in K
 | 
|---|
| 84 |                         self.meltingpoint=273.15;
 | 
|---|
| 85 | 
 | 
|---|
| 86 |                         %rate of change of melting point with pressure (K/Pa)
 | 
|---|
| 87 |                         self.beta=9.8*10^-8;
 | 
|---|
| 88 | 
 | 
|---|
| 89 |                         %mixed layer (ice-water interface) heat capacity (J/kg/K)
 | 
|---|
| 90 |                         self.mixed_layer_capacity=3974.;
 | 
|---|
| 91 | 
 | 
|---|
| 92 |                         %thermal exchange velocity (ice-water interface) (m/s)
 | 
|---|
| 93 |                         self.thermal_exchange_velocity=1.00*10^-4;
 | 
|---|
| 94 | 
 | 
|---|
| 95 |                         %Rheology law: what is the temperature dependence of B with T
 | 
|---|
| 96 |                         %available: none, paterson and arrhenius
 | 
|---|
| 97 |                         self.rheology_law='Paterson';
 | 
|---|
| 98 | 
 | 
|---|
| 99 |                         % GIA:
 | 
|---|
| 100 |                         self.lithosphere_shear_modulus  = 6.7*10^10;  % (Pa)
 | 
|---|
| 101 |                         self.lithosphere_density        = 3.32;       % (g/cm^-3)
 | 
|---|
| 102 |                         self.mantle_shear_modulus       = 1.45*10^11; % (Pa)
 | 
|---|
| 103 |                         self.mantle_density             = 3.34;       % (g/cm^-3)
 | 
|---|
| 104 | 
 | 
|---|
| 105 |                         %SLR
 | 
|---|
| 106 |                         self.earth_density= 5512;  % average density of the Earth, (kg/m^3)
 | 
|---|
| 107 | 
 | 
|---|
| 108 |                 end % }}}
 | 
|---|
| 109 |                 function md = checkconsistency(self,md,solution,analyses) % {{{
 | 
|---|
| 110 |                         md = checkfield(md,'fieldname','materials.rho_ice','>',0);
 | 
|---|
| 111 |                         md = checkfield(md,'fieldname','materials.rho_water','>',0);
 | 
|---|
| 112 |                         md = checkfield(md,'fieldname','materials.rho_freshwater','>',0);
 | 
|---|
| 113 |                         md = checkfield(md,'fieldname','materials.mu_water','>',0);
 | 
|---|
| 114 |                         md = checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1);
 | 
|---|
| 115 |                         md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]);
 | 
|---|
| 116 |                         md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'});
 | 
|---|
| 117 | 
 | 
|---|
| 118 |                         if ismember('GiaAnalysis',analyses),
 | 
|---|
| 119 |                                 md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1);
 | 
|---|
| 120 |                                 md = checkfield(md,'fieldname','materials.lithosphere_density','>',0,'numel',1);
 | 
|---|
| 121 |                                 md = checkfield(md,'fieldname','materials.mantle_shear_modulus','>',0,'numel',1);
 | 
|---|
| 122 |                                 md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',1);
 | 
|---|
| 123 |                         end
 | 
|---|
| 124 |                         if ismember('SealevelriseAnalysis',analyses),
 | 
|---|
| 125 |                                 md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1);
 | 
|---|
| 126 |                         end
 | 
|---|
| 127 | 
 | 
|---|
| 128 |                 end % }}}
 | 
|---|
| 129 |                 function disp(self) % {{{
 | 
|---|
| 130 |                         disp(sprintf('   Materials:'));
 | 
|---|
| 131 | 
 | 
|---|
| 132 |                         fielddisplay(self,'rho_ice','ice density [kg/m^3]');
 | 
|---|
| 133 |                         fielddisplay(self,'rho_water','ocean water density [kg/m^3]');
 | 
|---|
| 134 |                         fielddisplay(self,'rho_freshwater','fresh water density [kg/m^3]');
 | 
|---|
| 135 |                         fielddisplay(self,'mu_water','water viscosity [N s/m^2]');
 | 
|---|
| 136 |                         fielddisplay(self,'heatcapacity','heat capacity [J/kg/K]');
 | 
|---|
| 137 |                         fielddisplay(self,'thermalconductivity',['ice thermal conductivity [W/m/K]']);
 | 
|---|
| 138 |                         fielddisplay(self,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]');
 | 
|---|
| 139 |                         fielddisplay(self,'meltingpoint','melting point of ice at 1atm in K');
 | 
|---|
| 140 |                         fielddisplay(self,'latentheat','latent heat of fusion [J/kg]');
 | 
|---|
| 141 |                         fielddisplay(self,'beta','rate of change of melting point with pressure [K/Pa]');
 | 
|---|
| 142 |                         fielddisplay(self,'mixed_layer_capacity','mixed layer capacity [W/kg/K]');
 | 
|---|
| 143 |                         fielddisplay(self,'thermal_exchange_velocity','thermal exchange velocity [m/s]');
 | 
|---|
| 144 |                         fielddisplay(self,'rheology_B','flow law parameter [Pa/s^(1/n)]');
 | 
|---|
| 145 |                         fielddisplay(self,'rheology_n','Glen''s flow law exponent');
 | 
|---|
| 146 |                         fielddisplay(self,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']);
 | 
|---|
| 147 |                         fielddisplay(self,'lithosphere_shear_modulus','Lithosphere shear modulus [Pa]');
 | 
|---|
| 148 |                         fielddisplay(self,'lithosphere_density','Lithosphere density [g/cm^-3]');
 | 
|---|
| 149 |                         fielddisplay(self,'mantle_shear_modulus','Mantle shear modulus [Pa]');
 | 
|---|
| 150 |                         fielddisplay(self,'mantle_density','Mantle density [g/cm^-3]');
 | 
|---|
| 151 |                         fielddisplay(self,'earth_density','Mantle density [kg/m^-3]');
 | 
|---|
| 152 |                 end % }}}
 | 
|---|
| 153 |                 function marshall(self,prefix,md,fid) % {{{
 | 
|---|
| 154 |                         WriteData(fid,prefix,'name','md.materials.type','data',3,'format','Integer');
 | 
|---|
| 155 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double');
 | 
|---|
| 156 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double');
 | 
|---|
| 157 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double');
 | 
|---|
| 158 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','mu_water','format','Double');
 | 
|---|
| 159 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','heatcapacity','format','Double');
 | 
|---|
| 160 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','latentheat','format','Double');
 | 
|---|
| 161 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double');
 | 
|---|
| 162 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double');
 | 
|---|
| 163 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double');
 | 
|---|
| 164 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double');
 | 
|---|
| 165 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double');
 | 
|---|
| 166 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double');
 | 
|---|
| 167 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 | 
|---|
| 168 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2);
 | 
|---|
| 169 |                         WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String');
 | 
|---|
| 170 | 
 | 
|---|
| 171 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_shear_modulus','format','Double');
 | 
|---|
| 172 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_density','format','Double','scale',10^3);
 | 
|---|
| 173 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','mantle_shear_modulus','format','Double');
 | 
|---|
| 174 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','mantle_density','format','Double','scale',10^3);
 | 
|---|
| 175 |                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','earth_density','format','Double');
 | 
|---|
| 176 |                 end % }}}
 | 
|---|
| 177 |                 function savemodeljs(self,fid,modelname) % {{{
 | 
|---|
| 178 |                 
 | 
|---|
| 179 |                         writejsdouble(fid,[modelname '.materials.rho_ice'],self.rho_ice);
 | 
|---|
| 180 |                         writejsdouble(fid,[modelname '.materials.rho_water'],self.rho_water);
 | 
|---|
| 181 |                         writejsdouble(fid,[modelname '.materials.rho_freshwater'],self.rho_freshwater);
 | 
|---|
| 182 |                         writejsdouble(fid,[modelname '.materials.mu_water'],self.mu_water);
 | 
|---|
| 183 |                         writejsdouble(fid,[modelname '.materials.heatcapacity'],self.heatcapacity);
 | 
|---|
| 184 |                         writejsdouble(fid,[modelname '.materials.latentheat'],self.latentheat);
 | 
|---|
| 185 |                         writejsdouble(fid,[modelname '.materials.thermalconductivity'],self.thermalconductivity);
 | 
|---|
| 186 |                         writejsdouble(fid,[modelname '.materials.temperateiceconductivity'],self.temperateiceconductivity);
 | 
|---|
| 187 |                         writejsdouble(fid,[modelname '.materials.meltingpoint'],self.meltingpoint);
 | 
|---|
| 188 |                         writejsdouble(fid,[modelname '.materials.beta'],self.beta);
 | 
|---|
| 189 |                         writejsdouble(fid,[modelname '.materials.mixed_layer_capacity'],self.mixed_layer_capacity);
 | 
|---|
| 190 |                         writejsdouble(fid,[modelname '.materials.thermal_exchange_velocity'],self.thermal_exchange_velocity);
 | 
|---|
| 191 |                         writejsdouble(fid,[modelname '.materials.mixed_layer_capacity'],self.mixed_layer_capacity);
 | 
|---|
| 192 |                         writejs1Darray(fid,[modelname '.materials.rheology_B'],self.rheology_B);
 | 
|---|
| 193 |                         writejs1Darray(fid,[modelname '.materials.rheology_n'],self.rheology_n);
 | 
|---|
| 194 |                         writejsstring(fid,[modelname '.materials.rheology_law'],self.rheology_law);
 | 
|---|
| 195 |                         writejsdouble(fid,[modelname '.materials.lithosphere_shear_modulus'],self.lithosphere_shear_modulus);
 | 
|---|
| 196 |                         writejsdouble(fid,[modelname '.materials.lithosphere_density'],self.lithosphere_density);
 | 
|---|
| 197 |                         writejsdouble(fid,[modelname '.materials.mantle_shear_modulus'],self.mantle_shear_modulus);
 | 
|---|
| 198 |                         writejsdouble(fid,[modelname '.materials.mantle_density'],self.mantle_density);
 | 
|---|
| 199 |                         writejsdouble(fid,[modelname '.materials.earth_density'],self.earth_density);
 | 
|---|
| 200 | 
 | 
|---|
| 201 |                 end % }}}
 | 
|---|
| 202 |         end
 | 
|---|
| 203 | end
 | 
|---|