/[MITgcm]/MITgcm/utils/matlab/load_grid.m
ViewVC logotype

Contents of /MITgcm/utils/matlab/load_grid.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (show annotations) (download)
Tue Feb 12 16:06:32 2013 UTC (11 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64o, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint64n, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint65, checkpoint64j, checkpoint64m, checkpoint64l, HEAD
Changes since 1.1: +5 -5 lines
fix vertical coordinate output from MNC grid-files

1 function G = load_grid(varargin);
2 % function G = load_grid(dirName,option,nyAxis,nxAxis);
3 %
4 % load the MITgcm output grid-files into one structure array
5 % arguments:
6 % dirName = directory where grid-files are
7 % optional arguments:
8 % ncdf = last digit of integer "option":
9 % ncdf = 0,2 : read binary (MDSIO) files using rdmds (default: ncdf=0)
10 % ncdf = 1,3 : read NetCDF (MNC) files using rdmnc
11 % ncdf = 2,3 : as 0,1 + print elapsed-time for loading files
12 % option >= 10 : only read in Horiz.Grid spacing (without hFac)
13 % option >= 20 : only read in Verti.Grid spacing (without hFac)
14 % nxAxis = number of spacing in 1.D x-axis xAxC (default: nxAxis = grid-Nx)
15
16 % $Header: /u/gcmpack/MITgcm/utils/matlab/load_grid.m,v 1.1 2012/10/04 14:02:42 jmc Exp $
17 % $Name: $
18
19 if nargin == 0
20 rDir = '.';
21 else
22 rDir = varargin{1};
23 end
24 if nargin < 2,
25 ktyp=0;
26 ncdf=0;
27 else
28 ncdf= varargin{2};
29 ktyp=floor(ncdf/10); ncdf=rem(ncdf,10);
30 end
31
32 if ncdf > 1, TimeT0=clock; end
33
34 if rem(ncdf,2) == 0,
35 %- load MDSIO grid files :
36 if ktyp < 2,
37 xC=rdmds(fullfile(rDir,'XC'));
38 yC=rdmds(fullfile(rDir,'YC'));
39 xG=rdmds(fullfile(rDir,'XG'));
40 yG=rdmds(fullfile(rDir,'YG'));
41
42 dXc=rdmds(fullfile(rDir,'DXC'));
43 dYc=rdmds(fullfile(rDir,'DYC'));
44 dXg=rdmds(fullfile(rDir,'DXG'));
45 dYg=rdmds(fullfile(rDir,'DYG'));
46
47 rAc=rdmds(fullfile(rDir,'RAC'));
48 rAw=rdmds(fullfile(rDir,'RAW'));
49 rAs=rdmds(fullfile(rDir,'RAS'));
50 rAz=rdmds(fullfile(rDir,'RAZ'));
51
52 %- load grid orientation relative to West-East / South-North dir:
53 if isempty(dir(fullfile(rDir,'AngleCS.*'))),
54 fprintf(' no grid orientation (Cos & Sin) file to load ');
55 csAngle= ones(size(rAc));
56 snAngle=zeros(size(rAc));
57 fprintf(' => set COS=1,SIN=0\n');
58 else
59 csAngle=rdmds(fullfile(rDir,'AngleCS'));
60 snAngle=rdmds(fullfile(rDir,'AngleSN'));
61 end
62 dims = size(rAc);
63 end
64
65 if rem(ktyp,2) == 0,
66 rC=rdmds(fullfile(rDir,'RC')); rC=squeeze(rC);
67 rF=rdmds(fullfile(rDir,'RF')); rF=squeeze(rF);
68 dRc=rdmds(fullfile(rDir,'DRC')); dRc=squeeze(dRc);
69 dRf=rdmds(fullfile(rDir,'DRF')); dRf=squeeze(dRf);
70 dims = size(rC);
71 end
72
73 if ktyp == 0,
74 %- define domain:
75 hFacC=rdmds(fullfile(rDir,'hFacC'));
76 hFacW=rdmds(fullfile(rDir,'hFacW'));
77 hFacS=rdmds(fullfile(rDir,'hFacS'));
78 depth=rdmds(fullfile(rDir,'Depth'));
79 dims = size(hFacC);
80 end
81
82 if ncdf > 1, TimeT1=clock; end
83
84 else
85 %- load NetCDF grid files :
86 %S=rdmnc([rDir,'grid.*.nc']);
87 S=rdmnc(fullfile(rDir,'grid.*.nc'));
88 if ncdf > 1, TimeT1=clock; end
89
90 %--
91 xC=S.XC;
92 yC=S.YC;
93 xG=S.XG(1:end-1,1:end-1);
94 yG=S.YG(1:end-1,1:end-1);
95 rC=S.RC';
96 rF=S.RF';
97 dXc=S.dxC(1:end-1,:);
98 dYc=S.dyC(:,1:end-1);
99 dXg=S.dxG(:,1:end-1);
100 dYg=S.dyG(1:end-1,:);
101 dRc=S.drC';
102 dRf=S.drF';
103 rAc=S.rA;
104 rAw=S.rAw(1:end-1,:);
105 rAs=S.rAs(:,1:end-1);
106 rAz=S.rAz(1:end-1,1:end-1);
107 if isfield(S,'AngleCS') & isfield(S,'AngleSN'),
108 csAngle=S.AngleCS;
109 snAngle=S.AngleSN;
110 else
111 fprintf(' no grid orientation (Cos & Sin) in grid file ');
112 csAngle= ones(size(rAc));
113 snAngle=zeros(size(rAc));
114 fprintf(' => set COS=1,SIN=0\n');
115 end
116 hFacC=S.HFacC;
117 hFacW=S.HFacW(1:end-1,:,:);
118 hFacS=S.HFacS(:,1:end-1,:);
119 depth=S.Depth;
120 dims = size(hFacC);
121 end
122
123 if ncdf > 1, fprintf(' (took %6.4f s)\n',etime(TimeT1,TimeT0)); end
124
125 if length(dims) == 1, dims(2)=1; end
126 if length(dims) == 2, dims(3)=1; end
127
128 %- 1-D axis:
129 if nargin < 3,
130 %- yAxis
131 ny=dims(2);
132 yAxC=yC(1,:)';
133 if rem(ncdf,2) == 0,
134 yAxV=[yG(1,:) 2*yC(1,end)-yG(1,end)]';
135 else
136 yAxV=S.YG(1,:)';
137 end
138 else
139 ny=varargin{3};
140 dyAx=(max(yG(:))-min(yG(:)))/ny;
141 yAxV=min(yG(:))+[0:ny]'*dyAx;
142 yAxC=(yAxV(1:ny)+yAxV(2:ny+1))/2;
143 end
144 if nargin < 4,
145 %- xAxis
146 nx=dims(1);
147 xAxC=xC(:,1);
148 if rem(ncdf,2) == 0,
149 xAxU=[xG(:,1)' 2*xC(end,1)-xG(end,1)]';
150 else
151 xAxU=S.XG(:,1);
152 end
153 else
154 nx=varargin{4};
155 dxAx=max(max(xG(:)),max(xC(:)))-min(min(xG(:)),min(xC(:)));
156 if dxAx > 360*(1-1/nx), dxAx=360; end
157 dxAx=dxAx/nx;
158 xAxU=min(xG(:))+[0:nx]'*dxAx;
159 xAxC=(xAxU(1:nx)+xAxU(2:nx+1))/2;
160 end
161
162 %- volume:
163 if rem(ncdf,2) == 1 | ktyp == 0,
164 %fprintf(' rAc:'); fprintf(' %i',size(rAc)); fprintf('\n');
165 %fprintf('dims:'); fprintf(' %i',dims); fprintf('\n');
166 %fprintf(' dRf:'); fprintf(' %i',size(dRf)); fprintf('\n');
167 %vol=repmat(rAc,[1 1 dims(3)]).*hFacC;
168 vol=reshape(rAc,[dims(1)*dims(2) 1])*dRf'; vol=reshape(vol,dims).*hFacC;
169 end
170
171 % create the structure
172
173 if rem(ncdf,2) == 1 | ktyp == 0,
174 G = struct('dims',dims, ...
175 'nx',nx,'ny',ny,'xAxC',xAxC,'yAxC',yAxC,'xAxU',xAxU,'yAxV',yAxV, ...
176 'xC',xC,'yC',yC,'xG',xG,'yG',yG,'rC',rC,'rF',rF, ...
177 'dXc',dXc,'dYc',dYc,'dXg',dXg,'dYg',dYg,'dRc',dRc,'dRf',dRf, ...
178 'rAc',rAc,'rAw',rAw,'rAs',rAs,'rAz',rAz, ...
179 'csAngle',csAngle,'snAngle',snAngle, ...
180 'hFacC',hFacC,'hFacW',hFacW,'hFacS',hFacS,'depth',depth,'vol',vol);
181 elseif ktyp == 1,
182 G = struct('dims',dims, ...
183 'nx',nx,'ny',ny,'xAxC',xAxC,'yAxC',yAxC,'xAxU',xAxU,'yAxV',yAxV, ...
184 'xC',xC,'yC',yC,'xG',xG,'yG',yG, ...
185 'dXc',dXc,'dYc',dYc,'dXg',dXg,'dYg',dYg, ...
186 'rAc',rAc,'rAw',rAw,'rAs',rAs,'rAz',rAz, ...
187 'csAngle',csAngle,'snAngle',snAngle);
188 else
189 G = struct('dims',dims, ...
190 'rC',rC,'rF',rF, ...
191 'dRc',dRc,'dRf',dRf);
192 end
193
194 return

  ViewVC Help
Powered by ViewVC 1.1.22