MITgcm_contrib/mitgcm_tools/loadgrid.m

function [GRID] = loadgrid(varargin)
%loadgrid()
%loadgrid(DIRECTORY)
%
%Reads MITgcm output files and input "data" file to create a GRID structure
%If DIRECTORY is not specified the current working directory is used.
%
%The following files are expected to be in the directory:
% data XC.* YC.* XG.* YG.* DXG.* DYG.* RAC.* hFacC.* hFacW.* hFacS.*
%
%e.g.
%>> GRID=loadgrid
%GRID =
%     drf: [50 70 100 140 190 240 290 340 390 440 490 540 590 640 690]
%     drc: [25 60 85 120 165 215 265 315 365 415 465 515 565 615 665]
%      rf: [1x16 double]
%      rc: [1x15 double]
%      xc: [90x40 double]
%      yc: [90x40 double]
%      xg: [90x40 double]
%      yg: [90x40 double]
%     dxg: [90x40 double]
%     dyg: [90x40 double]
%     rac: [90x40 double]
%   hfacc: [90x40x15 double]
%   hfacw: [90x40x15 double]
%   hfacs: [90x40x15 double]
%    mskc: [90x40x15 double]
%    mAtl: [90x40 double]
%    matl: [90x40 double]
%    mPac: [90x40 double]
%    mpac: [90x40 double]
%    msoc: [90x40 double]
%>> GRID2=loadgrid('/scratch/john/run2/');
%
%Most elements of the structure are named to corresponding to MITgcm
%variables.
%
%Lower case regional masks (matl, mpac and msoc) denote Atlantic and Pacific
%Southern Ocean regions only and are exclusive (i.e. matl does not extend
%into the Southern Ocean).
%Upper case regional masks (mAtl and mPac) denotes sectors and do include the
%Southern Ocean.
%  sum(sum( GRID.mskc(:,:,1) )) = sum( GRID.mAtl(:) + GRID.mPac(:) )
%           = sum( GRID.matl(:) + GRID.mpac(:) + GRID.msoc(:) )
%
%Written by adcroft@mit.edu, 2001
%$Header:

if nargin==0
 Dir='./';
elseif nargin==1
 Dir=[varargin{1} '/'];
else
 error('I don''t know what to do with the second argument');
end

% Extract drF from "data" file
datafile=[Dir 'data'];
fid=fopen(datafile,'r');
if fid==-1
 error(['Could not open file:' datafile ' for reading']);
end
fclose(fid);
drf=evalc([ ...
      '! grep -v ''#'' ' datafile ...
      '| awk ''/[dD[eE][lL][rRzZpP]/,/XXX/ {printf "%s",$0}'' ' ...
      '| sed ''s/[dD][eE][lL][zZrRpP][ ]*=\([0-9Ee,\. +\-]*\).*/\1/''  ' ...
      '| sed ''s/[ ]//g''  ' ...
      '| sed ''s/,/ /g''  ' ...
      ';'
     ]);
eval(['drf=[' drf '];'])
drf=drf;
drc=(drf([1 1:end-1])+drf)/2; drc(1)=drc(1)/2;

%% % Extract drC from output file
%% outputfile='output.txt';
%% drc=evalc(['!head -1000 ' outputfile ...
%%            ' | awk ''/drC/,/;/ {print $3}'' -' ...
%%            ' | egrep "e|E"' ...
%%            ' | sed ''s/,//'' ']);
%% eval(['drc=[' drc '];'])
%% drc=drc';
%% % Extract drF from output file
%% drf=evalc('!head -1000 output.txt | awk "/drF/,/;/ {print \$3}" - | egrep "e|E" | sed s/,// ');
%% eval(['drf=[' drf '];'])
%% drf=drf';

rf=-cumsum([0 drf]);
rc=-cumsum([drc]);

GRID.drf=drf;
GRID.drc=drc;
GRID.rf=rf;
GRID.rc=rc;

msg1='Error: Grid data written by the model is needed. The following files are needed: XC.* YC.* XG.* YG.* DXG.* DYG.* RAC.* hFacC.* hFacW.* hFacS.*';
msg2='Error: The appropriate model output files are not present. If the following files were not written by the model check that you have version 1.11 or greater of ini_grid.F: DXG.* DYG.* RAC.*';

xc=locrdmds([Dir 'XC'],msg1); GRID.xc=xc;
yc=locrdmds([Dir 'YC'],msg1); GRID.yc=yc;
xg=locrdmds([Dir 'XG'],msg1); GRID.xg=xg;
yg=locrdmds([Dir 'YG'],msg1); GRID.yg=yg;
dxg=locrdmds([Dir 'DXG'],msg2); GRID.dxg=dxg;
dyg=locrdmds([Dir 'DYG'],msg2); GRID.dyg=dyg;
rac=locrdmds([Dir 'RAC'],msg2); GRID.rac=rac;
hfacc=locrdmds([Dir 'hFacC'],msg1); GRID.hfacc=hfacc;
hfacw=locrdmds([Dir 'hFacW'],msg1); GRID.hfacw=hfacw;
hfacs=locrdmds([Dir 'hFacS'],msg1); GRID.hfacs=hfacs;

mskc=hfacc; mskc(find(hfacc~=0))=1; GRID.mskc=mskc;

mskc=mskc(:,:,1);
j=[];
%j=[j find( yc>-32.5 & xc>290 )'];
%j=[j find( yc>-32.5 & xc<25 )'];
j=[j find( xc>290 )'];
j=[j find( xc<25 & xc>290-360 )'];
j=[j find( yc>9 & yc<60 & (yc-9)+(xc-276)>0 )'];
j=[j find( yc>9 & yc<60 & (yc-9)+(xc-276+360)>0 & xc<30 )'];
j=[j find( yc>17 & yc<60 & xc>261 )'];
j=[j find( yc>50 & (yc-70)-(xc-270)<0 )'];
j=[j find( yc>31 & xc<38 & xc>-90 )'];
j=[j find( yc>31 & xc>360-90 )'];
j=[j find( yc>64 )'];
matl=0*mskc;
matl(j)=1;
matl=matl.*mskc;
GRID.mAtl=matl;

j=[];
j=[j find( yc>-32.5 )'];
matl=0*mskc;
matl(j)=1;
matl=GRID.mAtl.*matl.*mskc;
GRID.matl=matl;

mpac=(1-matl).*mskc;
j=[];
j=[j find( yc<65 )'];
mpac=0*mskc;
mpac(j)=1;
mpac=mpac.*mskc.*(1-GRID.mAtl);
GRID.mPac=mpac;

j=[];
j=[j find( yc>-32.5 )'];
mpac=0*mskc;
mpac(j)=1;
mpac=GRID.mPac.*mpac.*mskc;
GRID.mpac=mpac;

msoc=(1-matl).*(1-mpac).*mskc;
j=[];
j=[j find( yc<0 )'];
GRID.msoc=msoc.*mskc.*(1-mpac).*(1-matl);

function [A] = locrdmds(fname,errmsg)
[fid,msg]=fopen([fname '.001.001.meta'],'r');
if fid == -1
 A=[];
 disp(errmsg);
 error(['The error occured while trying to open ' fname '.001.001.meta'])
else
 fclose(fid);
 A=rdmds(fname);
end
1	function [GRID] = loadgrid(varargin)
2	%loadgrid()
3	%loadgrid(DIRECTORY)
4	%
5	%Reads MITgcm output files and input "data" file to create a GRID structure
6	%If DIRECTORY is not specified the current working directory is used.
7	%
8	%The following files are expected to be in the directory:
9	% data XC.* YC.* XG.* YG.* DXG.* DYG.* RAC.* hFacC.* hFacW.* hFacS.*
10	%
11	%e.g.
12	%>> GRID=loadgrid
13	%GRID =
14	% drf: [50 70 100 140 190 240 290 340 390 440 490 540 590 640 690]
15	% drc: [25 60 85 120 165 215 265 315 365 415 465 515 565 615 665]
16	% rf: [1x16 double]
17	% rc: [1x15 double]
18	% xc: [90x40 double]
19	% yc: [90x40 double]
20	% xg: [90x40 double]
21	% yg: [90x40 double]
22	% dxg: [90x40 double]
23	% dyg: [90x40 double]
24	% rac: [90x40 double]
25	% hfacc: [90x40x15 double]
26	% hfacw: [90x40x15 double]
27	% hfacs: [90x40x15 double]
28	% mskc: [90x40x15 double]
29	% mAtl: [90x40 double]
30	% matl: [90x40 double]
31	% mPac: [90x40 double]
32	% mpac: [90x40 double]
33	% msoc: [90x40 double]
34	%>> GRID2=loadgrid('/scratch/john/run2/');
35	%
36	%Most elements of the structure are named to corresponding to MITgcm
37	%variables.
38	%
39	%Lower case regional masks (matl, mpac and msoc) denote Atlantic and Pacific
40	%Southern Ocean regions only and are exclusive (i.e. matl does not extend
41	%into the Southern Ocean).
42	%Upper case regional masks (mAtl and mPac) denotes sectors and do include the
43	%Southern Ocean.
44	% sum(sum( GRID.mskc(:,:,1) )) = sum( GRID.mAtl(:) + GRID.mPac(:) )
45	% = sum( GRID.matl(:) + GRID.mpac(:) + GRID.msoc(:) )
46	%
47	%Written by adcroft@mit.edu, 2001
48	%$Header:
49
50	if nargin==0
51	Dir='./';
52	elseif nargin==1
53	Dir=[varargin{1} '/'];
54	else
55	error('I don''t know what to do with the second argument');
56	end
57
58	% Extract drF from "data" file
59	datafile=[Dir 'data'];
60	fid=fopen(datafile,'r');
61	if fid==-1
62	error(['Could not open file:' datafile ' for reading']);
63	end
64	fclose(fid);
65	drf=evalc([ ...
66	'! grep -v ''#'' ' datafile ...
67	'\| awk ''/[dD[eE][lL][rRzZpP]/,/XXX/ {printf "%s",$0}'' ' ...
68	'\| sed ''s/[dD][eE][lL][zZrRpP][ ]=\([0-9Ee,\. +\-]\).*/\1/'' ' ...
69	'\| sed ''s/[ ]//g'' ' ...
70	'\| sed ''s/,/ /g'' ' ...
71	';'
72	]);
73	eval(['drf=[' drf '];'])
74	drf=drf;
75	drc=(drf([1 1:end-1])+drf)/2; drc(1)=drc(1)/2;
76
77	%% % Extract drC from output file
78	%% outputfile='output.txt';
79	%% drc=evalc(['!head -1000 ' outputfile ...
80	%% ' \| awk ''/drC/,/;/ {print $3}'' -' ...
81	%% ' \| egrep "e\|E"' ...
82	%% ' \| sed ''s/,//'' ']);
83	%% eval(['drc=[' drc '];'])
84	%% drc=drc';
85	%% % Extract drF from output file
86	%% drf=evalc('!head -1000 output.txt \| awk "/drF/,/;/ {print \$3}" - \| egrep "e\|E" \| sed s/,// ');
87	%% eval(['drf=[' drf '];'])
88	%% drf=drf';
89
90	rf=-cumsum([0 drf]);
91	rc=-cumsum([drc]);
92
93	GRID.drf=drf;
94	GRID.drc=drc;
95	GRID.rf=rf;
96	GRID.rc=rc;
97
98	msg1='Error: Grid data written by the model is needed. The following files are needed: XC.* YC.* XG.* YG.* DXG.* DYG.* RAC.* hFacC.* hFacW.* hFacS.*';
99	msg2='Error: The appropriate model output files are not present. If the following files were not written by the model check that you have version 1.11 or greater of ini_grid.F: DXG.* DYG.* RAC.*';
100
101	xc=locrdmds([Dir 'XC'],msg1); GRID.xc=xc;
102	yc=locrdmds([Dir 'YC'],msg1); GRID.yc=yc;
103	xg=locrdmds([Dir 'XG'],msg1); GRID.xg=xg;
104	yg=locrdmds([Dir 'YG'],msg1); GRID.yg=yg;
105	dxg=locrdmds([Dir 'DXG'],msg2); GRID.dxg=dxg;
106	dyg=locrdmds([Dir 'DYG'],msg2); GRID.dyg=dyg;
107	rac=locrdmds([Dir 'RAC'],msg2); GRID.rac=rac;
108	hfacc=locrdmds([Dir 'hFacC'],msg1); GRID.hfacc=hfacc;
109	hfacw=locrdmds([Dir 'hFacW'],msg1); GRID.hfacw=hfacw;
110	hfacs=locrdmds([Dir 'hFacS'],msg1); GRID.hfacs=hfacs;
111
112	mskc=hfacc; mskc(find(hfacc~=0))=1; GRID.mskc=mskc;
113
114	mskc=mskc(:,:,1);
115	j=[];
116	%j=[j find( yc>-32.5 & xc>290 )'];
117	%j=[j find( yc>-32.5 & xc<25 )'];
118	j=[j find( xc>290 )'];
119	j=[j find( xc<25 & xc>290-360 )'];
120	j=[j find( yc>9 & yc<60 & (yc-9)+(xc-276)>0 )'];
121	j=[j find( yc>9 & yc<60 & (yc-9)+(xc-276+360)>0 & xc<30 )'];
122	j=[j find( yc>17 & yc<60 & xc>261 )'];
123	j=[j find( yc>50 & (yc-70)-(xc-270)<0 )'];
124	j=[j find( yc>31 & xc<38 & xc>-90 )'];
125	j=[j find( yc>31 & xc>360-90 )'];
126	j=[j find( yc>64 )'];
127	matl=0*mskc;
128	matl(j)=1;
129	matl=matl.*mskc;
130	GRID.mAtl=matl;
131
132	j=[];
133	j=[j find( yc>-32.5 )'];
134	matl=0*mskc;
135	matl(j)=1;
136	matl=GRID.mAtl.matl.mskc;
137	GRID.matl=matl;
138
139	mpac=(1-matl).*mskc;
140	j=[];
141	j=[j find( yc<65 )'];
142	mpac=0*mskc;
143	mpac(j)=1;
144	mpac=mpac.mskc.(1-GRID.mAtl);
145	GRID.mPac=mpac;
146
147	j=[];
148	j=[j find( yc>-32.5 )'];
149	mpac=0*mskc;
150	mpac(j)=1;
151	mpac=GRID.mPac.mpac.mskc;
152	GRID.mpac=mpac;
153
154	msoc=(1-matl).(1-mpac).mskc;
155	j=[];
156	j=[j find( yc<0 )'];
157	GRID.msoc=msoc.mskc.(1-mpac).*(1-matl);
158
159	function [A] = locrdmds(fname,errmsg)
160	[fid,msg]=fopen([fname '.001.001.meta'],'r');
161	if fid == -1
162	A=[];
163	disp(errmsg);
164	error(['The error occured while trying to open ' fname '.001.001.meta'])
165	else
166	fclose(fid);
167	A=rdmds(fname);
168	end