% SPURS-2 region for Kyla Drushka and Luc Rainville
% extract Eta, PhiBot, Salt, and Theta
% for 13-Sep-2011 to 15-Nov-2012
% lats 5N to 15N, lons 230E to 240E
% (example extraction on face 4 and 5, i.e., rotated UV fields)

% define desired region
nx=4320;
prec='real*4';
region_name='SPURS2';
minlat=5;
maxlat=15;
minlon=230;
maxlon=240;
mints=dte2ts('13-Sep-2011',25,2011,9,10);
maxts=dte2ts('15-Nov-2012',25,2011,9,10);

% extract indices for desired region
gdir='~dmenemen/llc_4320/grid/';
fnam=[gdir 'Depth.data'];
[tmp fc ix jx] = ...
    quikread_llc(fnam,nx,1,prec,gdir,minlat,maxlat,minlon,maxlon);
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))),:)=tmp{fc(f)};
end
quikpcolor(fld')
load ~dmenemen/llc_4320/grid/thk90.mat
bot90=dpt90+thk90/2;
kx=1:min(find(bot90>mmax(fld)));

% Get and save grid information
close all
pin='~dmenemen/llc_4320/grid/';
pout=['~dmenemen/llc_4320/regions/' region_name '/grid/'];
eval(['mkdir ' pout])
eval(['cd ' pout])
suf1=['_' int2str(sum(m)) 'x' int2str(n)];
suf2=[suf1 'x' int2str(length(kx))];

% Grid cell center
for fnm={'Depth','RAC','XC','YC','hFacC'}
    fin=[pin fnm{1} '.data'];
    switch fnm{1}
      case{'hFacC'}
        fout=[fnm{1} suf2];
        fld=zeros(sum(m),n,length(kx));
        for f=1:length(fc)
            fld((sum(m(1:f))+1):sum(m(1:(f+1))),:,:) = ...
                read_llc_fkij(fin,nx,fc(f),kx,ix{fc(f)},jx{fc(f)});
        end
      otherwise
        fout=[fnm{1} suf1];
        fld=zeros(sum(m),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
    end
    writebin(fout,fld);
end

% Southwest corner (vorticity) points, no direction
fld=zeros(sum(m),n);
for fnm={'XG','YG','RAZ'}
    fin=[pin fnm{1} '.data'];
    fout=[fnm{1} suf1];
    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
    writebin(fout,fld);
end

% West edge points, no direction
fnx='DXC';
fny='DYC';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx suf1];
fouty=[fny suf1];
fldx=zeros(sum(m),n);
fldy=zeros(sum(m),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 corner (vorticity) points, no direction
fnx='DXV';
fny='DYU';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx suf1];
fouty=[fny suf1];
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)}-1); % <<<<<<<<
    end
end
writebin(foutx,fldx);
writebin(fouty,fldy);

% Southwest edge points, no direction
fnx='DXG';
fny='DYG';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx suf1];
fouty=[fny suf1];
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
writebin(foutx,fldx);
writebin(fouty,fldy);

% create commands for extracting model output fields
pout=['~dmenemen/llc_4320/regions/' region_name '/'];
extract='/home4/bcnelson/MITgcm/extract/v1.10/extract4320 -g ';
timesteps=[int2str(mints) '-' int2str(maxts) ' '];
startPoint=[int2str(2*nx+min(ix{4})) ',' int2str(min(jx{4})) ',1 '];

% get and save regional 2D tracer fields
extent=[int2str(sum(m)) ',' int2str(n) ',1'];
for fnm={'Eta','PhiBot','KPPhbl','oceFWflx','oceQnet','oceQsw','oceSflux'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    disp(['cd ' pout fnm{1}])
    fieldNames=[fnm{1} ' '];
    system([extract timesteps  fieldNames  startPoint  extent '  > joblist']);
    disp('parallel --slf $PBS_NODEFILE -j2 -a joblist')
end

% Get and save regional S/T/W
extent=[int2str(sum(m)) ',' int2str(n) ',' int2str(kx(end))];
for fnm={'Salt','Theta','W'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    disp(['cd ' pout fnm{1}])
    fieldNames=[fnm{1} ' '];
    system([extract timesteps  fieldNames  startPoint  extent '  > joblist']);
    disp('parallel --slf $PBS_NODEFILE -j2 -a joblist')
end

% 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
% there is a -1 index offset -U meridional velocity in faces 4/5
% example below is for face 4

% U
startPoint=[int2str(2*nx+min(ix{4})) ',' int2str(min(jx{4})) ',1 '];
extent=[int2str(sum(m)) ',' int2str(n) ',' int2str(kx(end))];
eval(['mkdir ' pout 'U'])
eval(['cd ' pout 'U'])
fieldNames=['V' ' '];
system([extract timesteps  fieldNames  startPoint  extent ' > joblist']);
system('sed -i ''s/_V_/_U_/'' joblist');
disp(['cd ' pout 'U'])
disp('parallel --slf $PBS_NODEFILE -j2 -a joblist')

% V
startPoint=[int2str(2*nx+min(ix{4})) ',' int2str(min(jx{4})-1) ',1 '];
extent=[int2str(sum(m)) ',' int2str(n) ',' int2str(kx(end))];
eval(['mkdir ' pout 'V'])
eval(['cd ' pout 'V'])
fieldNames=['U' ' '];
system([extract '-n ' timesteps  fieldNames  startPoint  extent ' > joblist']);
system('sed -i ''s/_U_/_V_/'' joblist');
system(['sed -i ''s/' int2str(length(kx)) '_Neg/' int2str(length(kx)) '/'' joblist']);
disp(['cd ' pout 'V'])
disp('parallel --slf $PBS_NODEFILE -j2 -a joblist')