% Purpose
% Read the QuikSCAT JPL vector wind product
% Convert ascending and descending pass winds to stresses
% average to daily stresses
% Read in AVISO currents and interpolate to daily values
% calculate wind power input

% History:
% Built from qscat_read_u10.m
% Oct. 19/2007
% switched to reading the matlab files. Only difference is that wspd
% is not saved to matlab, so must be calculated. 
% hopefully this accounts for the difference I find:
% std(wpi_wk(:)-wpi_wk243(:)) = 3e-5
% compare this too: std(wpi_wk(:)) = 2e-2
%			std(wpi_LP - wpi_Liu) = 1e-3

% Input data

% Aviso currents:

dir = '/disk/staff/rscott/Aviso/dt_upd_madt_UV';
ufilelist = 'ncfile_names_19991229_20060104.txt';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Setup parameters:
%%%%%%%%%%%%%%%%%%%

% Convert winds to stresses using 

% flag_algor = 1 % Liu algorithm
% flag_algor = 2 % Large and Pond

flag_algor = 1
flag_straub = 1

cd /disk/staff/rscott/QuikSCAT/Winds
filelist = 'file_u10_names_2000_2005.txt';
stress_files = textread(filelist,'%s');
disp('reading this many stress_files ...')

nstressfiles = size(stress_files,1)

lon = (360/1440)*([0:1439]+0.5);
lat = (180/720) *([0:719]+0.5)-90;

%------------------------------
% Bulk formulae coefficients:
%------------------------------
% for Liu algorithm only:
alpha_c = 0.011;
beta = 0.11;
max_it = 15;

% for both:
rho = 1.22;

% for Large and Pond:
a1 = 0.0027;
a2 = 0.000142;
a3 = 0.0000764;
%------------------------------
skipwk = 0;   % default=0, set to positive integer to skip first skipwk weeks.
                % must also adjust current data reading too.
part_wk = 4;
% i1 = 4; start on the 4th day of the first week of Dec.29, 1999, which is
% is Jan. 1, 2000. Note that i1 is set to i1=1 at the end of first week. 

dy = 0;	% start at dy = 0, so it's incremented to 1, reading Jan.2/2000 first

if skipwk > 0
   dy = (7-part_wk+1) + 7*(skipwk-1);	% start at dy = 1, so it's incremented to 2, reading Jan.2/2000 as first
elseif skipwk==0
   dy = 0;
else
   error('set skipwk >= 0')
end
n = 1 + skipwk;	% for uv files, so it's incremented to 2, reading Jan.5/2000 first

% initialize arrays:

unity = ones(720,1440);

taux_bar  = zeros(720,1440);
tauy_bar  = zeros(720,1440);
count_tot = zeros(720,1440);

% wind power input arrays:

utau_aviso  = zeros(915,1080);
ngood_aviso = zeros(915,1080);
unity_aviso =  ones(915,1080);

% all uv file names:
filenamelen = 37;
files = textread(strcat(dir,'/',ufilelist),'%s');
nfiles = size(files,1);

% get first file, set as the "next week", so that it's shifted to 
% "last week" on first step
disp('first weeks current file')
infile = char(files(n))
outfile = strcat(infile(1:filenamelen),'.mat');
%getnetcdffile(strcat(dir,'/',infile),strcat(dir,'/',outfile))
load(strcat(dir,'/',outfile))
% convert cm/s --> m/s: orient so first index corresponds to latitude:
u_next = Grid_0001'/100;          
u_next(u_next>1e17) = NaN;
v_next = Grid_0002'/100;          
v_next(v_next>1e17) = NaN;

% loop over all weeks

for wk = skipwk+1:314

%---------------------------
% step 1:	get currents 
%---------------------------
% get last weeks current file
   u_last = u_next;
   v_last = v_next;
% get next weeks current file
   n = n + 1;
   disp('next weeks current file')
   infile = char(files(n))
   outfile = strcat(infile(1:filenamelen),'.mat');
%   getnetcdffile(strcat(dir,'/',infile),strcat(dir,'/',outfile))
   load(strcat(dir,'/',outfile))
% convert cm/s --> m/s: orient so first index corresponds to latitude:
   u_next = Grid_0001'/100;          
   u_next(u_next>1e17) = NaN;
   v_next = Grid_0002'/100;          
   v_next(v_next>1e17) = NaN;
%---------------------------
% hold weekly mean WPI
%---------------------------
   wpi_wk   = zeros(915,1080);
   count_wk = zeros(915,1080);
%---------------------------
%---------------------------
% hold weekly mean tau
%---------------------------
   taux_wk   = zeros(915,1080);
   tauy_wk   = zeros(915,1080);
   n_tau_wk  = zeros(915,1080);
%---------------------------


% first and last week are partial

   switch wk
      case {1,314}
         i1 = part_wk;
      otherwise
         i1 = 1;
   end
   for i = i1:7

% interpolate current to this day
      u_today = (7+1-i)/7 * u_last + (i-1)/7 * u_next;
      v_today = (7+1-i)/7 * v_last + (i-1)/7 * v_next;
      dy = dy+1;

% get daily stress on aviso grid
      taux_dy  = zeros(915,1080);
      tauy_dy  = zeros(915,1080);
      count_dy = zeros(915,1080);

      infile = char(stress_files(dy))
      load(strcat(infile(1:22),'.mat'))
%%%%%%%%%%%%%%reading files %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%variables: u10a, v10a, counta
%           u10d, v10d, countd
% need:
% variables: wspda
%           wspdd

% asc
      u10  = NaN*zeros(720,1440);
      v10  = NaN*zeros(720,1440);
      u10(gooda) = u10a;
      v10(gooda) = v10a;

% hereafter, on aviso grid

      u10a = regrid_map(lon,lat,u10,NbLongitudes,NbLatitudes);
      v10a = regrid_map(lon,lat,v10,NbLongitudes,NbLatitudes);

      if flag_straub == 1
         u10a = u10a+u_today;
	 v10a = v10a+v_today;
      end 
      gooda = find(~isnan(u10a+v10a));
      u10a = u10a(gooda);
      v10a = v10a(gooda);
      wspda = sqrt(u10a.^2+v10a.^2);

% stress calcul
% 0<= wspd <= 0.1 acconts for about 2/1000 of points!
%      gooda = find(wspda>0.1 & wspda<50);

      if flag_algor == 1
        [ustar,flag,smalla,diffa]=charnock_2d(wspda,alpha_c,beta,max_it);
        if flag == max_it
           disp(strcat('charnock did not converge iterations',int2str(max_it)))
	   disp(strcat('for ascending pass of ',infile))
           gooda = gooda(abs(diffa) < 1e-6);
	   ustar = ustar(abs(diffa) < 1e-6);
        end
      elseif flag_algor == 2
         ustar = sqrt(a1*wspda+a2*wspda.^2+a3*wspda.^3);
      end

      tauxa = rho * ustar.^2 .* u10a ./wspda;
      tauya = rho * ustar.^2 .* v10a ./wspda;

%      goodd = find(wspdd>0.1 & wspdd<50);
% des
      u10  = NaN*zeros(720,1440);
      v10  = NaN*zeros(720,1440);
      u10(goodd) = u10d;
      v10(goodd) = v10d;

% hereafter on aviso grid
     u10d = regrid_map(lon,lat,u10,NbLongitudes,NbLatitudes);
     v10d = regrid_map(lon,lat,v10,NbLongitudes,NbLatitudes);

      if flag_straub == 1
         u10d = u10d+u_today;
	 v10d = v10d+v_today;
      end 

      goodd = find(~isnan(u10d+v10d));
      u10d = u10d(goodd);
      v10d = v10d(goodd);
      wspdd = sqrt(u10d.^2+v10d.^2);

      switch flag_algor
         case 1 
	   [ustar,flag,smalld,diffd]=charnock_2d(wspdd,alpha_c,beta,max_it);
           if flag == max_it
             disp(strcat('charnock did not converge iterations',int2str(max_it)))
	     disp(strcat('for ascending pass of ',infile))
             goodd = goodd(abs(diffd) < 1e-6);
             ustar = ustar(abs(diffd) < 1e-6);
           end
         case 2
           ustar = sqrt(a1*wspdd+a2*wspdd.^2+a3*wspdd.^3);
      end

      tauxd = rho * ustar.^2 .* u10d ./wspdd;
      tauyd = rho * ustar.^2 .* v10d ./wspdd;

% accumulate good pass data

      taux_dy(gooda) = taux_dy(gooda) + tauxa; 
      taux_dy(goodd) = taux_dy(goodd) + tauxd;
      tauy_dy(gooda) = tauy_dy(gooda) + tauya; 
      tauy_dy(goodd) = tauy_dy(goodd) + tauyd;
      count_dy(gooda) = count_dy(gooda) + unity_aviso(gooda);
      count_dy(goodd) = count_dy(goodd) + unity_aviso(goodd);

% daily average: average ascending and descending
      good = find(count_dy>0);
      taux_dy(good) = taux_dy(good) ./ count_dy(good);
      tauy_dy(good) = tauy_dy(good) ./ count_dy(good);   
      taux_dy(count_dy<1) = NaN;
      tauy_dy(count_dy<1) = NaN;
% accumulate good daily data
      taux_wk(good) = taux_wk(good) + taux_dy(good);
      tauy_wk(good) = tauy_wk(good) + tauy_dy(good);
      n_tau_wk(good) = n_tau_wk(good) + unity_aviso(good);
% interpolate to AVISO grid
%      taux_aviso = regrid_map(lon,lat,taux_dy,NbLongitudes,NbLatitudes);
%      tauy_aviso = regrid_map(lon,lat,tauy_dy,NbLongitudes,NbLatitudes);

% accumulate WPI
      wpi = u_today.*taux_dy + v_today.*tauy_dy;
      g = find(~isnan(wpi));

% for long term average
      ngood_aviso(g) = ngood_aviso(g) + unity_aviso(g);
      utau_aviso(g)  = utau_aviso(g)  + wpi(g);

% for weekly average
      wpi_wk(g)   = wpi_wk(g) + wpi(g);
      count_wk(g) = count_wk(g) + unity_aviso(g);
   end
% wpi wkly ave
   g = find(count_wk>0);
   wpi_wk(g) = wpi_wk(g)./count_wk(g);
   wpi_wk(count_wk<1) = NaN;

   g = find(n_tau_wk>0);
   taux_wk(g) = taux_wk(g) ./ n_tau_wk(g);
   tauy_wk(g) = tauy_wk(g) ./ n_tau_wk(g);
   b = find(n_tau_wk<1);
   taux_wk(b) = NaN;
   tauy_wk(b) = NaN;
   
   switch flag_algor
      case 1
         save(strcat('wpi_Liu_wk_',int2str(wk),'straub',int2str(flag_straub),...
	             '.mat'),'wpi_wk','count_wk','taux_wk','tauy_wk','n_tau_wk')
      case 2
         save(strcat( 'wpi_LP_wk_',int2str(wk),'straub',int2str(flag_straub),...
	             '.mat'),'wpi_wk','count_wk','taux_wk','tauy_wk','n_tau_wk')
   end
   disp(strcat('done week: ',int2str(wk)))
end

% 6-yr average WPI
g = find(ngood_aviso>0);
utau_aviso(g) = utau_aviso(g) ./ ngood_aviso(g);
b = find(ngood_aviso<1);
utau_aviso(b) = NaN;

switch flag_algor
   case 1
      outfile_fin = 'Pin_u_tau_20000101_20051231_QS_Liu_aviso_upd_daily.mat';
   case 2 
      outfile_fin = 'Pin_u_tau_20000101_20051231_QS_LP_aviso_upd_daily.mat';
end

save(outfile_fin,'utau_aviso','ngood_aviso','NbLatitudes','NbLongitudes')