%======================================================================

%                    L O A D C T D . M 

%                    doc: Sat Jun 26 15:56:43 2004

%                    dlm: Fri Mar 21 11:01:39 2014

%                    (c) 2004 M. Visbeck & A. Thurnherr

%                    uE-Info: 52 1 NIL 0 0 72 0 2 8 NIL ofnI

%======================================================================

function [d,p]=loadctd(f,d,p)

% function [d,p]=loadctd(f,d,p)

% This routine works for generic ASCII files, containing CTD time series,

% with fields pressure, temperature, salinity, and time.

%====================

% TWEAKABLES

%====================

% PRE-AVERAGING OF GPS TIME (IN DAYS)

p=setdefv(p,'navtime_av',2/60/24);

% set ctdmaxlag.=100 to the maximum pings that the ADCP data can be shifted to

% best match W calculated from CTD pressure time series (loadctd)

% If you have good times set it to 10... if your time base is questionable

% you can set it 100 or more

p=setdefv(p,'ctdmaxlag',150);

% NUMBER OF DATA CHUNKS TO USE FOR TIME LAGGING

p=setdefv(p,'ctdmaxlagnp',600);

% restrict time range to profile and disregard data close to surface

% p.cut = 0 dont restrict

% p.cut > 0 restrict time to adcp depth below a depth of p.cut

p=setdefv(p,'cut',10);

% INTERPOLATE IRREGULAR CTD TIME SERIES

p=setdefv(p,'interp_ctd_times',1);

% FILE LAYOUT

f = setdefv(f,'ctd_header_lines',0);

f = setdefv(f,'ctd_fields_per_line',4);

f = setdefv(f,'ctd_time_field',1);

f = setdefv(f,'ctd_pressure_field',2);

f = setdefv(f,'ctd_temperature_field',3);

f = setdefv(f,'ctd_salinity_field',4);

f = setdefv(f,'ctd_badvals',-9e99);

% TIME BASE

% 	0 for elapsed time in seconds

% 	1 for year-day (1.0 = Jan 1, 00:00)

% 	2 for Visbeck's Gregorian (see gregoria.m)

f = setdefv(f,'ctd_time_base',0);

%======================================================================

% MODIFICATIONS BY ANT:

%   Jun 26, 2004: - totally re-wrote the file-reading, which reduces

%		    run time by factor 6 or so for deep casts

%   Jun 29, 2004: - estimate time offset at 90% max depth during downcast

%		    instead of at max depth (which is not well defined

%		    if there is a bottle stop)

%   Jun 30, 2004: - added handling for different time bases

%   		  - added optional interpolation code

%   Jul  1, 2004: - removed ipos argument

%   Dec 18, 2004: - BUG: header_lines > 0 did not work

%   Jul 19, 2006: - removed NaN before call to interp1 to avoid Matlab 7.2 warning

%   Jan  5, 2007: - merged from 2 different old versions

%		  - added CTD file layout into p structure

%   Jan 26, 2007: - BUG: file layout default was in p (rather than f) structure

%   May 23, 2008: - added error message on imaginary soundspeed values

%		  - added f.ctd_badvals

%   Jun 26, 2008: - change default of p.interp_ctd_times to 1, because this

%		    improves time matching of irregularly spaced CTD time

%		    series without affecting regularly spaced ones

%   Jun 30, 2008: - BUG: typo in log output

%   Jul  2, 2008: - removed "not recomended" (sic) comment from log file when

%		    profile is created from time series

%   Jul 17, 2008: - moved some code from [loadnav.m] to fix bug associated with

%		    adjusting of start/end positions in case of significant

%		    ADCP vs GPS/CTD clock offset (unclear whether this happened

%		    only when elapsed time was used)

%   Sep 18, 2008: - BUG: moved code had assumed that nav data were loaded

%   Jan  7, 2009: - tightened use of exist()

%   Jun 16, 2009: - BUG: patching short nav time series did not work correctly

%   Mar 21, 2014: - BUG: f.ctd_time_base used p.ctd_time_base set as default

% read SEABIRD ctd timeseries file

disp(['LOADCTD: load CTD time series ',f.ctd])

if ~exist(f.ctd,'file')

 warn=([' LOADCTD can not find ',f.ctd]);

 p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;

 disp(warn)

 return

end

% construct input format

cur_field = 1; input_format = '';

for i=1:f.ctd_fields_per_line

  switch i,

    case f.ctd_pressure_field,

      i_press = cur_field; cur_field = cur_field + 1;

      input_format = [input_format ' %g'];

    case f.ctd_temperature_field

      i_temp = cur_field; cur_field = cur_field + 1;

      input_format = [input_format ' %g'];

    case f.ctd_salinity_field

      i_salin = cur_field; cur_field = cur_field + 1;

      input_format = [input_format ' %g'];

    case f.ctd_time_field

      i_time = cur_field; cur_field = cur_field + 1;

      input_format = [input_format ' %g'];

    otherwise

      input_format = [input_format ' %*g'];

  end

end

if cur_field ~= 5

  error('File format definition error');

end

% open input & skip header

fp=fopen(f.ctd);

header_lines = f.ctd_header_lines;

while header_lines > 0

  header_lines = header_lines - 1;

  fgets(fp);

end

% read time series

[A,nread] = fscanf(fp,input_format,[4,inf]);

fclose(fp);

% remove bad values

ibad = find(A == f.ctd_badvals);

if (length(ibad) > 0)

  A(ibad) = nan;

  disp(sprintf(' Warning: %d bad values (%g) removed from CTD data',length(ibad),f.ctd_badvals));

end

% CTD time

timctd=A(i_time,:)';

switch f.ctd_time_base

  case 0 % elapsed time in seconds

    timctd = timctd/24/3600 + julian(p.time_start);

  case 1 % year-day

    timctd = timctd + julian([p.time_start(1) 1 0 0 0 0]);

end

disp(sprintf(' read %d CTD scans; median delta_t = %.2f seconds',...

		length(timctd),median(diff(timctd))*24*3600));

% pressure, temperature, salinity

data=A([i_press i_temp i_salin],:)';

% interpolate to regular time series

if p.interp_ctd_times

  min_t = min(timctd);

  max_t = max(timctd);

  delta_t = median(diff(timctd));

  data = interp1q(timctd,A([i_press i_temp i_salin],:)',[min_t:delta_t:max_t]');

  timctd = [min_t:delta_t:max_t]';

  disp(sprintf(' interpolated to %d CTD scans; delta_t = %.2f seconds',...

  		length(timctd),median(diff(timctd))*24*3600));

end

% calc LADCP depth

w=meannan(d.rw);

ii=find(~isfinite(w));

w(ii)=0;

dt=diff(d.time_jul)*24*3600;

dt=mean([0,dt;dt,0]);

z=cumsum(w.*dt);

% guess time offset; in Martin's version, this was done based on bottom

% time, which does not work well when there is a bottle stop at the bottom, as

% the LADCP depth can drift slowly. Now, the guess is done based on the time

% when the package reaches 90% of the maximum depth/pressure.

[zmax,izmax]=max(z);

[pmax,ipmax]=max(data(:,1));

i_ctd_near_bottom   = min(find(data(:,1) > 0.9*pmax));

i_ladcp_near_bottom = min(find(z         > 0.9*zmax));

delta_t = (d.time_jul(i_ladcp_near_bottom) - timctd(i_ctd_near_bottom))*24*60;

disp([' 90% LADCP depth  ',int2str(0.9*zmax),' at ',datestrj(d.time_jul(i_ladcp_near_bottom))])

disp([' 90% CTD pressure ',int2str(0.9*pmax),' at ',datestrj(timctd(i_ctd_near_bottom))])

% in case elapsed time is given, estimate absolute time

if f.ctd_time_base == 0

  p.ctdtimoff = delta_t/60/24;

  timctd = timctd + p.ctdtimoff;

  delta_t = 0;

end

% now, handle warnings; NB: cannot happen if elapsed time is given

if abs(delta_t)>5

  disp('WARNING WARNING ');

  warn=(sprintf(' estimated time offset between CTD/LADCP is %.1f minutes',delta_t));

  disp(warn); p.warn(size(p.warn,1)+1,1:length(warn))=warn;

end

if abs(delta_t)>60

  p.ctdtimoff = round(delta_t/60)/24;

  p.ctdtimoff = delta_t/60/24;

  disp(sprintf(' maybe the hour was set wrong: adjusting CTD time by %d days',...

 	         p.ctdtimoff));

  timctd = timctd + p.ctdtimoff;

end

if p.navdata & f.nav_time_base == 0

  disp(sprintf(' adjusting GPS time to CTD time (%+d seconds)',floor(p.ctdtimoff*24*3600)))

  d.navtime_jul = d.navtime_jul + p.ctdtimoff;

end

% plot profile

figure(4);

clf;

orient tall;

subplot(325)

tim0=timctd(1);

plot((timctd-tim0)*24,-data(:,1))

axis tight

title('Cut CTD profile')

xlabel('Time in hours')

ylabel('CTD pressure [dBar]')

hold on

ax=axis;

%   start time

if p.cut~=0

 i10=find(abs(data(1:ipmax,1)-p.cut/2)<p.cut/2);

 if length(i10)>0

  istart=i10(end);

  p=setdefv(p,'ctdtime',1);

  plot(24*(timctd(i10)-tim0),-data(i10,1),'.g')

 else

  istart=1;

  if abs(data(1,1)-p.cut)<10

   p.ctdtime=1;

  else

   disp([' WARNING first CTD depth is: ',int2str(data(1,1)),' [dB]'])

   disp([' WARNING ignore pressure time series '])

   p.ctddepth=0;

   p.ctdtime=0;

   return

  end

 end

else

 istart=1;

 p.ctdtime=0;

end

p.ctd_starttime=timctd(istart);

plot((timctd(istart)-tim0)*24,-data(istart,1),'or')

plot([1 1]*24*(timctd(istart)-tim0),ax(3:4),'r--')

%   end time

if p.cut~=0

 i10=find(abs(data(ipmax:end,1)-p.cut/2)<p.cut/2)+ipmax-1;

 if length(i10)>0

  iend=i10(1);

  plot((timctd(i10)-tim0)*24,-data(i10,1),'.g')

 else

  iend=length(timctd);

  p.ctdtime=0;

 end

else

 iend=length(timctd);

end

p.ctd_endtime=timctd(iend);

plot((timctd(iend)-tim0)*24,-data(iend,1),'or')

plot([1 1]*24*(timctd(iend)-tim0),ax(3:4),'r--')

pause(0.01)

if p.ctdtime==1

 disp([' Changed Start Time : ',datestrj(d.time_jul(1)),...

      '  to ',datestrj(p.ctd_starttime)])

 p.time_start=gregoria(p.ctd_starttime);

 disp([' Changed End   Time : ',datestrj(d.time_jul(end)),...

      '  to ',datestrj(p.ctd_endtime)])

 p.time_end=gregoria(p.ctd_endtime);

end

% ===================================================================

% at this point you NEED to have read the raw CTD data

%

% timectd to be in the same time units as the ADCP data

% data a matrix starting with pressure (not depth) in dbar

%

% if you have synchronized CTD ADCP data jump below

% =================================================================

% fix time for interpolation

if min(timctd)>max(d.time_jul) | max(timctd)<min(d.time_jul)

  disp('CTD timeseries does not overlap WRONG STATION????')

  disp(['  CTD-data : ',datestrj(timctd(1)),'  to ',datestrj(timctd(end))])

  disp([' ADCP-data : ',datestrj(d.time_jul(1)),'  to ',datestrj(d.time_jul(end))])

  disp(' will try to use W-data to get depth')

  p.ctddepth=0;

  return

end

p.ctddepth=1;

% check for latitude

if existf(p,'poss')

  lat = p.poss(1);

elseif p.navdata

  lat = medianan(d.slat);

else

  disp('WARNING WARNING ');

  warn=(' do not have required latitude for accurate pressure->depth conversion; using 0 instead');

  disp(warn); p.warn(size(p.warn,1)+1,1:length(warn))=warn;

  lat = 0;

end

ii=find(diff(timctd)==0);

data(ii,:)=[];

timctd(ii)=[];

[ln,lv]=size(data);

% check for spikes in pressure data

%

dtadcp=medianan(diff(d.time_jul));

dtctd=medianan(diff(timctd),length(timctd)/5);

disp([' median CTD time difference ',num2str(dtctd*24*3600),' s'])

% fix jitter in CTD time base

timctd2=cumsum(round((diff(timctd))/dtctd)*dtctd);

timctd2=[timctd(1); timctd(1)+timctd2];

ibad=1:10;

sbad=0;

nn=0;

while length(ibad)>9 & nn<10

 z=-p2z(data(:,1),lat);

 wctd=[-diff(z)./(diff(timctd2)*24*3600)];

 ibad=find(abs(wctd)>(3*std(wctd)));

 ibad=[ibad; ibad+1];

 data(ibad,:)=[];

 timctd(ibad,:)=[];

 timctd2(ibad,:)=[];

 sbad=sbad+length(ibad);

 nn=nn+1;

end

disp([' removed ',int2str(sbad),' pressure spikes'])

z=-p2z(data(:,1),lat);

% compute wctd similar to what ADCP would do

nshift=max([1,round(dtadcp/dtctd)])+1;

nshift2=fix(nshift/2);

disp([' use ',int2str(nshift2),' time base for W_ctd'])

wctd=z+nan;

i2=(nshift):(length(z)-nshift2);

i1=1:length(i2);

% get W from CTD 

wctd(i1+nshift2)=-(z(i2)-z(i1))./((timctd2(i2)-timctd2(i1))*24*3600);

if sbad>10

  warn=[' removed ',int2str(sbad),' pressure spikes during: ',int2str(nn),' scans'];

  disp(warn)

  p.warn(size(p.warn,1)+1,1:length(warn))=warn;

end

% interpolate on ADCP time series

timctd(1)=0;

timctd(end)=1e32;

data_int=interp1(timctd,data,d.time_jul','nearest');

igood = find(isfinite(wctd));

d.wctd=interp1(timctd(igood),wctd(igood),d.time_jul','nearest')';

d.z=-p2z(data_int(:,1),lat)';

disp([' CTD max depth : ',int2str(-min(d.z))])

if p.ctdmaxlag>0

 % prepare for time lag check

 dt=mean(diff(d.time_jul));

 w=medianan(d.rw,2);

 iw=sum(isfinite(d.rw));

 ii=find(iw<4);

 w(ii)=nan;

 % check for timelag

 dtctd=medianan(diff(timctd));

 % make up array to check for lag

 ctdtw=[timctd, wctd];

 adcptw=[d.time_jul', w'];

 [lag,co]=besttlag(ctdtw,adcptw,p.ctdmaxlag,p.ctdmaxlagnp);

 lagdt=-lag*dtctd;

 disp([' best lag W: ',int2str(lag),' CTD scans ~ ',...

             int2str(lagdt*24*3600),' seconds  corr:',num2str(co)]);

 p=setdefv(p,'ctdmincorr',max(0.9-100/length(wctd),0.6));

 if co<p.ctdmincorr

 % try whole cast

  [lag,i1,i2,co]=bestlag(d.wctd,w,p.ctdmaxlag);

  lagdt=-lag*dt;

  disp([' best lag W: ',int2str(lag),' CTD scans ~ ',...

             int2str(lagdt*24*3600),' seconds  corr:',num2str(co)]);

 end

 % recompute W-ctd

 data_int=interp1(timctd,data,d.time_jul'+lagdt,'linear');

 igood = find(isfinite(wctd));

 d.wctd=interp1(timctd(igood),wctd(igood),d.time_jul'+lagdt,'nearest')';

 d.z=-p2z(data_int(:,1),lat)';

 i1=1:length(w);

% plot two small sections for visual control

 ii=[-12:12]+fix(length(i1)*0.7);

 ii=[ii, [-12:12]+fix(length(i1)*0.3)];

 ii=ii(find(ii>0 & ii<length(i1)));

 subplot(326)

 plot(w(i1(ii)),'-b')

 hold on

 plot(d.wctd(i1(ii)),'-r')

 axis tight

 ax=axis;

 plot([25 25],ax(3:4),'-k')

 text(5,mean(w(i1(ii(25:end)))),'down cast')

 text(30,mean(w(i1(ii(1:25)))),'up cast')

 title(['best lag W: ',int2str(lag),' scans ~ ',...

               int2str(lagdt*24*3600),' sec.  c:',num2str(co)]);

 ylabel('W used for lag correlation')

 xlabel('sample scans (b)ADCP (r)CTD')

 grid

 pause(0.1)

 if p.debug==9, keyboard, elseif p.debug>1, pause, end

 % decide if the correlation is good enought to shift ADCP time

 if abs(lag)<p.ctdmaxlag & co>p.ctdmincorr

  disp(' adjust ADCP time to CTD time and shift depth record ')

  d.time_jul=d.time_jul+lagdt;

  if lagdt*24*3600>10

   disp('WARNING WARNING WARNING')

   warn=[' shifted ADCP timeseries by ',int2str(lagdt*24*3600),' seconds '];

   disp(warn)

   p.warn(size(p.warn,1)+1,1:length(warn))=warn;

  end

 else

  error('Cannot determine time offset between CTD and LADCP time series --- aborting');

  % The following old code has been disabled for IX_10 because it is known not to work correctly

  % for combined NAV/CTD data, because the NAV time series *must* be matched to the LADCP data

  % for the GPS constraint to be applied.

  text(1,0,'LOW CORRELATION not used','color','r')

  disp('WARNING WARNING WARNING')

  warn=(' lag too big or correlation too weak: ignore CTD time series');

  disp(warn)

  p.warn(size(p.warn,1)+1,1:length(warn))=warn;

  p.ctddepth=0;

  if abs(lag)<p.ctdmaxlag 

   disp(' will still use max ctd depth to constrain int W')

   p.zpar(2)=-min(d.z);

  end

  rmfield(d,'z');

  rmfield(d,'wctd');

  return

 end

end

% save other CTD data 

d.ctd_data=data_int(:,2:end)';

% pressure data

d.ctd_p=data_int(:,1)';

% if you have a temperatue time series provide here

d.ctd_temp=data_int(:,2)';

% if you have sound speed time series provide here

% d.ctd_ss=data_int(:,iss-1)';

% have salinity time series provide here

d.ctd_s=data_int(:,3)';

[pmax,imax]=max(d.ctd_p);

pp=[max(d.ctd_p(1),0):pmax]';

ipp=1:imax;

[dum,is]=sort(d.ctd_p(ipp));

ipps=ipp(is);

if existf(d,'ctdprof_p')==0

 % make ctd-profile from data;

 d.ctdprof_p=pp;

 d.ctdprof_lat=lat;

 d.ctdprof_z=interp1q(d.ctd_p(ipps)',-d.z(ipps)',pp);

 d.ctdprof_t=interp1q(d.ctd_p(ipps)',d.ctd_temp(ipps)',pp);

 d.ctdprof_s=interp1q(d.ctd_p(ipps)',d.ctd_s(ipps)',pp);

% get N^2

 if exist('sw_bfrq','file')

  d.ctdprof_N2=sw_bfrq(d.ctdprof_s,d.ctdprof_t,d.ctdprof_p,d.ctdprof_lat);

  d.ctdprof_N2(end+1)=d.ctdprof_N2(end);

  d.ctdprof_ss=sw_svel(d.ctdprof_s,d.ctdprof_t,d.ctdprof_p);

  if ~isreal(d.ctdprof_ss)

    error('sound-speed profile has imaginary values --- check CTD data');

  end

 else

  disp(' download SW routines to get more accurate soundspeed and N^2 ')

 end

 disp(' made CTD profile from time series data. ')

end

if p.ctdtime==1 		% clear all velocity data outside time range

 ii=(d.time_jul>=p.ctd_starttime) & (d.time_jul<=p.ctd_endtime);

 if length(ii)~=sum(ii)

  d=cutstruct(d,ii);

%%% CODE DISABLED BECAUSE IT IS UNCLEAR WHETHER IT IS CONSISTENT WITH

%%% REMAINDER OF CHANGES FOR VERSION IX_6

%%%  if ~existf(d,'slon')

%%%   % don't have time variable navigation data, adjust end position

%%%   % to account for shorter cast time

%%%   slat=p.poss(1)+p.poss(2)/60;

%%%   slon=p.poss(3)+p.poss(4)/60;

%%%   elat=p.pose(1)+p.pose(2)/60;

%%%   elon=p.pose(3)+p.pose(4)/60;

%%%   elon=slon+(elon-slon)*sum(ii)/length(ii);

%%%   elat=slat+(elat-slat)*sum(ii)/length(ii);

%%%   p.pose=[fix(elat), (elat-fix(elat))*60, fix(elon), (elon-fix(elon))*60];

%%%   disp(' adjusted end position to account for short CTD time series')

%%%  end

 end

end