function [FORCING] = loadforcing(varargin)
%loadforcing()
%loadforcing(DIRECTORY)
%
%Reads MITgcm parameter input files ("data") file to create a FORCING structure
%If DIRECTORY is not specified the current working directory is used.
%
%e.g.
% >> FORCING=loadforcing
% FORCING = 
% >> FORCING2=loadforcing('/scratch/john/run2/');
%
%Written by adcroft@mit.edu, 2002
%

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 names of forcing files from "data" file
datafile=[Dir 'data'];

fid=fopen(datafile,'r');
if fid==-1
 error(['Could not open file:' datafile ' for reading']);
end

tinitfile=grepparameter(datafile,'hydrogthetafile');
sinitfile=grepparameter(datafile,'hydrogsaltfile');
sstfile=grepparameter(datafile,'thetaclimfile');
sssfile=grepparameter(datafile,'saltclimfile');
tauxfile=grepparameter(datafile,'zonalwindfile');
tauyfile=grepparameter(datafile,'meridwindfile');
qfile=grepparameter(datafile,'surfqfile');
empmrfile=grepparameter(datafile,'empmrfile');

GRID=loadgrid(Dir);
nxy=size(GRID.rac);
nxyr=size(GRID.hfacc);
clear GRID

FORCING.fmt=grepparameter(datafile,'readbinaryprec');
if isempty(FORCING.fmt)
 FORCING.fmt=32;
end
if FORCING.fmt==64
 fmt='real*8';
else
 fmt='real*4';
end
FORCING.fmt=fmt;

FORCING.tinit=readforcing(Dir,tinitfile,fmt,nxy,nxyr);
FORCING.sinit=readforcing(Dir,sinitfile,fmt,nxy,nxyr);
FORCING.sst=readforcing(Dir,sstfile,fmt,nxy,nxyr);
FORCING.sss=readforcing(Dir,sssfile,fmt,nxy,nxyr);
FORCING.taux=readforcing(Dir,tauxfile,fmt,nxy,nxyr);
FORCING.tauy=readforcing(Dir,tauyfile,fmt,nxy,nxyr);
FORCING.q=readforcing(Dir,qfile,fmt,nxy,nxyr);
FORCING.empmr=readforcing(Dir,empmrfile,fmt,nxy,nxyr);

function [forcingdata] = readforcing(Dir,filename,fmt,nxy,nxyr)

forcingdata=[];
if isempty(filename)
 return;
end

fid=fopen([Dir filename],'r','b');
if fid==-1
  error(['Could not open file:' filename ' for reading']);
end
forcingdata=fread(fid,inf,fmt);
N=size(forcingdata,1);
if N>prod(nxyr)
 nshape=[nxyr size(forcingdata,1)/prod(nxyr)];
else
 nshape=[nxy N/prod(nxy)];
end
forcingdata=reshape(forcingdata,nshape);