function ERAI_singv_MFD(outdir)

% Computes moisture flux divergence
%
% CEB Sept 2012

g = 9.81;

% RCM dates and times
dates = [];
dates = [dates, [datenum([2000 12 1 0 0 0]):datenum([0 0 1 0 0 0]):datenum([2000 12 31 0 0 0])]];
for yy=[2001:2009]
  dates = [dates, [datenum([yy 1 1 0 0 0]):datenum([0 0 1 0 0 0]):datenum([yy 2 28 0 0 0])]];
  dates = [dates, [datenum([yy 12 1 0 0 0]):datenum([0 0 1 0 0 0]):datenum([yy 12 31 0 0 0])]];
end
dates = [dates, [datenum([2010 1 1 0 0 0]):datenum([0 0 1 0 0 0]):datenum([2010 2 28 0 0 0])]];

times = [];
times = [times, [datenum([2000 12 1 0 0 0]):datenum([0 0 0 6 0 0]):datenum([2000 12 31 18 0 0])]];
for yy=[2001:2009]
  times = [times, [datenum([yy 1 1 0 0 0]):datenum([0 0 0 6 0 0]):datenum([yy 2 28 18 0 0])]];
  times = [times, [datenum([yy 12 1 0 0 0]):datenum([0 0 0 6 0 0]):datenum([yy 12 31 18 0 0])]];
end
times = [times, [datenum([2010 1 1 0 0 0]):datenum([0 0 0 6 0 0]):datenum([2010 2 28 18 0 0])]];

p = nc_varget('/nfs/a137/earceb/ERA-interim/2012/ap/ggap201208070000.nc','p'); 
p = p.*100; % Work in SI units!
lat = nc_varget('/nfs/a137/earceb/ERA-interim/2012/as/ggas201208070000.nc','latitude');
lon = nc_varget('/nfs/a137/earceb/ERA-interim/2012/as/ggas201208070000.nc','longitude');
lt = find(lat>=-6 & lat<=8);
ln = find(lon>=96 & lon<=109);
lat = lat(lt);
lon = lon(ln);

% Initialise daily variables
MFD = zeros(4*length(dates),length(lt),length(ln));

% For each date and time read in data and compute moisture flux divergence
count = 0;
for pp = 1:length(dates)
%for pp = 1
  disp(['Reading data for ', datestr(dates(pp))])  
  
  for ff=1:4  
    count = count + 1;  
    
    time = datestr(dates(pp),30);
    rawfiles = dir(['/nfs/a137/earceb/ERA-interim/',time(1:4),'/ap/ggap',time(1:8),'*.nc']); 
    
    U = nc_varget(['/nfs/a137/earceb/ERA-interim/',time(1:4),'/ap/',rawfiles(ff).name],...
        'U',[1-1 1-1 lt(1)-1 ln(1)-1],[1 37 length(lt) length(ln)]);    
    V = nc_varget(['/nfs/a137/earceb/ERA-interim/',time(1:4),'/ap/',rawfiles(ff).name],...
        'V',[1-1 1-1 lt(1)-1 ln(1)-1],[1 37 length(lt) length(ln)]); 
    q = nc_varget(['/nfs/a137/earceb/ERA-interim/',time(1:4),'/ap/',rawfiles(ff).name],...
        'Q',[1-1 1-1 lt(1)-1 ln(1)-1],[1 37 length(lt) length(ln)]);  
    LWC = nc_varget(['/nfs/a137/earceb/ERA-interim/',time(1:4),'/ap/',rawfiles(ff).name],...
        'CLWC',[1-1 1-1 lt(1)-1 ln(1)-1],[1 37 length(lt) length(ln)]);
    IWC = nc_varget(['/nfs/a137/earceb/ERA-interim/',time(1:4),'/ap/',rawfiles(ff).name],...
        'CIWC',[1-1 1-1 lt(1)-1 ln(1)-1],[1 37 length(lt) length(ln)]);  
    q = q + LWC + IWC;    
    
    % Compute moisture flux (MF) by integrating q*u and q*v
    MFu = zeros(length(lt),length(ln));
    MFv = zeros(length(lt),length(ln));
    for ltt=1:length(lt)
      for lnn=1:length(ln)      
        MFu(ltt,lnn) = trapz(p,q(:,ltt,lnn).*U(:,ltt,lnn)).*-1.*(1/g);
        MFv(ltt,lnn) = trapz(p,q(:,ltt,lnn).*V(:,ltt,lnn)).*-1.*(1/g);      
      end
    end  
       
    % Compute moisture flux divergence (MFD)
    MFDtemp = zeros(length(lt),length(ln));
    for ltt=2:length(lt)-1
      for lnn=2:length(ln)-1
        dU = MFu(ltt,lnn+1) - MFu(ltt,lnn-1);
        dV = MFv(ltt+1,lnn) - MFv(ltt-1,lnn);        
        [dx,bearing]=latlon2dist(lat(ltt),lon(lnn+1),lat(ltt),lon(lnn-1));
        dx = dx.*1000; % Change from km to m
        [dy,bearing]=latlon2dist(lat(ltt+1),lon(lnn),lat(ltt-1),lon(lnn));
        dy = dy.*1000;
        MFDtemp(ltt,lnn) = (dU/dx) + (dV/dy); 
      end
    end  
    MFD(count,:,:) = MFDtemp;
    
  end  
end    

% Compute mean over each t and t-1 to get data at half hourly intervals
diffs = diff(MFD,1,1);
MFD_ERAI.MFD = MFD(1:length(times)-1,:,:) + (diffs./2);  

% Sort out time series (between the 00,06,12,18 UTC times)
MFD_ERAI.time = times(1:length(times)-1)+datenum([0 0 0 3 0 0]);
    
savelist = ['MFD_ERAI'];
printname = ['MJO_dc_ERAI_MFD.mat'];
eval(['save ',outdir,'/',printname,' ',savelist]);