% script file to process ship board ADCP data from Georges Bank
% obtained with the transect program and reformated with tranmat.f
%
% a header file "*hd#" is first read
% this file has the following information in each column:
%
%
% 1 -- consecutive ensemble number
% 2 -- decimal julian day
% 3-5 -- year,month,day
% 6-9 -- hour,min,sec,sechundreth
% 10 -- recorded ensemble #
% 11 -- number of ensembles average in section
% 12 -- pitch
% 13 -- roll
% 14 -- heading
% 15 -- temperature
% 16 -- bt vel. u
% 17 -- bt vel. v
% 18 -- bt vel. ve
% 19 -- bt vel. err
% 20 -- depth reading beam 1
% 21 -- depth reading beam 2
% 22 -- depth reading beam 3
% 23 -- depth reading beam 4
% 24 -- Latitude (deg)
% 25 -- Longitude (deg)
% 26 -- ship's u vel (cm/s)
% 27 -- ship's v vel (cm/s)
% 28 -- elapse distance (m) from bt
% 29 -- elapse time (s)
% 30 -- number of valid depth bins in profile
% 
% a sequence of profile files "*pf#" are then read
% these files contain the following information by column:
%
% 1 -- consecutive ensemble #
% 2 -- bin depth (m)
% 3-6 -- Earth velocity E, N, V, Err
% 7-10 -- Intensity beam 1, 2, 3, 4
% 11 -- percent good 
%
%
% if when converting into ascii files with in TRANSECT or PLAYBACK
% the reference was set to NAV then inav=1 else if NONE was use inav=0;
inav=0;

% if ibottom = 1 then use only bottom track velocities to calculate
% averages.
ibottom=1;

iprt=0; vel=1; velmax=1.5; %maximum velocity admitted in output
 cruise='sj9508',
%
% flag = -2^15 in new transect program
% 
  flag = -2^15; pcgood = 95;
klh=0; % index for hourly average values
knav=0; kbt=0;
tree = '/vol/mar/d/georges/sjohnson/';

% constants
% define scales and labels for plots

if vel,
 scale= 5.e-1; stext='50 cm/s'; svec=.50; 
 seli= 5; stexeli='1.0 m/s'; sveli=.5; 
 sscale= 16; sstext='0.20 m/s'; ssvec=.2;
else 
 scale= 5.e-3; stext='50 m^2/s'; svec=50.; 
 seli= 5.e-3; stexeli='100 m^2/s'; sveli=50.; 
 sscale= .2; sstext='20 m^2/s'; ssvec=20.;
end

 nsec = input(' # of cruise sections to process ? ');

for k=1:nsec %begin big k loop for sections

filen = int2str(k);
 filehd = [tree cruise '/' cruise,'hd',filen],
 eval(['load ',filehd])
tic
%
% calculate transports for the whole cruise where bottom track is available.

en=h(:,1); djul=h(:,2); ubt=h(:,16); vbt=h(:,17); vebt=h(:,18); dep = h(:,20:23);
lat = h(:,24); lon = h(:,25); unav = h(:,26); vnav = h(:,27); et = h(:,29);

nbt = find(ubt>flag); lnbt=length(nbt)
nav = find(ubt==flag); lnav = length(nav)

if ibottom
 ubt=ubt(nbt); vbt=vbt(nbt); vebt=vebt(nbt);
 lat = lat(nbt); lon = lon(nbt); unav = unav(nbt); vnav = vnav(nbt);
 en = en(nbt); djul = djul(nbt); dep=dep(nbt,:);
end

 dep = dep'; dep = mean(dep); dep = dep';
 len = length(en);

% calculate indices for bt profiles
dbin = h(:,30); cdbin = cumsum(dbin); cdbin = [0;cdbin];
if ibottom
  dbin = dbin(nbt); cdbin = cdbin(nbt);
end

% PROFILES

  filepf = [tree cruise '/' cruise,'pf',filen];
  eval(['load ',filepf])

enp=p(:,1); uu=p(:,3); vv=p(:,4); vee=p(:,5); err=p(:,6); pcc=p(:,11);

% find profiles correct for ship's motion and calculate vertical averages

ut = zeros(len,1); vt = ut; vep = ut;

   for kp=1:len
  %  nen = find(enp==en(kp));
     nen = cdbin(kp)+1:cdbin(kp)+dbin(kp);
if ~isempty(nen)

% verify that profiles are correctly chosen
     pf = enp(nen) ;
 
     if pf(1)~=en(kp),
     'profile not chosen correctly pause' , kp, pause, end

     if ~isempty(find(diff(pf))), pf, 'mixed profiles pause', pause,end

     u=uu(nen); v=vv(nen); ve=vee(nen); pc=pcc(nen); er=err(nen);
     indpc = find(pc>=pcgood);
    if ~isempty(indpc)
     u=u(indpc); v=v(indpc); ve=ve(indpc); er=er(indpc);
    end
     nu=find(u~=flag); u=u(nu);v=v(nu); ve=ve(nu); er=er(nu);
  % only consider velocity components that exceed the error velocity
  % in magnitude at each depth bin
  % aer = abs(er); u=u(abs(u)>aer); v=v(abs(v)>aer); ve=ve(abs(ve)>aer);
  %   u = mean(u); v = mean(v); ve = mean(ve);
 if isempty(u),
               ut(kp) = nan; vt(kp) = nan;     vep(kp) = nan ;
 else
	if ubt(kp)>flag,
		if inav, u = u - unav(kp); v = v - vnav(kp);end
     		u = u - ubt(kp); v = v - vbt(kp); ve = ve - vebt(kp);
     ut(kp) = mean(u); vt(kp) = mean(v); vep(kp) = mean(ve);
	else
            if ~ibottom,
		if unav(kp)>flag,
	            if ~inav, u = u + unav(kp); v = v + vnav(kp); end
     ut(kp) = mean(u); vt(kp) = mean(v);     vep(kp) = mean(ve) ;
		else
     ut(kp) = nan; vt(kp) = nan;     vep(kp) = nan ;
		end
            else 
               ut(kp) = nan; vt(kp) = nan;     vep(kp) = nan ;
            end
        end
 end
end
   end
%
% eliminate nan values 
	na=find(~isnan(ut));
	ut=ut(na);vt=vt(na);vep=vep(na);
	djul=djul(na); lat=lat(na); lon=lon(na);dep=dep(na);
len=length(na); % redefine len

%convert to m/s and write files

ut=.01*ut; vt=.01*vt; vep=.01*vep;

eval(['save ' tree  cruise '/tran',cruise,num2str(k),' djul lat lon ut vt vep dep'])
toc
clear p h uu vv vee err pcc

%Calculate hourly averages

dh = djul*24; idh = fix(dh); indh=find(diff(idh)); indh=[0;indh;len];
lindh = length(indh);

 for kp=1:lindh-1
    klh = klh + 1; si=0;
       inh = indh(kp)+1:indh(kp+1); linh=length(inh);
 if ~ibottom,
     for kq=1:linh
	if ~isempty(find(inh(kq)==nav)), si=1; end
     end
 end
	if si, knav=knav+1; navk(knav)=klh;
	else, kbt=kbt+1; nbtk(kbt)=klh; end
	uth(klh) = mean(ut(inh));
	vth(klh) = mean(vt(inh));
	veph(klh) = mean(vep(inh));
	th(klh) = mean(dh(inh));
	lonh(klh) = mean(lon(inh));
	lath(klh) = mean(lat(inh));
	deph(klh) = mean(dep(inh));
 end

% plot partial results

if k==1,
 figure (1);clf, gbmap; 
 set(gca,'xlim',[-71 -64.5],'ylim',[39.5 43.2],'position',[.15 .13 .7 .6806]) 
 xlabel('Longitude'), ylabel('Latitude'),grid
 title(['ADCP data from Cruise ',cruise,' to Georges Bank.']),
 plot(lon,lat,'c.'), h=gca; set(h,'box','on'),pause(1)
else
 plot(lon,lat,'c.'),pause(1)
end

end %end of k loop for each section

hold off
if iprt, orient landscape, eval(['print -dpsc ',cruise,'track']), end

% plot results

% plot hourly average transport data vectors 
 
figure(2),clf, gbmap , %georbat;
if vel,
 title(['Observed hourly and vertically average currents from ',cruise,'.'])
else
 title(['Observed hourly average transport vectors from ',cruise,'.']), end
 xlabel('Longitude'), ylabel('Latitude'); h=gca;
 set(h,'xlim',[-71 -64.5],'ylim',[39.5 43.2],'position',[.15 .13 .7 .6806])
 xyaxis=[-71 -64.5 39.5 43.2]; gsca = [11 8.5 .7 .6806];
 pltvecgb(lonh(nbtk),lath(nbtk),uth(nbtk),vth(nbtk),scale,xyaxis,'g',gsca)
if ~ibottom
 pltvecgb(lonh(navk),lath(navk),uth(navk),vth(navk),scale,xyaxis,'r',gsca)
end
 scavecgb(-65.5,40,svec,stext,scale,xyaxis,'g',gsca)
 set(h,'box','on')
 hold off , pause(3)

'number of hourly currents with bottom tracking =', length(nbtk)

nu=find(abs(uth)<velmax);
th=th(nu); lath=lath(nu); lonh=lonh(nu); uth=uth(nu); 
vth=vth(nu); veph=veph(nu); deph=deph(nu); 
nu=find(abs(vth)<velmax);
th=th(nu); lath=lath(nu); lonh=lonh(nu); uth=uth(nu); 
vth=vth(nu); veph=veph(nu); deph=deph(nu);

eval(['save ' tree  cruise '/tranh',cruise,' uth vth deph lath lonh veph th'])

figure(3),clf, gbmap , %georbat;

 title(['Saved hourly and vertically average currents from ',cruise,'.'])
 xlabel('Longitude'), ylabel('Latitude'); h=gca;
 set(h,'xlim',[-71 -64.5],'ylim',[39.5 43.2],'position',[.15 .13 .7 .6806])
 xyaxis=[-71 -64.5 39.5 43.2]; gsca = [11 8.5 .7 .6806];
 pltvecgb(lonh,lath,uth,vth,scale,xyaxis,'g',gsca)
 scavecgb(-65.5,40,svec,stext,scale,xyaxis,'g',gsca)
 set(h,'box','on')
 hold off , pause(3)
if iprt, orient landscape,  eval(['print -dpsc ',cruise,'tranh']), end

'number of hourly currents saved in output file =', length(uth)
