% ti_2006_ctd_horsec.m
% Tioga-171, Tioga-177, ... CTD data: Plot horizontal maps 
% Olga Kosnyreva, AOPE/WHOI, 5/5/2006
%
% I N P U T:
% (1) CTD data
% Filenames: ti171***.cdn and ti177***.cdn, where *** is station #
% Variables:
% 1 - timeS: time [elapsed seconds]
% 2 - pr: pressure [db]
% 3 - t090: temperature, pri, ITS-90 [deg C]
% 4 - c0S/m: conductivity, primary [S/m]
% 5 -  fluorescence
% 6 - irradiance (PAR), 4550, 13 March 1998
% 7 - transmissometer, cst-758dr, 14 July 2004
% 8 -  altimeter
% 9 -  seapoint turbidity meter, 1660
% 10 - oxygen, SBE, primary, 430113, 11/18/2004
% 11 - surface irradiance (SPAR), degrees = 0.0
% 12 - salinity, PSS-78 [PSU]
% 13 - sigmaT, (kg/m3)
% 14 - flag:  0.000e+00
% (2) CTD Stations list data.
% Filenames: ti171_ctd_list.mat (or .txt), ti177_ctd_list.mat (or .txt), ...

clear;
flops(0);
% Save png files.
%saveOption='yes';
saveOption='no';
% ---- MENU: Tioga cruise
crNum=f_menu_of_numbers_title('Tioga cruise',[171 177]);
if isempty(crNum), return; end
if crNum==171
   path_name='~/tioga_2006/ti171/CTD/process/';
   ctdList=1:29;
   date_str='April (11 - 12), 2006'; 
   cruise_title='R/V Tioga TI171';
   cruise_id='ti171';
elseif crNum==177
   path_name='~/tioga_2006/ti177/ctd/process/';
   ctdList=1:29;
   date_str='April (25 - 26), 2006';
   cruise_title='R/V Tioga TI177';
   cruise_id='ti177';
end
%********************************************************
%**************         C O N S T A N T S        ********
%********************************************************
% MENU: Depth.
mapDep=f_menu_of_numbers_title('Depth (m)',[3 10 20 50]);
if isempty(mapDep), return; end
%
% Draw coast line.
drawCoast='yes';
%drawCoast='no';
if strcmp(drawCoast,'yes')
   %  load 'structure fields' for MARMAP1 mesh:
   %  mp - cartesian, mm - geographycal coordinate 
   %marm_ll_g2s_moody
   load land_data;
   mm=land_data;
   mm.x=-mm.x;
   mm.lnd(:,1)=-mm.lnd(:,1);
end
% Land, isobaths - color.
bnd_color=[0.71 0.66 0.43];
isobaths_color=[0.38 0.41 1];
color_lbl=[0 0 0.43];
% Option: colored or black and white map.
%optColor='inBlackWhite';
optColor='inColor';
if strcmp(optColor,'inBlackWhite')
   ocean_color=[1 1 1];
else
   ocean_color=[0.8 0.96 1];    % Ocean
end
%
% Print file type.
printFile='ps'; % and PNG
%printFile='jpg';
if strcmp(printFile,'jpg')
   fnt_lbl_isobath=8;
   fnt_ax=7;
   fnt_lbl=7;
   fnt_title=10;
   fnt_title_sub=fnt_title;
   fnt_text=8;
   marker_size=2;
   text_title_color=[0 0 0.7];
   text_min_max_color=[0.4 0.4 0.26];
   text_cb_color=text_min_max_color;
   marker_edge_color=[0.5 0.5 0.5];
   marker_face_color=[1 1 0.9];
   line_width=0.5;
   land_color=[0.8 0.8 0.66];
else
   fnt_lbl_isobath=8;
   fnt_ax=7;
   fnt_lbl=9;
   fnt_title=10;
   fnt_title_sub=fnt_title;
   fnt_text=6;
   marker_size=2;
   text_title_color=[0 0 0];
   text_min_max_color=[0 0 0];
   text_cb_color=text_min_max_color;
   marker_edge_color=[0.5 0.5 0.5];
   marker_face_color=[1 1 0.9];
   line_width=0.7;
   land_color=[0.855 0.783 0.663];
end
% Lower depth in depth range.
lowerDep=160;
% Load ColorMap.
load palet_12_section_v2;
cmap=palet_12_section_v2;
nband=size(cmap,1);
% CTD list file: lat, lon columns#.
colLat=8;
colLon=9;
% CTD datafile:
headerLastLineNum=95;
n_par = 14;
% cdn CTD datafile columns#.
colPres=2;
colTemp=3;
colCnd=4;
colFluor=5; % (v1) Fluor real units
colPar=6;   % (v2) Par  real units
colTrans=7; % (v3) Trans real units
colAltim=8; % (v4) Altim real units
colTurb=9;  % (v5) Turb real units
colOxy=10;  % (v6) Oxy ml/L real units
colSpar=11; % (v9) Spar real units
colSal=12;
colSig=13;
% Flag - 'bad data'.
bad_flag=-9.990e-29;
% Initiate the variables.
data=ones(300,12);
da=ones(3600,1);
DATAC=ones(1000,5);
% Mapping C O N S T A N T S:
   makeContourLabel='no';
   makeShadingFlat='yes';
   makeContours='no';
% C O N S T A N T S
% dependig on Mapping area:
%mapArea='full';
%mapArea='area1';
%mapArea='area2';
%mapArea='forJPG';
mapArea='spec1_Mass';
if strcmp(mapArea,'full')
   latmin = 41;
   latmax = 45.5;
   lonmax = 71;
   lonmin = 65;
elseif strcmp(mapArea,'area1')
   latmin = 43.1;
   latmax = 44.1;
   lonmin = 68.3;
   lonmax = 70;
elseif strcmp(mapArea,'area2')
   latmin = 43.1;
   latmax = 44.2;
   lonmin = 67.9;
   lonmax = 70.25;
elseif strcmp(mapArea,'forJPG')
   latmin = 42;
   latmax = 45.2;
   lonmax = 71;
   lonmin = 65.8-0.2;
elseif strcmp(mapArea,'spec1_Mass')
   latmin = 41.45;
   latmax = 42.8;
   lonmax = 71+0.1;
   lonmin = 69.7;
else
   latmin = 43.1;
   latmax = 44.2;
   lonmin = 67.9;
   lonmax = 70;
end
% Meshgrid limits.
meshLonMin=lonmin;
meshLonMax=lonmax;
meshLatMin=latmin;
meshLatMax=latmax;
% Meshgrid interval for mapping [deg].
meshInt=0.025;
% Calculate relation relyx = axisy/axisx (for PlotBoxAspectRatio) 
weg = 111*cos(pi*(latmin+(latmax-latmin)/2)/180);% Weight of one degree longitude
xc = (lonmax-lonmin)*weg;                        % Longitude, km
yc = (latmax-latmin)*111;                        % Latitude, km
relyx = yc/xc;
% For CTD data interpolation onto depth vector: depth limits.
z_min_int=1;
z_max_int= lowerDep;
z_int_int=1;
%**************************************************************************
%***************************  LOAD CTD LIST DATAFILE   ********************
%**************************************************************************
p=['load ti' num2str(crNum) '_ctd_list;'];
eval(p);
p=['c_list=ti' num2str(crNum) '_ctd_list;'];
eval(p);
%**************************************************************************
%*************************     LOOP# 2:    QUANTITY#     ******************
%**************************************************************************
% QUANTITIES LIST.
quant_list =[1 2 3]; % Temp, Sal, Sigmat 
countParNum=0;
for parNum=quant_list   % Always any 3 quantities
   countParNum=countParNum+1;
   % CONSTANTS depending on parameter name:
   switch parNum
   case 1
         colNum = colTemp;
         jN='temp';
         j_title='Temperature';
         j_label='Temperature (\circ C)';
         if crNum==171%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& TEMP &&&
            if mapDep==3
               vv=[4.68:0.14:6.36];
            elseif mapDep==10
               vv=[4.84:0.055:5.5];
            elseif mapDep==20
               vv=[4.33:0.075:5.23];
            elseif mapDep==50
               vv=[4.05:0.075:4.95];
            else
                return;
            end
        elseif crNum==177
            if mapDep==3
               vv=[6.14:0.14:7.82];
            elseif mapDep==10
               vv=[6.09:0.13:7.65];
            elseif mapDep==20
               vv=[5.58:0.134:7.188];%
            elseif mapDep==50
               vv=[4.25:0.115:5.63];
            else
               return;
            end
         else
            return;
         end %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& TEMP  E N D  &&&
    case 2
         colNum = colSal;
         N='sal';
         j_title='Salinity';
         j_label='Salinity (psu)';
         if crNum==171%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& SAL &&&
            if mapDep==3
               vv=[30.65:0.17:32.69];
            elseif mapDep==10
               vv=[31.80:0.075:32.7];
            elseif mapDep==20
               vv=[32.03:0.055:32.69];
            elseif mapDep==50
               vv=[32.38:0.039:32.848];
           else
                return;
            end
        elseif crNum==177
             if mapDep==3
               vv=[30.39:0.18:32.55];
            elseif mapDep==10
               vv=[31.40:0.095:32.54];
            elseif mapDep==20
               vv=[31.41:0.095:32.55];
            elseif mapDep==50
               vv=[32.07:0.052:32.694];
           else
                return;
            end
    end %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& SAL  E N D  &&&
   case 3
         colNum = colSig;
         jN='sigmat';
         j_title='Sigmat';
         j_label='Sigmat (kg/m^3)';
         if crNum==171%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& SIGT &&&
            if mapDep==3
                vv=[24.13:0.136:25.762];
            elseif mapDep==10
               vv=[25.12:0.058:25.816];
            elseif mapDep==20
               vv=[25.33:0.04:25.81];
            elseif mapDep==50
               vv=[25.67:0.025:25.97];
           else
                return;
            end
        elseif crNum==177
            if mapDep==3
                vv=[23.70:0.155:25.56];
            elseif mapDep==10
               vv=[24.54:0.085:25.56];
            elseif mapDep==20
               vv=[24.57:0.085:25.59];
            elseif mapDep==50
               vv=[25.28:0.048:25.856];
           else
                return;
            end
        end%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& SIGT E N D  &&&
   case 4
         colNum = colFluor;
         jN='Fluorescence';
         j_title='Fluorescence';
         j_label='Fluorescence  (mkg/L)';  
         vv=0:0.05:0.6; %
   case 5
         colNum = colTrans;
         jN='Transmissometer';
         j_title='Transmissometer';
         j_label='Transmissometer';
         vv=4:0.05:4.6;
   case 6
         colNum = colTurb;
         jN='Seapoint turbidity';
         j_title='Seapoint turbidity meter';
         j_label='Seapoint turbidity meter (FTU)';
         vv=0.04:0.01:0.16;
   case 7
         colNum = colPar;
         jN='Irradiance'; 
         j_title='Irradiance, PAR';
         j_label='Irradiance, PAR';
         vv=0.1:0.25:3.1; 
   case 8
         colNum = colAltim;
         jN='Altimeter';
         j_title='Altimeter';
         j_label='Altimeter  (m)';
   case 9
         colNum = colOxy;
         jN='Oxygen';
         j_title='Oxygen';
         j_label='Oxygen (ml/L)';
         %vv=0.1:0.25:3.1; % 12 colors
   case 10
         colNum = colSpar;
         jN='';
         j_title='Surface Irradiance, SPAR';
         j_label='Surface Irradiance, SPAR';
   case 11
         colNum = colOxy; 
         j_label='Oxygen (\muMol/kg)';
         j_title='Oxygen (\muMol/kg)';
         jN='Oxygen (\muMol/kg)';
   case 12
         colNum = colOxy; 
         j_label='Oxygen saturation (%)';
         j_title='Oxygen saturation (%)';
         jN='Oxygen saturation (%)';
   end
   vv_label=vv;
%*************************
% Set figure properties;
%*************************
   if countParNum==1 
      figure;
      set(gcf,'color',[1 1 1],'position',[526   229   235   625],...
        'paperPositionMode','auto');
      colormap(cmap);
%*************************
% Make a matrix with subplot positions
%*************************
      POS=[];
      plot_rows=3;
      plot_cols=1;
      for i=1:(plot_rows*plot_cols) 
         subplot(plot_rows,plot_cols,i);
         POS=[POS;get(gca,'Position')];
      end 
      clf;
      dx_ax=0;
      dy_ay=-0.1;
      ddx=0;
      ddy=-0.07;
      kx=1.13;
      ky=1.13;
   end
%*****************************************************************
%*************************    For each   STATION#    *************
%*****************************************************************
   % For each station from survey:
   % get one station data, calculate depth, interpolate data onto depth vector.
   DATAC=[];
   countNotUsed=0;
   count_bad=0;
   for stNum=ctdList
      % Write input file names to string.
      if stNum < 10
         j = ['ti' num2str(crNum) '00' num2str(stNum) '.cdn'];
      elseif (stNum<100) & (stNum>9) 
         j = ['ti' num2str(crNum) '0' num2str(stNum) '.cdn'];
      else 
         j = ['ti' num2str(crNum) num2str(stNum) '.cdn'];
      end
      % Open file for reading.
      pfid = fopen([path_name j]);
      % Read first lines (headers).
      for lineNum=1:1:headerLastLineNum
         headerLine=fgets(pfid);
      end
      % Read formated data from file.
      da=fscanf(pfid,'%f',[1 inf])';
      % Reshape data matrix.
      data = (reshape(da,n_par,size(da,1)/n_par))';
      % Close file.
      fclose(pfid);
      % Current station latitude, longitude.
      r_st=find(c_list(:,1)==stNum);
      latSt = c_list(r_st,colLat);
      lonSt = c_list(r_st,colLon);
      % Current station Pressure.
      presSt=data(:,colPres);
      % Current station depth.
      depSt=sw_dpth(presSt,latSt);
      % Current station quantity.
      if parNum==11
         % ----------Oxygen Concentration from ml/L to  (mkMol/kg).
         OX = data(:,colOxy);
         T=data(:,colTemp);
         S=data(:,colSal);
         P=zeros(size(T));
         % (2) Find and remove "bad" data.
         r_bad=find(T(:,1)==bad_flag | S(:,1)==bad_flag | OX(:,1)==bad_flag);
         if ~isempty(r_bad)
            str_disp=strvcat(str_disp,['Station No ' num2str(stNum) ', ' jN ': ' num2str(length(r_bad)) ' line(s) with "bad" data (flag=-9.990e-29)'],...
              'Pressure  Quantity',num2str([presSt(r_bad,:) parSt(r_bad,:)]));
            T(r_bad,:)=nan;
            S(r_bad,:)=nan;
            P(r_bad,:)=nan;
            OX(r_bad,:)=nan;
            presSt(r_bad,:)=nan;
            depSt(r_bad,:)=nan;
         end
         [swan,sigT]=swstate(S,T,P);
         parSt = OX * 44.66 ./ (1000 + sigT) * 1000; % Oxygen (mkMol/kg)
      elseif parNum==12
         % ----------Oxygen Saturation in %.
         OX = data(:,colOxy);
         T=data(:,colTemp);
         S=data(:,colSal);
         % Find and remove "bad" data.
         r_bad=find(T(:,1)==bad_flag | S(:,1)==bad_flag | OX(:,1)==bad_flag);
         if ~isempty(r_bad)
            str_disp=strvcat(str_disp,['Station No ' num2str(stNum) ', ' jN ': ' num2str(length(r_bad)) ' line(s) with "bad" data (flag=-9.990e-29)'],...
              'Pressure  Quantity',num2str([presSt(r_bad,:) parSt(r_bad,:)]));
            T(r_bad,:)=nan;
            S(r_bad,:)=nan;
            OX(r_bad,:)=nan;
            presSt(r_bad,:)=nan;
         end
         % METHOD 3 for Oxygen solubility calculation in ml/l
         % (Ocean Toolbox, Phil Morgan) 
         soly=sw_satO2(S,T);
         % Oxygen Saturation calculation in %.
         parSt = OX ./ soly *100;
      else
         parSt=data(:,colNum);
         r_bad=find(parSt(:,1)==bad_flag);
         if ~isempty(r_bad)
            str_disp=strvcat(str_disp,['Station No ' num2str(stNum) ', ' jN ': ' num2str(length(r_bad)) ' line(s) with "bad" data (flag=-9.990e-29)'],...
              'Pressure  Quantity',num2str([presSt(r_bad,:) parSt(r_bad,:)]));
            parSt(r_bad,:)=nan;
            presSt(r_bad,:)=nan;
            depSt(r_bad,:)=nan;
         end
      end
%********************
% CHECK: is depth vector monotonic. 
% If not then remove rows with nonmonotonic depth.
%********************
      a=diff(depSt);
      r_nonmon=find(a(:)<=0);
      if ~isempty(r_nonmon)
         countNotUsed = countNotUsed + length(r_nonmon);
         disp(['Station ' num2str(stNum) ': nonmonotone points = ' ...
                 num2str(length(r_nonmon)) ', first "bad" line# is ' num2str(r_nonmon(1))]);
         depSt(r_nonmon+1,:)=[];
         parSt(r_nonmon+1,:)=[];
         return;
      end
%******************
% Interpolate data onto the depth vector.
%******************
      % Depth vector.
      depSti = (z_min_int:z_int_int:z_max_int)';
      % Interpolate one station data onto depth vector.
      parSti = interp1(depSt,parSt,depSti,'linear');
      % Station number vector.
      sti = stNum*ones(size(depSti,1),1);
      % Latitude and longitude vectors.
      lati = latSt*ones(size(depSti,1),1);
      loni = lonSt*ones(size(depSti,1),1);
%********************
% Make a matrix with interpolated data.
%********************
      DATAC = [DATAC;sti depSti parSti lati loni];
      clear parSt depSt presSt sti depSti parSti lati loni;
   end % For for stNum=ctdList
%*****************************************************************
%****************************      PLOT MAPS         *************
%*****************************************************************
   % Loop beginning: enter depth level.
   plotCount = 0;
   countSt=0;
%********************
% Find data at depth.
%********************
   rDep = find(DATAC(:,2)==mapDep);
   DATAI = DATAC(rDep,:);
%********************
% Lon, lat, parameter values.
%********************
   r_nan = find(isnan(DATAI(:,3)));
   DATAI(r_nan,:)=[];
   if size(DATAI,1)>2
      countSt = countSt+size(rDep,1);
      LON = DATAI(:,5);
      LAT = DATAI(:,4);
      PAR = DATAI(:,3);
%********************
% Meshgrid and griddata. 
%********************
      [LONi,LATi] = meshgrid(meshLonMin:meshInt:meshLonMax,meshLatMin:meshInt:meshLatMax);
      PARi = griddata(LON,LAT,PAR,LONi,LATi,'linear');
%********************
% Qantity Min, Max at depth. 
%********************
      parMin = min(PAR(:));
      parMax = max(PAR(:));
%********************
% Qantity Min, Max after Griddata.
%**************************
      parMin_i = min(PARi(find(~isnan(PARi(:)))));
      parMax_i = max(PARi(find(~isnan(PARi(:)))));
      if parMin_i==parMax_i, return; end    
%*********************************************************************************
%**************************        SET AXES           ****************************
%*********************************************************************************
      axPos=POS(countParNum,:);  
      rNum=floor((countParNum-1)/plot_cols)+1;
      cNum=countParNum-(rNum-1)*plot_cols;
      d_x = ddx * (cNum-1); 
      d_y = ddy * (rNum-1);
      h(countParNum)=axes('position',[axPos(1)+axPos(3)*(dx_ax+d_x) axPos(2)+axPos(4)*(dy_ay+d_y) axPos(3)*kx axPos(4)*ky]);
%*************************
% Plot filled contours,
% use shading interpolation.
%*************************
      contourf(LONi,LATi,PARi,vv);
      %contourf3(LONi,LATi,PARi,vv);
      set(gca,'xdir','reverse');
%*************************
% Set pseudocolor axis scaling.
%*************************
      if parNum==2 | parNum==3 | parNum==1 % Salinity , Sigmat, Temp
         caxis([min(vv-0.000001) max(vv-0.000001)]);
         % NOTE: vv-0.00000001 useed to avoid repeating colors between contours
      else
         caxis([min(vv) max(vv)]);
      end
%*************************
% Set shading interpolation.
%*************************
      if strcmp(makeShadingFlat,'yes'), shading flat; end
%*************************
% Plot and label contour lines.
%*************************
         % if strcmp(makeContours,'yes')
         %   hold on;
         %    [cs,h2]=contour(LONi,LATi,PARi,vv_label,'k-');
         %    hold off;
         %    % Put contour labels.
         %    bandW = vv(2)-vv(1);
         %    % Put contour labels.
         %    if ~isempty(cs), clabel_2(cs,'manual'); end
         % end % For if strcmp(makeContours,'yes')
%*********************************************************************************
%**************************       PLOT  LAND          ****************************
%*********************************************************************************
      if strcmp(drawCoast,'yes')
         % Draw land fill.
         hold on;
         lndfill(mm,land_color);
         hold off;
      end
%*************************
% Plot station locations.
%*************************
      hold on;
      plot(LON,LAT,'o','markersize',marker_size,...
          'markerEdgeColor',marker_edge_color,'markerFaceColor',marker_face_color,...
          'linewidth',line_width);
      hold off;
%*****************************
% Text: min, max.
%*****************************
      str_txt(1)={['Min=' num2str(parMin)]};
      str_txt(2)={['Max=' num2str(parMax)]};
      text(lonmin+range([lonmin lonmax])*0.38,latmin+range([latmin latmax])*0.08,str_txt,...
          'FontSize',fnt_text,'horizontalAlignment','left','color',text_min_max_color);
%*************************
% Colorbar.
%*************************
      axPos=get(gca,'position');
      h1 = colorbar;
      set(h1,'ytick',vv,'yLim',[min(vv) max(vv)],...
           'FontSize',fnt_ax,'ticklength',[0.05*0.7 0.025*0.7],...
           'xcolor',text_cb_color,'ycolor',text_cb_color);
      if nband>10
         tlb=f_ticklabel(vv,2);
         set(h1,'ytickLabel',tlb); 
      else
         set(h1,'yticklabel',vv);
      end
      aPos = get(h1,'Position');
      set(h1,'Position',[aPos(1,1)-aPos(1,3)*1.9 aPos(1,2)-aPos(1,4)*0.005 aPos(1,3)*0.4 aPos(1,4)*1]);
%*****************************
% Set axes labels and properties.
%*****************************
      axis([lonmin lonmax latmin latmax]);
      set(gca,'FontSize',fnt_ax,...
        'box','on','PlotBoxAspectRatio',[1 relyx*1 1],...
        'ticklength',get(gca,'ticklength')*1,'tickdir','out','color',ocean_color);
      if countParNum==3, xlabel('Longitude (degW)','FontSize',fnt_lbl); end
      ylabel('Latitude (degN)','FontSize',fnt_lbl);
%*****************************
% Page title.
%*****************************
      if countParNum==1
         str_text={[cruise_title ' CTD, depth ' num2str(mapDep) 'm']};
         text(lonmax-range([lonmax lonmin])*0.5,latmax+range([latmax latmin])*0.2,str_text,...
             'FontSize',fnt_title,....
             'horizontalAlignment','center','color',[0 0 1]);
      end
%******************************
% Subplot title.
%******************************
      str_title={[j_title]};
      title(str_title,...
         'FontSize',fnt_title_sub,'verticalAlignment','bottom',...
         'color',[0 0 1]);
      clear str_title;
%**************************
% Plot axes box 
% (if it gets invisible).
%**************************
      hold on;
      plot([lonmax lonmax],[latmin latmax],'-k',...
        [lonmin lonmin],[latmin latmax],'-k',...
        [lonmax lonmin],[latmin latmin],'-k',...
        [lonmax lonmin],[latmax latmax],'-k');
      hold off;
   else % For if length(rDep)>2
      disp(['No data at depth ' num2str(mapDep) 'm']);
   end % For if length(rDep)>2
   clear DATAC DATAI;
end % for parNum=1:7 
%*******************   
% Get print file - png
%*******************
if strcmp(saveOption,'yes')
    pause(1);
    p=['print -r100 ' cruise_id '_ctd_' num2str(mapDep) 'm.png -dpng'];
    eval(p);
end 
