%compute Chelton's Ekman pumping with filter

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%try to replicate Chelton's color map
%cmap=flipud(colormap(hsv));
%cmap=cmap(12:end,:);
%colormap(cmap);
cmap=(colormap(jet));

%blue-red colormap
br_cmap(1:64,1:3)=1;
for i=1:31
  br_cmap(i,:)=[1-(31-i)/31,1-(31-i)/31,1];
end
for i=33:64
  br_cmap(i,:)=[1,1+(33-i)/31,1+(33-i)/31];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

m=1202;
n=1802; 


dx=500;  %from email, grid spacing is uniform 500m.
dx_filt_ss=5000;  %subsample grid spacing is uniform 5000m.
dy=dx;
dy_filt_ss=dx_filt_ss;

rho0=1020;  %can calculate this from model sst and salinity if needed, but
            %this should work for now.
rhoa=1.2;
Cd=1e-3;   %not clear how Chelton computed this...so make constant for now

for cnt = 1:22;
cnt
switch cnt	
	case 1;
		filt_wv = '1km' %filtering wavelength
		L = 3; %Width of Parzen window
	case 2;
		filt_wv = '2.5km'
		L = 5;
	case 3;
		filt_wv = '5km'
		L = 9;
	case 4;
		filt_wv = '10km'
		L = 19;
	case 5;
		filt_wv = '15km'
		L = 27;
	case 6;
		filt_wv = '20km'
		L = 37;
	case 7;
		filt_wv = '25km'
		L = 45;
	case 8;
		filt_wv = '30km'
		L = 55;
	case 9;
		filt_wv = '35km'
		L = 63;
	case 10;
		filt_wv = '40km'
		L = 73;
	case 11;
		filt_wv = '45km'
		L = 81;
	case 12;
		filt_wv = '50km'
		L = 91;
	case 13;
		filt_wv = '55km'
		L = 101;
	case 14;
		filt_wv = '60km'
		L = 109;
	case 15;
		filt_wv = '65km'
		L = 119;
	case 16;
		filt_wv = '70km'
		L = 127;
	case 17;
		filt_wv = '75km'
		L = 137;
	case 18;
		filt_wv = '80km'
		L = 145;
	case 19;
		filt_wv = '85km'
		L = 155;
	case 20;
		filt_wv = '90km'
		L = 163;
	case 21;
		filt_wv = '95km'
		L = 173;
	case 22;
		filt_wv = '100km'
		L = 181;
end

uv_file='/Net/woce/people/jelya/SMCC/ext_smcc_his_uv.dta';
str_file='/Net/woce/people/jelya/SMCC/smcc_frc_oro.dta';
grd_file='/home/morey/Codes/Scatterometer/WaCM/SMCC/smcc_grd_mod.dta';
tim_file='/Net/woce/people/jelya/SMCC/ext_smcc_his_time.txt';

grd_fid=fopen(grd_file,'r');
dum=fread(grd_fid,1,'int32');
msk_rho=fread(grd_fid,[m n],'double');
dum=fread(grd_fid,1,'int32');
dum=fread(grd_fid,1,'int32');
lat_rho=fread(grd_fid,[m n],'double');
dum=fread(grd_fid,1,'int32');
dum=fread(grd_fid,1,'int32');
lon_rho=fread(grd_fid,[m n],'double');
dum=fread(grd_fid,1,'int32');
dum=fread(grd_fid,1,'int32');
ang=fread(grd_fid,[m n],'double'); %(radians ccw from east, I think)
dum=fread(grd_fid,1,'int32');
fclose(grd_fid);

msk_rho(msk_rho<1)=nan;
f=2.*7.2921e-5.*sind(lat_rho); 


A=textread(tim_file);
his_times=mod(squeeze(A(:,2))/86400,360);

disp('reading wind stress file...');
str_fid=fopen(str_file,'r');
str_times=[15:30:345]; %360-day climatology stress?
ustr(length(str_times),m,n)=nan;
vstr(length(str_times),m,n)=nan;
for i=1:length(str_times)
  tmp=fread(grd_fid,[m n],'float');
  ustr(i,:,:)=tmp;
  tmp=fread(grd_fid,[m n],'float');
  vstr(i,:,:)=tmp;
end
fclose(str_fid);


uv_fid=fopen(uv_file,'r');

%i=1;

for i=1:length(his_times);
i
time = strcat('time',int2str(i));
  %interpolate wind field in time...this code needs to be modified if his_time
  %are earlier than 15 days or later than 345 days
  j1=max(find(str_times<=his_times(i)));
  j2=j1+1;
  w1=(str_times(j2)-his_times(i))/(str_times(j2)-str_times(j1));
  w2=1-w1;
  taux=squeeze(w1*ustr(j1,:,:)+w2*ustr(j2,:,:));
  tauy=squeeze(w1*vstr(j1,:,:)+w2*vstr(j2,:,:));

  u=fread(uv_fid,[m n],'double');
  v=fread(uv_fid,[m n],'double');
  u=u.*msk_rho;
  v=v.*msk_rho;

  disp('filtering...');

%Apply 2-d Parzen Filter
  w=parzenwin(L);
  w=w/sum(w);
  is=(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;

  for j=is+1:m-is
    for i=1:n;
      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;
      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

%Subsample grids
	lon_rho_ss = lon_rho(1:10:m,1:10:n);
	lat_rho_ss = lat_rho(1:10:m,1:10:n);
	u_filt_ss = u_filt(1:10:m,1:10:n);
	v_filt_ss = v_filt(1:10:m,1:10:n);
	taux_filt_ss = taux_filt(1:10:m,1:10:n);
	tauy_filt_ss = tauy_filt(1:10:m,1:10:n);

	f_ss=2.*7.2921e-5.*sind(lat_rho_ss); 

%calculate curl -- arrays are x by y, with 1,1 being lower left corner
  dvdx=v*nan;
  dvdx(2:end-1,:)=(v(3:end,:)-v(1:end-2,:))/(2*dx);
  dudy=u*nan;
  dudy(:,2:end-1)=(u(:,3:end)-u(:,1:end-2))/(2*dy);
  zeta=dvdx-dudy;

  dzeta_dx=zeta*nan;
  dzeta_dx(2:end-1,:)=(zeta(3:end,:)-zeta(1:end-2,:))/(2*dx);
  dzeta_dy=zeta*nan;
  dzeta_dy(:,2:end-1)=(zeta(:,3:end)-zeta(:,1:end-2))/(2*dy);

  w_zeta=1./(rho0.*f.*f).*(taux.*dzeta_dy-tauy.*dzeta_dx);

%calculate curl (filtered) -- arrays are x by y, with 1,1 being lower left corner
  dvdx_filt=v_filt*nan;
  dvdx_filt(2:end-1,:)=(v_filt(3:end,:)-v_filt(1:end-2,:))/(2*dx);
  dudy_filt=u_filt*nan;
  dudy_filt(:,2:end-1)=(u_filt(:,3:end)-u_filt(:,1:end-2))/(2*dy);
  zeta_filt=dvdx_filt-dudy_filt;

  dzeta_dx_filt=zeta_filt*nan;
  dzeta_dx_filt(2:end-1,:)=(zeta_filt(3:end,:)-zeta_filt(1:end-2,:))/(2*dx);
  dzeta_dy_filt=zeta_filt*nan;
  dzeta_dy_filt(:,2:end-1)=(zeta_filt(:,3:end)-zeta_filt(:,1:end-2))/(2*dy);

  w_zeta_filt=1./(rho0.*f.*f).*(taux_filt.*dzeta_dy_filt-tauy_filt.*dzeta_dx_filt);

%calculate curl (filtered and subsampled) -- arrays are x by y, with 1,1 being lower left corner
  dvdx_filt_ss=v_filt_ss*nan;
  dvdx_filt_ss(2:end-1,:)=(v_filt_ss(3:end,:)-v_filt_ss(1:end-2,:))/(2*dx_filt_ss);
  dudy_filt_ss=u_filt_ss*nan;
  dudy_filt_ss(:,2:end-1)=(u_filt_ss(:,3:end)-u_filt_ss(:,1:end-2))/(2*dy_filt_ss);
  zeta_filt_ss=dvdx_filt_ss-dudy_filt_ss;

  dzeta_dx_filt_ss=zeta_filt_ss*nan;
  dzeta_dx_filt_ss(2:end-1,:)=(zeta_filt_ss(3:end,:)-zeta_filt_ss(1:end-2,:))/(2*dx_filt_ss);
  dzeta_dy_filt_ss=zeta_filt_ss*nan;
  dzeta_dy_filt_ss(:,2:end-1)=(zeta_filt_ss(:,3:end)-zeta_filt_ss(:,1:end-2))/(2*dy_filt_ss);

  w_zeta_filt_ss=1./(rho0.*f_ss.*f_ss).*(taux_filt_ss.*dzeta_dy_filt_ss-tauy_filt_ss.*dzeta_dx_filt_ss);

	x_ind = zeros(m,n);
	y_ind = zeros(m,n);
	for i=1:m
		x_ind(i,:)=i;
	end
	for j=1:n
		y_ind(:,j)=j;
	end

%calculate Wc term.  
 uu=taux/(rhoa*Cd);
  vv=tauy/(rhoa*Cd);
  strmag=(uu.*uu+vv.*vv).^.25;
  ubg=uu./strmag;
  vbg=vv./strmag;

  ua=(ubg-u);
  va=(vbg-v);
  mag=sqrt(ua.*ua+va.*va);
  ua=ua.*mag;
  va=va.*mag;

  dvadx=va*nan;
  dvadx(2:end-1,:)=(va(3:end,:)-va(1:end-2,:))/(2*dx);
  duady=ua*nan;
  duady(:,2:end-1)=(ua(:,3:end)-ua(:,1:end-2))/(2*dy);
  crl=dvadx-duady;

  w_c=(rhoa.*Cd./(rho0.*f)) .*crl;

%calculate Wc term (filtered). 
  uu_filt=taux_filt/(rhoa*Cd);
  vv_filt=tauy_filt/(rhoa*Cd);
  strmag_filt=(uu_filt.*uu_filt+vv_filt.*vv_filt).^.25;
  ubg_filt=uu_filt./strmag_filt;
  vbg_filt=vv_filt./strmag_filt;

  ua_filt=(ubg_filt-u_filt);
  va_filt=(vbg_filt-v_filt);
  mag_filt=sqrt(ua_filt.*ua_filt+va_filt.*va_filt);
  ua_filt=ua_filt.*mag_filt;
  va_filt=va_filt.*mag_filt;

  dvadx_filt=va_filt*nan;
  dvadx_filt(2:end-1,:)=(va_filt(3:end,:)-va_filt(1:end-2,:))/(2*dx);
  duady_filt=ua_filt*nan;
  duady_filt(:,2:end-1)=(ua_filt(:,3:end)-ua_filt(:,1:end-2))/(2*dy);
  crl_filt=dvadx_filt-duady_filt;

  w_c_filt=(rhoa.*Cd./(rho0.*f)) .*crl_filt;

%calculate Wc term (filtered and subsampled).  
  uu_filt_ss=taux_filt_ss/(rhoa*Cd);
  vv_filt_ss=tauy_filt_ss/(rhoa*Cd);
  strmag_filt_ss=(uu_filt_ss.*uu_filt_ss+vv_filt_ss.*vv_filt_ss).^.25;
  ubg_filt_ss=uu_filt_ss./strmag_filt_ss;
  vbg_filt_ss=vv_filt_ss./strmag_filt_ss;

  ua_filt_ss=(ubg_filt_ss-u_filt_ss);
  va_filt_ss=(vbg_filt_ss-v_filt_ss);
  mag_filt_ss=sqrt(ua_filt_ss.*ua_filt_ss+va_filt_ss.*va_filt_ss);
  ua_filt_ss=ua_filt_ss.*mag_filt_ss;
  va_filt_ss=va_filt_ss.*mag_filt_ss;

  dvadx_filt_ss=va_filt_ss*nan;
  dvadx_filt_ss(2:end-1,:)=(va_filt_ss(3:end,:)-va_filt_ss(1:end-2,:))/(2*dx_filt_ss);
  duady_filt_ss=ua_filt_ss*nan;
  duady_filt_ss(:,2:end-1)=(ua_filt_ss(:,3:end)-ua_filt_ss(:,1:end-2))/(2*dy_filt_ss);
  crl_filt_ss=dvadx_filt_ss-duady_filt_ss;

  w_c_filt_ss=(rhoa.*Cd./(rho0.*f_ss)) .*crl_filt_ss;

%calculate Kinetic Energy Flux (u dot tau) (unfiltered)
	keflux = u.*taux+v.*tauy;

%calculate Kinetic Energy Flux (u dot tau) (filtered)
	keflux_filt = u_filt.*taux_filt+v_filt.*tauy_filt;

end
fclose(uv_fid);
end