% define desired region
nx=2160;
region_name='latlon';
minlat=-70;
maxlat=57;
mints=92160;
maxts=1586400;

% extract indices for desired region
prec='real*4';
gdir='~dmenemen/llc_2160/grid/';
fnam=[gdir 'YC.data'];
tmp=read_llc_fkij(fnam,nx,1,1,1);
tmp(~tmp)=nan;
fc=[1 2 4 5];
ix{1}=1:nx;
ix{2}=1:nx;
ix{3}=[];
ix{4}=1:nx;
ix{5}=1:nx;
jx{1}=find(tmp>=minlat&tmp<=maxlat);
jx{2}=find(tmp>=minlat&tmp<=maxlat);
jx{3}=[];
jx{4}=find(tmp>=minlat&tmp<=maxlat);
jx{5}=find(tmp>=minlat&tmp<=maxlat);
m(1)=0;
for f=1:length(fc)
    m(f+1)=length(ix{fc(f)});
end
n=length(jx{fc(1)});
fld=zeros(sum(m),n);
for f=1:length(fc)
    fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
        read_llc_fkij(fnam,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
end
%clf, quikpcolor(fld'), thincolorbar
%close all

% get and save grid information
pin='~dmenemen/llc_2160/grid/';
pout=['~dmenemen/llc_2160/regions/' region_name '/grid/'];
eval(['mkdir ' pout])
eval(['cd ' pout])
for fnm={'AngleCS','AngleSN','Depth','RAC', ...
         'U2zonDir','V2zonDir','XC','XG','YC','YG','hFacC'}
    fin=[pin fnm{1} '.data'];
    fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n)];
    for f=1:length(fc)
        fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(fin,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
    end
    switch fnm{1}
      case {'AngleCS','U2zonDir'}, fld(:)=1;
      case {'AngleSN','V2zonDir'}, fld(:)=0;
      case 'RAC'
        fld2=fld;
        for j=1:(nx*2)
            tmp=double(fld(:,j));
            fld2(:,j)=mean(tmp(find(tmp)));
        end
        fld=fld2;
        case {'XC','XG'}, fld=fld(:,500)*ones(1,2*nx);
        case {'YC','YG'}, fld=ones(4*nx,1)*fld(500,:);
    end
    writebin(fout,fld);
end

% Southwest center points, no direction, fill blank tiles
fnx='DXC';
fny='DYC';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx '_' int2str(sum(m)) 'x' int2str(n)];
fouty=[fny '_' int2str(sum(m)) 'x' int2str(n)];
for f=1:length(fc)
    switch fc(f)
      case {1,2}
        fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
        fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
      case {4,5}
        fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
        fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)}-1); % <<<<<<<<
    end
end
fldx=ones(4*nx,1)*fldx(1,:);
fldy=ones(4*nx,1)*fldy(1,:);
writebin(foutx,fldx);
writebin(fouty,fldy);

% Southwest center points, no direction, fill blank tiles
fnx='RAW';
fny='RAS';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx '_' int2str(sum(m)) 'x' int2str(n)];
fouty=[fny '_' int2str(sum(m)) 'x' int2str(n)];
for f=1:length(fc)
    switch fc(f)
      case {1,2}
        fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
        fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
      case {4,5}
        fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
        fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)}-1); % <<<<<<<<
    end
end
fldx2=fldx;
for j=1:(nx*2)
    tmp=double(fldx(:,j));
    fldx2(:,j)=mean(tmp(find(tmp)));
end
fldx=fldx2;
fldy2=fldy;
for j=1:(nx*2)
    tmp=double(fldy(:,j));
    fldy2(:,j)=mean(tmp(find(tmp)));
end
fldy=fldy2;
writebin(foutx,fldx);
writebin(fouty,fldy);

% Southwest center points, no direction
fnx='hFacW';
fny='hFacS';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx '_' int2str(sum(m)) 'x' int2str(n)];
fouty=[fny '_' int2str(sum(m)) 'x' int2str(n)];
for f=1:length(fc)
    switch fc(f)
      case {1,2}
        fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
        fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
      case {4,5}
        fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
        fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)}-1); % <<<<<<<<
    end
end
writebin(foutx,fldx);
writebin(fouty,fldy);

% Southwest edges, no direction, fill blank tiles
fnx='DXG';
fny='DYG';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx '_' int2str(sum(m)) 'x' int2str(n)];
fouty=[fny '_' int2str(sum(m)) 'x' int2str(n)];
for f=1:length(fc)
    switch fc(f)
      case {1,2}
        fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
        fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
      case {4,5}
        fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)}-1); % <<<<<<<<
        fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
    end
end
fldx=ones(4*nx,1)*fldx(1,:);
fldy=ones(4*nx,1)*fldy(1,:);
writebin(foutx,fldx);
writebin(fouty,fldy);

% Southwest corner (vorticity) points, no direction, fill blank tiles
fnm='RAZ';
fin=[pin fnm '.data'];
finy=[pin fny '.data'];
fout=[fnm '_' int2str(sum(m)) 'x' int2str(n)];
for f=1:length(fc)
    switch fc(f)
      case {1,2}
        fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(fin,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
      case {4,5}
        fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(fin,nx,fc(f),1,ix{fc(f)},jx{fc(f)}-1); % <<<<<<<<
    end
end
fld=ones(4*nx,1)*fld(1,:);
writebin(fout,fld);

% get and save regional fields of U and V
% note that zonal velocity is U in faces 1/2 and V in faces 4/5
% and meridional velocity is V in faces 1/2 and -U in faces 4/5
ts=mints:80:maxts;
pn='~dmenemen/llc_2160/MITgcm/run';
pout=['~dmenemen/llc_2160/regions/' region_name '/'];
kx=1:1;
eval(['mkdir ' pout 'U'])
eval(['mkdir ' pout 'V'])
eval(['cd ' pout])
for t=1:length(ts), mydisp(t)
    if (ts(t)<1198080)
        pin=[pn '_day49_624/'];
    else
        pin=[pn '/'];
    end
    finu=[pin 'U.' myint2str(ts(t),10) '.data'];
    finv=[pin 'V.' myint2str(ts(t),10) '.data'];
    dy=ts2dte(ts(t),45,2011,1,17,30);
    foutu=['U/U_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) '.' dy];
    foutv=['V/V_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) '.' dy];
    for k=1:length(kx);
        for f=1:length(fc)
            switch fc(f)
              case {1,2}
                fldu((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                    read_llc_fkij(finu,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
                fldv((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                    read_llc_fkij(finv,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
              case {4,5}
                fldu((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                    read_llc_fkij(finv,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
                fldv((sum(m(1:f))+1):sum(m(1:(f+1))),:) = - ...               % <<<<<<<<
                    read_llc_fkij(finu,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)}-1); % <<<<<<<<
            end
        end
        writebin(foutu,fldu,1,'real*4',k-1);
        writebin(foutv,fldv,1,'real*4',k-1);
    end
end

% create commands for extracting model output fields
pout=['~dmenemen/llc_2160/regions/' region_name '/'];
extract='/home4/bcnelson/MITgcm/extract/v1.10/extract2160 -g ';
timesteps=[int2str(mints) '-' int2str(maxts) ' '];
startPoint='1,1441,1 ';
extent='8640,4320,1';
for fnm={'Salt','Theta'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    fieldNames=[fnm{1} ' '];
    system([extract timesteps  fieldNames  startPoint  extent '  > joblist']);
    disp(['cd ' pout fnm{1}])
    disp('parallel --slf $PBS_NODEFILE -j2 -a joblist')
end