function wb_4km_40day_MFD(outdir)

% Computes moisture flux divergence
%
% CEB Sept 2012

load /nfs/see-fs-02_users/earceb/AMMA/Cascade/cascade_4km_grid_coords.mat

g = 9.81;

filenames_qcf = dir(['/nfs/a90/earceb/CASCADE/netcdf/4km/','*d*_12.nc']);
filenames_qcl = dir(['/nfs/a90/earceb/CASCADE/netcdf/4km/','*d*_254.nc']);
filenames_q = dir(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pbfiles/','pbfiles*d1.nc']);
filenames_p = dir(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pbfiles/','pbfiles*c2.nc']);
filenames_u = dir(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pffiles/','pffiles*a1.nc']);
filenames_v = dir(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pffiles/','pffiles*a2.nc']);

% p index and qc index
ht37 = nc_varget('/nfs/see-fs-02_users/earceb/AMMA/Cascade/4km_exp_40day/pbfiles/pbfiles12d1.nc','hybrid_ht');
ht43 = nc_varget('/nfs/see-fs-02_users/earceb/AMMA/Cascade/4km_exp_40day/pbfiles/pbfiles38c2.nc','hybrid_ht');
ht70 = nc_varget('/nfs/a90/earceb/CASCADE/netcdf/4km/xfjeaa_d01_12.nc','hybrid_ht');

pindex = zeros(37,1);
for yy=1:length(ht37)
  ii = find(ht43==ht37(yy));
  pindex(yy) = ii;
end  
qcindex = zeros(37,1);
for yy=1:length(ht37)
  ii = find(ht70==ht37(yy));
  qcindex(yy) = ii;
end 

theta_levs = nc_varget('/nfs/cascade3/Cascade40/UM71km4/netcdf/pbfiles/pbfiles27c2.nc','hybrid_ht');
theta_levs = theta_levs(pindex);
rho_levs = nc_varget('/nfs/cascade3/Cascade40/UM71km4/netcdf/pbfiles/pbfiles27c1.nc','hybrid_ht');
rho_levs = rho_levs(pindex);

MFDall = zeros(240*4,776,1110);
MFDall_monsoon = zeros(240*4,776,1110);
MFDall_AEJ = zeros(240*4,776,1110);
time = zeros(240*4,1);
tcount = 0;

% Read in data ------------------------------------------------------------
for ff = 1:length(filenames_q)
    
%   % This makes sure a continuous data set is picked out where job IDs change  
%   t2 = nc_varget(['/nfs/a90/earceb/CASCADE/netcdf/12km_exp_40day_SM/',filenames_b(ff).name],'t'); 
%   [m] = length(t2);
%   if m==5
%     st = 2; fn = 5;    
%   elseif m==4
%     st = 1; fn = 4;  
%   else
%     disp('Warning!!! time series if not 4 or 5 long!!')
%     size(t2)
%     return
%   end  
    
  for tt=1:24   

    tcount = tcount+1  
        
    % Read in data
    disp('Reading data')
    t = nc_varget(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pbfiles/',filenames_q(ff).name],'t',[tt-1],[1]); 
    if ff<=10
      tm2add = datenum([2006 7 25 0 0 0]);
    elseif ff>10 & ff<=20
      tm2add = datenum([2006 8 4 0 0 0]);   
    elseif ff>20 & ff<=30
      tm2add = datenum([2006 8 14 0 0 0]);      
    elseif ff>30 & ff<=40  
      tm2add = datenum([2006 8 24 0 0 0]);
    else
      disp(['Warning ',num2str(ff),' ',num2str(tt)])    
    end
    
    time(tcount) = t + tm2add;
    disp('p')    
    p = nc_varget(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pbfiles/',filenames_p(ff).name],'air_pressure',[tt-1 0 0 0],[1 43 776 1110]);  
    p = p(pindex,:,:);
    %p = p./100; % Work in Pascals! SI units!      
    disp('q')     
    q = nc_varget(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pbfiles/',filenames_q(ff).name],'specific_humidity',[tt-1 0 0 0],[1 37 776 1110]); 
    disp('qcl and qcf')    
    qcl = nc_varget(['/nfs/a90/earceb/CASCADE/netcdf/4km/',filenames_qcl(ff).name],'mass_fraction_of_cloud_liquid_water_in_air',[tt-1 0 0 0],[1 70 776 1110]);
    qcl = qcl(qcindex,:,:);
    qcf = nc_varget(['/nfs/a90/earceb/CASCADE/netcdf/4km/',filenames_qcf(ff).name],'mass_fraction_of_cloud_ice_in_air',[tt-1 0 0 0],[1 70 776 1110]);
    qcf = qcf(qcindex,:,:);
    q = q + qcl + qcf; % Use total water, not just vapour
    disp('winds')        
    u = nc_varget(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pffiles/',filenames_u(ff).name],'x-wind',[tt-1 0 0 0],[1 37 776 1110]);
    v = nc_varget(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pffiles/',filenames_v(ff).name],'y-wind',[tt-1 0 0 0],[1 37 775 1110]);
    
    % Interpolate v onto same grid as u, p and q
    disp('Interpolate v onto same grid as u, p and q')
    x = nc_varget(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pffiles/',filenames_u(ff).name],'x');
    x_1 = nc_varget(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pffiles/',filenames_v(ff).name],'x_1');
    y = nc_varget(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pffiles/',filenames_u(ff).name],'y');
    y_1 = nc_varget(['/nfs/cascade3/Cascade40/UM71km4/netcdf/pffiles/',filenames_v(ff).name],'y_1');  
    % Interpolation (table lookup)
    v_int = zeros(37,774,1109);
    for zz=1:37  
      v_int(zz,:,:) = interp2(x_1',y_1,squeeze(v(zz,:,:)),x(1:1109)',y(2:775));
    end
   
    % Put into matrix the same size as u, p and q (missing values at the
    % edge for zeros)
    v = zeros(37,776,1110);
    v(:,2:775,1:1109) = v_int;
    
    % Interpolate u and v onto theta model levels and
    % Compute moisture flux (MF) by integrating q*u and q*v
    disp('Interpolate u and v onto theta model levels and compute MF')
    u_theta = zeros(37,776,1110);
    v_theta = zeros(37,776,1110); 
    MFu = zeros(776,1110);
    MFv = zeros(776,1110);
    MFu_monsoon = zeros(776,1110);
    MFv_monsoon = zeros(776,1110);    
    MFu_AEJ = zeros(776,1110);
    MFv_AEJ = zeros(776,1110);     
    for lt=1:776
      for ln=1:1110    
        % Interp u and v  
        u_theta(:,lt,ln) = interp1(rho_levs,u(:,lt,ln),theta_levs,'linear','extrap');
        v_theta(:,lt,ln) = interp1(rho_levs,v(:,lt,ln),theta_levs,'linear','extrap');
        % Compute MF  
        MFu(lt,ln) = trapz(p(:,lt,ln),q(:,lt,ln).*u_theta(:,lt,ln)).*-1.*(1/g);
        MFv(lt,ln) = trapz(p(:,lt,ln),q(:,lt,ln).*v_theta(:,lt,ln)).*-1.*(1/g);   
        MFu_monsoon(lt,ln) = trapz(p(1:18,lt,ln),q(1:18,lt,ln).*u_theta(1:18,lt,ln)).*-1.*(1/g);
        MFv_monsoon(lt,ln) = trapz(p(1:18,lt,ln),q(1:18,lt,ln).*v_theta(1:18,lt,ln)).*-1.*(1/g);
        MFu_AEJ(lt,ln) = trapz(p(20:24,lt,ln),q(20:24,lt,ln).*u_theta(20:24,lt,ln)).*-1.*(1/g);
        MFv_AEJ(lt,ln) = trapz(p(20:24,lt,ln),q(20:24,lt,ln).*v_theta(20:24,lt,ln)).*-1.*(1/g);        
      end
    end  
       
    % Compute moisture flux divergence (MFD)
    disp('Compute moisture flux divergence (MFD)') 
    MFDtemp = zeros(776,1110);
    MFDtemp_monsoon = zeros(776,1110);  
    MFDtemp_AEJ = zeros(776,1110);    
    for lt=2:775
      for ln=2:1109
        dU = MFu(lt,ln+1) - MFu(lt,ln-1);
        dV = MFv(lt+1,ln) - MFv(lt-1,ln); 
        dU_monsoon = MFu_monsoon(lt,ln+1) - MFu_monsoon(lt,ln-1);
        dV_monsoon = MFv_monsoon(lt+1,ln) - MFv_monsoon(lt-1,ln);        
        dU_AEJ = MFu_AEJ(lt,ln+1) - MFu_AEJ(lt,ln-1);
        dV_AEJ = MFv_AEJ(lt+1,ln) - MFv_AEJ(lt-1,ln);        
        [dx,bearing]=latlon2dist(lat2d(lt,ln+1),lon2d(lt,ln+1),lat2d(lt,ln-1),lon2d(lt,ln-1));
        dx = dx.*1000; % Change from km to m
        [dy,bearing]=latlon2dist(lat2d(lt+1,ln),lon2d(lt+1,ln),lat2d(lt-1,ln),lon2d(lt-1,ln));
        dy = dy.*1000;
        MFDtemp(lt,ln) = (dU/dx) + (dV/dy);
        MFDtemp_monsoon(lt,ln) = (dU_monsoon/dx) + (dV_monsoon/dy);
        MFDtemp_AEJ(lt,ln) = (dU_AEJ/dx) + (dV_AEJ/dy);        
      end
    end  
    MFDall(tcount,:,:) = MFDtemp;
    MFDall_monsoon(tcount,:,:) = MFDtemp_monsoon;
    MFDall_AEJ(tcount,:,:) = MFDtemp_AEJ;    
    
   % Save out file for each timestep 
   timee = time(tcount); 
   savelist = ['MFDtemp MFDtemp_monsoon MFDtemp_AEJ lat2d lon2d timee'];
   printname = [['cascade_4km_MFD_',datestr(timee,30),'.mat']];
   eval(['save ',outdir,'/',printname,' ',savelist]);
    
  end  

end  
% -------------------------------------------------------------------------


% Compute mean over each t and t-1 to get data at half hourly intervals
diffs = diff(MFDall,1,1);
MFD_40day_4km.MFD = MFDall(1:959,:,:) + (diffs./2);  
diffs = diff(MFDall_monsoon,1,1);
MFD_40day_4km.MFD_monsoon = MFDall_monsoon(1:959,:,:) + (diffs./2); 
diffs = diff(MFDall_AEJ,1,1);
MFD_40day_4km.MFD_AEJ = MFDall_AEJ(1:959,:,:) + (diffs./2); 

% Sort out time series (on the half hour to match PWVt, E and P)
MFD_40day_4km.time = time(1:959)+datenum([0 0 0 0 30 0]);

MFD_40day_4km.lat2d = lat2d;
MFD_40day_4km.lon2d = lon2d;
MFD_40day_4km.lat = lat2d(:,555);
MFD_40day_4km.lon = lon2d(388,:);

savelist = ['MFD_40day_4km'];
printname = ['cascade_4km_40day_MFD.mat'];
eval(['save ',outdir,'/',printname,' -v7.3 ',savelist]);