% Average 3d model velocity fields


clear
%define desired region
nx=2160;
% region_name='DopplerScat/subtropic';
% minlat=13.5;
% maxlat=16.5;
% minlon=-131.5;
% maxlon=-128.5;
% kx=1:1;
%
% nx=2160;
% region_name='DopplerScat/subpolar';
% minlat=45-1.5;
% maxlat=45+1.5;
% minlon=180-1.5;
% maxlon=180+1.5;
% kx=1:1;
% %
nx=2160;
region_name='DopplerScat/southern';
minlat=-50-1.5;
maxlat=-50+1.5;
minlon=90-1.5;
maxlon=90+1.5;
kx=1:1;

nx=2160;

prec='real*4';
gdir='/nobackupp5/dmenemen/llc_2160/grid/';
fnam=[gdir 'Depth.data'];
[fld fc ix jx]=quikread_llc(fnam,nx,1,prec,gdir,minlat,maxlat,minlon,maxlon);
%quikpcolor(fld')
%thincolorbar

eval('cd ~dmenemen/llc_2160/regions/DopplerScat/latlon/')
eval(['mkdir ' region_name])
% 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} '_' int2str(length(ix)) 'x' int2str(length(jx))];
    switch fnm{1}
        case{'hFacC','hFacS','hFacW'}
            fld=read_llc_fkij(fin,nx,fc,kx,ix,jx);
            fout=[fout 'x' int2str(length(kx))];
        otherwise
            fld=read_llc_fkij(fin,nx,fc,1,ix,jx);
    end
    writebin(fout,fld);
end

% get and save regional S/T/U/V/W
pn='/nobackupp9/dmenemen/llc_2160/MITgcm/';
pout=['/nobackupp9/dmenemen/llc_2160/regions/' region_name '/'];
for fnm={'Eta','U','V'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=92160:(80*24):1586400;
        
        if (ts<1198080)
            pin=[pn 'run_day49_624/'];
        else
            pin=[pn 'run/'];
        end
        
        fin=[pin fnm{1} '.' myint2str(ts,10) '.data'];
        dy=ts2dte(ts,45,2011,1,17,30);
        fout=[fnm{1} '_' int2str(length(ix)) 'x' int2str(length(jx)) ...
            'x' int2str(length(kx)) '.' dy];
        if strcmp('Eta',fnm{1})
            fld=read_llc_fkij(fin,nx,fc,1,ix,jx);
            writebin(fout,fld,1,'real*4');
        else
            for k=1:length(kx);
                fld=read_llc_fkij(fin,nx,fc,kx(k),ix,jx);
                writebin(fout,fld,1,'real*4',k-1);
            end
        end
    end
end

% get and save regional surface fluxes
for fnm={'oceTAUX','oceTAUY'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=92160:(80*24):1586400;
        
        
            if (ts<1198080)
                pin=[pn 'run_day49_624/'];
            else
                pin=[pn 'run/'];
            end
            
            fin=[pin fnm{1} '.' myint2str(ts,10) '.data'];
            dy=ts2dte(ts,45,2011,1,17,30);
            fout=[fnm{1} '_' int2str(length(ix)) 'x' int2str(length(jx)) '.' dy];
            fld=read_llc_fkij(fin,nx,fc,1,ix,jx);
            writebin(fout,fld);
        
    end
end