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  'dopplerSCAT2/KE_U'])
eval(['mkdir ' pout  'dopplerSCAT2/KE_V'])
%eval(['cd ' pout 'dopplerSCAT'])
[dx,dy] = quikreadpcolor_dxdy_llc('~dmenemen/llc_2160/grid/DXC.data','~dmenemen/llc_2160/grid/DYC.data',nx);

for  ts=92160:(80*24):92160+10*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])
            
        else
            
            % do computations
            % Appy 10km parzen filter
            
            minRes = max(max(dy(~isnan(dy)&dy>0)));
 
            filt_wv = '10km';
                            
            
            % Apply 2-d Parzen Filter
     
            
            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);

            for ii = 1 : n
                
                    res = dx(:,ii);
                    r(ii) = mean(res(~isnan(res)&res>0&abs(res)<inf));
                    if isnan(r(ii))
                        r(ii) = 2000;
                    end
                    
                    L(ii) = 2*round(10000/r(ii)/2)+1;
            end
            
            disp('line 78')
            for ii = (L(1)-1)/2+1 : n-(L(end)+1)/2-1
                
                w = fspecial('gaussian',L(ii),.55);
                is = (L(ii)-1)/2;
                
                for jj = (L(ii)-1)/2+1 : m-(L(ii)+1)/2-1
                    
                    u_filt(jj,ii) = sum(sum(w.*u(jj-is:jj+is,ii-is:ii+is)));
                    v_filt(jj,ii) = sum(sum(w.*v(jj-is:jj+is,ii-is:ii+is)));
                    taux_filt(jj,ii) = sum(sum(w.*taux(jj-is:jj+is,ii-is:ii+is)));
                    tauy_filt(jj,ii) = sum(sum(w.*tauy(jj-is:jj+is,ii-is:ii+is)));
                    
                end
            end
            
            disp('line 92')
            % convert filtered data to neutral wind
            
            cd = .001;
            rho = 1.292;
            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
            for ii = 1:n
                dd = dy(1,ii);
                LL = ceil(10/dd);
                jj = LL;
                while jj<=m-LL
                    noiseu(ii:ii+LL,jj:jj+LL) =  .5.*randn(m,n);
                    noisev(ii:ii+LL,jj:jj+LL) =  .5.*randn(m,n);
                    noiseu10(ii:ii+LL,jj:jj+LL) =  1.1.*randn(m,n);
                    noisev10(ii:ii+LL,jj:jj+LL) =  1.1.*randn(m,n);
                    jj = jj + LL;
                end
            end
                    
            u_filt_n = u_filt + .5.*randn(m,n);
            v_filt_n = v_filt + .5.*randn(m,n);
            
            u10_filt_n = u10_filt + 1.1.*randn(m,n);
            v10_filt_n = v10_filt + 1.1.*randn(m,n);
            %{
              curr_tot = sqrt(u_filt.^2 + v_filt.^2);
            wind_tot = sqrt(u10_filt.^2 + v10_filt.^2);
            
            curr_ang = atan2(v_filt,u_filt);
            wind_ang = atan2(v10_filt,u10_filt);
            
            curr_tot_n = curr_tot + .5.*randn(m,n);
            wind_tot_n = wind_tot + 1.1.*randn(m,n);
            
            
            u_filt_n = curr_tot_n.*cos(curr_ang);
            v_filt_n = curr_tot_n.*sin(curr_ang);
            
            u10_filt_n = wind_tot_n.*cos(wind_ang);
            v10_filt_n = wind_tot_n.*sin(wind_ang);
            %}
            % 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 'dopplerSCAT3/KE_U' '/KE_U_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
            writebin(fileName,keU);
            
            fileName = [pout 'dopplerSCAT3/KE_V' '/KE_V_' int2str((nx)) 'x' int2str((nx*13)) '.' dy];
            writebin(fileName,keV);
        end
    else
        warning(['MISSING FILE ' fin])
        break
    end
    
end