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	adcroft	1.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	adcroft	1.2	%The following files are expected to be in the directory:
9			% data XC.* YC.* XG.* YG.* DXG.* DYG.* RAC.* hFacC.* hFacW.* hFacS.*
10			%
11	adcroft	1.1	%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	adcroft	1.2	'! grep -v ''#'' ' datafile ...
67			'\| awk ''/[dD[eE][lL][rRzZpP]/,/XXX/ {printf "%s",$0}'' ' ...
68			'\| sed ''s/[dD][eE][lL][zZrRpP][ ]=\([0-9Ee,\. +\-]\).*/\1/'' ' ...
69	adcroft	1.1	'\| 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	adcroft	1.2	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	adcroft	1.1
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	adcroft	1.2	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