% 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:27;
% 
% nx=2160;
% region_name='DopplerScat/southern_BIG';
% minlat=-70;
% maxlat=57;
% minlon=0;
% maxlon=360;
% kx=1:27;

%
% nx=2160;
% region_name='DopplerScat/subpolar';
% minlat=45-1.5;
% maxlat=45+1.5;
% minlon=180-1.5;
% maxlon=180+1.5;
% kx=1:27;
%
% nx=2160;
% region_name='DopplerScat/southern';
% minlat=-50-1.5;
% maxlat=-50+1.5;
% minlon=90-1.5;
% maxlon=90+1.5;
% kx=1:27;

for lat = -68.5:3:48.5
    for lon = 1.5:3:358.5
        % extract indices for desired region
        
        nx=2160;
        if lat<0
            region_name=['DopplerScat/latlon/n_' num2str(abs(lat)) '_' num2str(lon)];
        else
            region_name=['DopplerScat/latlon/p_' num2str(lat) '_' num2str(lon)];
        end
        
        minlat=lat-1.5;
        maxlat=lat+1.5;
        minlon=lon-1.5;
        maxlon=lon+1.5;
        kx=1:22;
        
        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
        if ~iscell(fld)
            if fld > 200 
            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+((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
                    
                    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+((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
                    
                    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
        end
        end
    end
end