% Equatorial 140W Line for Frank O. Bryan
% extract 06-Mar-2011 to 22-Apr-2013
% lats 8S to 8N, lons 140W, 0-500m
% (example extraction on face 4 and 5, i.e., rotated UV fields)

% define desired region
nx=2160;
prec='real*4';
region_name='Eq140W';
minlat=-8;
maxlat=8;
minlon=220;
maxlon=220+60/nx;
mints=dte2ts('06-Mar-2011',45,2011,1,17);
maxts=dte2ts('22-Apr-2013',45,2011,1,17);

% extract indices for desired region
gdir=['~dmenemen/llc_' int2str(nx) '/grid/'];
fnam=[gdir 'Depth.data'];
[fld fc ix jx] = ...
    quikread_llc(fnam,nx,1,prec,gdir,minlat,maxlat,minlon,maxlon);
quikpcolor(fld')
load ~dmenemen/llc_4320/grid/thk90.mat
kx=find(dpt90<=500);

% Get and save grid information
close all
pin=['~dmenemen/llc_' int2str(nx) '/grid/'];
pout=['~dmenemen/llc_' int2str(nx) '/regions/' region_name '/grid/'];
eval(['mkdir ' pout])
eval(['cd ' pout])
suf1=['_' int2str(length(ix)) 'x' int2str(length(jx))];
suf2=[suf1 'x' int2str(length(kx))];
for fnm={'Depth','RAC','XC','YC','hFacC'}
    fin=[pin fnm{1} '.data'];
    switch fnm{1}
      case{'hFacC'}
        fld=read_llc_fkij(fin,nx,fc,kx,ix,jx);
        fout=[fnm{1} suf2];
      otherwise
        fld=read_llc_fkij(fin,nx,fc,1,ix,jx);
        fout=[fnm{1} suf1];
    end
    writebin(fout,fld);
end

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

% get and save regional 2D fields
extent=[int2str(length(ix)) ',' int2str(length(jx)) ',1'];
for fnm={'Eta','PhiBot','KPPhbl','oceFWflx','oceQnet','oceQsw','oceSflux'}
    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

%  Get and save regional S/T/W
extent=[int2str(length(ix)) ',' int2str(length(jx)) ',' int2str(kx(end))];
for fnm={'Salt','Theta','W'}
    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

% 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((fc-2)*nx+min(ix)) ',' int2str(min(jx)) ',1 '];
extent=[int2str(length(ix)) ',' int2str(length(jx)) ',' 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((fc-2)*nx+min(ix)) ',' int2str(min(jx)-1) ',1 '];
extent=[int2str(length(ix)) ',' int2str(length(jx)) ',' 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(kx(end)) '_Neg/'  int2str(kx(end)) '/'' joblist']);
disp(['cd ' pout 'V'])
disp('parallel --slf $PBS_NODEFILE -j2 -a joblist')