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.
%
%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
%

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([ ...
      '!awk ''/[dD[eE][lL][rRzZpP]/,/XXX/ {printf "%s",$0}'' ' datafile ...
      '| sed ''s/[dD][eE][lL][zZrRpP][ ]*=\([0-9,. ]*\).*/\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;

xc=rdmds([Dir 'XC']); GRID.xc=xc;
yc=rdmds([Dir 'YC']); GRID.yc=yc;
xg=rdmds([Dir 'XG']); GRID.xg=xg;
yg=rdmds([Dir 'YG']); GRID.yg=yg;

dxg=rdmds([Dir 'DXG']); GRID.dxg=dxg;
dyg=rdmds([Dir 'DYG']); GRID.dyg=dyg;

rac=rdmds([Dir 'RAC']); GRID.rac=rac;

hfacc=rdmds([Dir 'hFacC']); GRID.hfacc=hfacc;
hfacw=rdmds([Dir 'hFacW']); GRID.hfacw=hfacw;
hfacs=rdmds([Dir 'hFacS']); 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);