clc
clear

region_name='global';
nx=2160;


% get and save daily-averaged regional fields
pn='~dmenemen/llc_2160/MITgcm/';
pout=['~dmenemen/llc_2160/regions/' region_name '/'];

eval(['mkdir ' pout  'dopplerSCAT/'])
%eval(['cd ' pout 'dopplerSCAT'])
[dx,dy] = quikreadpcolor_dxdy_llc('~dmenemen/llc_4320/grid/DXC.data','~dmenemen/llc_4320/grid/DYC.data',nx);

for  ts=92160+(80*24):1586400;
    
    if (ts<1198080)
        pin=[pn 'run_day49_624/'];
    else
        pin=[pn 'run/'];
    end
    
    dy=ts2dte(ts,45,2011,1,17,30);
     
    fld=0; nn=0;
    
    finU=[pin 'U' '.' myint2str(ts,10) '.data'];
    finV=[pin 'V' '.' myint2str(ts,10) '.data'];
    finTx=[pin 'oceTAUX' '.' myint2str(ts,10) '.data'];
    finTy=[pin 'oceTAUY' '.' myint2str(ts,10) '.data'];
    
    if exist(finU) && exist(finV) && exist(finTx) && exist(finTy)
        [u, v] = quikreadpcolor_uv_llc(finU,finV,nx);
        [taux, tauy] = quikreadpcolor_uv_llc(finTx,finTy,nx);
        
        if sum(sum(u(~isnan(u)))) == 0 || sum(sum(v(~isnan(v)))) == 0 || sum(sum(taux(~isnan(taux)))) == 0 || sum(sum(tauy(~isnan(tauy)))) == 0
            warning(['ZEROs FILE ' finU])
            break
        else
            
            % do computations
            % Appy 10m parzen filter
            
            minRes = max(max(dy(~isnan(dy)&dy>0)));
 
            filt_wv = '10km';
            L = round(9500/3000); % NOTE: this isn't a very good way to do this.
                                    % numbers found by plotting L vs
                                    % filt_wv, multiplying L by the
                                    % resolution originally used (500m) to
                                    % generate the L and wv numbers. Then,
                                    % 10km was found to match with 9500m.
                                    % This is then divided by the mean
                                    % value for resolution in our work, 3000m.
                           
            
            % Apply 2-d Parzen Filter
            w=parzenwin(L);
            w=w/sum(w); 
            
            
            % TEST PART
          %  d1 = dx(~isnan(dx)&abs(dx)~=inf&dx~=0);
           % d2 = dy(~isnan(dy)&abs(dy)~=inf&dy~=0);
          %  L = round(9500/((min(d2)+ min(d2))/2));  
            %%%%%
            
            is=10;%(L-1)/2;
            u_filt=u.*nan;
            v_filt=v.*nan;
            u_filt1=u_filt;
            v_filt1=v_filt;
            
            taux_filt=taux.*nan;
            tauy_filt=tauy.*nan;
            taux_filt1=taux_filt;
            tauy_filt1=tauy_filt;
            
            [m,n] = size(u);
            %{
            %w1 = repmat(w,1,n);
%                 
%             for j=is+1:m-is
%                 %for i=1:n;
%                   %  L1 = round(9500./(((dx(j,:))+ (dy(j,:)))./2));
%                  %   for kk = 1 : length(L1);
%                  %   w1(:,kk = parzenwin(L1(kk));
%                 %    end
%                     u_filt1(j,:) = sum(w1.*u(j-is:j+is,:));
%                     v_filt1(j,:) = sum(w1.*v(j-is:j+is,:));
%                     taux_filt1(j,:) = sum(w1.*taux(j-is:j+is,:));
%                     tauy_filt1(j,:) = sum(w1.*tauy(j-is:j+is,:));
%                     
%                     %u_filt1(j,i)=sum(w1.*u(j-is:j+is,i));
%                     %v_filt1(j,i)=sum(w1.*v(j-is:j+is,i));
%                     %taux_filt1(j,i)=sum(w1.*taux(j-is:j+is,i));
%                    % tauy_filt1(j,i)=sum(w1.*tauy(j-is:j+is,i));
%                 %end
%             end
%             
%             w2 = repmat(w',m,1);
%             
%            % for j=1:m
%                 for i=is+1:n-is;
%                     
%                     
%                     u_filt(:,i) = sum(w2.*u_filt1(:,i-is:i+is),2);
%                     v_filt(:,i) = sum(w2.*v_filt1(:,i-is:i+is),2);
%                     taux_filt(:,i) = sum(w2.*taux_filt1(:,i-is:i+is),2);
%                     tauy_filt(:,i) = sum(w2.*tauy_filt1(:,i-is:i+is),2);
%                     
%                    % u_filt(j,i)=sum(w.*u_filt1(j,i-is:i+is)');
%                    % v_filt(j,i)=sum(w.*v_filt1(j,i-is:i+is)');
%                    % taux_filt(j,i)=sum(w.*taux_filt1(j,i-is:i+is)');
%                    % tauy_filt(j,i)=sum(w.*tauy_filt1(j,i-is:i+is)');
%                 end
           % end
            %}
            
           
               
            %for j=is+1:m-is
                for i=1:n;
                    
                    res = dx(:,i);
                    r = mean(res(~isnan(res)&res>0&abs(res)<inf));
                    if isnan(r)
                        r = 2000;
                    end
                    
                    L = round(10000/r);
                    w = parzenwin(L);
                    w = w/sum(w);
                    is = (L-1)/2;
                    
                    w1 = repmat(w,m/length(w),1);
                    % filter by row
                    
                    u_filt1(:,i) = sum(w1.*u(:,i));
                    v_filt1(:,i) = sum(w1.*v(:,i));
                    taux_filt1(:,i) = sum(w1.*taux(:,i));
                    tauy_filt1(:,i) = sum(w1.*taux(:,i));
                    
                    
                   % u_filt1(j,i)=sum(w.*u(j-is:j+is,i));
                    %v_filt1(j,i)=sum(w.*v(j-is:j+is,i));
                    %taux_filt1(j,i)=sum(w.*taux(j-is:j+is,i));
                  %  tauy_filt1(j,i)=sum(w.*tauy(j-is:j+is,i));
                end
           % end
            
       
           
            for j=1:m
                for i=is+1:n-is;
                    
                    res = dx(:,i);
                    r = mean(res(~isnan(res)&res>0&abs(res)<inf));
                    if isnan(r)
                        r = 2000;
                    end
                    L = round(10000/r);
                    w = parzenwin(L);
                    w = w/sum(w);
                    
                    % filter by column
                    u_filt(j,i)=sum(w.*u_filt1(j,i-is:i+is)');
                    v_filt(j,i)=sum(w.*v_filt1(j,i-is:i+is)');
                    taux_filt(j,i)=sum(w.*taux_filt1(j,i-is:i+is)');
                    tauy_filt(j,i)=sum(w.*tauy_filt1(j,i-is:i+is)');
                end
            end
            
            
            
            % convert filtered data to neutral wind
            
            cd = .002;
            rho = 1020;
            u10_filt = sqrt(abs(taux_filt)/rho/cd).*taux_filt./abs(taux_filt);
            v10_filt = sqrt(abs(tauy_filt)/rho/cd).*tauy_filt./abs(tauy_filt);
            
            % add zero mean white gaussian noise
            u_filt_n = u_filt + 1.1.*randn(m,n);
            v_filt_n = v_filt + 1.1.*randn(m,n);
            
            u10_filt_n = u10_filt + .5.*randn(m,n);
            v10_filt_n = v10_filt + .5.*randn(m,n);
            
            % convert wind back to stress
            taux_filt_n = rho*cd*u10_filt_n.*abs(u10_filt_n);
            tauy_filt_n = rho*cd*v10_filt_n.*abs(v10_filt_n);
            
            % compute kinetic energy flux
            
            keU = taux_filt_n.*u_filt_n;
            keV = tauy_filt_n.*v_filt_n;
            
            fileName = [pout 'dopplerSCAT/KE_U' '/KE_U_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
            writebin(fileName,keU);
            
            fileName = [pout 'dopplerSCAT/KE_V' '/KE_V_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
            writebin(fileName,keV);
        end
    else
        warning(['MISSING FILE ' fin])
        break
    end
    
end