% 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/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;

%cancel at ts =
%
%       643200
% lat
%
% lat =
%
%   -59.5000
%
% lon
%
% lon =
%
%    67.5000

for lat = -41.5:3:48.5
    for lon = 1.5:3:358.5
        % extract indices for desired region
        disp(num2str(lat));
        disp(num2str(lon));
        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) % this makes sure I don't cross faces.. just leave those out
            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','oceTAUX','oceTAUY'}
                    fnm{1}
                    eval(['mkdir ' pout fnm{1}])
                    eval(['cd ' pout fnm{1}])
                    n = 0;
                    fld = 0;
                    fout=[fnm{1} '_winter_average_' int2str(length(ix)) 'x' int2str(length(jx)) ...
                        'x' int2str(length(kx))];
                    for ts=92160+((90+6*30)*80*24):(80*24):92160+((90+6*30+90)*80*24);
                        ts
                        if (ts<1198080)
                            pin=[pn 'run_day49_624/'];
                        else
                            pin=[pn 'run/'];
                        end
                        
                        fin=[pin fnm{1} '.' myint2str(ts,10) '.data'];
                        if exist(fin)
                            n = n + 1;
                            no = 0;
                            %dy=ts2dte(ts,45,2011,1,17,30);
                            
                            if strcmp('Eta',fnm{1}) || strcmp('oceTAUY',fnm{1}) || strcmp('oceTAUX',fnm{1})
                                tmp1 = read_llc_fkij(fin,nx,fc,1,ix,jx);
                                if sum(sum(tmp1)) == 0
                                    n = n-1;
                                else
                                    fld= fld + tmp1;
                                end
                            else
                                for k=1:length(kx);
                                    tmp = read_llc_fkij(fin,nx,fc,kx(k),ix,jx);
                                    if n == 1 && k == 1
                                        fld = zeros([size(tmp),22]);
                                        tmpT = fld;
                                    end
                                    if sum(sum(tmp)) == 0;
                                        no = 1;
                                        n = n-1;
                                        break;
                                    else
                                        tmpT(:,:,k) = tmp;
                                    end

                                end
                                if ~no
                                 fld= fld + tmpT;
                                end
                            end
                        end
                        
                    end
                    n
                    writebin(fout,fld./n,1,'real*4');
                end
                
                
            end
        end
    end
end