% define desired region
clear
nx=2160;
region_name='DopplerScat/latlon';
minlat=-70;
maxlat=57;

% extract indices for desired region
prec='real*4';
gdir='/nobackupp9/dmenemen/llc_2160/grid/';
fnam=[gdir 'Depth.data'];
[tmp fc ix jx]=quikread_llc(fnam,nx,1,prec,gdir,minlat,maxlat);
jx{1} = 1:6480;
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),6480); 
tmp{1} = [tmp{1},zeros(2160,6480-length(tmp{1}))];
for f=1:length(fc)
    %fld((sum(m(1:f))+1):sum(m(1:(f+1))),1:length(tmp{f}))=tmp{fc(f)};
   
    fld = [tmp{1};tmp{2};tmp{4};tmp{5}];
end
quikpcolor(fld')

% get and save grid information
pin='/nobackupp9/dmenemen/llc_2160/grid/';
pout=['/nobackupp9/dmenemen/llc_2160/regions/' region_name '/grid/'];
eval(['mkdir ' pout])
eval(['cd ' pout])
for fnm={'AngleCS','AngleSN','DXC','DXG','DYC','DYG','Depth', ...
         'RAC','RAS','RAW','RAZ','U2zonDir','V2zonDir', ...
         'XC','XG','YC','YG','hFacC','hFacS','hFacW'}
    fin=[pin fnm{1} '.data'];
    fout=[fnm{1} '_' 10800 'x' 6480];
    for f=1:length(fc)
        if f ~=3
            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

% 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
pn='/nobackupp9/dmenemen/llc_2160/MITgcm/';
pout=['/nobackupp9/dmenemen/llc_2160/regions/' region_name '/'];
kx=1:22;
eval(['mkdir ' pout 'U'])
eval(['mkdir ' pout 'V'])
eval(['cd ' pout])

for ts=92160+((90+6*30)*80*24):(80*24):92160+((90+6*30+90)*80*24);
    
    if (ts<1198080)
        pin=[pn 'run_day49_624/'];
    else
        pin=[pn 'run/'];
    end
    
    
    dy=ts2dte(ts,45,2011,1,17,30);

    finu=[pin 'U.' myint2str(ts,10) '.data'];
    finv=[pin 'V.' myint2str(ts,10) '.data'];
    dy=ts2dte(ts,25,2011,9,10,30);
    foutu=['U/U_' 10800 'x' 6480 'x' int2str(length(kx)) '.' dy];
    foutv=['V/V_' 10800 'x' 6480 'x' int2str(length(kx)) '.' dy];
    for k=1:length(kx); mydisp(k)
        for f=1:length(fc)
            switch fc(f)
              case {1,2}
                fldu(1+(f-1)*2160:(f-1)*2160+2160,:,k) = ...
                    read_llc_fkij(finu,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
                fldv(1+(f-1)*2160:(f-1)*2160+2160,:,k) = ...
                    read_llc_fkij(finv,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
              case {4,5}
                fldu(1+(f-2)*2160:(f-2)*2160+2160,:,k) = ...
                    read_llc_fkij(finv,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
                fldv(1+(f-2)*2160:(f-2)*2160+2160,:,k) = - ...
                    read_llc_fkij(finu,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
            end
        end
        writebin(foutu,fldu,1,'real*4',k-1);
        writebin(foutv,fldv,1,'real*4',k-1);
    end
end

pn='/nobackupp9/dmenemen/llc_2160/MITgcm/';
pout=['/nobackupp9/dmenemen/llc_2160/regions/' region_name '/'];
kx=1:1;
eval(['mkdir ' pout 'Eta'])
eval(['mkdir ' pout 'oceTAUX'])
eval(['mkdir ' pout 'oceTAUY'])
eval(['cd ' pout])

for ts=92160+((90+6*30)*80*24):(80*24):92160+((90+6*30+90)*80*24);
    
    if (ts<1198080)
        pin=[pn 'run_day49_624/'];
    else
        pin=[pn 'run/'];
    end
    
    
    dy=ts2dte(ts,45,2011,1,17,30);

    finu=[pin 'oceTAUX.' myint2str(ts,10) '.data'];
    finv=[pin 'oceTAUY.' myint2str(ts,10) '.data'];
    dy=ts2dte(ts,25,2011,9,10,30);
    foutu=['oceTAUX/oceTAUX_' 10800 'x' 6480 'x' int2str(length(kx)) '.' dy];
    foutv=['oceTAUY/oceTAUY_' 10800 'x' 6480 'x' int2str(length(kx)) '.' dy];
    for k=1:length(kx); mydisp(k)
        for f=1:length(fc)
            switch fc(f)
              case {1,2}
                fldu(1+(f-1)*2160:(f-1)*2160+2160,:,k) = ...
                    read_llc_fkij(finu,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
                fldv(1+(f-1)*2160:(f-1)*2160+2160,:,k) = ...
                    read_llc_fkij(finv,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
              case {4,5}
                fldu(1+(f-2)*2160:(f-2)*2160+2160,:,k) = ...
                    read_llc_fkij(finv,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
                fldv(1+(f-2)*2160:(f-2)*2160+2160,:,k) = - ...
                    read_llc_fkij(finu,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
            end
        end
        writebin(foutu,fldu,1,'real*4',k-1);
        writebin(foutv,fldv,1,'real*4',k-1);
    end
end

pn='/nobackupp9/dmenemen/llc_2160/MITgcm/';
pout=['/nobackupp9/dmenemen/llc_2160/regions/' region_name '/'];
kx=1:1;
eval(['mkdir ' pout 'Eta'])

for ts=92160+((90+6*30)*80*24):(80*24):92160+((90+6*30+90)*80*24);
    
    if (ts<1198080)
        pin=[pn 'run_day49_624/'];
    else
        pin=[pn 'run/'];
    end
    
    
    dy=ts2dte(ts,45,2011,1,17,30);

    fin=[pin 'Eta.' myint2str(ts,10) '.data'];
   
    dy=ts2dte(ts,25,2011,9,10,30);
    fout=['Eta/Eta_' 10800 'x' 6480 'x' int2str(length(kx)) '.' dy];
    for k=1:length(kx); mydisp(k)
        for f=1:length(fc)
            switch fc(f)
              case {1,2}
                fld(1+(f-1)*2160:(f-1)*2160+2160,:,k) = ...
                    read_llc_fkij(fin,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
             
              case {4,5}
                fld(1+(f-2)*2160:(f-2)*2160+2160,:,k) = ...
                    read_llc_fkij(fin,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
                
            end
        end
        writebin(fout,fld,1,'real*4',k-1);
      
    end
end