source: issm/oecreview/Archive/24684-25833/ISSM-24847-24848.diff

Last change on this file was 25834, checked in by Mathieu Morlighem, 4 years ago

CHG: added 24684-25833
File size: 6.0 KB
RevLine 
[25834]1Index: ../trunk-jpl/src/m/contrib/larour/glaciermip.m
2===================================================================
3--- ../trunk-jpl/src/m/contrib/larour/glaciermip.m (nonexistent)
4+++ ../trunk-jpl/src/m/contrib/larour/glaciermip.m (revision 24848)
5@@ -0,0 +1,152 @@
6+%GLACIERMIP class definition
7+%
8+% Usage:
9+% gmip = glaciermip('ncfile',glaciermipncfile,'version',1)
10+%
11+%
12+
13+classdef glaciermip < handle
14+ properties (SetAccess=public) %Model fields
15+ ncfile = '';
16+ version = NaN;
17+
18+ %version 2:
19+ area =[];
20+ mass =[];
21+ region =[];
22+ time =[];
23+ glaciermodel =[];
24+ climatemodel =[];
25+ scenario =[];
26+
27+ %from version 1:
28+ run=[];
29+ forcingmodel=[];
30+ realization=[];
31+ volume=[];
32+ end
33+ methods
34+ function self = glaciermip(varargin) % {{{
35+
36+ if nargin==0,
37+ self=setdefaultparameters(self);
38+ else
39+ self=setdefaultparameters(self);
40+
41+ options=pairoptions(varargin{:});
42+
43+ self.ncfile=getfieldvalue(options,'ncfile');
44+ self.version=getfieldvalue(options,'version');
45+
46+
47+ %read variables:
48+ if self.version==1,
49+ self.region=ncread(self.ncfile,'region');
50+ self.time=ncread(self.ncfile,'time');
51+ self.run=ncread(self.ncfile,'run');
52+ self.glaciermodel=ncread(self.ncfile,'glaciermodel');
53+ self.forcingmodel=ncread(self.ncfile,'forcingmodel');
54+ self.scenario=ncread(self.ncfile,'scenario');
55+ self.realization=ncread(self.ncfile,'realization');
56+ self.area=ncread(self.ncfile,'area');
57+ self.volume=ncread(self.ncfile,'volume');
58+ elseif self.version==2,
59+ self.area=ncread(self.ncfile,'Area');
60+ self.mass=ncread(self.ncfile,'Mass');
61+ self.region=ncread(self.ncfile,'Region');
62+ self.time=ncread(self.ncfile,'Time');
63+ self.glaciermodel=ncread(self.ncfile,'Glacier_Model');
64+ self.climatemodel=ncread(self.ncfile,'Climate_Model');
65+ self.scenario=ncread(self.ncfile,'Scenario');
66+ %mass(region,time,climatemodel,glaciermodel,scenario)
67+ else
68+ error(sprintf('glaciermipfile constructor error message: version %i for MIP not supported!'),self.version);
69+ end
70+ end
71+ end
72+ %}}}
73+ function inv = setdefaultparameters(inv) % {{{
74+ end
75+ %}}}
76+ function t=gettime(self,varargin) % {{{
77+
78+ options=pairoptions(varargin{:});
79+ t=self.time;
80+
81+ end % }}}
82+ function masses=getmass(self,varargin) % {{{
83+
84+ options=pairoptions(varargin{:});
85+
86+ rg=getfieldvalue(options,'region',1:length(self.region));
87+ cm=getfieldvalue(options,'climatemodel',1:length(self.climatemodel));
88+ gm=getfieldvalue(options,'glaciermodel',1:length(self.glaciermodel));
89+ sc=getfieldvalue(options,'scenario',1:length(self.scenario));
90+ zerostonan=getfieldvalue(options,'zerostonan',0);
91+
92+ if self.version==1,
93+ error(sprintf('getmass not supported yet for Glacier MIP version %i',self.version));
94+ end
95+
96+ %serialize:
97+ masses=[];
98+ for i=rg,
99+ for j=cm,
100+ for k=gm,
101+ for l=sc,
102+ masses=[masses; squeeze(self.mass(i,:,j,k,l))];
103+ end
104+ end
105+ end
106+ end
107+ if zerostonan,
108+ masses(find(masses==0))=NaN;
109+ end
110+
111+ end % }}}
112+ function massrates=getmassrates(self,varargin) % {{{
113+
114+ options=pairoptions(varargin{:});
115+
116+ %get mass first:
117+ masses=self.getmass(varargin{:});
118+
119+ %compute mass rates:
120+ dt=diff(self.time);
121+ dm=diff(masses,1,2);
122+
123+ massrates=dm;
124+ for i=1:size(massrates,1),
125+ massrates(i,:)= massrates(i,:)./dt';
126+ end
127+
128+
129+ end % }}}
130+ function disp(self) % {{{
131+ disp(sprintf(' Glacier MIP (version %i):',self.version));
132+
133+ if self.version==1,
134+ fielddisplay(self,'ncfile','netcdf file for GlacierMIP results');
135+ fielddisplay(self,'time','time scale in yr');
136+ fielddisplay(self,'region','region');
137+ fielddisplay(self,'run','run');
138+ fielddisplay(self,'glaciermodel','glacier model');
139+ fielddisplay(self,'forcingmodel','forcing model');
140+ fielddisplay(self,'scenario','scenario');
141+ fielddisplay(self,'realization','realization');
142+ fielddisplay(self,'area','area');
143+ fielddisplay(self,'volume','volume');
144+ end
145+ if self.version==2,
146+ fielddisplay(self,'ncfile','netcdf file for GlacierMIP results');
147+ fielddisplay(self,'time','time');
148+ fielddisplay(self,'region','region');
149+ fielddisplay(self,'glaciermodel','glaciermodel');
150+ fielddisplay(self,'climatemodel','climatemodel');
151+ fielddisplay(self,'scenario','scenario');
152+ fielddisplay(self,'area','area');
153+ fielddisplay(self,'mass','mass');
154+ end
155+ end % }}}
156+ end
157+end
158Index: ../trunk-jpl/src/m/contrib/larour/glacier_inventory.m
159===================================================================
160--- ../trunk-jpl/src/m/contrib/larour/glacier_inventory.m (revision 24847)
161+++ ../trunk-jpl/src/m/contrib/larour/glacier_inventory.m (revision 24848)
162@@ -92,6 +92,22 @@
163 end
164 end
165 %}}}
166+ function totalarea=area(self,varargin) % {{{
167+ region=-1;
168+ totalarea=0;
169+ if nargin==2,
170+ region=varargin{1};
171+ end
172+ if region==-1,
173+ %figure out the areas of everybody:
174+ for i=1:self.nregions(),
175+ totalarea=totalarea+sum(self.regions(i).Area);
176+ end
177+ else
178+ totalarea=totalarea+sum(self.regions(region).Area);
179+ end
180+
181+ end % }}}
182 function [mpartition,npartition]=partition(self,varargin) % {{{
183 mpartition=zeros(self.nregions(),1);
184 npartition=zeros(self.nregions(),1);
185@@ -239,7 +255,7 @@
186 %}}}
187 function mesh_connectivity(self,mesh) % {{{
188
189- %initialize glacier_connectivity:
190+ %initialize glacier and element connectivity:
191 glacier_connectivity=zeros(self.nglaciers,1);
192
193 %assume maximum width for connectivity table and initialize
194@@ -325,14 +341,11 @@
195 %For each element, this will tell us which glaciers belong to this element.
196 for j=1:length(glacier_connectivity),
197 el=glacier_connectivity(j);
198- if el==0,
199- continue;
200- end
201 count=element_connectivity(el,ny);
202 if count>ny,
203 error('need to enlarge connectivity table');
204 end
205- element_connectivity(el,count+1)=partition_counter+j;
206+ element_connectivity(el,count+1)=j;
207 element_connectivity(el,ny)=count+1;
208 end
209
Note: See TracBrowser for help on using the repository browser.