diff --git a/@tsg_nc/tsg_platform_info.mat b/@tsg_nc/tsg_platform_info.mat
index a424c2f3209f1634c1ea2f6ba88841a870eee1d8..5a747398c5b22dccf2faa18cc587c3c494b1849b 100644
Binary files a/@tsg_nc/tsg_platform_info.mat and b/@tsg_nc/tsg_platform_info.mat differ
diff --git a/M_Csiro/Contents.m b/M_Csiro/Contents.m
new file mode 100644
index 0000000000000000000000000000000000000000..8aa83aec52b2e938b1bc7b9a18f0bad7c1faa815
--- /dev/null
+++ b/M_Csiro/Contents.m
@@ -0,0 +1,63 @@
+% SEAWATER Library 
+% Version 2.0.1   22-Apr-1998
+%
+%                  ******************************* 
+%                  *      SEAWATER Library       * 
+%                  *                             * 
+%                  *       Version 2.0.1         * 
+%                  *      (for Matlab 5.x)       * 
+%                  *                             * 
+%                  *                             *
+%                  *     Phillip P. Morgan       * 
+%                  *           CSIRO             * 
+%                  *                             *
+%                  * Phil.Morgan@marine.csiro.au *
+%                  ******************************* 
+%
+% LIST OF ROUTINES:
+%
+%     SW_NEW     What's new in this version of seawater.
+%
+%     SW_ADTG    Adiabatic temperature gradient 
+%     SW_ALPHA   Thermal expansion coefficient (alpha) 
+%     SW_AONB    Calculate alpha/beta (a on b) 
+%     SW_BETA    Saline contraction coefficient (beta) 
+%     SW_BFRQ    Brunt-Vaisala Frequency Squared (N^2)
+%     SW_COPY    Copyright and Licence file
+%     SW_CP      Heat Capacity (Cp) of Sea Water 
+%     SW_DENS    Density of sea water 
+%     SW_DENS0   Denisty of sea water at atmospheric pressure 
+%     SW_DIST    Distance between two lat, lon coordinates
+%     SW_DPTH    Depth from pressure 
+%     SW_F       Coriolis factor "f" 
+%     SW_FP      Freezing Point of sea water 
+%     SW_G       Gravitational acceleration 
+%     SW_GPAN    Geopotential anomaly  
+%     SW_GVEL    Geostrophic velocity 
+%     SW_INFO    Information on the SEAWATER library. 
+%     SW_PDEN    Potential Density 
+%     SW_PRES    Pressure from depth 
+%     SW_PTMP    Potential temperature 
+%     SW_SALS    Salinity of sea water 
+%     SW_SALT    Salinity from cndr, T, P 
+%     SW_SATAr   Solubility (saturation) of Ar in seawater
+%     SW_SATN2   Solubility (saturation) of N2 in seawater
+%     SW_SATO2   Solubility (saturation) of O2 in seawater
+%     SW_SVAN    Specific volume anomaly 
+%     SW_SVEL    Sound velocity of sea water 
+%     SW_SMOW    Denisty of standard mean ocean water (pure water) 
+%     SW_TEMP    Temperature from potential temperature 
+%     SW_TEST    Run test suite on library
+%     SW_VER     Version number of SEAWATER library
+%
+% LOW LEVEL ROUTINES CALLED BY ABOVE: (also available for you to use)
+%
+%     SW_C3515   Conductivity at (35,15,0) 
+%     SW_CNDR    Conductivity ratio   R = C(S,T,P)/C(35,15,0) 
+%     SW_SALDS   Differiential dS/d(sqrt(Rt)) at constant T. 
+%     SW_SALRP   Conductivity ratio   Rp(S,T,P) = C(S,T,P)/C(S,T,0) 
+%     SW_SALRT   Conductivity ratio   rt(T)     = C(35,T,0)/C(35,15,0) 
+%     SW_SECK    Secant bulk modulus (K) of sea water 
+%=======================================================================
+
+% Contents.m $Revision: 1.6 $  $Date: 1998/04/22 02:12:17 $
diff --git a/M_Csiro/README b/M_Csiro/README
new file mode 100644
index 0000000000000000000000000000000000000000..a2881aa87bc3549c5a322ef4a71c4e8953abe4b2
--- /dev/null
+++ b/M_Csiro/README
@@ -0,0 +1,7 @@
+% README file  $Revision: 1.1 $  $Date: 1994/10/11 02:54:21 $
+%
+%    SEAWATER is a toolkit of MATLAB routines for calculating the
+%    properties of sea water. They are a self contained library and
+%    are extremely easy to use.  
+%
+% See the file sw_info.m for info on installation and usage
diff --git a/M_Csiro/bm_bvfrq.m b/M_Csiro/bm_bvfrq.m
new file mode 100644
index 0000000000000000000000000000000000000000..2de700ed39d8cac1b80f9c849695ee9c61868ec1
--- /dev/null
+++ b/M_Csiro/bm_bvfrq.m
@@ -0,0 +1,95 @@
+function [n2, e] = bm_bvfrq( S, T, P, LAT, DP)
+%*************************************************************
+%                ***** brunt-vaisala freq *****
+%
+%  uses 1980 equation of state
+%
+% in :
+%       P   pressure         decibars
+%       T   temperature      deg celsius (ipts-68)
+%       S   salinity         psu (pss-78)
+%       LAT latitude
+%       DP  intervalle de pression (db) sur lequel est calcule n2
+%           
+% out :
+%       n2  bouyancy freq    cph
+%       e   stability        1/m
+% 
+% r. millard   feb. 1991
+%
+%*************************************************************
+%------
+% BEGIN
+%------
+
+% change call list xlat to glat & pass gravity m/s^2
+% after formulation of n. p. fofonoff & breck owen's
+%
+
+% Calcul de la gravité
+% note that sw_g expects height as argument
+Z  = sw_dpth(P,LAT);
+gp = sw_g(LAT,-Z);
+g2 = gp .* gp;
+g2 = g2 * 1.e-4;
+
+% compute least squares estimate of specific volume anamoly gradient
+
+[m,n] = size(P);
+n2    = zeros(size(P));
+
+% Calcul du nombre d'observation sur lequel est calcule n2
+
+dp = P(2)-P(1);
+nobs = DP / dp;
+
+% Debut du profil
+
+for i = 1:(nobs/2)
+  j = 1:i+nobs/2;
+  [n2(i), e(i)] = cal_bvfrq( S, T, P, g2(i), gp(i), nobs, j);
+end
+
+% Fin du profil
+
+for i = m-(nobs/2):m
+  j = i-nobs/2:m;
+  [n2(i), e(i)] = cal_bvfrq( S, T, P, g2(i), gp(i), nobs, j);
+end
+
+% Profil
+
+for i = (nobs/2+1):m-(nobs/2+1)
+  j = i-nobs/2:i+nobs/2;
+  [n2(i), e(i)] = cal_bvfrq( S, T, P, g2(i), gp(i), nobs, j);
+end
+
+% Corps du calcul
+
+function [n2, e] = cal_bvfrq( S, T, P, g2, gp, nobs, j)
+  
+%  default gravity & rad/sec to cph conversion
+%
+radsec = 572.9578;
+
+  pav  = mean(P(j));
+  data = sw_svan( S(j), sw_ptmp( S(j), T(j), P(j), pav), pav);
+  cxy  = sum( data .* (P(j)-pav) );
+  cy   = sum( data );
+  cxx  = sum( (P(j)-pav).*(P(j)-pav) );
+  
+  sigma = sw_pden(S(j), sw_ptmp( S(j), T(j), P(j), pav), P(j), pav);
+
+  a0    = cxy/cxx;
+
+  v350p = 1./sigma - data;
+  vbar  = v350p(end) + cy/(nobs+1);
+  dvdp  = a0;
+
+  e     = -g2*dvdp/(vbar)^2;
+  n2    = radsec*sqrt(abs(e))*sign(e);
+
+  %  define stability parameter units (1/m)
+  e = e/gp;
+
+  return
\ No newline at end of file
diff --git a/M_Csiro/ctd.txt b/M_Csiro/ctd.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3167e6f8211dad1f583fc054da1e8b6e539039dd
--- /dev/null
+++ b/M_Csiro/ctd.txt
@@ -0,0 +1,31 @@
+   0 25.698 35.221
+  10 26.673 36.106
+  20 26.678 36.106
+  30 26.676 36.107
+  50 24.528 36.561
+  75 22.753 36.614
+ 100 21.427 36.637
+ 125 20.633 36.627
+ 150 19.522 36.558
+ 200 18.798 36.555
+ 250 18.431 36.537
+ 300 18.189 36.526
+ 400 17.726 36.477
+ 500 17.165 36.381
+ 600 15.592 36.105
+ 700 13.458 35.766
+ 800 11.109 35.437
+ 900  8.798 35.178
+1000  6.292 35.044
+1100  5.249 35.004
+1200  4.813 34.995
+1300  4.554 34.986
+1400  4.357 34.977
+1500  4.245 34.975
+1750  4.028 34.973
+2000  3.852 34.975
+2500  3.424 34.968
+3000  2.963 34.946
+3500  2.462 34.920
+4000  2.259 34.904
+4327  2.221 34.896
\ No newline at end of file
diff --git a/M_Csiro/sw_CSIRO.zip b/M_Csiro/sw_CSIRO.zip
new file mode 100644
index 0000000000000000000000000000000000000000..70b24762db14d6795222e1fac3379bb8302ed29d
Binary files /dev/null and b/M_Csiro/sw_CSIRO.zip differ
diff --git a/M_Csiro/sw_adtg.m b/M_Csiro/sw_adtg.m
new file mode 100644
index 0000000000000000000000000000000000000000..9badeac75b134c422d659773581f1d80c3c20311
--- /dev/null
+++ b/M_Csiro/sw_adtg.m
@@ -0,0 +1,119 @@
+
+function ADTG = sw_adtg(S,T,P)
+
+% SW_ADTG    Adiabatic temperature gradient
+%===========================================================================
+% SW_ADTG   $Revision: 1.4 $  $Date: 1994/10/10 04:16:37 $
+%           Copyright (C) CSIRO, Phil Morgan  1992.
+%
+% adtg = sw_adtg(S,T,P)
+%
+% DESCRIPTION:
+%    Calculates adiabatic temperature gradient as per UNESCO 1983 routines.
+%
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78) ]
+%   T = temperature [degree C (IPTS-68)]
+%   P = pressure    [db]
+%       (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+% OUTPUT:
+%   ADTG = adiabatic temperature gradient [degree_C/db]
+%
+% AUTHOR:  Phil Morgan 92-04-03  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonoff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater. Unesco Tech. Pap. in Mar. Sci., No. 44, 53 pp.  Eqn.(31) p.39
+%
+%    Bryden, H. 1973.
+%    "New Polynomials for thermal expansion, adiabatic temperature gradient
+%    and potential temperature of sea water."
+%    DEEP-SEA RES., 1973, Vol20,401-408.
+%=========================================================================
+
+%-------------
+% CHECK INPUTS
+%-------------
+if nargin ~= 3
+   error('sw_adtg.m: Must pass 3 parameters ')
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+%-------------
+% BEGIN
+%-------------
+a0 =  3.5803E-5;
+a1 = +8.5258E-6;
+a2 = -6.836E-8;
+a3 =  6.6228E-10;
+
+b0 = +1.8932E-6;
+b1 = -4.2393E-8;
+
+c0 = +1.8741E-8;
+c1 = -6.7795E-10;
+c2 = +8.733E-12;
+c3 = -5.4481E-14;
+
+d0 = -1.1351E-10;
+d1 =  2.7759E-12;
+
+e0 = -4.6206E-13;
+e1 = +1.8676E-14;
+e2 = -2.1687E-16;
+
+ADTG =      a0 + (a1 + (a2 + a3.*T).*T).*T ... 
+         + (b0 + b1.*T).*(S-35)  ...
+	 + ( (c0 + (c1 + (c2 + c3.*T).*T).*T) + (d0 + d1.*T).*(S-35) ).*P ...
+         + (  e0 + (e1 + e2.*T).*T ).*P.*P;
+
+if Transpose
+   ADTG = ADTG';
+end %if
+
+return
+%==========================================================================
+
diff --git a/M_Csiro/sw_alpha.m b/M_Csiro/sw_alpha.m
new file mode 100644
index 0000000000000000000000000000000000000000..828bf3bbfb4273aa71ce86296890cdb22a8a8972
--- /dev/null
+++ b/M_Csiro/sw_alpha.m
@@ -0,0 +1,105 @@
+
+function [ALPHA] = sw_alpha(S, T, P, keyword)
+
+% SW_ALPHA   Thermal expansion coefficient (alpha)
+%================================================================
+% SW_ALPHA  $Revision: 1.6 $   $Date: 1998/04/21 05:42:10 $
+%           Copyright (C) CSIRO, Nathan Bindoff 1993.
+%
+% USAGE:  [ALPHA] = alpha(S, T, P, keyword)
+%
+%         [ALPHA] = alpha(S, T,    P, 'temp')    %default
+%         [ALPHA] = alpha(S, PTMP, P, 'ptmp')
+%
+% DESCRIPTION:
+%    A function to calculate the thermal expansion coefficient.
+%
+% INPUT:
+%   S       = salinity              [psu      (PSS-78) ]
+% * PTMP    = potential temperature [degree C (IPTS-68)]
+% * T       = temperature           [degree C (IPTS-68)]
+%   P       = pressure              [db]
+%             (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+%   keyword = optional string to identify if temp or ptmp passed.
+%           = No argument defaults to 'temp'
+%           = 'temp' assumes (S,T,P) passed.    Will execute slower
+%                    as ptmp will be calculated internally.
+%           = 'ptmp' assumes (S,PTMP,P) passed. Will execute faster.
+%
+% OUTPUT:
+%   ALPHA = Thermal expansion coeff (alpha) [degree_C.^-1]
+%
+% AUTHOR:   N.L. Bindoff  1993
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCE: 
+%    McDougall, T.J. 1987.  "Neutral Surfaces"
+%    Journal of Physical Oceanography vol 17 pages 1950-1964, 
+%
+% CHECK VALUE:
+%    See sw_beta.m amd sw_aonb.m
+%================================================================
+
+% Modifications
+% 93-04-22. Phil Morgan,  Help display modified to suit library
+% 93-04-23. Phil Morgan,  Input argument checking
+% 94-10-15. Phil Morgan,  Pass S,T,P and keyword for 'ptmp'
+
+  
+% CHECK INPUT ARGUMENTS
+if ~(nargin==3 | nargin==4)
+  error('sw_alpha.m: requires 3 or 4 input arguments')
+end %if
+if nargin == 3
+  keyword = 'temp';
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+% BEGIN
+
+ALPHA = sw_aonb(S,T,P,keyword).*sw_beta(S,T,P,keyword);
+
+return
+%------------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_aonb.m b/M_Csiro/sw_aonb.m
new file mode 100644
index 0000000000000000000000000000000000000000..e1f3370ca58fba7b79639ec5dec67ab87506c373
--- /dev/null
+++ b/M_Csiro/sw_aonb.m
@@ -0,0 +1,130 @@
+
+function [AONB]= aonb(S, T, P, keyword)
+
+% SW_AONB    Calculate alpha/beta (a on b)
+%================================================================
+% SW_AONB   $Revision: 1.5 $   $Date: 1994/11/15 04:10:45 $
+%           Copyright (C) CSIRO, Nathan Bindoff 1993
+%
+% USAGE: [AONB] = aonb(S, T, P, {keyword} )
+%
+%        [AONB] = aonb(S, T,    P, 'temp' )      %default
+%        [AONB] = aonb(S, PTMP, P, 'ptmp' )
+%
+% DESCRIPTION
+%    Calculate alpha/beta.  See sw_alpha.m and sw_beta.m
+%
+% INPUT:  (all must have same dimensions)
+%   S       = salinity              [psu      (PSS-78) ]
+% * PTMP    = potential temperature [degree C (IPTS-68)]
+% * T       = temperature           [degree C (IPTS-68)]
+%   P       = pressure              [db]
+%             (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+%   keyword = optional string to identify if temp or ptmp passed.
+%           = No argument defaults to 'temp'
+%           = 'temp' assumes (S,T,P) passed.    Will execute slower
+%                    as ptmp will be calculated internally.
+%           = 'ptmp' assumes (S,PTMP,P) passed. Will execute faster.
+%
+% OUTPUT
+%   AONB  = alpha/beta [psu/degree_C]
+%
+% AUTHOR:   N.L. Bindoff  1993
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCE: 
+%    McDougall, T.J. 1987. "Neutral Surfaces"
+%    Journal of Physical Oceanography vol 17 pages 1950-1964, 
+%
+% CHECK VALUE: 
+%    aonb=0.34763 psu C^-1 at S=40.0 psu, ptmp=10.0 C, p=4000 db
+%================================================================
+
+% Modifications
+% 93-04-22. Phil Morgan,  Help display modified to suit library
+% 93-04-23. Phil Morgan,  Input argument checking
+% 94-10-15. Phil Morgan,  Pass S,T,P and keyword for 'ptmp'
+
+% CHECK INPUT ARGUMENTS
+if ~(nargin==3 | nargin==4)
+  error('sw_aonb.m: requires 3 input arguments')
+end %if
+if nargin == 3
+  keyword = 'temp';
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+% ENSURE WE USE PTMP IN CALCULATIONS
+if strcmp(lower(keyword),'ptmp')
+  % already have ptmp
+else
+  T = sw_ptmp(S,T,P,0); % now have ptmp
+end %if
+
+% BEGIN
+	 c1=fliplr([ 0.665157e-1, 0.170907e-1, ...
+	    -0.203814e-3, 0.298357e-5, ...
+            -0.255019e-7]);
+         c2=fliplr([ 0.378110e-2, ...
+            -0.846960e-4]);
+         c2a=fliplr([0.0 -0.164759e-6, ...
+            -0.251520e-11]);
+         c3=[-0.678662e-5];
+         c4=fliplr([+0.380374e-4, -0.933746e-6, ...
+            +0.791325e-8]);
+         c5=[0.512857e-12];
+         c6=[-0.302285e-13];
+%
+% Now calaculate the thermal expansion saline contraction ratio adb
+%
+        [m,n] = size(S);
+        sm35  = S-35.0*ones(m,n);
+        AONB  = polyval(c1,T) + sm35.*(polyval(c2,T)...
+               + polyval(c2a,P)) ...
+               + sm35.^2*c3 + P.*polyval(c4,T) ...
+               + c5*(P.^2).*(T.^2) + c6*P.^3;
+
+return
+%----------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_beta.m b/M_Csiro/sw_beta.m
new file mode 100644
index 0000000000000000000000000000000000000000..3be9894ad8686b439f3771bddb689471d53d2cc3
--- /dev/null
+++ b/M_Csiro/sw_beta.m
@@ -0,0 +1,126 @@
+
+function [BETA] = sw_beta(S, T, P, keyword)
+
+% SW_BETA    Saline contraction coefficient (beta)
+%========================================================================
+% SW_BETA  $Revision: 1.4 $   $Date: 1994/11/15 04:10:05 $
+%   %      Copyright (C) CSIRO, Nathan Bindoff 1993.
+%
+% USAGE:  [BETA] = sw_beta(S, T, P, {keyword} )
+%
+%         [BETA] = sw_beta(S, T,    P, 'temp')     %default
+%         [BETA] = sw_beta(S, PTMP, P, 'ptmp')
+%
+% DESCRIPTION
+%   The saline contraction coefficient as defined by T.J. McDougall.
+%
+% INPUT:  (all must have same dimensions)
+%   S       = salinity              [psu      (PSS-78) ]
+% * PTMP    = potential temperature [degree C (IPTS-68)]
+% * T       = temperature           [degree C (IPTS-68)]
+%   P       = pressure              [db]
+%             (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+%   keyword = optional string to identify if temp or ptmp passed.
+%           = No argument defaults to 'temp'
+%           = 'temp' assumes (S,T,P) passed.    Will execute slower
+%                    as ptmp will be calculated internally.
+%           = 'ptmp' assumes (S,PTMP,P) passed. Will execute faster.
+%
+% OUTPUT
+%   BETA = Saline Contraction Coefficient  [psu.^-1]
+%
+% AUTHOR:   N.L. Bindoff  1993
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCE: 
+%    McDougall, T.J. 1987. "Neutral Surfaces"
+%    Journal of Physical Oceanography vol 17 pages 1950-1964, 
+%
+% CHECK VALUE: 
+%    beta=0.72088e-3 psu.^-1 at S=40.0 psu, ptmp = 10.0 C, p=4000 db
+%========================================================================
+
+% Modifications
+% 93-04-22. Phil Morgan,  Help display modified to suit library
+% 93-04-23. Phil Morgan,  Input argument checking
+% 94-10-15. Phil Morgan,  Pass S,T,P and keyword for 'ptmp'
+
+% CHECK INPUT ARGUMENTS
+if ~(nargin==3 | nargin==4)
+  error('sw_beta.m: requires 3 or 4 input arguments')
+end %if
+if nargin == 3
+  keyword = 'temp';
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+% ENSURE WE USE PTMP IN CALCULATIONS
+if strcmp(lower(keyword),'ptmp')
+  % already have ptmp
+else
+  T = sw_ptmp(S,T,P,0); % now have ptmp
+end %if
+
+% BEGIN
+
+        c1=fliplr([ 0.785567e-3, -0.301985e-5 ...
+	     0.555579e-7, -0.415613e-9]);
+	c2=fliplr([ -0.356603e-6, 0.788212e-8]);
+	c3=fliplr([0.0 0.408195e-10, -0.602281e-15]);
+	c4=[0.515032e-8];
+	c5=fliplr([-0.121555e-7, 0.192867e-9, -0.213127e-11]);
+        c6=fliplr([0.176621e-12 -0.175379e-14]);
+	c7=[0.121551e-17];
+%
+% Now calaculate the thermal expansion saline contraction ratio adb
+%
+	[m,n] = size(S);
+	sm35  = S-35*ones(m,n);
+	BETA  = polyval(c1,T) + sm35.*(polyval(c2,T) + ...
+	        polyval(c3,P)) + c4*(sm35.^2) + ...
+	        P.*polyval(c5,T) + (P.^2).*polyval(c6,T) ...
+                +c7*( P.^3);
+
+return
+%------------------------------------------------------------------------
diff --git a/M_Csiro/sw_bfrq.m b/M_Csiro/sw_bfrq.m
new file mode 100644
index 0000000000000000000000000000000000000000..1c3d26210d866dbae192202448bb64731ca782b8
--- /dev/null
+++ b/M_Csiro/sw_bfrq.m
@@ -0,0 +1,164 @@
+
+function [n2,q,p_ave] = sw_bfrq(S,T,P,LAT)
+
+% SW_BFRQ    Brunt-Vaisala Frequency Squared (N^2)
+%===========================================================================
+% SW_BFRQ  $Revision: 1.12 $   $Date: 1994/11/15 04:13:34 $
+%          Copyright (C) CSIRO, Phil Morgan  1993. 
+%
+% USAGE:  [bfrq,vort,p_ave] = sw_bfrq(S,T,P,{LAT}) 
+%
+% DESCRIPTION:
+%    Calculates Brunt-Vaisala Frequency squared (N^2) at the mid depths
+%    from the equation,
+%
+%               -g      d(pdens)
+%         N2 =  ----- x --------
+%               pdens     d(z)
+%
+%    Also returns Potential Vorticity from q = f*N2/g.  
+%
+% INPUT:  (all must have same dimensions MxN)
+%   S   = salinity    [psu      (PSS-78) ]
+%   T   = temperature [degree C (IPTS-68)]
+%   P   = pressure    [db]   
+%
+%   OPTIONAL:
+%      LAT     = Latitude in decimal degrees north [-90..+90]
+%                May have dimensions 1x1 or 1xn where S(mxn).
+%                (Will use sw_g instead of the default g=9.8 m^2/s)
+%                (Will also calc d(z) instead of d(p) in numerator)
+% OUTPUT:
+%   bfrq  = Brunt-Vaisala Frequency squared (M-1xN)  [s^-2]
+%   vort  = Planetary Potential Vorticity   (M-1xN)  [(ms)^-1]
+%           (if isempty(LAT) vort=NaN ) 
+%   p_ave = Mid pressure between P grid     (M-1xN)  [db]
+%
+% AUTHOR:  Phil Morgan 93-06-24  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%   A.E. Gill 1982. p.54  eqn 3.7.15
+%   "Atmosphere-Ocean Dynamics"
+%   Academic Press: New York.  ISBN: 0-12-283522-0
+%
+%   Jackett, D.R. and McDougall, T.J. 1994.
+%   Minimal adjustment of hydrographic properties to achieve static
+%   stability.  submitted J.Atmos.Ocean.Tech.
+%
+%   Greg Johnson (gjohnson@pmel.noaa.gov) 
+%                added potential vorticity calcuation
+%=========================================================================
+
+% CALLER:  general purpose
+% CALLEE:  sw_dens.m sw_pden.m
+
+%$Id: sw_bfrq.M,v 1.12 1994/11/15 04:13:34 morgan Exp $
+
+%-------------
+% CHECK INPUTS
+%-------------
+if ~(nargin==3 | nargin==4) 
+   error('sw_bfrq.m: Must pass 3 or 4 parameters ')
+end %if
+
+if nargin == 3
+  LAT = [];
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+% IF LAT PASSED THEN VERIFY DIMENSIONS
+if ~isempty(LAT)
+   [mL,nL] = size(LAT);
+   if mL==1 & nL==1
+      LAT = LAT*ones(size(S));
+   %end  % Je commente le end et remplace le if suivant par elseif
+
+    elseif (ms~=mL) | (ns~=nL)              % S & LAT are not the same shape
+       if (ns==nL) & (mL==1)           % copy LATS down each column
+          LAT = LAT( ones(1,ms), : );  % s.t. dim(S)==dim(LAT)    
+       else
+          error('sw_bfrq.m:  Inputs arguments have wrong dimensions')
+       end %if
+   end %if
+end %if
+
+   
+
+%------
+% BEGIN
+%------
+if ~isempty(LAT)
+   % note that sw_g expects height as argument
+   Z = sw_dpth(P,LAT);
+   g = sw_g(LAT,-Z);
+   f = sw_f(LAT);
+else
+   Z = P;
+   g = 9.8*ones(size(P));
+   f = NaN*ones(size(P));
+end %if
+
+[m,n] = size(P);
+iup   = 1:m-1;
+ilo   = 2:m;
+p_ave = (P(iup,:)+P(ilo,:) )/2;
+pden_up = sw_pden(S(iup,:),T(iup,:),P(iup,:),p_ave);
+pden_lo = sw_pden(S(ilo,:),T(ilo,:),P(ilo,:),p_ave);
+ 
+mid_pden = (pden_up + pden_lo )/2;
+dif_pden = pden_up - pden_lo;
+mid_g    = (g(iup,:)+g(ilo,:))/2;
+dif_z    = diff(Z);
+n2       = -mid_g .* dif_pden ./ (dif_z .* mid_pden);
+
+mid_f    = f(iup,:);
+q        = mid_f .* dif_pden ./  (dif_z .* mid_pden);
+  
+if Transpose
+  n2    = n2';
+  q     = q';
+  p_ave = p_ave';
+end %if
+return
+%-------------------------------------------------------------------
diff --git a/M_Csiro/sw_bfrq_mean.m b/M_Csiro/sw_bfrq_mean.m
new file mode 100644
index 0000000000000000000000000000000000000000..816a3879ee092b115b434e67032e6899d599de73
--- /dev/null
+++ b/M_Csiro/sw_bfrq_mean.m
@@ -0,0 +1,173 @@
+
+function [n2,q,p_ave] = sw_bfrq_mean(S,T,P,LAT)
+
+% SW_BFRQ    Brunt-Vaisala Frequency Squared (N^2)
+%===========================================================================
+% SW_BFRQ  $Revision: 1.12 $   $Date: 1994/11/15 04:13:34 $
+%          Copyright (C) CSIRO, Phil Morgan  1993. 
+%
+% USAGE:  [bfrq,vort,p_ave] = sw_bfrq(S,T,P,{LAT}) 
+%
+% DESCRIPTION:
+%    Calculates Brunt-Vaisala Frequency squared (N^2) at the mid depths
+%    from the equation,
+%
+%               -g      d(pdens)
+%         N2 =  ----- x --------
+%               pdens     d(z)
+%
+%    Also returns Potential Vorticity from q = f*N2/g.  
+%
+% INPUT:  (all must have same dimensions MxN)
+%   S   = salinity    [psu      (PSS-78) ]
+%   T   = temperature [degree C (IPTS-68)]
+%   P   = pressure    [db]   
+%
+%   OPTIONAL:
+%      LAT     = Latitude in decimal degrees north [-90..+90]
+%                May have dimensions 1x1 or 1xn where S(mxn).
+%                (Will use sw_g instead of the default g=9.8 m^2/s)
+%                (Will also calc d(z) instead of d(p) in numerator)
+% OUTPUT:
+%   bfrq  = Brunt-Vaisala Frequency squared (M-1xN)  [s^-2]
+%   vort  = Planetary Potential Vorticity   (M-1xN)  [(ms)^-1]
+%           (if isempty(LAT) vort=NaN ) 
+%   p_ave = Mid pressure between P grid     (M-1xN)  [db]
+%
+% AUTHOR:  Phil Morgan 93-06-24  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%   A.E. Gill 1982. p.54  eqn 3.7.15
+%   "Atmosphere-Ocean Dynamics"
+%   Academic Press: New York.  ISBN: 0-12-283522-0
+%
+%   Jackett, D.R. and McDougall, T.J. 1994.
+%   Minimal adjustment of hydrographic properties to achieve static
+%   stability.  submitted J.Atmos.Ocean.Tech.
+%
+%   Greg Johnson (gjohnson@pmel.noaa.gov) 
+%                added potential vorticity calcuation
+%=========================================================================
+
+% CALLER:  general purpose
+% CALLEE:  sw_dens.m sw_pden.m
+
+%$Id: sw_bfrq.M,v 1.12 1994/11/15 04:13:34 morgan Exp $
+
+%-------------
+% CHECK INPUTS
+%-------------
+if ~(nargin==3 | nargin==4) 
+   error('sw_bfrq.m: Must pass 3 or 4 parameters ')
+end %if
+
+if nargin == 3
+  LAT = [];
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+% IF LAT PASSED THEN VERIFY DIMENSIONS
+if ~isempty(LAT)
+   [mL,nL] = size(LAT);
+   if mL==1 & nL==1
+      LAT = LAT*ones(size(S));
+   %end  % Je commente le end et remplace le if suivant par elseif
+
+    elseif (ms~=mL) | (ns~=nL)              % S & LAT are not the same shape
+       if (ns==nL) & (mL==1)           % copy LATS down each column
+          LAT = LAT( ones(1,ms), : );  % s.t. dim(S)==dim(LAT)    
+       else
+          error('sw_bfrq.m:  Inputs arguments have wrong dimensions')
+       end %if
+   end %if
+end %if
+
+   
+
+%------
+% BEGIN
+%------
+if ~isempty(LAT)
+   % note that sw_g expects height as argument
+   Z = sw_dpth(P,LAT);
+   g = sw_g(LAT,-Z);
+   f = sw_f(LAT);
+else
+   Z = P;
+   g = 9.8*ones(size(P));
+   f = NaN*ones(size(P));
+end %if
+
+[m,n] = size(P);
+iup   = 1:m-1;
+ilo   = 2:m;
+p_ave = (P(iup,:)+P(ilo,:) )/2;
+pden_up = sw_pden(S(iup,:),T(iup,:),P(iup,:),p_ave);
+pden_lo = sw_pden(S(ilo,:),T(ilo,:),P(ilo,:),p_ave);
+ 
+mid_pden = (pden_up + pden_lo )/2;
+dif_pden = pden_up - pden_lo;
+mid_g    = (g(iup,:)+g(ilo,:))/2;
+dif_z    = diff(Z);
+
+% Calcul des moyennes
+% -------------------
+Mmid_pden = mean(mid_pden,2);
+Mdif_pden = mean(dif_pden,2);
+Mmid_g = mid_g(:,1);
+Mdif_z = dif_z(:,1);
+
+n2       = -Mmid_g .* Mdif_pden ./ (Mdif_z .* Mmid_pden);
+
+mid_f    = f(iup,:);
+Mmid_f   = mid_f(:,1);
+q        = Mmid_f .* Mdif_pden ./  (Mdif_z .* Mmid_pden);
+  
+if Transpose
+  n2    = n2';
+  q     = q';
+  p_ave = p_ave';
+end %if
+return
+%-------------------------------------------------------------------
diff --git a/M_Csiro/sw_bfrq_old.m b/M_Csiro/sw_bfrq_old.m
new file mode 100644
index 0000000000000000000000000000000000000000..72714c5148f03a0c6970a98e8ea4e5d39ae5e9d2
--- /dev/null
+++ b/M_Csiro/sw_bfrq_old.m
@@ -0,0 +1,164 @@
+
+function [n2,q,p_ave] = sw_bfrq(S,T,P,LAT)
+
+% SW_BFRQ    Brunt-Vaisala Frequency Squared (N^2)
+%===========================================================================
+% SW_BFRQ  $Revision: 1.12 $   $Date: 1994/11/15 04:13:34 $
+%          Copyright (C) CSIRO, Phil Morgan  1993. 
+%
+% USAGE:  [bfrq,vort,p_ave] = sw_bfrq(S,T,P,{LAT}) 
+%
+% DESCRIPTION:
+%    Calculates Brunt-Vaisala Frequency squared (N^2) at the mid depths
+%    from the equation,
+%
+%               -g      d(pdens)
+%         N2 =  ----- x --------
+%               pdens     d(z)
+%
+%    Also returns Potential Vorticity from q = f*N2/g.  
+%
+% INPUT:  (all must have same dimensions MxN)
+%   S   = salinity    [psu      (PSS-78) ]
+%   T   = temperature [degree C (IPTS-68)]
+%   P   = pressure    [db]   
+%
+%   OPTIONAL:
+%      LAT     = Latitude in decimal degrees north [-90..+90]
+%                May have dimensions 1x1 or 1xn where S(mxn).
+%                (Will use sw_g instead of the default g=9.8 m^2/s)
+%                (Will also calc d(z) instead of d(p) in numerator)
+% OUTPUT:
+%   bfrq  = Brunt-Vaisala Frequency squared (M-1xN)  [s^-2]
+%   vort  = Planetary Potential Vorticity   (M-1xN)  [(ms)^-1]
+%           (if isempty(LAT) vort=NaN ) 
+%   p_ave = Mid pressure between P grid     (M-1xN)  [db]
+%
+% AUTHOR:  Phil Morgan 93-06-24  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%   A.E. Gill 1982. p.54  eqn 3.7.15
+%   "Atmosphere-Ocean Dynamics"
+%   Academic Press: New York.  ISBN: 0-12-283522-0
+%
+%   Jackett, D.R. and McDougall, T.J. 1994.
+%   Minimal adjustment of hydrographic properties to achieve static
+%   stability.  submitted J.Atmos.Ocean.Tech.
+%
+%   Greg Johnson (gjohnson@pmel.noaa.gov) 
+%                added potential vorticity calcuation
+%=========================================================================
+
+% CALLER:  general purpose
+% CALLEE:  sw_dens.m sw_pden.m
+
+%$Id: sw_bfrq.M,v 1.12 1994/11/15 04:13:34 morgan Exp $
+
+%-------------
+% CHECK INPUTS
+%-------------
+if ~(nargin==3 | nargin==4) 
+   error('sw_bfrq.m: Must pass 3 or 4 parameters ')
+end %if
+
+if nargin == 3
+  LAT = [];
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+% IF LAT PASSED THEN VERIFY DIMENSIONS
+if ~isempty(LAT)
+   [mL,nL] = size(LAT);
+   if mL==1 & nL==1
+      LAT = LAT*ones(size(S));
+   end %if 
+
+   if (ms~=mL) | (ns~=nL)              % S & LAT are not the same shape
+       if (ns==nL) & (mL==1)           % copy LATS down each column
+          LAT = LAT( ones(1,ms), : );  % s.t. dim(S)==dim(LAT)    
+       else
+          error('sw_bfrq.m:  Inputs arguments have wrong dimensions')
+       end %if
+   end %if
+end %if
+
+   
+
+%------
+% BEGIN
+%------
+if ~isempty(LAT)
+   % note that sw_g expects height as argument
+   Z = sw_dpth(P,LAT);
+   g = sw_g(LAT,-Z);
+   f = sw_f(LAT);
+else
+   Z = P;
+   g = 9.8*ones(size(P));
+   f = NaN*ones(size(P));
+end %if
+
+[m,n] = size(P);
+iup   = 1:m-1;
+ilo   = 2:m;
+p_ave = (P(iup,:)+P(ilo,:) )/2;
+pden_up = sw_pden(S(iup,:),T(iup,:),P(iup,:),p_ave);
+pden_lo = sw_pden(S(ilo,:),T(ilo,:),P(ilo,:),p_ave);
+ 
+mid_pden = (pden_up + pden_lo )/2;
+dif_pden = pden_up - pden_lo;
+mid_g    = (g(iup,:)+g(ilo,:))/2;
+dif_z    = diff(Z);
+n2       = -mid_g .* dif_pden ./ (dif_z .* mid_pden);
+
+mid_f    = f(iup,:);
+q        = mid_f .* dif_pden ./  (dif_z .* mid_pden);
+  
+if Transpose
+  n2    = n2';
+  q     = q';
+  p_ave = p_ave';
+end %if
+return
+%-------------------------------------------------------------------
diff --git a/M_Csiro/sw_c3515.m b/M_Csiro/sw_c3515.m
new file mode 100644
index 0000000000000000000000000000000000000000..0dd8a57a6123d6fbe4ac9548d1e71bcfbda97d9d
--- /dev/null
+++ b/M_Csiro/sw_c3515.m
@@ -0,0 +1,39 @@
+
+function c3515 = sw_c3515()
+
+% SW_C3515   Conductivity at (35,15,0)
+%=========================================================================
+% SW_c3515  $Revision: 1.4 $   $Date: 1994/10/10 04:35:22 $
+%       %   Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE:  c3515 = sw_c3515
+%
+% DESCRIPTION:
+%   Returns conductivity at S=35 psu , T=15 C [ITPS 68] and P=0 db).
+%
+% INPUT: (none)
+%
+% OUTPUT:
+%   c3515  = Conductivity   [mmho/cm == mS/cm] 
+% 
+% AUTHOR:  Phil Morgan 93-04-17  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    R.C. Millard and K. Yang 1992.
+%    "CTD Calibration and Processing Methods used by Woods Hole 
+%     Oceanographic Institution"  Draft April 14, 1992
+%    (Personal communication)
+%=========================================================================
+
+% CALLER: none
+% CALLEE: none
+%
+
+c3515 = 42.914;
+
+return
+%-------------------------------------------------------------------------
diff --git a/M_Csiro/sw_cndr.m b/M_Csiro/sw_cndr.m
new file mode 100644
index 0000000000000000000000000000000000000000..0060af820af5c7c26a7637a5bdbdcef7f0f478b5
--- /dev/null
+++ b/M_Csiro/sw_cndr.m
@@ -0,0 +1,145 @@
+
+function R = sw_cndr(S,T,P)
+
+% SW_CNDR    Conductivity ratio   R = C(S,T,P)/C(35,15,0)
+%=========================================================================
+% SW_CNDR  $Revision: 1.3 $  $Date: 1994/10/10 04:36:58 $
+%          Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE:  cndr = sw_cndr(S,T,P)
+%
+% DESCRIPTION:
+%   Calculates conductivity ratio from S,T,P.
+%   
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78) ]
+%   T = temperature [degree C (IPTS-68)]
+%   P = pressure    [db]
+%       (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+% OUTPUT:
+%   cndr = Conductivity ratio     R =  C(S,T,P)/C(35,15,0) [no units]
+% 
+% AUTHOR:  Phil Morgan 93-04-21  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonoff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%=========================================================================
+
+% CALLER: general purpose
+% CALLEE: sw_salds.m sw_sals.m sw_salrt.m
+
+%--------------
+% check inputs
+%-------------
+if nargin~=3
+  error('sw_cndr.m: must have 3 input arguments')
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+%-------
+% BEGIN
+%-------
+del_T = T - 15;
+
+for i = 1:ms
+  for j = 1:ns
+    %---------------------------------------------------------------------
+    % DO A NEWTON-RAPHSON ITERATION FOR INVERSE INTERPOLATION OF Rt FROM S.
+    %---------------------------------------------------------------------
+    S_loop  = S(i,j);  % S in the loop
+    T_loop  = T(i,j);  % T in the loop
+    Rx_loop = sqrt(S_loop/35.0);                % first guess at Rx = sqrt(Rt)
+    SInc    = sw_sals(Rx_loop.*Rx_loop,T_loop); % S INCrement (guess) from Rx
+    iloop    = 0;
+    end_loop = 0;
+    while ~end_loop
+       Rx_loop = Rx_loop + (S_loop - SInc)./sw_salds(Rx_loop,del_T(i,j));
+       SInc    = sw_sals(Rx_loop.*Rx_loop,T_loop);
+       iloop   = iloop + 1;
+       dels    = abs(SInc-S_loop);
+       if (dels>1.0e-4 & iloop<10) 
+          end_loop = 0;
+       else
+          end_loop = 1;
+       end %if
+    end %while  
+
+    Rx(i,j) = Rx_loop;
+
+  end %for j
+end %for i
+
+%------------------------------------------------------
+% ONCE Rt FOUND, CORRESPONDING TO EACH (S,T) EVALUATE R
+%------------------------------------------------------
+% eqn(4) p.8 Unesco 1983
+
+d1 =  3.426e-2;
+d2 =  4.464e-4;
+d3 =  4.215e-1;
+d4 = -3.107e-3;
+
+e1 =  2.070e-5;
+e2 = -6.370e-10;
+e3 =  3.989e-15;
+
+A  = (d3 + d4.*T);
+B  = 1 + d1.*T + d2.*T.^2;
+C  = P.*(e1 + e2.*P + e3.*P.^2); 
+
+% eqn(6) p.9 UNESCO 1983.
+Rt    = Rx.*Rx;
+rt    = sw_salrt(T);
+Rtrt  = rt.*Rt;
+D     = B - A.*rt.*Rt;
+E     = rt.*Rt.*A.*(B+C);
+R     = sqrt(abs(D.^2+4*E)) - D;
+R     = 0.5*R./A;
+
+return
+%-----------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_copy.m b/M_Csiro/sw_copy.m
new file mode 100644
index 0000000000000000000000000000000000000000..69a0229a7ddb964aeac58bd55b9b98b96b7b6ce3
--- /dev/null
+++ b/M_Csiro/sw_copy.m
@@ -0,0 +1,165 @@
+% SW_COPY    Copyright and licence information on SEAWATER library.
+% =================================================================
+% SW_COPY  $Revision: 1.6 $  $Date: 1994/12/06 03:05:54 $
+%          Copyright (C) Phil Morgan 1993
+%
+%      SOFTWARE LICENCE AGREEMENT
+%
+%      1.0 Grant of Licence
+%
+%      1.1  The CSIRO Division of  Oceanography (herein referred to as
+%           "CSIRO") hereby grants you (hereinafter  referred to  as 
+%           the  "Licencee"),  subject  to  the  Licencee agreeing  to
+%           comply with the terms and  conditions of this Agreement, a
+%           non-transferable,   non-exclusive  licence   to  use   the
+%           computer  programs described in this document (hereinafter
+%           referred to  as the   "Software")  for  the  purpose   of  
+%           the  Licencee's computing activity.
+%
+%      1.2  CSIRO hereby grants the Licencee the right to  make copies
+%           of  the  Software  for   the  purpose  of  the  Licencee's
+%           computing activity only.
+%
+%      1.3  The benefit of  the rights granted to the Licencee  by the
+%           Licence and this Agreement  generally shall be personal to
+%           the Licencee and the  Licencee shall not mortgage, charge,
+%           assign,  rent,  lease,  sell or  otherwise  dispose  of or
+%           transfer the same or any part to any third party.
+%
+%      1.4  Unless  otherwise agreed  in  writing or  provided  for in
+%           this  Agreement, CSIRO  shall  be under  no  obligation or
+%           responsibility  to provide the Licencee with any training,
+%           maintenance  services,  enhancements  or  updates  of  the
+%           Software or any services whatsoever.
+%
+%      2.0 Acknowledgment by the Licencee
+%
+%      2.1  The Licencee acknowledges and agrees that it shall not:
+%
+%           (i)   sell,  let for  hire or  by way  of  trade, offer  or
+%                 exhibit  or  expose  for sale  or  hire  or otherwise
+%  	          distribute the Software for  the purposes of trade or
+%                 any other purpose;
+%  
+%           (ii)  authorise  or assist any third person  to do any
+%                 of the acts set out in (i)  above;  
+%
+%           (iii) modify  the  Software source  code  without  advising
+%                 CSIRO. 
+%
+%      2.2 The Licencee agrees that:
+%
+%           (a)  CSIRO  is  the  owner  of  all  copyright  and  other
+%                Intellectual  Property   Rights  subsisting  in   the
+%                Software;
+%       
+%           (b)  this   document  must  be   properly   cited  in  any
+%                publication  reporting  results  derived   from  this
+%                document or obtained from application and use of this
+%                software. Any  of   the  Licencee's  documentation   
+%                describing  results  generated  by  the  Licencee's  
+%                use  of  the Software will contain an acknowledgement 
+%                of CSIRO's  ownership of the Software;
+%
+%           (c)  CSIRO reserves all rights  in the Software other than
+%                the  rights   granted  to   the   Licencee  by   this
+%                Agreement;
+%
+%           (d)  each  item  of  the Software  will  display  a  banner
+%                summarising   the  terms   of   this   Agreement  and
+%                acknowledging  the source  of  the Software,  and the
+%                contents of  a banner  will not be  modified and  its
+%                display  will  not  be  inactivated  by  the Licencee
+%                without the approval of CSIRO.
+%
+%      3.0 Indemnity
+%
+%      3.1  To the full  extent permitted by  law, CSIRO  excludes any
+%           and  all liability  in  respect  of any  loss  or  damage,
+%           whether  personal  (includes  death  or  illness)  or  of
+%           property  and  whether direct,  consequential  or  special
+%           (including  consequential  financial  loss or  damage) of
+%           the Licencee, its  officers, agents and  employees or  any
+%           third party  howsoever caused,  which may  be suffered  or
+%           incurred  or which  may  arise directly  or  indirectly in
+%           respect of  or  arising  out  of  the  Licencee's  use  or
+%           inability to use  the Software or the failure or  omission
+%           on the  part of CSIRO  to comply  with the conditions  and
+%           warranties under  this  Licence  Agreement.    Insofar  as
+%           liability for  loss or damages  under or  pursuant to such
+%           legislation cannot  be  excluded,  CSIRO's  liability  for
+%           loss or  damages shall  be limited  to the  amount of  One
+%           Dollar ($1.00).
+%
+%      3.2  CSIRO  make  no  warranties,  expressed  or  implied,  and
+%           excludes all  other warranties  representations, terms  or
+%           conditions, whether  express or implied,  oral or written,
+%           statutory  or  otherwise,  relating  in  any  way  to  the
+%           Software, or  to  this  Agreement, including  any  implied
+%           warranty of merchantability  or of fitness for  particular
+%           purpose.   To the full extent permitted  by the law of the
+%           Commonwealth of  Australia  or the  laws of  any State  or
+%           Territory  of  Australia,  any  conditions  or  warranties
+%           imposed by such  legislation are hereby  excluded.   In so
+%           far as  liability under  or pursuant  to such  legislation
+%           may not  be excluded,  CSIRO's liability  to the  Licencee
+%           pursuant to this Agreement shall be limited as  set out in
+%           clause 3.1 hereof.
+%
+%      3.3  The  Licencee acknowledges  and agrees  that  the Software
+%           was developed  for CSIRO  research purposes  and may  have
+%           inherent  defects, errors or deficiencies, and  that it is
+%           the  responsibility  of  the   Licencee  to  make  its  own
+%           assessment  of the  suitability  of the  Software  for the
+%           purpose  of  the   Licencee's  computing  activity.    The
+%           Licencee will use  the Software, and  advice, opinions  or
+%           information supplied by CSIRO,  its officers, employees or
+%           agents  concerning  the  Software  at the  Licencee's  own
+%           risk.
+%
+%      3.4  The  Licencee hereby  releases and  indemnifies and  shall
+%           continue to  release and  indemnify  CSIRO, its  officers,
+%           employees  and  agents  from  and  against  all   actions,
+%           claims, proceedings  or demands  (including those  brought
+%           by third parties) which may be bought against  it or them,
+%           whether  on their  own or  jointly  with the  Licencee and
+%           whether at common law, in equity or pursuant to statute or
+%           otherwise, in respect of  any loss, death, injury, illness
+%           or  damage  (whether  personal  or  property,  and whether
+%           direct   or    consequential,   including    consequential
+%           financial  loss)   and  any  infringement  of   copyright,
+%           patents,   trade  marks,  designs  or  other  Intellectual
+%           Property Rights, howsoever  arising out of the  Licencee's
+%           exercise of its  rights under this  Agreement and from and
+%           against  all  damages,  costs  and  expenses  incurred  in
+%           defending  or  settling  any  such  claim,  proceeding  or
+%           demand.
+%
+%      3.5  The  Licencee's  obligation to  indemnify  CSIRO  and  its
+%           officers,  employees  and  agents set  out  in  clause 3.4
+%           hereof   is   a   continuing   obligation   separate  from
+%           and independent of the Licencee's other obligations  under
+%           this  Agreement,  and  shall  survive  all  expiration  or
+%           termination of this Agreement.
+%
+%      4.0 Termination 
+%
+%      4.1  The Licence shall terminate  immediately upon the Licencee
+%           breaching any term or  condition of this Agreement whether
+%           or  not CSIRO is aware of  the occurrence of the breach at
+%           the time that it happens.
+%
+%      4.2  CSIRO  may terminate the Licence on  reasonable grounds by
+%           notice in  writing  to the  Licencee, and  such notice  of
+%           termination shall  be effective  immediately upon  receipt
+%           by  the  Licencee. 
+%
+%=========================================================================
+
+%$Id: sw_copy.M,v 1.6 1994/12/06 03:05:54 morgan Exp $
+
+more on
+help sw_copy
+more off
+return
+%--------------------------------------------------------------------
diff --git a/M_Csiro/sw_cp.m b/M_Csiro/sw_cp.m
new file mode 100644
index 0000000000000000000000000000000000000000..8094e061d38532980fbc489960fe3b89be6b72a7
--- /dev/null
+++ b/M_Csiro/sw_cp.m
@@ -0,0 +1,176 @@
+
+function cp = sw_cp(S,T,P)
+
+% SW_CP      Heat Capacity (Cp) of sea water
+%=========================================================================
+% SW_CP  $Revision: 1.3 $  $Date: 1994/10/10 04:38:05 $
+%         Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE: cp = sw_cp(S,T,P)
+%
+% DESCRIPTION:
+%    Heat Capacity of Sea Water using UNESCO 1983 polynomial.
+%
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78)]
+%   T = temperature [degree C (IPTS-68)]
+%   P = pressure    [db]
+%       (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+% OUTPUT:
+%   cp = Specific Heat Capacity  [J kg^-1 C^-1] 
+% 
+% AUTHOR:  Phil Morgan 93-04-20  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%=========================================================================
+
+% CALLER: general purpose
+% CALLEE: none
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+if nargin ~=3
+   error('sw_cp.m: Must pass 3 parameters')
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+%------
+% BEGIN
+%------
+P = P/10; % to convert db to Bar as used in Unesco routines
+
+%------------
+% eqn 26 p.32
+%------------
+c0 = 4217.4;
+c1 =   -3.720283;
+c2 =    0.1412855;
+c3 =   -2.654387e-3;
+c4 =    2.093236e-5;
+
+a0 = -7.64357;
+a1 =  0.1072763;
+a2 = -1.38385e-3;
+
+b0 =  0.1770383;
+b1 = -4.07718e-3;
+b2 =  5.148e-5;
+
+Cpst0 =  c0 + c1.*T + c2.*T.^2 + c3.*T.^3 + c4.*T.^4 + ...
+        (a0 + a1.*T + a2.*T.^2).*S + ...
+	(b0 + b1.*T + b2.*T.^2).*S.*sqrt(S);
+    
+%------------
+% eqn 28 p.33
+%------------
+a0 = -4.9592e-1;
+a1 =  1.45747e-2;
+a2 = -3.13885e-4;
+a3 =  2.0357e-6;
+a4 =  1.7168e-8;
+
+b0 =  2.4931e-4;
+b1 = -1.08645e-5;
+b2 =  2.87533e-7;
+b3 = -4.0027e-9;
+b4 =  2.2956e-11;
+
+c0 = -5.422e-8;
+c1 =  2.6380e-9;
+c2 = -6.5637e-11;
+c3 =  6.136e-13;
+
+del_Cp0t0 =  (a0 + a1.*T + a2.*T.^2 + a3.*T.^3 + a4.*T.^4).*P +    ...
+	     (b0 + b1.*T + b2.*T.^2 + b3.*T.^3 + b4.*T.^4).*P.^2 + ...   
+             (c0 + c1.*T + c2.*T.^2 + c3.*T.^3).*P.^3;
+
+%------------	 
+% eqn 29 p.34
+%------------
+d0 =  4.9247e-3;
+d1 = -1.28315e-4;
+d2 =  9.802e-7;
+d3 =  2.5941e-8;
+d4 = -2.9179e-10;
+
+e0 = -1.2331e-4;
+e1 = -1.517e-6;
+e2 =  3.122e-8;
+
+f0 = -2.9558e-6;
+f1 =  1.17054e-7;
+f2 = -2.3905e-9;
+f3 =  1.8448e-11;
+
+g0 =  9.971e-8;
+
+h0 =  5.540e-10;
+h1 = -1.7682e-11;
+h2 =  3.513e-13;
+
+j1 = -1.4300e-12;
+S3_2  = S.*sqrt(S);
+
+del_Cpstp = [(d0 + d1.*T + d2.*T.^2 + d3.*T.^3 + d4.*T.^4).*S + ...
+             (e0 + e1.*T + e2.*T.^2).*S3_2].*P                + ...
+	    [(f0 + f1.*T + f2.*T.^2 + f3.*T.^3).*S            + ...
+	     g0.*S3_2].*P.^2                                  + ...
+	     [(h0 + h1.*T + h2.*T.^2).*S                      + ...
+	     j1.*T.*S3_2].*P.^3;
+     
+
+cp = Cpst0 + del_Cp0t0 + del_Cpstp;
+
+if Transpose
+   cp = cp';
+end %if
+
+return
+%--------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_data.mat b/M_Csiro/sw_data.mat
new file mode 100644
index 0000000000000000000000000000000000000000..c12ce2d917c580ea3790aed4f2ae55421b45fc31
Binary files /dev/null and b/M_Csiro/sw_data.mat differ
diff --git a/M_Csiro/sw_dens.m b/M_Csiro/sw_dens.m
new file mode 100644
index 0000000000000000000000000000000000000000..a1613d52e7b0c60a537e76ecf252c7c753244cd8
--- /dev/null
+++ b/M_Csiro/sw_dens.m
@@ -0,0 +1,103 @@
+
+function dens = sw_dens(S,T,P)
+
+% SW_DENS    Density of sea water
+%=========================================================================
+% SW_DENS  $Revision: 1.3 $  $Date: 1994/10/10 04:39:16 $
+%          Copyright (C) CSIRO, Phil Morgan 1992.
+%
+% USAGE:  dens = sw_dens(S,T,P)
+%
+% DESCRIPTION:
+%    Density of Sea Water using UNESCO 1983 (EOS 80) polynomial.
+%
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78)]
+%   T = temperature [degree C (IPTS-68)]
+%   P = pressure    [db]
+%       (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+% OUTPUT:
+%   dens = density  [kg/m^3] 
+% 
+% AUTHOR:  Phil Morgan 92-11-05  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonoff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%
+%    Millero, F.J., Chen, C.T., Bradshaw, A., and Schleicher, K.
+%    " A new high pressure equation of state for seawater"
+%    Deap-Sea Research., 1980, Vol27A, pp255-264.
+%=========================================================================
+
+% CALLER: general purpose
+% CALLEE: sw_dens0.m sw_seck.m
+
+% UNESCO 1983. eqn.7  p.15
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+if nargin ~=3
+   error('sw_dens.m: Must pass 3 parameters')
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+%------
+% BEGIN
+%------
+densP0 = sw_dens0(S,T);
+K      = sw_seck(S,T,P);
+P      = P/10;  % convert from db to atm pressure units
+dens   = densP0./(1-P./K);
+
+if Transpose
+   dens = dens';
+end %if
+
+return
+%--------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_dens0.m b/M_Csiro/sw_dens0.m
new file mode 100644
index 0000000000000000000000000000000000000000..0d63ca8078d24d5f733e6089ce53d6b1dfaebff4
--- /dev/null
+++ b/M_Csiro/sw_dens0.m
@@ -0,0 +1,91 @@
+
+function dens = sw_dens0(S,T)
+
+% SW_DENS0   Denisty of sea water at atmospheric pressure
+%=========================================================================
+% SW_DENS0  $Revision: 1.3 $  $Date: 1994/10/10 04:54:09 $
+%           Copyright (C) CSIRO, Phil Morgan 1992
+%
+% USAGE:  dens0 = sw_dens0(S,T)
+%
+% DESCRIPTION:
+%    Density of Sea Water at atmospheric pressure using
+%    UNESCO 1983 (EOS 1980) polynomial.
+%
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78)]
+%   T = temperature [degree C (IPTS-68)]
+%
+% OUTPUT:
+%   dens0 = density  [kg/m^3] of salt water with properties S,T,
+%           P=0 (0 db gauge pressure)
+% 
+% AUTHOR:  Phil Morgan 92-11-05  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%     Unesco 1983. Algorithms for computation of fundamental properties of 
+%     seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%
+%     Millero, F.J. and  Poisson, A.
+%     International one-atmosphere equation of state of seawater.
+%     Deep-Sea Res. 1981. Vol28A(6) pp625-629.
+%=========================================================================
+
+% CALLER: general purpose, sw_dens.m
+% CALLEE: sw_smow.m
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+if nargin ~=2
+   error('sw_dens0.m: Must pass 2 parameters')
+end %if
+
+[mS,nS] = size(S);
+[mT,nT] = size(T);
+
+if (mS~=mT) | (nS~=nT)
+   error('sw_dens0.m: S,T inputs must have the same dimensions')
+end %if
+
+Transpose = 0;
+if mS == 1  % a row vector
+  S = S(:);
+  T = T(:);
+  Transpose = 1;
+end %if
+
+%----------------------
+% DEFINE CONSTANTS
+%----------------------
+%     UNESCO 1983 eqn(13) p17.
+
+b0 =  8.24493e-1;
+b1 = -4.0899e-3;
+b2 =  7.6438e-5;
+b3 = -8.2467e-7;
+b4 =  5.3875e-9;
+
+c0 = -5.72466e-3;
+c1 = +1.0227e-4;
+c2 = -1.6546e-6;
+
+d0 = 4.8314e-4;
+
+%$$$ dens = sw_smow(T) + (b0 + b1*T + b2*T.^2 + b3*T.^3 + b4*T.^4).*S  ...
+%$$$                    + (c0 + c1*T + c2*T.^2).*S.*sqrt(S) + d0*S.^2;
+
+dens = sw_smow(T) + (b0 + (b1 + (b2 + (b3 + b4*T).*T).*T).*T).*S  ...
+                   + (c0 + (c1 + c2*T).*T).*S.*sqrt(S) + d0*S.^2;	       
+
+if Transpose
+  dens = dens';
+end %if
+
+return
+%--------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_dist.m b/M_Csiro/sw_dist.m
new file mode 100644
index 0000000000000000000000000000000000000000..329ca194a5a2ae888ee03a301cbc036ec1009883
--- /dev/null
+++ b/M_Csiro/sw_dist.m
@@ -0,0 +1,112 @@
+
+function [dist,phaseangle] = distance(lat,lon,units)
+
+% SW_DIST    Distance between two lat,lon coordinates
+%===================================================================
+% SW_DIST  $Revision: 1.4 $  $Date: 1994/10/10 04:55:23 $
+%          Copyright (C) CSIRO, Phil Morgan & Steve Rintoul 1992. 
+%
+% USAGE:  [dist,phaseangle] = distance(lat,lon {,units} )
+%
+% DESCRIPTION:
+%   Calculate distance between two positions on glode using the "Plane
+%   Sailing" method.  Also uses simple geometry to calculate the bearing of
+%   the path between position pairs.
+% 
+% INPUT:
+%    lat      = decimal degrees (+ve N, -ve S) [- 90.. +90]
+%    lon      = decimal degrees (+ve E, -ve W) [-180..+180]
+%    units    = optional string specifing units of distance
+%               'nm'  = nautical miles (default)
+%               'km'  = kilometres
+%
+% OUTPUT:
+%    dist        = distance between positions in units
+%    phaseangle  = angle of line between stations with x axis (East).
+%                  Range of values are -180..+180. (E=0, N=90, S=-90)
+%
+% AUTHOR:   Phil Morgan and Steve Rintoul 92-02-10
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+% 
+% REFERENCE:
+%    The PLANE SAILING method as descriibed in "CELESTIAL NAVIGATION" 1989 by
+%    Dr. P. Gormley. The Australian Antartic Division.
+%==================================================================
+
+% CALLER:   general purpose
+% CALLEE:   none
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+if nargin > 3
+  error('sw_dist.m: No more than 3 arguments allowed')
+elseif nargin==3 
+  if ~isstr(units)
+      error('sw_dist.m: units argument must be string')
+  end %if
+elseif nargin==2
+  units = 'nm';  % default units
+else
+  error('sw_dist.m: wrong number of arguments')
+end%if
+
+[mlat,nlat] = size(lat);
+if mlat~=1 & nlat~=1
+   error('sw_dist.m: lat, lon must be vectors.  No matrices allowed')
+else
+  if mlat == 1
+    Transpose = 1;  % row vector passed in
+  else
+    Transpose = 0;  % accept column vector
+  end%if
+end%if
+lat = lat(:); %force to column vectors
+lon = lon(:);
+if length(lat)~=length(lon)
+   error('sw_dist.m: lat and lon must have same number of elements')
+end%if
+
+%-----------------
+% DEFINE CONSTANTS
+%-----------------
+DEG2RAD = (2*pi/360);
+RAD2DEG = 1/DEG2RAD;
+DEG2MIN = 60;
+DEG2NM  = 60;
+NM2KM   = 1.8520;    % Defined in Pond & Pickard p303.
+
+% BEGIN
+npositions = length(lat);
+ind=1:npositions-1;     % index to first of position pairs
+
+dlon = diff(lon);
+if any(abs(dlon)>180)
+   flag = find(abs(dlon)>180);
+   for ii=1:length(flag)
+     dlon(flag(ii))= -sign(dlon(flag(ii))) * (360 - abs(dlon(flag(ii))) );
+   end %for
+end %if
+latrad = abs(lat*DEG2RAD);
+dep    = cos( (latrad(ind+1)+latrad(ind))./2 ) .* dlon;
+dlat   = diff(lat);
+dist   = DEG2NM*sqrt(dlat.^2 + dep.^2);  % in n.miles
+
+if strcmp(units,'km')   % defaults to n.miles
+    dist = dist * NM2KM;
+end %if
+
+% CALCUALTE ANGLE TO X AXIS
+phaseangle  = angle(dep+dlat*sqrt(-1))*RAD2DEG;
+ 
+if Transpose
+  dist = dist';
+  phaseangle = phaseangle';
+end %if
+
+return
+%--------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_dpth.m b/M_Csiro/sw_dpth.m
new file mode 100644
index 0000000000000000000000000000000000000000..91f3ab295b910853462335172f2cc04ef835620d
--- /dev/null
+++ b/M_Csiro/sw_dpth.m
@@ -0,0 +1,87 @@
+
+function DEPTHM = sw_dpth(P,LAT)
+
+% SW_DPTH    Depth from pressure
+%===========================================================================
+% SW_DPTH   $Revision: 1.3 $  $Date: 1994/10/10 04:56:32 $
+%           Copyright (C) CSIRO, Phil Morgan 1992.
+%
+% USAGE:  dpth = sw_dpth(P,LAT)
+%
+% DESCRIPTION:
+%    Calculates depth in metres from pressure in dbars.
+%
+% INPUT:  (all must have same dimensions)
+%   P   = Pressure    [db]
+%   LAT = Latitude in decimal degress north [-90..+90]
+%         (lat may have dimensions 1x1 or 1xn where P(mxn).
+%
+% OUTPUT:
+%  dpth = depth [metres]
+%
+% AUTHOR:  Phil Morgan 92-04-06  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%=========================================================================
+
+% CALLER:  general purpose
+% CALLEE:  none
+
+%-------------
+% CHECK INPUTS
+%-------------
+[mP,nP] = size(P);
+[mL,nL] = size(LAT);
+if mL==1 & nL==1
+  LAT = LAT*ones(size(P));
+  [mL,nL] = size(LAT);
+end %if  
+
+if (mP~=mL) | (nP~=nL)              % P & LAT are not the same shape
+     if (nP==nL) & (mL==1)          % LAT for each column of P
+        LAT = LAT( ones(1,mP), : ); %     copy LATS down each column
+                                    %     s.t. dim(P)==dim(LAT)
+     else
+        error('sw_depth.m:  Inputs arguments have wrong dimensions')
+     end %if
+end %if
+
+Transpose = 0;
+if mP == 1  % row vector
+   P         =  P(:);
+   LAT       =  LAT(:);
+   Transpose = 1;
+end %if
+
+%-------------
+% BEGIN
+%-------------
+% Eqn 25, p26.  Unesco 1983.
+
+DEG2RAD = pi/180;
+c1 = +9.72659;
+c2 = -2.2512E-5;
+c3 = +2.279E-10;
+c4 = -1.82E-15;
+gam_dash = 2.184e-6;
+
+LAT = abs(LAT);
+X   = sin(LAT*DEG2RAD);  % convert to radians
+X   = X.*X;
+bot_line = 9.780318*(1.0+(5.2788E-3+2.36E-5*X).*X) + gam_dash*0.5*P;
+top_line = (((c4*P+c3).*P+c2).*P+c1).*P;
+DEPTHM   = top_line./bot_line;
+
+if Transpose
+   DEPTHM = DEPTHM';
+end %if
+
+return
+%===========================================================================
+%
diff --git a/M_Csiro/sw_f.m b/M_Csiro/sw_f.m
new file mode 100644
index 0000000000000000000000000000000000000000..caa36eb60201a23f9790cb70a23d50888d5d4aa7
--- /dev/null
+++ b/M_Csiro/sw_f.m
@@ -0,0 +1,57 @@
+
+function f = sw_f(lat)
+
+% SW_F       Coriolis factor "f"
+%===========================================================================
+% SW_F   $Revision: 1.3 $  $Date: 1994/10/10 04:57:08 $
+%        Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE:  f = sw_f(lat)
+%
+% DESCRIPTION:
+%    Calculates the Coriolis factor "f" defined by
+%       f = 2*Omega*Sin(lat)  where Omega = 7.292e-5 radians/sec
+%
+% INPUT:  
+%   lat = Latitude in decimal degress north [-90..+90]
+%
+% OUTPUT:
+%  f    = Coriolis Factor "f" [s-1]
+%
+% AUTHOR:  Phil Morgan 93-04-20  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCE: 
+%   S. Pond & G.Pickard  2nd Edition 1986
+%   Introductory Dynamical Oceanogrpahy
+%   Pergamon Press Sydney.  ISBN 0-08-028728-X
+%   
+%   A.E. Gill 1982. p.597
+%   "Atmosphere-Ocean Dynamics"
+%   Academic Press: New York.  ISBN: 0-12-283522-0
+%=========================================================================
+
+% CALLER:  general purpose
+% CALLEE:  none
+
+%-------------
+% CHECK INPUTS
+%-------------
+if nargin ~= 1
+   error('sw_f.m:  Requires one input argument')
+end %if  
+
+%-------------
+% BEGIN
+%-------------
+% Eqn p27.  Unesco 1983.
+DEG2RAD = pi/180;
+OMEGA   = 7.292e-5;     %s-1   A.E.Gill p.597
+f       = 2*OMEGA*sin(lat*DEG2RAD);
+
+return
+%===========================================================================
+
diff --git a/M_Csiro/sw_fp.m b/M_Csiro/sw_fp.m
new file mode 100644
index 0000000000000000000000000000000000000000..acd959c0d2f234100f54b6f550dea7f0fe1731f7
--- /dev/null
+++ b/M_Csiro/sw_fp.m
@@ -0,0 +1,90 @@
+
+function fp = sw_fp(S,P)
+
+% SW_FP      Freezing point of sea water
+%=========================================================================
+% SW_FP % $Revision: 1.3 $  $Date: 1994/10/10 04:57:50 $
+%         Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE:  fp = sw_fp(S,P)
+%
+% DESCRIPTION:
+%    Heat Capacity of Sea Water using UNESCO 1983 polynomial.
+%
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78)]
+%   P = pressure    [db]
+%       (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+% OUTPUT:
+%   fp = Freezing Point temperature [degree C (IPTS-68)]
+% 
+% AUTHOR:  Phil Morgan 93-04-20  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%=========================================================================
+
+% CALLER: general purpose
+% CALLEE: none
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+if nargin ~=2
+   error('sw_fp.m: Must pass 3 parameters')
+end %if
+
+[ms,ns] = size(S);
+[mp,np] = size(P);
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+  error('sw_fp.m: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   S       =  S(:);   
+   Transpose = 1;
+end %if
+
+%------
+% BEGIN
+%------
+%P = P/10; % to convert db to Bar as used in Unesco routines
+
+%------------
+% eqn  p.29
+%------------
+a0 = -0.0575;
+a1 = 1.710523e-3;
+a2 = -2.154996e-4;
+b  = -7.53e-4;
+
+fp = a0.*S + a1.*S.*sqrt(S) + a2.*S.^2 + b.*P;
+
+if Transpose
+   fp = fp';
+end %if
+
+return
+%--------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_g.m b/M_Csiro/sw_g.m
new file mode 100644
index 0000000000000000000000000000000000000000..c2a3c90218887dfd4466e6734e4f7b3e9490df22
--- /dev/null
+++ b/M_Csiro/sw_g.m
@@ -0,0 +1,70 @@
+
+function g = sw_g(LAT,z)
+
+% SW_G       Gravitational acceleration
+%===========================================================================
+% SW_G   $Revision: 1.4 $  $Date: 1994/10/11 00:00:54 $
+%        Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE:  g = sw_g(lat,z)
+%
+% DESCRIPTION:
+%    Calculates acceleration due to gravity as function of latitude.
+%
+% INPUT:  (all must have same dimensions)
+%   lat = Latitude in decimal degress north [-90..+90]
+%   z   = height in metres (+ve above sea surface, -ve below)
+%
+% OUTPUT:
+%  g    = gravity [m/s^2]
+%
+% AUTHOR:  Phil Morgan 93-04-20  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%   Unesco 1983. Algorithms for computation of fundamental properties of 
+%   seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%
+%   A.E. Gill 1982. p.597
+%   "Atmosphere-Ocean Dynamics"
+%   Academic Press: New York.  ISBN: 0-12-283522-0
+%=========================================================================
+
+% CALLER:  general purpose
+% CALLEE:  none
+
+%-------------
+% CHECK INPUTS
+%-------------
+if ~(nargin==1 | nargin==2)
+   error('sw_g.m:  Requires one or two input arguments')
+end %if  
+if nargin == 1
+  z = zeros(size(LAT));
+end %if
+
+[mL,nL] = size(LAT);
+[mz,nz] = size(z);
+if ~(mL==mz | nL==nz)
+   error('sw_g.m:  Input arguments should have same dimensions')
+end %if
+
+%-------------
+% BEGIN
+%-------------
+% Eqn p27.  Unesco 1983.
+a       = 6371000;    % mean radius of earth  A.E.Gill
+DEG2RAD = pi/180;
+LAT     = abs(LAT);
+X       = sin(LAT*DEG2RAD);  % convert to radians
+sin2    = X.*X;
+g       = 9.780318*(1.0+(5.2788E-3+2.36E-5*sin2).*sin2);
+if any(any(z))
+   g    = g./((1+z/a).^2);    % from A.E.Gill p.597
+end %if
+return
+%===========================================================================
+
diff --git a/M_Csiro/sw_gpan.m b/M_Csiro/sw_gpan.m
new file mode 100644
index 0000000000000000000000000000000000000000..2f5c1ad4b1174e69b9fdf4a0d34eb66d9a2cf86d
--- /dev/null
+++ b/M_Csiro/sw_gpan.m
@@ -0,0 +1,126 @@
+
+function [ga, pe] = sw_gpan(S,T,P,LAT)
+
+% SW_GPAN    Geopotential anomaly
+%=========================================================================
+% SW_GPAN  $Revision: 1.3 $  $Date: 1994/10/10 05:01:00 $
+%          Copyright (C) CSIRO, Phil Morgan 1992.
+%
+% USAGE:  [gpan, pe]= sw_gpan(S,T,P,LAT)
+%
+% DESCRIPTION:
+%   Geopotential Anomaly calculated as the integral of svan from the
+%   the sea surface to the bottom.  Thus RELATIVE TO SEA SURFACE.
+%
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78)]
+%   T = temperature [degree C (ITP-68)]
+%   P = Pressure    [db]
+%   LAT = latitude
+%       (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+% OUTPUT:
+%  gpan = Geopotential Anomaly  [m^3 kg^-1 Pa == m^2 s^-2 == J kg^-1]
+%  pe   = anomalie d'energie potentielle en kg/s^2 * 10
+%
+% AUTHOR:  Phil Morgan 92-11-05  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCE: S. Pond & G.Pickard  2nd Edition 1986
+%            Introductory Dynamical Oceanogrpahy
+%            Pergamon Press Sydney.  ISBN 0-08-028728-X
+%
+% Note that older literature may use units of "dynamic decimeter' for above.
+%
+% Adapted method from Pond and Pickard (p76) to calc gpan rel to sea 
+% surface whereas P&P calculated relative to the deepest common depth.
+%=========================================================================
+
+%
+% CALLER: general purpose
+% CALLEE: sw_svan.m 
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+if nargin ~= 4
+   error('sw_gpan.m: Must pass 4 parameters')
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+%------
+% BEGIN
+%------
+db2Pascal  = 1e4;
+[m,n]      = size(P);
+svan       = sw_svan(S,T,P);
+mean_svan  = 0.5*(svan(2:m,:) + svan(1:m-1,:) );
+
+if n==1
+   top = svan(1,1).*P(1,1)*db2Pascal;
+else
+   top = svan(1,:).*P(1,:)*db2Pascal;
+end %if
+
+%press_diff = diff(P);
+
+delta_ga   = (mean_svan.*diff(P))*db2Pascal;
+ga         = cumsum([ top; delta_ga]);
+
+% Rajoute par Y. Gouriou en se basant sur la routine de Millard
+%  pe[ref] = hdy[ref] * (ctd[ref].pres + ctd[ind].pres)*0.5F / 
+%		          (float)grav( (double)ctd[ref].pres,lat ); */
+%  pe[ref] =  (float) ( fabs( (double)( pres[ref] - pres[ind]) )
+%                     * (anvs[ref]*pres[ref] + anvs[ind]*pres[ind])
+%                     * 0.5e-5 / (float)grav( (double)pres[ref],lat ));
+mean_P     = 0.5*( P(2:m,:) + P(1:m-1,:) );
+delta_pe   = delta_ga .* mean_P ./  sw_g( LAT, -sw_dpth(P(2:end),LAT));
+pe         = cumsum([ top/sw_g( LAT, 0); delta_pe]); 
+
+if Transpose
+   ga = ga';
+   pe = pe';
+end %if
+
+return
+%--------------------------------------------------------------------
diff --git a/M_Csiro/sw_gvel.m b/M_Csiro/sw_gvel.m
new file mode 100644
index 0000000000000000000000000000000000000000..31ae07df662d3cc57db51f9de82daa8ea08ead0a
--- /dev/null
+++ b/M_Csiro/sw_gvel.m
@@ -0,0 +1,62 @@
+
+function vel = sw_gvel(ga,lat,lon)
+
+% SW_GVEL    Geostrophic velocity
+%===================================================================
+% GEOVEL   $Revision: 1.5 $  $Date: 1994/11/15 04:00:36 $
+%          Copyright (C) CSIRO, Phil Morgan 1992
+%
+% USAGE:  vel = sw_gvel(ga,lat,lon)
+%
+% DESCRIPTION:
+%    Calculates geostrophic velocity given the geopotential anomaly
+%    and position of each station.
+% 
+% INPUT:
+%    ga   = geopotential anomoly relative to the sea surface.
+%           dim(mxnstations)
+%    lat  = latitude  of each station (+ve = N, -ve = S) [ -90.. +90]
+%    lon  = longitude of each station (+ve = E, -ve = W) [-180..+180]
+%
+% OUTPUT:
+%    vel  = geostrophic velocity RELATIVE to the sea surface.
+%           dim(m,nstations-1)
+%
+% AUTHOR:   Phil Morgan   1992/03/26  (morgan@ml.csiro.au)
+%           
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCE: S. Pond & G.Pickard  2nd Edition 1986
+%            Introductory Dynamical Oceanogrpahy
+%            Pergamon Press Sydney.  ISBN 0-08-028728-X
+%            Equation 8.9A p73  Pond & Pickard
+%
+% NOTE: This calls sw_dist.m.  You can replace the call to this
+%       routine if you have a more appropraite distance routine.
+%==================================================================
+
+% CALLER:   general purpose
+% CALLEE:   sw_dist.m
+%
+
+  
+DEG2RAD = pi/180;
+RAD2DEG = 180/pi;
+OMEGA   = 7.292e-5;  % Angular velocity of Earth  [radians/sec]
+
+% You may replace the call to sw_dist if you have
+% a more appropriate distance routine.
+distm = 1000*sw_dist(lat,lon,'km');
+[m,n] = size(ga);
+f     = 2*OMEGA*sin( (lat(1:n-1)+lat(2:n))*DEG2RAD/2 );
+lf    = f.*distm;
+
+LF = lf(ones(m,1),:);
+
+vel   = -( ga(:,2:n)-ga(:,1:n-1) ) ./ LF;  
+
+return
+%--------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_info.m b/M_Csiro/sw_info.m
new file mode 100644
index 0000000000000000000000000000000000000000..21e57fcccbf24ff0049e28dae6b063c4deb6fd00
--- /dev/null
+++ b/M_Csiro/sw_info.m
@@ -0,0 +1,30 @@
+
+% SW_INFO    Computational routines for the properties of sea water
+%
+% SEAWATER - devloped by Phil Morgan, CSIRO
+%       
+% DESCRIPTION:
+%    SEAWATER is a toolkit of MATLAB routines for calculating the
+%    properties of sea water. They are a self contained library and
+%    are extremely easy to use and will run on all computers that 
+%    support MATLAB.  
+%
+% MATLAB:
+%    For information on MATLAB contact info@mathworks.com
+%       
+% DISCLAIMER
+%   This software is provided "as is" without warranty of any kind.
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% MORE INFORMATION:
+%   http://www.marine.csiro.au/~morgan/seawater
+%
+% $Revision: 1.10 $  $Date: 1998/04/21 05:50:33 $
+% Copyright (C) CSIRO, Phil Morgan 1993.
+%=========================================================================
+
+more on
+help sw_info
+more off
+return
+%--------------------------------------------------------------------
diff --git a/M_Csiro/sw_new.m b/M_Csiro/sw_new.m
new file mode 100644
index 0000000000000000000000000000000000000000..0034026316f4a07874dacc8c85922cec66809177
--- /dev/null
+++ b/M_Csiro/sw_new.m
@@ -0,0 +1,64 @@
+
+% SW_NEW    What's new in this version of seawater.
+%
+% 22 April 1998  release 2.0.1 (For version 5.x of Matlab)
+% ********************************************************
+% This version is not optimised but will run under Matlab 5.x
+% sw_satAr    New routine. Solubility of Ar in seawater
+% sw_satN2    New routine. Solubility of Ar in seawater
+% sw_satO2    New routine. Solubility of Ar in seawater
+% sw_test     Updated to include tests for above
+% 
+% April 1998  release 1.2e (For version 4.x of Matlab)
+% ************************
+% sw_alpha    Fixed bug where temp used in calculations regardless of
+%             whether 'temp' or 'pmpt' was passed as keyword.
+%
+% sw_info     Shorter version. Refer users to web pages
+%             http://www.marine.csiro.au
+%
+% sw_ver      New routine. Returns version number of SEAWATER
+%
+% sw_test     New Routine. Run a test on the SEAWATER routines
+%             and compare results with literature values
+%
+% 94/11/15 release 1.2d
+% **********************
+% sw_bfrq.m   Now also returns potential vorticity.
+%             Thanks to Greg Johnson (gjohnson@pmel.noaa.gov)
+%
+% sw_gvel.m   OMEGA=7.29e-5 changed to OMEGA=7.292e-5 to be
+%             consistent with sw_f.m
+%
+%             IMPORTANT CHANGE: The usage of the following 
+%             routines has changed!
+%
+% sw_alpha.m |    All these routines expect (S,T,P) to
+% sw_beta.m  |--  be passed instead of (S,PTMP,P) as in 
+% sw_aonb.m  |    previous releases of seawater.
+%                 Fast execution can still be obtained by passing
+%                 ptmp with a string flag 'ptmp' see help.
+%
+% 94/10/19 release 1.2c
+% **********************
+% Added routine sw_new.m to inform of updates and new features.
+% sw_bfrq.m   Fixed bug where LAT = [] was needed as argument when
+%             no latitude values are being passed.
+%             Now pass PRESSURE instead of DEPTH -> more consistent
+%             though only a negligible change is answers.
+%
+% sw_info.m   Updated to include a registration section.
+%             Noted that software is FREE.  
+%             Noted best email address is seawater@ml.csiro.au
+%             Requests for Report also via email to library@ml.csiro.au
+%
+% 94/10/12 release 1.2b
+% ********************
+% First official release and announcement on the networks.
+%
+
+more on
+help sw_new
+more off
+
+%-------------
diff --git a/M_Csiro/sw_pden.m b/M_Csiro/sw_pden.m
new file mode 100644
index 0000000000000000000000000000000000000000..9f47a8bbd2bd277dafa5dacbc0a169d96058e414
--- /dev/null
+++ b/M_Csiro/sw_pden.m
@@ -0,0 +1,57 @@
+
+function pden = sw_pden(S,T,P,PR)
+
+% SW_PDEN    Potential density
+%===========================================================================
+% SW_PDEN  $Revision: 1.3 $  $Date: 1994/10/10 05:05:21 $
+%          Copyright (C) CSIRO, Phil Morgan  1992. 
+%
+% USAGE:  pden = sw_pden(S,T,P,PR) 
+%
+% DESCRIPTION:
+%    Calculates potential density of water mass relative to the specified
+%    reference pressure by pden = sw_dens(S,ptmp,PR).
+%   
+% INPUT:  (all must have same dimensions)
+%   S  = salinity    [psu      (PSS-78) ]
+%   T  = temperature [degree C (IPTS-68)]
+%   P  = pressure    [db]
+%   PR = Reference pressure  [db]
+%       (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+% OUTPUT:
+%   pden = Potential denisty relative to the ref. pressure [kg/m^3] 
+%
+% AUTHOR:  Phil Morgan 1992/04/06  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%   A.E. Gill 1982. p.54
+%   "Atmosphere-Ocean Dynamics"
+%   Academic Press: New York.  ISBN: 0-12-283522-0
+%=========================================================================
+
+% CALLER:  general purpose
+% CALLEE:  sw_ptmp.m sw_dens.m
+
+%-------------
+% CHECK INPUTS
+%-------------
+if nargin ~= 4
+   error('sw_pden.m: Must pass 4 parameters ')
+end %if
+
+% LET sw_ptmp.m DO DIMENSION CHECKING
+
+%------
+% BEGIN
+%------
+ptmp = sw_ptmp(S,T,P,PR);
+pden = sw_dens(S,ptmp,PR);
+
+return      
+%=========================================================================
+
diff --git a/M_Csiro/sw_pres.m b/M_Csiro/sw_pres.m
new file mode 100644
index 0000000000000000000000000000000000000000..1d2d7ac7a396db622eac9f9a29e0a65a414384e5
--- /dev/null
+++ b/M_Csiro/sw_pres.m
@@ -0,0 +1,81 @@
+
+function pres = sw_pres(DEPTH,LAT)
+
+% SW_PRES    Pressure from depth
+%===========================================================================
+% SW_PRES   $Revision: 1.5 $  $Date: 1994/10/11 01:23:32 $
+%           Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE:  pres = sw_pres(depth,lat)
+%
+% DESCRIPTION:
+%    Calculates pressure in dbars from depth in meters.
+%
+% INPUT:  (all must have same dimensions)
+%   depth = depth [metres]  
+%   lat   = Latitude in decimal degress north [-90..+90]
+%           (LAT may have dimensions 1x1 or 1xn where depth(mxn) )
+%
+% OUTPUT:
+%  pres   = Pressure    [db]
+%
+% AUTHOR:  Phil Morgan 93-06-25  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Saunders, P.M. 1981
+%    "Practical conversion of Pressure to Depth"
+%    Journal of Physical Oceanography, 11, 573-574
+%
+% CHECK VALUE: 
+%    P=7500.00 db for LAT=30 deg, depth=7321.45 meters
+%=========================================================================
+
+% CALLER:  general purpose
+% CALLEE:  none
+
+%-------------
+% CHECK INPUTS
+%-------------
+
+[mD,nD] = size(DEPTH);
+[mL,nL] = size(LAT);
+if mL==1 & nL==1
+  LAT = LAT*ones(size(DEPTH));
+  [mL,nL] = size(LAT);
+end %if  
+
+if (mD~=mL) | (nD~=nL)              % DEPTH & LAT are not the same shape
+     if (nD==nL) & (mL==1)          % LAT for each column of DEPTH
+        LAT = LAT( ones(1,mD), : ); %     copy LATS down each column
+                                    %     s.t. dim(DEPTH)==dim(LAT)
+     else
+        error('sw_pres.m:  Inputs arguments have wrong dimensions')
+     end %if
+end %if
+
+Transpose = 0;
+if mD == 1  % row vector
+   DEPTH   =  DEPTH(:);
+   LAT     =  LAT(:);
+   Transpose = 1;
+end %if
+
+%-------------
+% BEGIN
+%-------------
+
+DEG2RAD = pi/180;
+X       = sin(abs(LAT)*DEG2RAD);  % convert to radians
+C1      = 5.92E-3+X.^2*5.25E-3;
+pres    = ((1-C1)-sqrt(((1-C1).^2)-(8.84E-6*DEPTH)))/4.42E-6;
+
+if Transpose
+   pres = pres';
+end %if
+
+return
+%===========================================================================
diff --git a/M_Csiro/sw_ptmp.m b/M_Csiro/sw_ptmp.m
new file mode 100644
index 0000000000000000000000000000000000000000..b0c37ebab7c7ab702f84b3c44a3e8ac53fffe817
--- /dev/null
+++ b/M_Csiro/sw_ptmp.m
@@ -0,0 +1,137 @@
+
+function PT = sw_ptmp(S,T,P,PR)
+
+% SW_PTMP    Potential temperature
+%===========================================================================
+% SW_PTMP  $Revision: 1.3 $  $Date: 1994/10/10 05:45:13 $
+%          Copyright (C) CSIRO, Phil Morgan 1992. 
+%
+% USAGE:  ptmp = sw_ptmp(S,T,P,PR) 
+%
+% DESCRIPTION:
+%    Calculates potential temperature as per UNESCO 1983 report.
+%   
+% INPUT:  (all must have same dimensions)
+%   S  = salinity    [psu      (PSS-78) ]
+%   T  = temperature [degree C (IPTS-68)]
+%   P  = pressure    [db]
+%   PR = Reference pressure  [db]
+%        (P & PR may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+% OUTPUT:
+%   ptmp = Potential temperature relative to PR [degree C (IPTS-68)]
+%
+% AUTHOR:  Phil Morgan 92-04-06  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonoff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%    Eqn.(31) p.39
+%
+%    Bryden, H. 1973.
+%    "New Polynomials for thermal expansion, adiabatic temperature gradient
+%    and potential temperature of sea water."
+%    DEEP-SEA RES., 1973, Vol20,401-408.
+%=========================================================================
+
+% CALLER:  general purpose
+% CALLEE:  sw_adtg.m
+
+%-------------
+% CHECK INPUTS
+%-------------
+if nargin ~= 4
+   error('sw_ptmp.m: Must pass 4 parameters ')
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+[mpr,npr] = size(PR);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+% CHECK OPTIONAL SHAPES FOR PR
+if     mpr==1  & npr==1      % PR is a scalar.  Fill to size of S
+   PR = PR(1)*ones(ms,ns);
+elseif npr==ns & mpr==1      % PR is row vector with same cols as S
+   PR = PR( ones(1,ms), : ); %   Copy down each column.
+elseif mpr==ms & npr==1      % P is column vector
+   PR = PR( :, ones(1,ns) ); %   Copy across each row
+elseif mpr==ms & npr==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+  error('check_stp: PR has wrong dimensions')
+end %if
+[mpr,npr] = size(PR);
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   PR      = PR(:);
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+%------
+% BEGIN
+%------
+
+% theta1
+del_P  = PR - P;
+del_th = del_P.*sw_adtg(S,T,P);
+th     = T + 0.5*del_th;
+q      = del_th;
+
+% theta2
+del_th = del_P.*sw_adtg(S,th,P+0.5*del_P);
+th     = th + (1 - 1/sqrt(2))*(del_th - q);
+q      = (2-sqrt(2))*del_th + (-2+3/sqrt(2))*q;
+
+% theta3
+del_th = del_P.*sw_adtg(S,th,P+0.5*del_P);
+th     = th + (1 + 1/sqrt(2))*(del_th - q);
+q      = (2 + sqrt(2))*del_th + (-2-3/sqrt(2))*q;
+
+% theta4
+del_th = del_P.*sw_adtg(S,th,P+del_P);
+PT     = th + (del_th - 2*q)/6;
+
+if Transpose
+  PT = PT';
+end %if
+
+return      
+%=========================================================================
diff --git a/M_Csiro/sw_salds.m b/M_Csiro/sw_salds.m
new file mode 100644
index 0000000000000000000000000000000000000000..134e970134d9b26b3b95bad079e84b383074f4a2
--- /dev/null
+++ b/M_Csiro/sw_salds.m
@@ -0,0 +1,74 @@
+
+function dS = sw_salds(Rtx,delT)
+
+% SW_SALDS   Differiential dS/d(sqrt(Rt)) at constant T.
+%=========================================================================
+% SW_SALDS   $Revision: 1.3 $  $Date: 1994/10/10 05:46:08 $
+%            Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE:  dS = sw_salds(Rtx,delT)
+%
+% DESCRIPTION:
+%   Calculates Salinity differential dS/d(sqrt(Rt)) at constant T.
+%   UNESCO 1983 polynomial.
+%
+% INPUT: (all must have same dimensions)
+%   Rtx   = sqrt(Rt) where Rt defined in sw_salt.m
+%   delT  = T-15     [degree C (IPTS-68)]
+%
+% OUTPUT:
+%   dS = S differential dS/d(sqrt(Rt)) at constant T.
+% 
+% AUTHOR:  Phil Morgan 93-04-21  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonoff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%=========================================================================
+
+% CALLER: sw_cndr.m
+% CALLEE: none
+
+%-------------
+% CHECK INPUTS
+%-------------
+if nargin~=2
+   error('sw_salds.m: must have 2 input arguments')
+end %if
+
+[m1,n1] = size(Rtx);
+[m2,n2] = size(delT);
+if ~(m1==m2 | n1==n2)
+  error('sw_salds.m: Rtx and delT must have the same shape')
+end %if
+
+%-------
+% BEGIN
+%-------
+a0 =  0.0080;
+a1 = -0.1692;
+a2 = 25.3851;
+a3 = 14.0941;
+a4 = -7.0261;
+a5 =  2.7081;
+
+b0 =  0.0005;
+b1 = -0.0056;
+b2 = -0.0066;
+b3 = -0.0375;
+b4 =  0.0636;
+b5 = -0.0144;
+
+k  =  0.0162;
+
+dS =  a1 + (2*a2 + (3*a3 + (4*a4 + 5*a5.*Rtx).*Rtx).*Rtx).*Rtx + ...
+     (delT./(1+k*delT))* ...
+     (b1 + (2*b2 + (3*b3 + (4*b4 + 5*b5.*Rtx).*Rtx).*Rtx).*Rtx);
+
+return
+%-----------------------------------------------------------------------
diff --git a/M_Csiro/sw_salrp.m b/M_Csiro/sw_salrp.m
new file mode 100644
index 0000000000000000000000000000000000000000..c795f22290dec2653b6cba528415ef3fb8c1dfd5
--- /dev/null
+++ b/M_Csiro/sw_salrp.m
@@ -0,0 +1,69 @@
+
+function Rp = sw_salrp(R,T,P)
+
+% SW_SALRP   Conductivity ratio   Rp(S,T,P) = C(S,T,P)/C(S,T,0)
+%=========================================================================
+% SW_SALRP   $Revision: 1.3 $  $Date: 1994/10/10 05:47:27 $
+%            Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE:  Rp = sw_salrp(R,T,P)
+%
+% DESCRIPTION:
+%    Equation Rp(S,T,P) = C(S,T,P)/C(S,T,0) used in calculating salinity.
+%    UNESCO 1983 polynomial.
+%
+% INPUT: (All must have same shape)
+%   R = Conductivity ratio  R =  C(S,T,P)/C(35,15,0) [no units]
+%   T = temperature [degree C (IPTS-68)]
+%   P = pressure    [db]
+%
+% OUTPUT:
+%   Rp = conductivity ratio  Rp(S,T,P) = C(S,T,P)/C(S,T,0)  [no units] 
+% 
+% AUTHOR:  Phil Morgan 93-04-17  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonoff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%=========================================================================
+
+% CALLER: sw_salt
+% CALLEE: none
+
+%-------------------
+% CHECK INPUTS
+%-------------------
+if nargin~=3
+  error('sw_salrp.m: requires 3 input arguments')
+end %if
+
+[mr,nr] = size(R);
+[mp,np] = size(P);
+[mt,nt] = size(T);
+if ~(mr==mp | mr==mt | nr==np | nr==nt)
+   error('sw_salrp.m: R,T,P must all have the same shape')
+end %if   
+
+%-------------------
+% eqn (4) p.8 unesco.
+%-------------------
+d1 =  3.426e-2;
+d2 =  4.464e-4;
+d3 =  4.215e-1;
+d4 = -3.107e-3;
+
+e1 =  2.070e-5;
+e2 = -6.370e-10;
+e3 =  3.989e-15;
+
+Rp = 1 + ( P.*(e1 + e2.*P + e3.*P.^2) ) ...
+     ./ (1 + d1.*T + d2.*T.^2 +(d3 + d4.*T).*R);
+ 
+return
+%-----------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_salrt.m b/M_Csiro/sw_salrt.m
new file mode 100644
index 0000000000000000000000000000000000000000..db985e594d3d1543519f9a27c1d843f5237fb0bc
--- /dev/null
+++ b/M_Csiro/sw_salrt.m
@@ -0,0 +1,49 @@
+
+function rt = sw_salrt(T)
+
+% SW_SALRT   Conductivity ratio   rt(T)     = C(35,T,0)/C(35,15,0)
+%=========================================================================
+% SW_SALRT  $Revision: 1.3 $  $Date: 1994/10/10 05:48:34 $
+%           Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE:  rt = sw_salrt(T)
+%
+% DESCRIPTION:
+%    Equation rt(T) = C(35,T,0)/C(35,15,0) used in calculating salinity.
+%    UNESCO 1983 polynomial.
+%
+% INPUT: 
+%   T = temperature [degree C (IPTS-68)]
+%
+% OUTPUT:
+%   rt = conductivity ratio  [no units] 
+% 
+% AUTHOR:  Phil Morgan 93-04-17  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonoff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%=========================================================================
+
+% CALLER: sw_salt
+% CALLEE: none
+
+% rt = rt(T) = C(35,T,0)/C(35,15,0)
+% Eqn (3) p.7 Unesco.
+
+c0 =  0.6766097;
+c1 =  2.00564e-2;
+c2 =  1.104259e-4;
+c3 = -6.9698e-7;
+% c4 =  1.0031e-9;
+c4 =  1.e-9;
+
+rt = c0 + (c1 + (c2 + (c3 + c4.*T).*T).*T).*T;
+
+return
+%--------------------------------------------------------------------
diff --git a/M_Csiro/sw_sals.m b/M_Csiro/sw_sals.m
new file mode 100644
index 0000000000000000000000000000000000000000..fb7bdbd1f3cd822298eea1eeda0faa0cdd0b9c66
--- /dev/null
+++ b/M_Csiro/sw_sals.m
@@ -0,0 +1,79 @@
+
+function S = sw_sals(Rt,T)
+
+% SW_SALS    Salinity of sea water
+%=========================================================================
+% SW_SALS  $Revision: 1.3 $  $Date: 1994/10/10 05:49:13 $
+%          Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE:  S = sw_sals(Rt,T)
+%
+% DESCRIPTION:
+%    Salinity of sea water as a function of Rt and T.  
+%    UNESCO 1983 polynomial.
+%
+% INPUT:
+%   Rt = Rt(S,T) = C(S,T,0)/C(35,T,0)
+%   T  = temperature [degree C (IPTS-68)]
+%
+% OUTPUT:
+%   S  = salinity    [psu      (PSS-78)]
+% 
+% AUTHOR:  Phil Morgan 93-04-17  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonoff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%=========================================================================
+
+% CALLER: sw_salt
+% CALLEE: none
+
+%--------------------------
+% CHECK INPUTS
+%--------------------------
+if nargin~=2
+  error('sw_sals.m: requires 2 input arguments')
+end %if
+
+[mrt,nrt] = size(Rt);
+[mT,nT]   = size(T);
+if ~(mrt==mT | nrt==nT)
+   error('sw_sals.m: Rt and T must have the same shape')
+end %if
+
+%--------------------------
+% eqn (1) & (2) p6,7 unesco
+%--------------------------
+a0 =  0.0080;
+a1 = -0.1692;
+a2 = 25.3851;
+a3 = 14.0941;
+a4 = -7.0261;
+a5 =  2.7081;
+
+b0 =  0.0005;
+b1 = -0.0056;
+b2 = -0.0066;
+b3 = -0.0375;
+b4 =  0.0636;
+b5 = -0.0144;
+
+k  =  0.0162;
+
+Rtx   = sqrt(Rt);
+del_T = T - 15;
+del_S = (del_T ./ (1+k*del_T) ) .* ...
+        ( b0 + (b1 + (b2+ (b3 + (b4 + b5.*Rtx).*Rtx).*Rtx).*Rtx).*Rtx);
+	
+S = a0 + (a1 + (a2 + (a3 + (a4 + a5.*Rtx).*Rtx).*Rtx).*Rtx).*Rtx;
+
+S = S + del_S;
+
+return
+%----------------------------------------------------------------------
diff --git a/M_Csiro/sw_salt.m b/M_Csiro/sw_salt.m
new file mode 100644
index 0000000000000000000000000000000000000000..d0102b828cdad92201111d489763c52029cb45a4
--- /dev/null
+++ b/M_Csiro/sw_salt.m
@@ -0,0 +1,60 @@
+
+function S = sw_salt(cndr,T,P)
+
+% SW_SALT    Salinity from cndr, T, P
+%=========================================================================
+% SW_SALT  $Revision: 1.3 $  $Date: 1994/10/10 05:49:53 $
+%          Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE: S = sw_salt(cndr,T,P)
+%
+% DESCRIPTION:
+%   Calculates Salinity from conductivity ratio. UNESCO 1983 polynomial.
+%
+% INPUT:
+%   cndr = Conductivity ratio     R =  C(S,T,P)/C(35,15,0) [no units]
+%   T    = temperature [degree C (IPTS-68)]
+%   P    = pressure    [db]
+%
+% OUTPUT:
+%   S    = salinity    [psu      (PSS-78)]
+% 
+% AUTHOR:  Phil Morgan 93-04-17  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonoff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%=========================================================================
+
+% CALLER: general purpose
+% CALLEE: sw_sals.m sw_salrt.m sw_salrp.m
+
+  
+%----------------------------------
+% CHECK INPUTS ARE SAME DIMENSIONS
+%----------------------------------
+[mc,nc] = size(cndr);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+if ~(mc==mt | mc==mp | nc==nt | nc==np)
+  error('sw_salt.m: cndr,T,P must all have the same dimensions')
+end %if
+
+%-------
+% BEGIN
+%-------
+R  = cndr;
+rt = sw_salrt(T);
+Rp = sw_salrp(R,T,P);
+Rt = R./(Rp.*rt);
+S  = sw_sals(Rt,T);
+
+return
+%--------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_satAr.m b/M_Csiro/sw_satAr.m
new file mode 100644
index 0000000000000000000000000000000000000000..481271b49716147440fe16372fcbef2bc9515187
--- /dev/null
+++ b/M_Csiro/sw_satAr.m
@@ -0,0 +1,95 @@
+%$$$ 
+%$$$ #undef __PR
+%$$$ #include "VARIANT.h"
+
+function c = sw_satAr(S,T)
+
+% SW_SATAr   Satuaration of Ar in sea water
+%=========================================================================
+% sw_satAr $Revision: 1.1 $  $Date: 1998/04/22 02:15:56 $
+%          Copyright (C) CSIRO, Phil Morgan 1998.
+%
+% USAGE:  satAr = sw_satAr(S,T,P)
+%
+% DESCRIPTION:
+%    Solubility (satuaration) of Argon (Ar) in sea water
+%
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78)]
+%   T = temperature [degree C (IPTS-68)]
+%
+% OUTPUT:
+%   satAr = solubility of Ar  [ml/l] 
+% 
+% AUTHOR:  Phil Morgan 97-11-05  (morgan@ml.csiro.au)
+%
+%$$$ #include "disclaimer_in_code.inc"
+%
+% REFERENCES:
+%    Weiss, R. F. 1970
+%    "The solubility of nitrogen, oxygen and argon in water and seawater."
+%    Deap-Sea Research., 1970, Vol 17, pp721-735.
+%=========================================================================
+
+% CALLER: general purpose
+% CALLEE: 
+
+%$$$ #ifdef VARIANT_PRIVATE
+%$$$ %***********************************************************
+%$$$ %$Id: sw_satAr.M,v 1.1 1998/04/22 02:15:56 morgan Exp $
+%$$$ %
+%$$$ %$Log: sw_satAr.M,v $
+
+%$$$ %
+%$$$ %***********************************************************
+%$$$ #endif
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+if nargin ~=2
+   error('sw_satAr.m: Must pass 2 parameters')
+end %if
+
+% CHECK S,T dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('sw_satAr: S & T must have same dimensions')
+end %if
+
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if ms == 1  % row vector
+   T       =  T(:);
+   S       =  S(:);   
+   Transpose = 1;
+end %if
+
+%------
+% BEGIN
+%------
+
+% convert T to Kelvin
+T = 273.15 + T; 
+
+% constants for Eqn (4) of Weiss 1970
+a1 = -173.5146;
+a2 =  245.4510;
+a3 =  141.8222;
+a4 =  -21.8020;
+b1 =   -0.034474;
+b2 =    0.014934;
+b3 =   -0.0017729;
+
+% Eqn (4) of Weiss 1970
+lnC = a1 + a2.*(100./T) + a3.*log(T./100) + a4.*(T./100) + ...
+      S.*( b1 + b2.*(T./100) + b3.*((T./100).^2) );
+
+c = exp(lnC);
+
+return
+
diff --git a/M_Csiro/sw_satN2.m b/M_Csiro/sw_satN2.m
new file mode 100644
index 0000000000000000000000000000000000000000..f1a0ca8f74be7416a9d3b4d785d4a749cff21324
--- /dev/null
+++ b/M_Csiro/sw_satN2.m
@@ -0,0 +1,95 @@
+%$$$ 
+%$$$ #undef __PR
+%$$$ #include "VARIANT.h"
+
+function c = sw_satN2(S,T)
+
+% SW_SATN2   Satuaration of N2 in sea water
+%=========================================================================
+% sw_satN2 $Revision: 1.1 $  $Date: 1998/04/22 02:15:56 $
+%          Copyright (C) CSIRO, Phil Morgan 1998.
+%
+% USAGE:  satN2 = sw_satN2(S,T,P)
+%
+% DESCRIPTION:
+%    Solubility (satuaration) of Nitrogen (N2) in sea water
+%
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78)]
+%   T = temperature [degree C (IPTS-68)]
+%
+% OUTPUT:
+%   satN2 = solubility of N2  [ml/l] 
+% 
+% AUTHOR:  Phil Morgan 97-11-05  (morgan@ml.csiro.au)
+%
+%$$$ #include "disclaimer_in_code.inc"
+%
+% REFERENCES:
+%    Weiss, R. F. 1970
+%    "The solubility of nitrogen, oxygen and argon in water and seawater."
+%    Deap-Sea Research., 1970, Vol 17, pp721-735.
+%=========================================================================
+
+% CALLER: general purpose
+% CALLEE: 
+
+%$$$ #ifdef VARIANT_PRIVATE
+%$$$ %***********************************************************
+%$$$ %$Id: sw_satN2.M,v 1.1 1998/04/22 02:15:56 morgan Exp $
+%$$$ %
+%$$$ %$Log: sw_satN2.M,v $
+
+%$$$ %
+%$$$ %***********************************************************
+%$$$ #endif
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+if nargin ~=2
+   error('sw_satN2.m: Must pass 2 parameters')
+end %if
+
+% CHECK S,T dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('sw_satN2: S & T must have same dimensions')
+end %if
+
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if ms == 1  % row vector
+   T       =  T(:);
+   S       =  S(:);   
+   Transpose = 1;
+end %if
+
+%------
+% BEGIN
+%------
+
+% convert T to Kelvin
+T = 273.15 + T; 
+
+% constants for Eqn (4) of Weiss 1970
+a1 = -172.4965;
+a2 =  248.4262;
+a3 =  143.0738;
+a4 =  -21.7120;
+b1 =   -0.049781;
+b2 =    0.025018;
+b3 =   -0.0034861;
+
+% Eqn (4) of Weiss 1970
+lnC = a1 + a2.*(100./T) + a3.*log(T./100) + a4.*(T./100) + ...
+      S.*( b1 + b2.*(T./100) + b3.*((T./100).^2) );
+
+c = exp(lnC);
+
+return
+
diff --git a/M_Csiro/sw_satO2.m b/M_Csiro/sw_satO2.m
new file mode 100644
index 0000000000000000000000000000000000000000..ea1b5c73ea44be799c440d45ddb7c6f5d715fd20
--- /dev/null
+++ b/M_Csiro/sw_satO2.m
@@ -0,0 +1,95 @@
+%$$$ 
+%$$$ #undef __PR
+%$$$ #include "VARIANT.h"
+
+function c = sw_satO2(S,T)
+
+% SW_SATO2   Satuaration of O2 in sea water
+%=========================================================================
+% sw_satO2 $Revision: 1.1 $  $Date: 1998/04/22 02:15:56 $
+%          Copyright (C) CSIRO, Phil Morgan 1998.
+%
+% USAGE:  satO2 = sw_satO2(S,T,P)
+%
+% DESCRIPTION:
+%    Solubility (satuaration) of Oxygen (O2) in sea water
+%
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78)]
+%   T = temperature [degree C (IPTS-68)]
+%
+% OUTPUT:
+%   satO2 = solubility of O2  [ml/l] 
+% 
+% AUTHOR:  Phil Morgan 97-11-05  (morgan@ml.csiro.au)
+%
+%$$$ #include "disclaimer_in_code.inc"
+%
+% REFERENCES:
+%    Weiss, R. F. 1970
+%    "The solubility of nitrogen, oxygen and argon in water and seawater."
+%    Deap-Sea Research., 1970, Vol 17, pp721-735.
+%=========================================================================
+
+% CALLER: general purpose
+% CALLEE: 
+
+%$$$ #ifdef VARIANT_PRIVATE
+%$$$ %***********************************************************
+%$$$ %$Id: sw_satO2.M,v 1.1 1998/04/22 02:15:56 morgan Exp $
+%$$$ %
+%$$$ %$Log: sw_satO2.M,v $
+
+%$$$ %
+%$$$ %***********************************************************
+%$$$ #endif
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+if nargin ~=2
+   error('sw_satO2.m: Must pass 2 parameters')
+end %if
+
+% CHECK S,T dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('sw_satO2: S & T must have same dimensions')
+end %if
+
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if ms == 1  % row vector
+   T       =  T(:);
+   S       =  S(:);   
+   Transpose = 1;
+end %if
+
+%------
+% BEGIN
+%------
+
+% convert T to Kelvin
+T = 273.15 + T; 
+
+% constants for Eqn (4) of Weiss 1970
+a1 = -173.4292;
+a2 =  249.6339;
+a3 =  143.3483;
+a4 =  -21.8492;
+b1 =   -0.033096;
+b2 =    0.014259;
+b3 =   -0.0017000;
+
+% Eqn (4) of Weiss 1970
+lnC = a1 + a2.*(100./T) + a3.*log(T./100) + a4.*(T./100) + ...
+      S.*( b1 + b2.*(T./100) + b3.*((T./100).^2) );
+
+c = exp(lnC);
+
+return
+
diff --git a/M_Csiro/sw_seck.m b/M_Csiro/sw_seck.m
new file mode 100644
index 0000000000000000000000000000000000000000..ae1dc17696c03b0a2bf949e42527e0d3d254863d
--- /dev/null
+++ b/M_Csiro/sw_seck.m
@@ -0,0 +1,167 @@
+
+function K = sw_seck(S,T,P)
+
+% SW_SECK    Secant bulk modulus (K) of sea water
+%=========================================================================
+% SW_SECK  $Revision: 1.3 $  $Date: 1994/10/10 05:50:45 $
+%          Copyright (C) CSIRO, Phil Morgan 1992.
+%
+% USAGE:  dens = sw_seck(S,T,P)
+%
+% DESCRIPTION:
+%    Secant Bulk Modulus (K) of Sea Water using Equation of state 1980. 
+%    UNESCO polynomial implementation.
+%
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78) ]
+%   T = temperature [degree C (IPTS-68)]
+%   P = pressure    [db]
+%       (alternatively, may have dimensions 1*1 or 1*n where n is columns in S)
+%
+% OUTPUT:
+%   K = Secant Bulk Modulus  [bars]
+% 
+% AUTHOR:  Phil Morgan 92-11-05  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonoff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%    Eqn.(15) p.18
+%
+%    Millero, F.J. and  Poisson, A.
+%    International one-atmosphere equation of state of seawater.
+%    Deep-Sea Res. 1981. Vol28A(6) pp625-629.
+%=========================================================================
+
+% CALLER: sw_dens.m
+% CALLEE: none
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+if nargin ~=3
+   error('sw_seck.m: Must pass 3 parameters')
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+%--------------------------------------------------------------------
+% COMPUTE COMPRESSION TERMS
+%--------------------------------------------------------------------
+P = P/10;  %convert from db to atmospheric pressure units
+
+% Pure water terms of the secant bulk modulus at atmos pressure.
+% UNESCO eqn 19 p 18
+
+h3 = -5.77905E-7;
+h2 = +1.16092E-4;
+h1 = +1.43713E-3;
+h0 = +3.239908;   %[-0.1194975];
+
+%AW = h0 + h1*T + h2*T.^2 + h3*T.^3;
+AW  = h0 + (h1 + (h2 + h3.*T).*T).*T;
+
+k2 =  5.2787E-8;
+k1 = -6.12293E-6;
+k0 =  +8.50935E-5;   %[+3.47718E-5];
+
+BW  = k0 + (k1 + k2*T).*T;
+%BW = k0 + k1*T + k2*T.^2;
+
+e4 = -5.155288E-5;
+e3 = +1.360477E-2;
+e2 = -2.327105;
+e1 = +148.4206;
+e0 = 19652.21;    %[-1930.06];
+
+KW  = e0 + (e1 + (e2 + (e3 + e4*T).*T).*T).*T;   % eqn 19
+%KW = e0 + e1*T + e2*T.^2 + e3*T.^3 + e4*T.^4;
+
+%--------------------------------------------------------------------
+% SEA WATER TERMS OF SECANT BULK MODULUS AT ATMOS PRESSURE.
+%--------------------------------------------------------------------
+j0 = 1.91075E-4;
+       
+i2 = -1.6078E-6;
+i1 = -1.0981E-5;
+i0 =  2.2838E-3;
+
+SR = sqrt(S);
+
+A  = AW + (i0 + (i1 + i2*T).*T + j0*SR).*S; 
+%A = AW + (i0 + i1*T + i2*T.^2 + j0*SR).*S;  % eqn 17
+  
+
+m2 =  9.1697E-10;
+m1 = +2.0816E-8;
+m0 = -9.9348E-7;
+
+B = BW + (m0 + (m1 + m2*T).*T).*S;   % eqn 18
+%B  = BW + (m0 + m1*T + m2*T.^2).*S;   % eqn 18
+
+    
+f3 =  -6.1670E-5;
+f2 =  +1.09987E-2;
+f1 =  -0.603459;
+f0 = +54.6746;
+
+g2 = -5.3009E-4;
+g1 = +1.6483E-2;
+g0 = +7.944E-2;
+
+K0 = KW + (  f0 + (f1 + (f2 + f3*T).*T).*T ...
+        +   (g0 + (g1 + g2*T).*T).*SR         ).*S;      % eqn 16
+
+%K0  = KW + (f0 + f1*T + f2*T.^2 + f3*T.^3).*S ...
+%         + (g0 + g1*T + g2*T.^2).*SR.*S;             % eqn 16
+
+K = K0 + (A + B.*P).*P;  % eqn 15
+
+if Transpose
+   K = K';
+end %if
+
+return
+%----------------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_smow.m b/M_Csiro/sw_smow.m
new file mode 100644
index 0000000000000000000000000000000000000000..725e5e2155de287f99f259f764f82492d14c4d89
--- /dev/null
+++ b/M_Csiro/sw_smow.m
@@ -0,0 +1,69 @@
+
+function dens = sw_smow(T)
+
+% SW_SMOW    Denisty of standard mean ocean water (pure water)
+%=========================================================================
+% SW_SMOW  $Revision: 1.3 $  $Date: 1994/10/10 05:51:46 $
+%          Copyright (C) CSIRO, Phil Morgan 1992.
+%
+% USAGE:  dens = sw_smow(T)
+%
+% DESCRIPTION:
+%    Denisty of Standard Mean Ocean Water (Pure Water) using EOS 1980. 
+%
+% INPUT: 
+%   T = temperature [degree C (IPTS-68)]
+%
+% OUTPUT:
+%   dens = density  [kg/m^3] 
+% 
+% AUTHOR:  Phil Morgan 92-11-05  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%     Unesco 1983. Algorithms for computation of fundamental properties of 
+%     seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%     UNESCO 1983 p17  Eqn(14)
+%
+%     Millero, F.J & Poisson, A.
+%     INternational one-atmosphere equation of state for seawater.
+%     Deep-Sea Research Vol28A No.6. 1981 625-629.    Eqn (6)
+%=========================================================================
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+% TEST INPUTS
+if nargin ~= 1
+   error('sw_smow.m: Only one input argument allowed')
+end %if
+
+Transpose = 0;
+[mT,nT] = size(T);
+if mT == 1 % a row vector
+   T = T(:);
+   Tranpose = 1;
+end %if
+
+%----------------------
+% DEFINE CONSTANTS
+%----------------------
+a0 = 999.842594;
+a1 =   6.793952e-2;
+a2 =  -9.095290e-3;
+a3 =   1.001685e-4;
+a4 =  -1.120083e-6;
+a5 =   6.536332e-9;
+
+dens = a0 + (a1 + (a2 + (a3 + (a4 + a5*T).*T).*T).*T).*T;
+
+if Transpose
+  dens = dens';
+end %if
+
+return
+%--------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_svan.m b/M_Csiro/sw_svan.m
new file mode 100644
index 0000000000000000000000000000000000000000..c714431ee11d21d32668da223e5c698cdc29253c
--- /dev/null
+++ b/M_Csiro/sw_svan.m
@@ -0,0 +1,104 @@
+
+function svan = sw_svan(S,T,P)
+
+% SW_SVAN    Specific volume anomaly
+%=========================================================================
+% SW_SVAN  $Revision: 1.3 $  $Date: 1994/10/10 05:52:20 $
+%          Copyright (C) CSIRO,  Phil Morgan 1992.
+%
+% USAGE:  svan = sw_svan(S,T,P)
+%
+% DESCRIPTION:
+%   Specific Volume Anomaly calculated as 
+%        svan = 1/sw_dens(s,t,p) - 1/sw_dens(35,0,p)
+%   Note that it is often quoted in literature as 1e8*units
+%
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78) ]
+%   T = temperature [degree C (IPTS-68)]
+%   P = Pressure    [db]
+%       (alternatively, may have dimensions 1*1 or 1*n where n is columns in S)
+%
+% OUTPUT:
+%  svan = Specific Volume Anomaly  [m^3 kg^-1]
+%
+% AUTHOR:  Phil Morgan 92-11-05  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCE:
+%     Fofonoff, N.P. and Millard, R.C. Jr 
+%     Unesco 1983. Algorithms for computation of fundamental properties of 
+%     seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%     Eqn (9) p.15. 
+%
+%     S. Pond & G.Pickard  2nd Edition 1986
+%     Introductory Dynamical Oceanogrpahy
+%     Pergamon Press Sydney.  ISBN 0-08-028728-X
+%=========================================================================
+
+%
+% CALLER: general purpose
+% CALLEE: sw_dens.m
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+if nargin ~=3
+   error('sw_svan.m: Must pass 3 parameters')
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+ 
+% -----
+% BEGIN
+% -----
+svan = (  ones(size(S)) ./ sw_dens(S,T,P)) - ...
+         (ones(size(S)) ./ sw_dens(35*ones(size(S)),zeros(size(S)),P) );
+
+if Transpose
+   svan = svan';
+end %if
+
+return
+%--------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_svel.m b/M_Csiro/sw_svel.m
new file mode 100644
index 0000000000000000000000000000000000000000..f3dff2ff58a98cc0d0a1b75e451e3817354d5a9b
--- /dev/null
+++ b/M_Csiro/sw_svel.m
@@ -0,0 +1,181 @@
+
+function svel = sw_svel(S,T,P)
+
+% SW_SVEL    Sound velocity of sea water
+%=========================================================================
+% SW_SVEL  $Revision: 1.3 $  $Date: 1994/10/10 05:53:00 $
+%          Copyright (C) CSIRO, Phil Morgan 1993.
+%
+% USAGE:  svel = sw_svel(S,T,P)
+%
+% DESCRIPTION:
+%    Sound Velocity in sea water using UNESCO 1983 polynomial.
+%
+% INPUT:  (all must have same dimensions)
+%   S = salinity    [psu      (PSS-78)]
+%   T = temperature [degree C (IPTS-68)]
+%   P = pressure    [db]
+%       (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+% OUTPUT:
+%   svel = sound velocity  [m/s] 
+% 
+% AUTHOR:  Phil Morgan 93-04-20  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonoff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%=========================================================================
+
+% CALLER: general purpose
+% CALLEE: none
+
+% UNESCO 1983. eqn.33  p.46
+
+%----------------------
+% CHECK INPUT ARGUMENTS
+%----------------------
+if nargin ~=3
+   error('sw_svel.m: Must pass 3 parameters')
+end %if
+
+% CHECK S,T,P dimensions and verify consistent
+[ms,ns] = size(S);
+[mt,nt] = size(T);
+[mp,np] = size(P);
+
+  
+% CHECK THAT S & T HAVE SAME SHAPE
+if (ms~=mt) | (ns~=nt)
+   error('check_stp: S & T must have same dimensions')
+end %if
+
+% CHECK OPTIONAL SHAPES FOR P
+if     mp==1  & np==1      % P is a scalar.  Fill to size of S
+   P = P(1)*ones(ms,ns);
+elseif np==ns & mp==1      % P is row vector with same cols as S
+   P = P( ones(1,ms), : ); %   Copy down each column.
+elseif mp==ms & np==1      % P is column vector
+   P = P( :, ones(1,ns) ); %   Copy across each row
+elseif mp==ms & np==ns     % PR is a matrix size(S)
+   % shape ok 
+else
+   error('check_stp: P has wrong dimensions')
+end %if
+[mp,np] = size(P);
+ 
+
+  
+% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
+Transpose = 0;
+if mp == 1  % row vector
+   P       =  P(:);
+   T       =  T(:);
+   S       =  S(:);   
+
+   Transpose = 1;
+end %if
+%***check_stp
+
+%---------
+% BEGIN
+%--------
+
+P = P/10;  % convert db to bars as used in UNESCO routines
+
+%------------
+% eqn 34 p.46
+%------------
+c00 = 1402.388;
+c01 =    5.03711;
+c02 =   -5.80852e-2;
+c03 =    3.3420e-4;
+c04 =   -1.47800e-6;
+c05 =    3.1464e-9;
+
+c10 =  0.153563;
+c11 =  6.8982e-4;
+c12 = -8.1788e-6;
+c13 =  1.3621e-7;
+c14 = -6.1185e-10;
+
+c20 =  3.1260e-5;
+c21 = -1.7107e-6;
+c22 =  2.5974e-8;
+c23 = -2.5335e-10;
+c24 =  1.0405e-12;
+
+c30 = -9.7729e-9;
+c31 =  3.8504e-10;
+c32 = -2.3643e-12;
+
+Cw =    c00 + c01.*T + c02.*T.^2 + c03.*T.^3 + c04.*T.^4 + c05.*T.^5   ...
+     + (c10 + c11.*T + c12.*T.^2 + c13.*T.^3 + c14.*T.^4).*P          ...
+     + (c20 + c21.*T + c22.*T.^2 + c23.*T.^3 + c24.*T.^4).*P.^2        ...
+     + (c30 + c31.*T + c32.*T.^2).*P.^3;
+ 
+%-------------
+% eqn 35. p.47
+%-------------
+a00 =  1.389;
+a01 = -1.262e-2;
+a02 =  7.164e-5;
+a03 =  2.006e-6;
+a04 = -3.21e-8;
+
+a10 =  9.4742e-5;
+a11 = -1.2580e-5;
+a12 = -6.4885e-8;
+a13 =  1.0507e-8;
+a14 = -2.0122e-10;
+
+a20 = -3.9064e-7;
+a21 =  9.1041e-9;
+a22 = -1.6002e-10;
+a23 =  7.988e-12;
+
+a30 =  1.100e-10;
+a31 =  6.649e-12;
+a32 = -3.389e-13;
+
+A =     a00 + a01.*T + a02.*T.^2 + a03.*T.^3 + a04.*T.^4       ...
+     + (a10 + a11.*T + a12.*T.^2 + a13.*T.^3 + a14.*T.^4).*P   ...
+     + (a20 + a21.*T + a22.*T.^2 + a23.*T.^3).*P.^2            ...
+     + (a30 + a31.*T + a32.*T.^2).*P.^3;
+
+ 
+%------------ 
+% eqn 36 p.47
+%------------ 
+b00 = -1.922e-2;
+b01 = -4.42e-5;
+b10 =  7.3637e-5;
+b11 =  1.7945e-7;
+
+B = b00 + b01.*T + (b10 + b11.*T).*P;
+
+%------------ 
+% eqn 37 p.47
+%------------ 
+d00 =  1.727e-3;
+d10 = -7.9836e-6;
+
+D = d00 + d10.*P;
+
+%------------
+% eqn 33 p.46
+%------------
+svel = Cw + A.*S + B.*S.*sqrt(S) + D.*S.^2;
+
+if Transpose
+   svel = svel';
+end %if
+
+return
+%--------------------------------------------------------------------------
+
diff --git a/M_Csiro/sw_temp.m b/M_Csiro/sw_temp.m
new file mode 100644
index 0000000000000000000000000000000000000000..c62e5690b29ba847d83405a84295660f28de74dd
--- /dev/null
+++ b/M_Csiro/sw_temp.m
@@ -0,0 +1,59 @@
+
+function PT = sw_temp(S,T,P,PR)
+
+% SW_TEMP    Temperature from potential temperature
+%===========================================================================
+% TEMP  $Revision: 1.3 $  $Date: 1994/10/10 05:53:39 $
+%       Copyright (C) CSIRO, Phil Morgan  1992. 
+%
+% USAGE:  temp = sw_temp(S,PTMP,P,PR) 
+%
+% DESCRIPTION:
+%    Calculates temperature from potential temperature at the reference
+%    pressure PR and in-situ pressure P.
+%   
+% INPUT:  (all must have same dimensions)
+%   S     = salinity              [psu      (PSS-78) ]
+%   PTMP  = potential temperature [degree C (IPTS-68)]
+%   P     = pressure              [db]
+%   PR    = Reference pressure    [db]
+%           (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
+%
+% OUTPUT:
+%   temp = temperature [degree C (IPTS-68)]
+%
+% AUTHOR:  Phil Morgan 92-04-06  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+% REFERENCES:
+%    Fofonoff, P. and Millard, R.C. Jr
+%    Unesco 1983. Algorithms for computation of fundamental properties of 
+%    seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
+%    Eqn.(31) p.39
+%
+%    Bryden, H. 1973.
+%    "New Polynomials for thermal expansion, adiabatic temperature gradient
+%    and potential temperature of sea water."
+%    DEEP-SEA RES., 1973, Vol20,401-408.
+%=========================================================================
+
+% CALLER:  general purpose
+% CALLEE:  sw_ptmp.m
+
+%-------------
+% CHECK INPUTS
+%-------------
+if nargin ~= 4
+   error('sw_temp.m: Must pass 4 parameters ')
+end %if
+% LET sw_ptmp.m DO DIMENSION CHECKING
+
+% CARRY OUT INVERSE CALCULATION BY SWAPPING P0 & PR.
+PT = sw_ptmp(S,T,PR,P);
+
+return      
+%=========================================================================
+
diff --git a/M_Csiro/sw_test.m b/M_Csiro/sw_test.m
new file mode 100644
index 0000000000000000000000000000000000000000..26f235063b9702ecd361c71aba3642c067097a1f
--- /dev/null
+++ b/M_Csiro/sw_test.m
@@ -0,0 +1,600 @@
+
+function [] = sw_test()
+
+% SW_TEST    Test SEAWATER Library Routines
+%=========================================================================
+% SW_TEST   $Revision: 1.4 $  $Date: 1998/04/22 02:08:04 $
+%           Copyright (C) CSIRO, Phil Morgan 1994
+%
+% sw_test
+%
+% DESCRIPTION:
+%    Execute test routines to test and verify SEAWATER Library routines
+%    for your platform.  Prints output to screen and to file sw_test.txt
+%
+%    Use the "more" command to scroll results to screen
+%
+% OUTPUT:
+%   file sw_test.txt
+%
+% AUTHOR:  Phil Morgan  (morgan@ml.csiro.au)
+%
+% DISCLAIMER:
+%   This software is provided "as is" without warranty of any kind.  
+%   See the file sw_copy.m for conditions of use and licence.
+%
+%=========================================================================
+
+delete sw_test.txt
+disp('OUTPUT FROM THIS TEST WILL ALSO BE SAVED IN FILE sw_test.txt')
+disp(' <enter> to continue...')
+pause
+reply        = input('Full listing of help for each routine (y/n) ? ','s');
+display_help = strcmp(reply,'y') | strcmp(reply,'Y');
+
+format compact
+echo off
+diary sw_test.txt
+
+disp( '***********************')
+disp( '    TEST REPORT    ')
+disp( ' ')
+disp( ' SEA WATER LIBRARY ')
+disp( ' ')
+sw_ver
+disp( ' ')
+disp(['Matlab Version ' version ])
+disp( ' ')
+disp(['   ' date ''])
+disp( '***********************')
+
+disp(' ')
+
+%--------------------------------
+% TEST MAIN MODULE  sw_ptmp.m
+%      SUB-MODULES  sw_atg.m
+%--------------------------------
+module     = 'sw_ptmp';
+submodules = 'sw_adtg.m';
+disp('*************************************')
+disp(['**  TESTING MODULE: ' module])
+disp(['**  and SUB-MODULE: ' submodules])
+disp('*************************************')
+if display_help
+   eval(['help ' module])
+   eval(['help ' submodules])
+end %if
+
+% TEST 1 - data from Unesco 1983 p45
+
+T    = [ 0  0  0  0  0  0; 
+        10 10 10 10 10 10;
+	20 20 20 20 20 20;
+	30 30 30 30 30 30;
+	40 40 40 40 40 40 ];
+    
+S    = [25 25 25 35 35 35;
+        25 25 25 35 35 35;
+	25 25 25 35 35 35 ;
+	25 25 25 35 35 35;
+	25 25 25 35 35 35 ];
+
+P    = [0 5000 10000 0 5000 10000;
+        0 5000 10000 0 5000 10000;
+	0 5000 10000 0 5000 10000 ;
+	0 5000 10000 0 5000 10000;
+	0 5000 10000 0 5000 10000 ];
+
+Pr = [0 0 0 0 0 0];
+
+UN_ptmp =      [ 0  -0.3061  -0.9667   0  -0.3856 -1.0974;
+                10   9.3531   8.4684  10   9.2906  8.3643;
+	        20  19.0438  17.9426  20  18.9985 17.8654;
+	        30  28.7512  27.4353  30  28.7231 27.3851;
+	        40  38.4607  36.9254  40  38.4498 36.9023];
+
+ptmp    = sw_ptmp(S,T,P,Pr);
+
+%----------------
+% DISPLAY RESULTS
+%----------------
+disp(' ')
+disp   ('********************************************************')
+disp   ('Comparison of accepted values from UNESCO 1983 ')
+disp   (' (Unesco Tech. Paper in Marine Sci. No. 44, p45)')
+disp   (['with computed results from ' module ' on ' computer ' computer'])
+disp   ('********************************************************')
+
+for icol = 1:length(S(1,:))
+disp(' ')
+disp   ('   Sal  Temp  Press     PTMP       sw_ptmp')
+disp   ('  (psu)  (C)   (db)     (C)          (C)')
+  fprintf(1,' %4.0f  %4.0f   %5.0f   %8.4f  %11.5f\n', ...
+  [S(:,icol) T(:,icol) P(:,icol) UN_ptmp(:,icol) ptmp(:,icol)]');
+end %for
+
+%-------------------------------------------------------------------------------
+% TEST MAIN MODULE  sw_svan.m
+%      SUB-MODULES  sw_dens.m sw_dens0.m sw_smow.m sw_seck.m sw_pden.m sw_ptmp.m
+%------------------------------------------------------------------------------
+module     = 'sw_svan.m';
+submodules = 'sw_dens.m sw_dens0.m sw_smow.m sw_seck.m sw_pden.m sw_ptmp.m';
+disp(' ')
+disp('************************************************************************')
+disp(['**  TESTING MODULE: ' module])
+disp(['**  and SUB-MODULE: ' submodules])
+disp('************************************************************************')
+if display_help
+   eval(['help ' module])
+   eval(['help ' submodules])
+end %if
+
+% TEST DATA FROM 
+% Unesco Tech. Paper in Marine Sci. No. 44, p22
+
+s = [0     0   0      0  35    35  35    35]';
+p = [0 10000   0  10000   0 10000   0 10000]';
+t = [0     0  30     30   0     0  30    30]';
+
+UN_svan = [2749.54 2288.61 3170.58 3147.85 ...
+              0.0     0.00  607.14  916.34]';
+	   
+svan    = sw_svan(s,t,p);
+
+%----------------
+% DISPLAY RESULTS
+%----------------
+disp(' ')
+disp   ('********************************************************')
+disp   ('Comparison of accepted values from UNESCO 1983')
+disp   (' (Unesco Tech. Paper in Marine Sci. No. 44, p22)')
+disp   (['with computed results from ' module ' on ' computer ' computer'])
+disp   ('********************************************************')
+disp(' ')
+disp   ('   Sal  Temp  Press        SVAN        sw_svan')
+disp   ('  (psu)  (C)   (db)    (1e-8*m3/kg)  (1e-8*m3/kg)')
+fprintf(1,' %4.0f  %4.0f   %5.0f   %11.2f    %11.3f\n',[s t p UN_svan 1e+8*svan]');
+
+%-------------------------------------------------------------------------------
+% TEST MAIN MODULE  
+%      SUB-MODULES  
+%------------------------------------------------------------------------------
+module     = 'sw_salt.m';
+submodules = 'sw_salrt.m sw_salrp.m sw_sals.m';
+disp(' ')
+disp('************************************************************************')
+disp(['**  TESTING MODULE: ' module])
+disp(['**  and SUB-MODULE: ' submodules])
+disp('************************************************************************')
+if display_help
+   eval(['help ' module])
+   eval(['help ' submodules])
+end %if
+
+% TEST 1 - data from Unesco 1983 p9
+%***************************************************************************
+
+R     = [1 1.2 0.65]';   % cndr = R
+T     = [15 20 5]';
+P     = [0 2000 1500]';
+
+Rt    = [1 1.0568875 0.81705885]';
+UN_S  = [35 37.245628 27.995347]';
+S     = sw_salt(R,T,P);
+
+%----------------
+% DISPLAY RESULTS
+%----------------
+disp(' ')
+disp   ('********************************************************')
+disp   ('Comparison of accepted values from UNESCO 1983 ')
+disp   (' (Unesco Tech. Paper in Marine Sci. No. 44, p9)')
+disp   (['with computed results from ' module ' on ' computer ' computer'])
+disp   ('********************************************************')
+disp(' ')
+disp   ('   Temp    Press       R              S           sw_salt')
+disp   ('   (C)     (db)    (no units)       (psu)          (psu) ')
+table = [T P R UN_S S]';
+fprintf(1,' %4.0f       %4.0f  %8.2f      %11.6f  %14.7f\n', table);
+
+%-------------------------------------------------------------------------------
+% TEST MAIN MODULE  
+%      SUB-MODULES  
+%------------------------------------------------------------------------------
+module     = 'sw_cndr.m';
+submodules = 'sw_salds.m';
+disp(' ')
+disp('************************************************************************')
+disp(['**  TESTING MODULE: ' module])
+disp(['**  and SUB-MODULE: ' submodules])
+disp('************************************************************************')
+if display_help
+   eval(['help ' module])
+   eval(['help ' submodules])
+end %if
+
+% TEST 1 - data from Unesco 1983 p9
+
+T    = [0   10     0   10  10  30]';
+P    = [0    0  1000 1000   0   0]';
+S    = [25  25    25   25  40  40]';
+UN_R = [ 0.498088 0.654990 0.506244 0.662975 1.000073 1.529967]';
+R    = sw_cndr(S,T,P);
+
+%----------------
+% DISPLAY RESULTS
+%----------------
+disp(' ')
+disp   ('********************************************************')
+disp   ('Comparison of accepted values from UNESCO 1983 ')
+disp   (' (Unesco Tech. Paper in Marine Sci. No. 44, p14)')
+disp   (['with computed results from ' module ' on ' computer ' computer'])
+disp   ('********************************************************')
+disp(' ')
+disp   ('   Temp    Press       S            cndr         sw_cndr')
+disp   ('   (C)     (db)      (psu)        (no units)    (no units) ')
+table = [T P S UN_R R]';
+fprintf(1,' %4.0f       %4.0f   %8.6f   %11.6f  %14.8f\n', table);
+
+%-------------------------------------------------------------------------------
+% TEST MAIN MODULE  
+%      SUB-MODULES  
+%------------------------------------------------------------------------------
+module     = 'sw_dpth.m';
+disp(' ')
+disp('************************************************************************')
+disp(['**  TESTING MODULE: ' module])
+disp('************************************************************************')
+if display_help
+   eval(['help ' module])
+end %if
+
+% TEST DATA - matrix "pressure", vector "lat"  Unesco 1983 data p30.
+
+lat = [0 30 45 90];
+P   = [  500   500   500   500;
+        5000  5000  5000  5000;
+       10000 10000 10000 10000];
+
+UN_dpth = [   496.65   496.00   495.34   494.03;
+             4915.04  4908.56  4902.08  4889.13;
+	     9725.47  9712.65  9699.84  9674.23];
+
+dpth = sw_dpth(P,lat);
+
+%----------------
+% DISPLAY RESULTS
+%----------------
+disp(' ')
+disp   ('********************************************************')
+disp   ('Comparison of accepted values from Unesco 1983 ')
+disp   (' (Unesco Tech. Paper in Marine Sci. No. 44, p28)')
+disp   (['with computed results from ' module ' on ' computer ' computer'])
+disp   ('********************************************************')
+
+for irow = 1:3
+   disp(' ')
+   disp   ('    Lat       Press     DPTH      sw_dpth')
+   disp   ('  (degree)    (db)     (meter)    (meter)')
+   table = [lat' P(irow,:)' UN_dpth(irow,:)' dpth(irow,:)'];
+   fprintf(1,'  %6.3f     %6.0f   %8.2f   %8.3f\n', table')
+end %for
+
+%-------------------------------------------------------------------------------
+% TEST MAIN MODULE  
+%      SUB-MODULES  
+%------------------------------------------------------------------------------
+module     = 'sw_fp.m';
+disp(' ')
+disp('************************************************************************')
+disp(['**  TESTING MODULE: ' module])
+disp('************************************************************************')
+if display_help
+   eval(['help ' module])
+end %if
+
+% TEST 1 - 
+% UNESCO DATA p.30
+%***************************************************************************
+S    = [ 5   10  15  20  25  30  35  40;
+         5   10  15  20  25  30  35  40];
+
+P    = [  0   0   0   0   0  0    0   0;
+        500 500 500 500 500 500 500 500];
+    
+
+UN_fp = [-0.274 -0.542 -0.812 -1.083 -1.358 -1.638 -1.922 -2.212;
+         -0.650 -0.919 -1.188 -1.460 -1.735 -2.014 -2.299 -2.589];
+
+fp    = sw_fp(S,P);
+
+%----------------
+% DISPLAY RESULTS
+%----------------
+disp(' ')
+disp   ('********************************************************')
+disp   ('Comparison of accepted values from UNESCO 1983 ')
+disp   (' (Unesco Tech. Paper in Marine Sci. No. 44, p30)')
+disp   (['with computed results from ' module ' on ' computer ' computer'])
+disp   ('********************************************************')
+
+for irow = 1:2
+  disp(' ')
+  disp   ('   Sal   Press      fp        sw_fp')
+  disp   ('  (psu)   (db)      (C)        (C)')
+  table = [S(irow,:); P(irow,:); UN_fp(irow,:); fp(irow,:)];
+  fprintf(1,' %4.0f   %5.0f   %8.3f  %11.4f\n', table)
+end %for
+
+%-------------------------------------------------------------------------------
+% TEST MAIN MODULE  
+%      SUB-MODULES  
+%------------------------------------------------------------------------------
+module     = 'sw_cp.m';
+disp(' ')
+disp('************************************************************************')
+disp(['**  TESTING MODULE: ' module])
+disp('************************************************************************')
+if display_help
+   eval(['help ' module])
+end %if
+
+% TEST 1 - 
+% DATA FROM POND AND PICKARD INTRO. DYNAMICAL OCEANOGRAPHY 2ND ED. 1986
+%***************************************************************************
+
+T    = [ 0  0  0  0  0  0;
+        10 10 10 10 10 10;
+	20 20 20 20 20 20;
+	30 30 30 30 30 30;
+	40 40 40 40 40 40 ];
+    
+S    = [25 25 25 35 35 35;
+        25 25 25 35 35 35;
+	25 25 25 35 35 35 ;
+	25 25 25 35 35 35;
+	25 25 25 35 35 35 ];
+
+P    = [0 5000 10000 0 5000 10000;
+        0 5000 10000 0 5000 10000;
+	0 5000 10000 0 5000 10000 ;
+	0 5000 10000 0 5000 10000;
+	0 5000 10000 0 5000 10000 ];
+
+UN_cp =      [  4048.4  3896.3  3807.7  3986.5  3849.3  3769.1;
+                4041.8  3919.6  3842.3  3986.3  3874.7  3804.4;
+	        4044.8  3938.6  3866.7  3993.9  3895.0  3828.3;
+	        4049.1  3952.0  3883.0  4000.7  3909.2  3844.3;
+	        4051.2  3966.1  3905.9  4003.5  3923.9  3868.3 ];
+
+cp    = sw_cp(S,T,P);
+
+%----------------
+% DISPLAY RESULTS
+%----------------
+disp(' ')
+disp   ('********************************************************')
+disp   ('Comparison of accepted values from UNESCO 1983 ')
+disp   (' (Unesco Tech. Paper in Marine Sci. No. 44, p37)')
+disp   (['with computed results from ' module ' on ' computer ' computer'])
+disp   ('********************************************************')
+
+for icol = 1:length(S(1,:))
+  disp(' ')
+  disp   ('   Sal  Temp  Press      Cp        sw_cp')
+  disp   ('  (psu)  (C)   (db)    (J/kg.C)   (J/kg.C)')
+  fprintf(1,' %4.0f  %4.0f   %5.0f   %8.1f  %11.2f\n', ...
+  [S(:,icol) T(:,icol) P(:,icol) UN_cp(:,icol) cp(:,icol)]');
+end %for
+
+%-------------------------------------------------------------------------------
+% TEST MAIN MODULE  
+%      SUB-MODULES  
+%------------------------------------------------------------------------------
+module     = 'sw_svel.m';
+disp(' ')
+disp('************************************************************************')
+disp(['**  TESTING MODULE: ' module])
+disp('************************************************************************')
+if display_help
+   eval(['help ' module])
+end %if
+
+% TEST 1 - 
+% DATA FROM POND AND PICKARD INTRO. DYNAMICAL OCEANOGRAPHY 2ND ED. 1986
+%***************************************************************************
+
+T    = [ 0  0  0  0  0  0;
+        10 10 10 10 10 10;
+	20 20 20 20 20 20;
+	30 30 30 30 30 30;
+	40 40 40 40 40 40 ];
+    
+S    = [25 25 25 35 35 35;
+        25 25 25 35 35 35;
+	25 25 25 35 35 35 ;
+	25 25 25 35 35 35;
+	25 25 25 35 35 35 ];
+
+P    = [0 5000 10000 0 5000 10000;
+        0 5000 10000 0 5000 10000;
+	0 5000 10000 0 5000 10000 ;
+	0 5000 10000 0 5000 10000;
+	0 5000 10000 0 5000 10000 ];
+
+UN_svel =      [1435.8  1520.4  1610.4  1449.1  1534.0  1623.2;
+                1477.7  1561.3  1647.4  1489.8  1573.4  1659.0;
+	        1510.3  1593.6  1676.8  1521.5  1604.5  1687.2;
+	        1535.2  1619.0  1700.6  1545.6  1629.0  1710.1;
+	        1553.4  1638.0  1719.2  1563.2  1647.3  1727.8 ];
+
+svel    = sw_svel(S,T,P);
+
+%----------------
+% DISPLAY RESULTS
+%----------------
+disp(' ')
+disp   ('********************************************************')
+disp   ('Comparison of accepted values from UNESCO 1983 ')
+disp   (' (Unesco Tech. Paper in Marine Sci. No. 44, p50)')
+disp   (['with computed results from ' module ' on ' computer ' computer'])
+disp   ('********************************************************')
+
+for icol = 1:length(S(1,:))
+  disp(' ')
+  disp   ('   Sal  Temp  Press     SVEL       sw_svel')
+  disp   ('  (psu)  (C)   (db)     (m/s)       (m/s)')
+  fprintf(1,' %4.0f  %4.0f   %5.0f   %8.1f  %11.3f\n', ...
+  [S(:,icol) T(:,icol) P(:,icol) UN_svel(:,icol) svel(:,icol)]');
+end %for
+
+%----------------------------------------------------------------------------
+% TEST MAIN MODULE  
+%      SUB-MODULES  
+%---------------------------------------------------------------------------
+submodules     = 'sw_alpha.m sw_beta.m sw_aonb.m';
+disp(' ')
+disp('**********************************************************************')
+disp(['**  TESTING MODULE: ' submodules])
+disp('**********************************************************************')
+if display_help
+   eval(['help ' submodules])
+end %if
+
+% DATA FROM MCDOUOGALL 1987
+s    = 40;
+ptmp = 10;
+p    = 4000;
+beta_lit  = 0.72088e-03;
+aonb_lit  = 0.34763;
+alpha_lit = aonb_lit*beta_lit;
+
+%$$$ % TEST ARGUMENT PASSING
+%$$$ beta = sw_beta(s,ptmp,p,'ptmp')
+%$$$ beta = sw_beta(s,ptmp,p,'temp')
+%$$$ beta = sw_beta(s,ptmp,p)
+%$$$ 
+%$$$ alpha = sw_alpha(s,ptmp,p,'ptmp')
+%$$$ alpha = sw_alpha(s,ptmp,p,'temp')
+%$$$ alpha = sw_alpha(s,ptmp,p)
+%$$$ 
+%$$$ aonb  = sw_aonb( s,ptmp,p,'ptmp')
+%$$$ aonb  = sw_aonb( s,ptmp,p,'temp')
+%$$$ aonb  = sw_aonb( s,ptmp,p)
+%$$$ 
+beta  = sw_beta( s,ptmp,p,'ptmp');
+alpha = sw_alpha(s,ptmp,p,'ptmp');
+aonb  = sw_aonb( s,ptmp,p,'ptmp');
+
+%----------------
+% DISPLAY RESULTS
+%----------------
+disp(' ')
+disp   ('********************************************************')
+disp   ('Comparison of accepted values from MCDOUGALL 1987 ')
+disp   (['with computed results on ' computer ' computer'])
+disp   ('********************************************************')
+
+  disp(' ')
+  disp   ('   Sal  Temp  Press     BETA       sw_beta')
+  disp   ('  (psu)  (C)   (db)   (psu^-1)     (psu^-1)')
+  fprintf(1,' %4.0f  %4.0f   %5.0f   %11.4e  %11.5e\n', ...
+  [s ptmp p beta_lit beta]');
+
+  disp(' ')
+  disp   ('   Sal  Temp  Press     AONB       sw_aonb')
+  disp   ('  (psu)  (C)   (db)   (psu C^-1)   (psu C^-1)')
+  fprintf(1,' %4.0f  %4.0f   %5.0f   %8.5f  %11.6f\n', ...
+  [s ptmp p aonb_lit aonb]');
+
+  disp(' ')
+  disp   ('   Sal  Temp  Press     ALPHA       sw_alpha')
+  disp   ('  (psu)  (C)   (db)    (psu^-1)     (psu^-1)')
+  fprintf(1,' %4.0f  %4.0f   %5.0f   %11.4e  %11.4e\n', ...
+  [s ptmp p alpha_lit alpha]');
+
+%--------------------------------
+% TEST MAIN MODULE  sw_satO2.m
+%      SUB-MODULES  
+%--------------------------------
+module     = 'sw_satO2 sw_satN2 sw_satAr';
+disp(' ')
+disp('*************************************')
+disp(['**  TESTING MODULE: ' module])
+disp(['**  and SUB-MODULE: ' submodules])
+disp('*************************************')
+if display_help
+   eval(['help ' module])
+end %if
+
+% Data from Weiss 1970
+
+T    = [-1 -1; 
+        10 10;
+	20 20 ;
+	40 40 ];
+    
+S    = [20 40;
+        20 40;
+	20 40 ;
+	20 40];
+
+lit_O2=  [ 9.162   7.984;
+           6.950   6.121;
+	   5.644   5.015;
+	   4.050   3.656];
+       
+lit_N2=  [16.28   14.01;
+          12.64   11.01;
+	  10.47    9.21;
+	   7.78    6.95];
+       
+lit_Ar=  [ 0.4456 0.3877;
+           0.3397 0.2989;
+	   0.2766 0.2457;
+	   0.1986 0.1794];       
+
+       
+satO2    = sw_satO2(S,T);
+satN2    = sw_satN2(S,T);
+satAr    = sw_satAr(S,T);
+
+%----------------
+% DISPLAY RESULTS
+%----------------
+disp(' ')
+disp   ('********************************************************')
+disp   ('Comparison of accepted values from Weiss, R.F. 1979 ')
+disp   ('"The solubility of nitrogen, oxygen and argon in water and seawater."')
+disp   (' Deap-Sea Research., 1970, Vol 17, pp721-735.')
+disp   (['with computed results from ' module ' on ' computer ' computer'])
+disp   ('********************************************************')
+
+for icol = 1:length(S(1,:))
+disp(' ')
+disp   ('   Sal  Temp      O2         sw_satO2')
+disp   ('  (psu)  (C)      (ml/l)     (ml/l)')
+  fprintf(1,' %4.0f  %4.0f    %8.2f   %9.3f\n', ...
+  [S(:,icol) T(:,icol)  lit_O2(:,icol) satO2(:,icol)]');
+end %for
+
+for icol = 1:length(S(1,:))
+disp(' ')
+disp   ('   Sal  Temp      N2         sw_satN2')
+disp   ('  (psu)  (C)      (ml/l)     (ml/l)')
+  fprintf(1,' %4.0f  %4.0f    %8.2f  %9.3f\n', ...
+  [S(:,icol) T(:,icol)  lit_N2(:,icol) satN2(:,icol)]');
+end %for
+
+for icol = 1:length(S(1,:))
+disp(' ')
+disp   ('   Sal  Temp      Ar         sw_satAr')
+disp   ('  (psu)  (C)      (ml/l)     (ml/l)')
+  fprintf(1,' %4.0f  %4.0f     %8.4f  %9.4f\n', ...
+  [S(:,icol) T(:,icol)  lit_Ar(:,icol) satAr(:,icol)]');
+end %for
+
+diary off
+return
+
+%--------------------------------------------------------------------
diff --git a/tsg_io/readTsgDataLabview.m b/tsg_io/readTsgDataLabview.m
index 4778fef782715c831babbc00d8ca51b3ce9e8f01..7a44159709e5a38cfb84de891f0c5dce1b076acb 100644
--- a/tsg_io/readTsgDataLabview.m
+++ b/tsg_io/readTsgDataLabview.m
@@ -165,6 +165,7 @@ if fid ~= -1
       % --------------------------------------------
       %tsg.SSPS_QC = tsg.qc.active.Code * ones(length(noNaN),1);
       tsg.SSPS_QC = castByteQC( tsg.qc.active.Code, noNaN );
+      tsg.SSJT_QC = castByteQC( tsg.qc.active.Code, noNaN );
 
       % populate tsg.file structure
       % ---------------------------
diff --git a/tsg_util/calibration.m b/tsg_util/calibration.m
new file mode 100644
index 0000000000000000000000000000000000000000..8ff58104e14e39de5448318a17665024da3c4278
--- /dev/null
+++ b/tsg_util/calibration.m
@@ -0,0 +1,39 @@
+function calibration( hMainFig )
+%
+% Compute salinity from calibrated conductivity and jacket temperature 
+%
+
+% Get tsg application data
+% ------------------------
+tsg = getappdata( hMainFig, 'tsg_data' );
+
+if ~isempty( tsg.CNDC )
+  cndcCal = tsg.CNDC_LINCOEF(1) * tsg.CNDC + tsg.CNDC_LINCOEF(2);
+else
+  msgbox( 'Conductivity not loaded',...
+          'Function ''Calibration''',...
+          'warn', 'modal');
+end
+
+if ~isempty( tsg.SSJT )
+  ssjtCal = tsg.SSJT_LINCOEF(1) * tsg.SSJT + tsg.SSJT_LINCOEF(2);
+else
+  msgbox( 'Jacket temperature not loaded',...
+          'Function ''Calibration''',...
+          'warn', 'modal');
+end
+
+% Compute salinity - Use CSIRO functions
+% --------------------------------------
+tsg.SSPS_CAL = sw_salt( ...
+               cndcCal/sw_c3515(), t90TOt68(ssjtCal), zeros(size(cndcCal)));
+            
+% Keep SSJT calibrated
+% --------------------
+tsg.SSJT_CAL = ssjtCal;
+
+% Save tsg application data
+% --------------------------
+setappdata( hMainFig, 'tsg_data', tsg );
+
+end
diff --git a/tsg_util/corTsgLinear.m b/tsg_util/corTsgLinear.m
index b5ce9e0ff1c3cd696d32bbf65e1fb58b5b246df0..f4667961d3613b5c342533f908c7ae8fdba65fa3 100644
--- a/tsg_util/corTsgLinear.m
+++ b/tsg_util/corTsgLinear.m
@@ -26,9 +26,9 @@ if dateMax > dateMin
   % intialisation
   % -------------
   if isempty( tsg.SSPS_ADJUSTED )
-    tsg.SSPS_ADJUSTED       = tsg.SSPS;
-    tsg.SSPS_ADJUSTED_ERROR = NaN * ones( size( tsg.SSPS ) );
-    tsg.SSPS_ADJUSTED_QC    = tsg.SSPS_QC;
+    msgbox( 'Variable SSPS_ADJUSTED should be initialise in updateTsgStruct',...
+            'Function ''corTsgLinear''',...
+            'warn', 'modal');
   end
  
   % Find samples within TIME_WINDOWS with Good and probably Good QC
diff --git a/tsg_util/corTsgMedian.m b/tsg_util/corTsgMedian.m
index 8245a2113f2e1a197193eee440d31ca5dd37273b..316feb0a36e326ce1a2efdde064aa1e511ebebb7 100644
--- a/tsg_util/corTsgMedian.m
+++ b/tsg_util/corTsgMedian.m
@@ -43,9 +43,9 @@ if dateMax > dateMin
   % intialisation
   % -------------
   if isempty( tsg.SSPS_ADJUSTED )
-    tsg.SSPS_ADJUSTED       = tsg.SSPS;
-    tsg.SSPS_ADJUSTED_ERROR = NaN * ones( size( tsg.SSPS ) );
-    tsg.SSPS_ADJUSTED_QC    = tsg.SSPS_QC;
+    msgbox( 'Variable SSPS_ADJUSTED should be initialise in updateTsgStruct',...
+            'Function ''corTsgMedian''',...
+            'warn', 'modal');
   end
   
   % Find the indices of samples within the time limits.
diff --git a/tsg_util/plot_Calibration.m b/tsg_util/plot_Calibration.m
index 2638031df4dbca3038b4237807a0473bd1b27a53..930574630d3da7b22f40a7a5305713fb504af0ac 100644
--- a/tsg_util/plot_Calibration.m
+++ b/tsg_util/plot_Calibration.m
@@ -7,7 +7,7 @@ tsg  = getappdata( hMainFig, 'tsg_data');
 switch nPlot
   
   % ---------------------------------------------------------------------
-  % Plot SSPS and SSPS_CAL
+  % Plot SSPS in black and SSPS_CAL in red
   case 1
     
     erase_Line( hPlotAxes, 1 );
@@ -18,23 +18,22 @@ switch nPlot
     
     if ~isempty( tsg.SSPS_CAL )
       plot_Tsg( hMainFig, hPlotAxes, 1, tsg.DAYD, tsg.SSPS_CAL, [],...
-                'SSPS_CAL','r','none','*',2);
+                'SSPS_CAL','r','none','o',2);
     end
 
   % ---------------------------------------------------------------------
-  % Plot SSPS in Black and SSPS with no position value in Red
+  % Plot SSJT in black and SSJT_CAL in red
   case 2
     
     erase_Line( hPlotAxes, 2 );
-    if ~isempty( tsg.SSPS )
-      plot_Tsg( hMainFig, hPlotAxes, 2, tsg.DAYD, tsg.SSPS,[],...
-                'SSPS','k','none','*',2);
+    if ~isempty( tsg.SSJT )
+      plot_Tsg( hMainFig, hPlotAxes, 2, tsg.DAYD, tsg.SSJT,[],...
+                'SSJT','k','none','*',2);
     end
     
-    ind = find( isnan(tsg.LATX) == 1 );
-    if ~isempty( ind )
-      plot_Tsg( hMainFig, hPlotAxes, 2, tsg.DAYD(ind), tsg.SSPS(ind),[],...
-        'SSPS_NaN','r','none','*',2);
+    if ~isempty( tsg.SSJT_CAL )
+      plot_Tsg( hMainFig, hPlotAxes, 2, tsg.DAYD, tsg.SSJT_CAL,[],...
+        'SSJT_CAL','r','none','*',2);
     end
 
     
diff --git a/tsg_util/t90TOt68.m b/tsg_util/t90TOt68.m
new file mode 100644
index 0000000000000000000000000000000000000000..15ade3f1d4d6e51867dc8e1aa963354c5293daa3
--- /dev/null
+++ b/tsg_util/t90TOt68.m
@@ -0,0 +1,7 @@
+function [t68] = t90TOt68( t90 )
+
+if ~isempty( t90 )
+  t68 = 1.00024 * t90;
+end
+
+end
diff --git a/tsg_util/updateAdjustedVariable.m b/tsg_util/updateAdjustedVariable.m
new file mode 100644
index 0000000000000000000000000000000000000000..63c1a7b736917af6d1dd64ffcb59c8d0c707745f
--- /dev/null
+++ b/tsg_util/updateAdjustedVariable.m
@@ -0,0 +1,59 @@
+function updateAdjustedVariable( hMainFig )
+%
+% Update adjusted SSPS and SSJT arrays once calibration has been
+% applied to these variable.
+%
+% The TRICK :  
+% Test ADJUSTED_ERROR - Only records corrected compared to water sample 
+%                       get error values
+%
+
+% Get tsg application data
+% ------------------------
+tsg = getappdata( hMainFig, 'tsg_data' );
+
+% Get VALUE_CHANGED
+% -------------------------------------------------------
+VALUE_CHANGED = get(tsg.qc.hash, 'VALUE_CHANGED', 'code');
+
+% SSPS
+% Get Adjusted records that were not corrected using Water samples
+% Only records corrected get error values
+% ----------------------------------------------------------------
+ind = find( isnan(tsg.SSPS_ADJUSTED_ERROR ) == 1);
+
+if ~isempty(tsg.SSPS_CAL) 
+  tsg.SSPS_ADJUSTED(ind)    = tsg.SSPS_CAL(ind);
+  tsg.SSPS_ADJUSTED_QC(ind) = VALUE_CHANGED;
+else
+  
+  % If the calibration has been canceled the ADJUSTED value is equal to
+  % the raw value
+  % -------------------------------------------------------------------
+  tsg.SSPS_ADJUSTED(ind)    = tsg.SSPS(ind);
+  tsg.SSPS_ADJUSTED_QC(ind) = tsg.SSPS_QC(ind);
+end
+
+% SSJT
+% Get Adjusted records that were not corrected using Water samples
+% Only records corrected get error values
+% ----------------------------------------------------------------
+ind = find( isnan(tsg.SSJT_ADJUSTED_ERROR ) == 1);
+
+if ~isempty(tsg.SSJT_CAL)
+  tsg.SSJT_ADJUSTED(ind)    = tsg.SSJT_CAL(ind);
+  tsg.SSJT_ADJUSTED_QC(ind) = VALUE_CHANGED;
+else
+  
+  % If the calibration has been canceled the ADJUSTED value is equal to
+  % the raw value
+  % -------------------------------------------------------------------
+  tsg.SSJT_ADJUSTED(ind)    = tsg.SSJT(ind);
+  tsg.SSJT_ADJUSTED_QC(ind) = tsg.SSJT_QC(ind);
+end
+
+% Save tsg application data
+% --------------------------
+setappdata( hMainFig, 'tsg_data', tsg );
+
+end
\ No newline at end of file
diff --git a/tsg_util/updateTsgStruct.m b/tsg_util/updateTsgStruct.m
index af236b07668f7cb12f843765d2a37e15c1c83b48..5a67ab54e1c73fa723296de6c0f97d92eb59660c 100644
--- a/tsg_util/updateTsgStruct.m
+++ b/tsg_util/updateTsgStruct.m
@@ -22,23 +22,56 @@ tsg.EAST_LONX  = max(tsg.LONX);
 
 % get date start and end value and set to globals attributes
 % -----------------------------------------------------------------
-date = datestr(min(tsg.DAYD),30);
+date           = datestr(min(tsg.DAYD),30);
 tsg.DATE_START = [date(1:8) date(10:15)];
-date = datestr(max(tsg.DAYD),30);
+date           = datestr(max(tsg.DAYD),30);
 tsg.DATE_END   = [date(1:8) date(10:15)];
 
 % Compute ship velocity from positions if sog not available
 % ---------------------------------------------------------
 if isempty(tsg.SPDC)
-  range = m_lldist(tsg.LONX,tsg.LATX);
-  ind = size(tsg.DAYD);
+  range    = m_lldist(tsg.LONX,tsg.LATX);
+  ind      = size(tsg.DAYD);
   tsg.SPDC = zeros(size(ind));
+  
   for i=1:length(tsg.DAYD)-1
     tsg.SPDC(i) = range(i) / ((tsg.DAYD(i+1)-tsg.DAYD(i)) * 24 * 1.854);
   end
   tsg.SPDC = [tsg.SPDC';0];
 end
 
+% Initialise ADJUSTED variables if empty
+% --------------------------------------
+if isempty( tsg.SSPS_ADJUSTED ) && ~isempty( tsg.SSPS )
+  tsg.SSPS_ADJUSTED       = tsg.SSPS;
+  tsg.SSPS_ADJUSTED_QC    = tsg.SSPS_QC;
+  tsg.SSPS_ADJUSTED_ERROR = NaN * ones( size( tsg.SSPS ));
+
+elseif isempty( tsg.SSPS )
+
+  msgbox('You must initialise the tsg.SSPS variable',...
+    'function ''updateTsgStruct''', ....
+    'warn', 'modal');
+end
+
+if isempty( tsg.SSJT_ADJUSTED ) && ~isempty( tsg.SSJT )
+  tsg.SSJT_ADJUSTED       = tsg.SSJT;
+  tsg.SSJT_ADJUSTED_QC    = tsg.SSJT_QC;
+  tsg.SSJT_ADJUSTED_ERROR = NaN * ones( size( tsg.SSJT ));
+
+elseif isempty( tsg.SSJT )
+
+  msgbox('You must initialise the tsg.SSJT variable',...
+    'function ''updateTsgStruct''', ....
+    'warn', 'modal');
+end
+
+if isempty( tsg.SSPS_ADJUSTED ) && ~isempty( tsg.SSPS )
+  tsg.SSPT_ADJUSTED       = tsg.SSPT;
+  tsg.SSPT_ADJUSTED_QC    = tsg.SSPT_QC;
+  tsg.SSPT_ADJUSTED_ERROR = NaN * ones( size( tsg.SSPS ));
+end
+  
 % Save tsg structure 
 % ------------------
 setappdata( hTsgGUI, 'tsg_data', tsg);
diff --git a/tsgcor_GUI.m b/tsgcor_GUI.m
deleted file mode 100644
index 6220fb15c97f5e8ac3cb2011b6034f54d9a44efb..0000000000000000000000000000000000000000
--- a/tsgcor_GUI.m
+++ /dev/null
@@ -1,665 +0,0 @@
-function tsgcor_GUI( hTsgGUI )
-% tsgcor_GUI 
-%
-% GUI for correction of TSG data by comparison to samples
-% this GUI is a children of tsgqc_GUI
-%
-%
-
-%% COPYRIGHT & LICENSE
-%  Copyright 2007 - IRD US191, all rights reserved.
-%
-%  This file is part of tsgqc_GUI.
-%
-%    Datagui is free software; you can redistribute it and/or modify
-%    it under the terms of the GNU General Public License as published by
-%    the Free Software Foundation; either version 2 of the License, or
-%    (at your option) any later version.
-%
-%    tsgqc_GUI is distributed in the hope that it will be useful,
-%    but WITHOUT ANY WARRANTY; without even the implied warranty of
-%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-%    GNU General Public License for more details.
-%
-%    You should have received a copy of the GNU General Public License
-%    along with Datagui; if not, write to the Free Software
-%    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
-
-
-%%  Initialization tasks
-%   ********************
-clc
-
-% Screen limits for the GUI
-% -------------------------
-set(0,'Units','normalized');
-guiLimits = get(0,'ScreenSize');
-guiLimits(1) = guiLimits(1) + 0.05;
-guiLimits(2) = guiLimits(2) + 0.05;
-guiLimits(3) = guiLimits(3) - 0.1;
-guiLimits(4) = guiLimits(4) - 0.2;
-
-% Create and then hide the GUI as it is being constructed.
-% --------------------------------------------------------
-hTsgCorGUI = figure(...
-  'Name', 'TSG Correction', ...
-  'NumberTitle', 'off', ...
-  'Resize', 'on', ...
-  'Menubar','none', ...
-  'Toolbar', 'none', ...
-  'Tag', 'ButtonMotionOff', ...
-  'WindowButtonMotionFcn', @MouseMotion, ...
-  'CloseRequestFcn', @QuitProgram,...
-  'HandleVisibility','callback',...
-  'Visible','on',...
-  'Units', 'normalized',...
-  'Position',guiLimits);
- 
-% Get the CONSTANTS
-% -----------------
-cst    = getappdata( hTsgGUI, 'constante');
-
-%  Construct the Menu
-%   -----------------
-hFileMenu = uimenu(...
-  'Parent', hTsgCorGUI,...
-  'HandleVisibility','callback',...
-  'Label', 'File');
-
-hQuitMenu = uimenu(...
-  'Parent',hFileMenu,...
-  'Label','Quit',...
-  'Separator','on',...
-  'Accelerator','Q',...
-  'HandleVisibility','callback',...
-  'Callback',@QuitMenuCallback);
-
-%%  Construct the Toolbar
-%   --------------------
-hToolbar       =   uitoolbar(...   % Toolbar for Open and Print buttons
-  'Parent',hTsgCorGUI, ...
-  'HandleVisibility','callback');
-hZoomToggletool  =   uitoggletool(...   % Open Zoom toolbar button
-  'Parent',hToolbar,...
-  'Separator', 'on', ...
-  'TooltipString','Zoom',...
-  'CData', iconRead(fullfile(matlabroot, ...
-  '/toolbox/matlab/icons/zoom.mat')),...
-  'HandleVisibility','callback', ...
-  'Tag','PUSHTOOL_ZOOM',...
-  'Enable', 'off',...
-  'OffCallback', @Zoom_OffMenuCallback,...
-  'ONCallback',  @Zoom_OnMenuCallback);
-hPanToggletool  =   uitoggletool(...   % Open Pan toolbar button
-  'Parent',hToolbar,...
-  'TooltipString','Pan',...
-  'CData',iconRead(fullfile(matlabroot, ...
-  '/toolbox/matlab/icons/pan.mat')),...
-  'HandleVisibility','callback', ...
-  'Tag','PUSHTOOL_PAN',...
-  'Enable', 'off',...
-  'OffCallback', @Pan_OffMenuCallback,...
-  'ONCallback',  @Pan_OnMenuCallback);
-  %'ClickedCallback', @PanMenuCallback);
-
-hBottlePushtool  = uipushtool(...   % Open toolbar button
-  'Parent',hToolbar,...
-  'TooltipString','Plot the Samples',...
-  'Separator', 'on', ...
-  'Tag', 'off', ...
-  'CData',iconRead(...
-  [DEFAULT_PATH_FILE 'tsg_icon' filesep 'bottleicon.mat']),...
-  'HandleVisibility','callback', ...
-  'ClickedCallback', @BottleMenuCallback);
-hStartlimitPushtool  = uipushtool(...   % Open toolbar button
-  'Parent',hToolbar,...
-  'TooltipString','Start time',...
-  'Separator', 'on', ...
-  'Tag', 'off', ...
-  'CData',iconRead(...
-  [DEFAULT_PATH_FILE 'tsg_icon' filesep 'startlimit.mat']),...
-  'HandleVisibility','callback', ...
-  'ClickedCallback', @TimeLimitCallback);
-hEndlimitPushtool  = uipushtool(...   % Open toolbar button
-  'Parent',hToolbar,...
-  'TooltipString','End time',...
-  'Tag', 'off', ...
-  'CData',iconRead(...
-  [DEFAULT_PATH_FILE 'tsg_icon' filesep 'endlimit.mat']),...
-  'HandleVisibility','callback', ...
-  'ClickedCallback', @TimeLimitCallback);
-
-% Static text that displays the position, salinity and temperature
-% ----------------------------------------------------------------
-hInfoText = uicontrol(...
-  'Parent', hTsgCorGUI, ...
-  'Style', 'text', ...
-  'Fontsize', 12, ...
-  'Fontweight', 'bold', ...
-  'Visible','on',...
-  'Units', 'normalized',...
-  'String', 'Information sur la position du curseur', ...
-  'Position', [.05, .95, .9, .03]);
-
-% Construct the plot axes
-% -----------------------
-hPlotAxes(1) = axes(...     % the axes for plotting Salinity
-  'Parent', hTsgCorGUI, ...
-  'Units', 'normalized', ...
-  'Visible', 'off', ...
-  'HandleVisibility','callback', ...
-  'Position',[.25, .6, .7, .32]);
-hPlotAxes(2) = axes(...     % the axes for plotting sample
-  'Parent', hTsgCorGUI, ...
-  'Units', 'normalized', ...
-  'Visible', 'off', ...
-  'HandleVisibility','callback', ...
-  'Position',[.25, .3, .7, .25]);
-hPlotAxes(3) = axes(...     % the axes for plotting ship velocity
-                'Parent', hTsgCorGUI, ...
-                'Units', 'normalized', ...
-                'Visible', 'off', ...
-                'HandleVisibility','callback', ...
-                'Position',[.25, .05, .7, .2]); 
-
-    % Choose the date limits for the correction
-    % --------------------------------------------------
-    % Create the uipanel
-    hpDateLimit = uipanel( ...
-                'Parent', hTsgCorGUI, ...
-                'Title', 'Date Limits', ...
-                'Units', 'normalized', ...
-                'FontSize', 11, ...
-                'Fontweight', 'bold', ...
-                'BackgroundColor', 'white',...
-                'Position',[.01 .75 .2 .19]);
-
-    htDateMin = uicontrol( ...
-                'Parent', hpDateLimit, ...
-                'Style', 'Text', ...
-                'String', 'Minimum  (yyyy-mm-dd hh:mm:ss)', ...
-                'HorizontalAlignment', 'left', ...
-                'Units', 'normalized', ...
-                'backgroundcolor', 'white', ...
-                'FontSize', 10, ...
-                'Position',[.05 .8 .9 .15]);
-
-    hetDateMin = uicontrol( ...
-                'Parent', hpDateLimit, ...
-                'Style', 'edit', ...
-                'Units', 'normalized', ...
-                'FontSize', 10, ...
-                'Position',[.05 .6 .9 .18]);
-
-    htDateMax = uicontrol( ...
-                'Parent', hpDateLimit, ...
-                'Style', 'Text', ...
-                'String', 'Maximum  (yyyy-mm-dd hh:mm:ss)', ...
-                'HorizontalAlignment', 'left', ...
-                'Units', 'normalized', ...
-                'backgroundcolor', 'white', ...
-                'FontSize', 10, ...
-                'Position',[.05 .3 .9 .15]);
-
-    hetDateMax = uicontrol( ...
-                'Parent', hpDateLimit, ...
-                'Style', 'edit', ...
-                'Units', 'normalized', ...
-                'FontSize', 10, ...
-                'Position',[.05 .1 .9 .18]);
-                      
-            
-    % Choose the correction method
-    % --------------------------------------------------
-    % Create the button group
-    hbgMethod = uibuttongroup( ...
-                'Parent',hTsgCorGUI, ...
-                'Title','Correction Method', ...
-                'Units', 'normalized', ...
-                'FontSize',11, ...
-                'Fontweight', 'bold', ...
-                'BackgroundColor','white',...
-                'Position',[.01 .6 .2 .14]);
-             
-    % Create 2 radio buttons in the button group
-    hrbLinear = uicontrol( ...
-                'Style','pushbutton', ...
-                'Parent',hbgMethod, ...
-                'Units', 'normalized', ...
-                'String','Linear adjustment',...
-                'pos',[.05 .55 .9 .4], ...
-                'HandleVisibility','callback', ...
-                'Callback', @CorLinearCallback);
-
-    hrbMedian = uicontrol( ...
-                'Style','pushbutton', ...
-                'Parent',hbgMethod, ...
-                'Units', 'normalized', ...
-                'String','Running median filter',...
-                'pos',[.05 .05 .9 .4],...
-                'HandleVisibility','callback', ...
-                'Callback', @CorMedianCallback);
-
-               
-    % Initialize some button group properties 
-    set(hbgMethod,'SelectionChangeFcn',@selcbk);
-    set(hbgMethod,'SelectedObject',[]);  % No selection
-    set(hbgMethod,'Visible','on');
-            
-    % Construct the context menu to delete samples
-    % --------------------------------------------
-    hSampleCmenu       = uicontextmenu(...
-                            'Parent', hTsgCorGUI, ...
-                            'HandleVisibility','callback' );
-    hQcSampleNocontrol = uimenu(...
-                            'Parent', hSampleCmenu,... 
-                            'HandleVisibility','off', ...
-                            'Label', 'No control',...
-                            'ForegroundColor', 'k',...
-                            'Callback', @Qc);
-    hQcSampleGood      = uimenu(...
-                            'Parent', hSampleCmenu,... 
-                            'HandleVisibility','off', ...
-                            'Label', 'Good',...
-                            'ForegroundColor', 'b',...
-                            'Callback', @Qc);
-    hQcSampleBad       = uimenu(...
-                            'Parent', hSampleCmenu,... 
-                            'HandleVisibility','off', ...
-                            'Label', 'Bad',...
-                            'ForegroundColor', 'r',...
-                            'Callback', @Qc);
-
-                          
-                          
-                          
-% Pointer set to watch during reading and plotting
-% ------------------------------------------------
-set( hTsgCorGUI, 'Pointer', 'watch' );
-
-% Get the data useful to the correction GUI
-% -----------------------------------------
-tsg    = getappdata( hTsgGUI, 'tsg_data' );
-
-% Merge bucket and external samples
-% ---------------------------------
-tsg_mergesample( hTsgGUI );
-
-
-
-dateMin = datestr(tsg.TIME(1), 31);
-dateMax = datestr(tsg.TIME(end), 31);
-set( hetDateMin, 'String', dateMin);
-set( hetDateMax, 'String', dateMax); 
-
-% dt between 2 tsg measurements
-TSG_SAMPLING_TIME = tsg.TIME(2) - tsg.TIME(1);
-
-% Running average of TSG time series over TSG_DT_SMOOTH hour.
-[psal_smooth, nval] = ...
-      dev_moveaverage(tsg.TIME, tsg.PSAL, cst.TSG_DT_SMOOTH, cst.TSG_STDMAX);
-
-% Compute the sample-TSG differences
-sample = dev_diffTsgSample(tsg, psal_smooth, sample, TSG_SAMPLING_TIME);
-
-% Tracé
-% -----
-tsg_plot_SalTsgSample( hTsgGUI, hPlotAxes );
-
-% Pointer reset to arrow
-% ----------------------
-set( hTsgCorGUI, 'Pointer', 'arrow' );
-
-% The callback to detect the mouse motion can be set to on
-% --------------------------------------------------------
-set( hTsgCorGUI, 'Tag', 'ButtonMotionOn');
-    
-%  Callbacks for tsgcqc_GUI
-%  ************************
-
-    %----------------------------------------------------------------------
-    function ZoomMenuCallback(hObject, eventdata)
-    % Callback function run when the ...    
-        
-        % Returns a zoom mode object for the figure handle
-        % ------------------------------------------------
-        hZoom = zoom(hTsgCorGUI);
-        
-        % Specifies whether this mode is currently enabled on the figure
-        % --------------------------------------------------------------
-        zoomOnOff = get(hZoom, 'Enable' );
-        switch zoomOnOff
-            case 'on'
-                zoom off
-                zoomAdaptiveDateTicks('off');
-            case 'off'
-                pan off
-                set(hTsgCorGUI,'Pointer','arrow');
-                set(hEndlimitPushtool, 'Tag', 'off' );
-                set(hStartlimitPushtool, 'Tag', 'off' );
-
-                zoom on
-                zoomAdaptiveDateTicks('on');
-        end
-    end
-
-    %----------------------------------------------------------------------
-    function PanMenuCallback(hObject, eventdata)
-    % Callback function run when the ....    
-        
-        % Returns a pan mode object for the figure handle
-        % -----------------------------------------------
-        hPan = pan(hTsgCorGUI);
-        
-        % Specifies whether this mode is currently enabled on the figure
-        % --------------------------------------------------------------
-        panOnOff = get(hPan, 'Enable' );
-        switch panOnOff
-            case 'on'
-                pan off
-                panAdaptiveDateTicks('off');
-            case 'off'
-                zoom off
-                set(hTsgCorGUI,'Pointer','arrow');
-                set(hEndlimitPushtool, 'Tag', 'off' );
-                set(hStartlimitPushtool, 'Tag', 'off' );
-
-                pan on
-                panAdaptiveDateTicks('on');                
-        end
-    end
-
-    %----------------------------------------------------------------------
-    function TimeLimitCallback(hObject, eventdata)
-    % Callback function run when the ....    
-
-        % Desactivate the Zoom and Pan functions.
-        % ---------------------------------------
-        zoom off; pan off
-        panAdaptiveDateTicks('off');zoomAdaptiveDateTicks('off');
-
-        % Retrieve named application data
-        % -------------------------------
-        tsg = getappdata( hTsgGUI, 'tsg_data');
-        
-        % Toggle the tag of the Qc pushbutton to 'on' or 'off'
-        % ----------------------------------------------------
-        switch get(hObject, 'Tag'); 
-            case 'off'
-                set(hObject, 'Tag', 'on' );
-                set( hTsgCorGUI, 'Pointer', 'crosshair');
-            case 'on'
-                set(hObject, 'Tag', 'off' );
-                set(hTsgCorGUI,'Pointer','arrow');
-        end
-
-        % Wait for key press
-        % ------------------
-        [x, y] = ginput(1);
-        
-        % Write the date in the Editable uicontrol
-        % ----------------------------------------
-        if hObject == hEndlimitPushtool
-            set( hetDateMax, 'String', datestr(x, 31));
-        else
-            set( hetDateMin, 'String', datestr(x, 31));
-        end
-        
-        set(hTsgCorGUI,'Pointer','arrow');
-
-    end
-
-%----------------------------------------------------------------------
-  function BottleMenuCallback(gcbo, eventdata,handles)
-    % Callback function run when the QC pushbutton is selected
-    
-    % Desactivate the Zoom and Pan functions.
-    % ---------------------------------------
-    zoom off; pan off
-    panAdaptiveDateTicks('off');zoomAdaptiveDateTicks('off');
-    
-    % Toggle the tag of the Qc pushbutton to 'on' or 'off'
-    % ----------------------------------------------------
-    switch get(hBottlePushtool, 'Tag');
-      case 'off'
-        set(hBottlePushtool, 'Tag', 'on' );
-        set(hPlotAxes(2),'UIContextMenu', hSampleCmenu);
-        set(hTsgCorGUI, 'Pointer', 'crosshair');
-      case 'on'
-        set(hBottlePushtool, 'Tag', 'off' );
-        set(hPlotAxes(2),'UIContextMenu', []);
-        set(hTsgCorGUI,'Pointer','arrow');
-    end
-
-    qualityCode = -1;
-    ind = [];
-    while strcmp( get(hBottlePushtool, 'Tag'),'on')
-      
-      k = waitforbuttonpress;
-      
-      % If the QC pushbutton is pressed we quit the callback
-      % ----------------------------------------------------
-      if strcmp( get(hBottlePushtool, 'Tag'),'off')
-        
-        % Desactivate the context menu use to choose the Quality Codes
-        % ----------------------------------------------
-        set(hSampleCmenu, 'Visible', 'off');
-        break
-      end
-
-      % Test if the right mouse button is clicked
-      % -----------------------------------------
-      if strcmp(get(hTsgCorGUI, 'SelectionType'), 'alt') && ~isempty(ind)
-        
-        % Wait for a QC Context Menu choice : The user choose the
-        % quality code
-        % -------------------------------------------------------
-        uiwait
-        qualityCode = 1;
-      
-      else
-
-        % Mouse motion callback desactivated when a selection is made.
-        % Otherwise there is a conflict with the map if it is activated
-        % -------------------------------------------------------
-        set( hTsgCorGUI, 'Tag', 'ButtonMotionOff');
-        
-        % Selection of the data within the figure
-        % ---------------------------------------
-        point1    = get(gca,'CurrentPoint');    % button down detected
-        finalRect = rbbox;                      % return figure units
-        point2    = get(gca,'CurrentPoint');    % button up detected
-        
-        point1    = point1(1,1:2);              % extract x and y
-        point2    = point2(1,1:2);
-        
-        p1 = min(point1,point2);
-        p2 = max(point1,point2);             % calculate locations
-        
-        % Retrieve named application data
-        % -------------------------------
-        sample = getappdata( hTsgGUI, 'sample');
-        
-        ind = find(sample.DAYD > p1(1,1) & sample.DAYD < p2(1,1) & ...
-                   sample.SSPS_DIF > p1(1,2) & sample.SSPS_DIF < p2(1,2));
-                        
-        % Selection made : Mouse motion callback re-activated
-        % ---------------------------------------------------
-        set( hTsgCorGUI, 'Tag', 'ButtonMotionOn');
-      
-      end
-
-      % Plot the data with the color of the chosen quality Code.
-      % Is it the right place for this source code ????
-      % --------------------------------------------------------
-      if qualityCode ~= -1
-        
-        tsg = getappdata( hTsgGUI, 'tsg_data');
-        
-        sample.SSPS_QC(ind) =  tsg.qc.Code.ACTIVE;
-        
-        % Save the modifications
-        % ----------------------
-        setappdata( hTsgGUI, 'sample', sample);
-        
-        axes(hPlotAxes(2));
-        hold on
-        color = ['o' tsg.qc.Color.ACTIVE];
-        plot(sample.DAYD(ind), sample.SSPS_DIF(ind), color );
-        hold off
-      end
-    end
-  end
-
-%---------------------------------------------------------------------
-function Qc(hObject, eventdata)
-  % Callback function run when the QC context menu is selected
-  
-  % Retrieve Default Quality Code and Color
-  % ---------------------------------------
-  tsg = getappdata( hTsgGUI, 'qcColor');
-  
-  switch hObject
-    case hQcSampleNocontrol
-      tsg.qc.Code.ACTIVE   = tsg.qc.Code.NO_CONTROL;
-      tsg.qc.Color.ACTIVE  = tsg.qc.Color.NO_CONTROL;
-    case hQcSampleGood
-      tsg.qc.Code.ACTIVE   = tsg.qc.Code.GOOD;
-      tsg.qc.Color.ACTIVE  = tsg.qc.Color.GOOD;
-    case hQcSampleBad
-      tsg.qc.Code.ACTIVE   = tsg.qc.Code.PROBABLY_BAD;
-      tsg.qc.Color.ACTIVE  = tsg.qc.Color.PROBABLY_BAD;
-  end
-
-  setappdata( hMainFig, 'tsg_data', tsg );
-  
-  % uiwait in the QCMenuCallback function
-  % -------------------------------------
-  uiresume
-end
-
-    %---------------------------------------------------------------------
-    function MouseMotion(hObject, eventdata)
-     
-      % Test if the callback can be activated
-      % -------------------------------------
-      if strcmp( get( hTsgCorGUI, 'Tag'), 'ButtonMotionOn')
-          
-        % Retrieve named application data
-        % -------------------------------
-        tsg = getappdata( hTsgGUI, 'tsg_data');
-        
-        % Get the mouse position
-        % ----------------------
-        point = get(gcf,'CurrentPoint');
-        
-        if point(1) > .05 && point(2) > .6 && point(1) < .95 && point(2) < .92
-    
-            % Get current position of cusor and return its coordinates in
-            % axes with handle h_axes
-            % -----------------------------------------------------------
-            [x, y] = gpos(hPlotAxes(1));
-                    
-            if x > tsg.TIME(1) && x < tsg.TIME(end)
-            
-                indCursor = find( tsg.TIME > x);
-                % use sprintf with format instead strcat & num2str but flag
-                % - don't work with 0, eg %+07.4f
-                set( hInfoText, 'String',...
-                  sprintf(['%s   -   Latitude = %s   -   Longitude = %s '...
-                           '  -   Salinity = %07.4f   -   Temperature = %07.4f'],...
-                    datestr(tsg.TIME(indCursor(1)),'dd/mm/yyyy HH:MM'),...
-                    dd2dm(tsg.LATITUDE(indCursor(1)),0), ...
-                    dd2dm(tsg.LONGITUDE(indCursor(1)),1), ...
-                    tsg.PSAL(indCursor(1)), ...
-                    tsg.TEMP_TSG(indCursor(1))...
-                ));
-            
-            end
-        end
-      end
-    end
-
-    % -----------------------------------------------------------------
-    function SaveMenuCallback(hObject, eventdata)
-    % Callback function run when the Save menu item is selected
-        
-%        [fileName, pathName, filterIndex] = uiputfile('*.txt', ...
-%                                            'Save file name');
-        
-%        fileName = [pathName fileName];
-%        error = tsg_writeTsgData( hTsgGUI, fileName );
-%        if ~error
-            % 
-%        end
-    end
-
-    % -----------------------------------------------------------------
-    function CorLinearCallback(hObject, eventdata)
-    % Callback function run when 
-    
-        msgbox('Method not yet implemented', 'modal' );
-    
-    end
-
-    % -----------------------------------------------------------------
-    function CorMedianCallback(hObject, eventdata)
-    % Callback function run when 
-    
-        msgbox('Method not yet implemented', 'modal' );
-        
-        % Get the time limits for the correction A TESTER
-        % --------------------------------------
-        %dateMin = get( hetDateMin, 'String');
-        %dateMax = get( hetDateMax, 'String'); 
-
-        % Correction
-        % ----------
-        % tsg = dev_corMethod1(tsg, sample, dateMin, dateMax, cst.COR_TIME_WINDOWS);
-        
-        % Update application data 'tsg'
-        % ----------------------------
-        %tsg    = getappdata( hTsgGUI, 'tsg_data' );
-   
-    end
-
-
-    % -----------------------------------------------------------------
-    function QuitMenuCallback(hObject, eventdata)
-    % Callback function run when the Quit menu item is selected
-            
-        % If the data have been modified and not save, the program
-        % propose to save the data
-        % --------------------------------------------------------
-        if  strcmp( get( hSaveMenu, 'Tag' ), 'on')
-            selection = ...
-            questdlg('The file has been modified.  Do you want to save it ?',...
-                      'Save before Quit?',...
-                      'Yes', 'No', 'Yes');
-            if strcmp(selection, 'Yes')
-                return;
-            else
-                QuitProgram;
-            end
-        else
-            selection = ...
-                    questdlg(['Quit ' get(hTsgCorGUI, 'Name') '?'],...
-                             ['Quit ' get(hTsgCorGUI, 'Name') '?'],...
-                              'Yes', 'No', 'Yes');
-            if strcmp(selection, 'No')
-                return;
-            else    
-                QuitProgram;
-            end
-        end
-
-    end
-  
-  % ----------------------------------------------------------------
-  function QuitProgram(hObject, eventdata)
-
-    delete(hTsgCorGUI);
-
-  end
-      
-end
diff --git a/tsgqc_GUI.m b/tsgqc_GUI.m
index ff48a8ed009e4c9a03bdf96cb11eeee5816d9d04..ed4e03ac11af35ca83e877bf96577a2ab50e57b1 100644
--- a/tsgqc_GUI.m
+++ b/tsgqc_GUI.m
@@ -753,94 +753,87 @@ hetDateMax = uicontrol( ...
 hpCalCoef = uipanel( ...
   'Parent', hMainFig, ...
   'Title', 'Calibration', ...
-  'Units', 'normalized', ...
   'FontSize', tsg.fontSize-1, 'Fontweight', 'bold', ...
   'Visible', 'off', ...
-  'Position', [.0, .68, .15, .28]);
+  'Units', 'normalized','Position', [.0, .64, .15, .32]);
 
  htCalCNDC1 = uicontrol( ...
    'Parent', hpCalCoef, ...
-   'Style', 'Text', 'String', 'Conductivity : C = A*C +B', ...
-   'Units', 'normalized', ...
+   'Style', 'Text', 'String', 'Conductivity : C = A*C + B', ...
    'HorizontalAlignment', 'left', 'FontSize', tsg.fontSize-1, ...
-   'Position',[.01 .85 .95 .1]);
+   'Units', 'normalized', 'Position',[.01 .89 .95 .1]);
  htCalCNDC2 = uicontrol( ...
    'Parent', hpCalCoef, ...
    'Style', 'Text', 'String', 'A', ...
-   'Units', 'normalized', ...
    'HorizontalAlignment', 'left', 'FontSize', tsg.fontSize-1, ...
-   'Position',[.01 .75 .08 .1]);
+   'Units', 'normalized', 'Position',[.01 .8 .08 .1]);
  hetCalCNDCSlope = uicontrol( ...
    'Parent', hpCalCoef, ...
    'Style', 'edit', ...
-   'Units', 'normalized', ...
-   'BackgroundColor', 'white',...
-   'FontSize', tsg.fontSize, ...
+   'FontSize', tsg.fontSize, 'BackgroundColor', 'white',...
    'Tag', 'CORRECT_CAL_CNDC_A',...
    'HandleVisibility','on', ...
-   'Position',[.1 .75 .85 .10]);
+   'Units', 'normalized', 'Position',[.1 .8 .85 .10]);
  htCalCNDC3 = uicontrol( ...
    'Parent', hpCalCoef, ...
    'Style', 'Text', 'String', 'B', ...
-   'Units', 'normalized', ...
    'HorizontalAlignment', 'left',  'FontSize', tsg.fontSize-1, ...
-   'Position',[.01 .6 .08 .1]);
+   'Units', 'normalized', 'Position',[.01 .68 .08 .1]);
  hetCalCNDCOffset = uicontrol( ...
    'Parent', hpCalCoef, ...
    'Style', 'edit', ...
-   'Units', 'normalized', ...
-   'BackgroundColor', 'white',...
-   'FontSize', tsg.fontSize, ...
+   'FontSize', tsg.fontSize, 'BackgroundColor', 'white',...
    'HandleVisibility','on', ...
    'Tag', 'CORRECT_CAL_CNDC_B',...
-   'Position',[.1 .6 .85 .10]);
+   'Units', 'normalized', 'Position',[.1 .68 .85 .10]);
 
  htCalSSJT1 = uicontrol( ...
    'Parent', hpCalCoef, ...
-   'Style', 'Text', 'String', 'SSJT : T = A*T +B', ...
-   'Units', 'normalized', ...
+   'Style', 'Text', 'String', 'SSJT : T = A*T + B', ...
    'HorizontalAlignment', 'left', 'FontSize', tsg.fontSize-1, ...
-   'Position',[.01 .45 .95 .1]);
+   'Units', 'normalized', 'Position',[.01 .55 .95 .1]);
  htCalSSJT2 = uicontrol( ...
    'Parent', hpCalCoef, ...
    'Style', 'Text', 'String', 'A', ...
-   'Units', 'normalized', ...
    'HorizontalAlignment', 'left', 'FontSize', tsg.fontSize-1, ...
-   'Position',[.01 .35 .08 .1]);
+   'Units', 'normalized', 'Position',[.01 .45 .08 .1]);
  hetCalSSJTSlope = uicontrol( ...
    'Parent', hpCalCoef, ...
    'Style', 'edit', ...
-   'Units', 'normalized', ...
    'BackgroundColor', 'white',...
    'FontSize', tsg.fontSize, ...
    'HandleVisibility','on', ...
    'Tag', 'CORRECT_CAL_SSJT_A',...
-   'Position',[.1 .35 .85 .10]);
+   'Units', 'normalized', 'Position',[.1 .45 .85 .10]);
  htCalSSJT3 = uicontrol( ...
    'Parent', hpCalCoef, ...
    'Style', 'Text', 'String', 'B', ...
-   'Units', 'normalized', ...
    'HorizontalAlignment', 'left', 'FontSize', tsg.fontSize-1, ...
-   'Position',[.01 .2 .08 .1]);
+   'Units', 'normalized', 'Position',[.01 .32 .08 .1]);
  hetCalSSJTOffset = uicontrol( ...
    'Parent', hpCalCoef, ...
    'Style', 'edit', ...
-   'Units', 'normalized', ...
    'FontSize', tsg.fontSize, 'BackgroundColor', 'white',...
    'HandleVisibility','on', ...
    'Tag', 'CORRECT_CAL_SSJT_B',...
-   'Position',[.1 .2 .85 .10]);
+   'Units', 'normalized', 'Position',[.1 .32 .85 .10]);
 
  hrbCal = uicontrol( ...
-  'Style','pushbutton', ...
-  'Parent',hpCalCoef, ...
-  'Units', 'normalized', ...
+  'Style','pushbutton', 'Parent',hpCalCoef, ...
   'String','Calibrate',...
   'FontSize',tsg.fontSize-1,...
   'Tag', 'CORRECT_CAL_PUSH', ...
-  'pos',[.05 .04 .9 .13], ...
+  'Units', 'normalized','pos',[.05 .17 .9 .1], ...
   'HandleVisibility','callback', ...
   'Callback', @CalibrateCallback);
+ hrbCancelCal = uicontrol( ...
+  'Style','pushbutton', 'Parent',hpCalCoef, ...
+  'String','Cancel calibration',...
+  'FontSize',tsg.fontSize-1,...
+  'Tag', 'CORRECT_CAL_PUSH', ...
+  'Units', 'normalized','pos',[.05 .04 .9 .1], ...
+  'HandleVisibility','callback', ...
+  'Callback', @CancelCalibrationCallback);
 
 
 
@@ -1151,6 +1144,73 @@ end
    
   end
 
+%% CalibrateCallback .......................................... Calibration
+  function CalibrateCallback(hObject, eventdata)
+    % Callback function run when
+    
+    % Get tsg application data
+    % ------------------------
+    tsg = getappdata( hMainFig, 'tsg_data' );
+    
+    % Get the calibration coefficients
+    % --------------------------------
+    tsg.CNDC_LINCOEF(1) = str2num(get( hetCalCNDCSlope, 'String'));
+    tsg.CNDC_LINCOEF(2) = str2num(get( hetCalCNDCOffset, 'String'));
+    tsg.SSJT_LINCOEF(1) = str2num(get( hetCalSSJTSlope, 'String'));
+    tsg.SSJT_LINCOEF(2) = str2num(get( hetCalSSJTOffset, 'String'));
+        
+    % Save tsg application data
+    % --------------------------
+    setappdata( hMainFig, 'tsg_data', tsg );
+
+    % Calibrate the sensors
+    % ---------------------
+    calibration( hMainFig );
+    
+    % Update the Adjusted variables (SSPS - SSJT) with calibrated records
+    % -------------------------------------------------------------------
+    updateAdjustedVariable( hMainFig );
+    
+    % Refresh plot #1
+    % ---------------
+    plot_Calibration( hMainFig, hPlotAxes, 1 );
+    plot_Calibration( hMainFig, hPlotAxes, 2 );
+    
+    % As soon as a modification took place the data should be saved
+    % -------------------------------------------------------------
+    set( hSaveMenu, 'UserData', 'on' );
+
+  end
+
+%% CancelCalibrationCallback .................................. Calibration
+  function CancelCalibrationCallback(hObject, eventdata)
+    % Callback function run when
+    
+    % Get tsg application data
+    % ------------------------
+    tsg = getappdata( hMainFig, 'tsg_data' );
+        
+    tsg.SSPS_CAL = [];
+    tsg.SSJT_CAL = [];
+    
+    % Update the Adjusted variables (SSPS - SSJT) with calibrated records
+    % -------------------------------------------------------------------
+    updateAdjustedVariable( hMainFig );
+    
+    % Save tsg application data
+    % --------------------------
+    setappdata( hMainFig, 'tsg_data', tsg );
+    
+    % Refresh plot #1
+    % ---------------
+    plot_Calibration( hMainFig, hPlotAxes, 1 );
+    plot_Calibration( hMainFig, hPlotAxes, 2 );
+    
+    % As soon as a modification took place the data should be saved
+    % -------------------------------------------------------------
+    set( hSaveMenu, 'UserData', 'on' );
+
+  end
 
 %% Zoom_OffMenuCallback
   %----------------------------------------------------------------------
@@ -1801,91 +1861,6 @@ end
 
   end
 
-%% SelectTime_OnMenuCallback ............................ Correction module
-%----------------------------------------------------------------------
-  function SelectTime_OnMenuCallback(hObject, eventdata)
-    % Callback function run when the ....
-            
-    % Desactivate Zoom and Pan functions.
-    % ----------------------------------
-    set( hZoomToggletool, 'state', 'off' );
-    set( hQCToggletool,   'state', 'off' );
-    set( hPanToggletool,  'state', 'off' );
-    
-    % Create a pointer to select the time limits
-    % ------------------------------------------
-    selTimePointer = ones(16)+1;
-    selTimePointer(1,:)       = 1; selTimePointer(16,:)      = 1;
-    selTimePointer(:,1)       = 1; selTimePointer(:,16)      = 1;
-    selTimePointer(1:4,8:9)   = 1; selTimePointer(13:16,8:9) = 1;
-    selTimePointer(8:9,1:4)   = 1; selTimePointer(8:9,13:16) = 1;
-    selTimePointer(5:12,5:12) = NaN; % Create a transparent region in the center
-    
-    % Activate clic mouse menu on second axes (salinity) for next rbbox
-    % ----------------------------------------------------------------
-    set(hMainFig,'WindowButtonDownFcn', @Time_SelectCallback);
-
-    % change cursor
-    % ---------------
-    set( hMainFig, 'Pointer', 'custom',...
-      'PointerShapeCData', selTimePointer, 'PointerShapeHotSpot',[9 9]);
-
-    % ----------------------------------------------------------------------
-    % nested function on mouse clic when Select Time toggle tool is selected
-    % ----------------------------------------------------------------------
-    function Time_SelectCallback(gcbo, eventdata)
-
-      % disable ButtonMotion on main fig during select
-      % prevent drawing to map
-      % ----------------------------------------------
-      set( hMainFig, 'WindowButtonMotionFcn', []);
-
-      % Retrieve named application data
-      % -------------------------------
-      tsg = getappdata( hMainFig, 'tsg_data');
-
-      % Selection of the data within the figure
-      % ---------------------------------------
-      point1    = get(gca,'CurrentPoint');    % button down detected
-      finalRect = rbbox;                      % return figure units
-      point2    = get(gca,'CurrentPoint');    % button up detected
-
-      point1 = point1(1,1:2);                 % extract x and y
-      point2 = point2(1,1:2);
-
-      p1 = min(point1,point2);
-      p2 = max(point1,point2);                % calculate locations
-
-      % get index on selected zone - Only on X axes (time)
-      % --------------------------------------------------
-      ind = find(tsg.DAYD >= p1(1,1) & tsg.DAYD <= p2(1,1));
-
-      % Write the date in the Editable uicontrol
-      % ----------------------------------------
-      set( hetDateMin, 'String', datestr(tsg.DAYD(ind(1)),   31));
-      set( hetDateMax, 'String', datestr(tsg.DAYD(ind(end)), 31));
-
-      % enable ButtonMotion on main fig after select QC area
-      % ----------------------------------------------------
-      set( hMainFig, 'WindowButtonMotionFcn', @MouseMotion);
-    end
-  end
-
-%% SelectTime_OffMenuCallback ........................... Correction module
-%--------------------------------------------------------------------------
-  function SelectTime_OffMenuCallback(hObject, eventdata)
-  % Callback function run when the ....
-  
-  % Desactivate time limit buttons
-  % ------------------------------
-  set( hTimelimitToggletool, 'state', 'off');
-
-  set( hMainFig,'WindowButtonDownFcn', []);
-
-  set(hMainFig,'Pointer','arrow');
-
-  end
-
 %% CorCancelCallback .................................... Correction Module
   function CorCancelCallback(hObject, eventdata)
     % Callback function run when ...
@@ -1989,6 +1964,91 @@ end
 
   end
 
+%% SelectTime_OnMenuCallback
+%---------------------------
+  function SelectTime_OnMenuCallback(hObject, eventdata)
+    % Callback function run when the ....
+            
+    % Desactivate Zoom and Pan functions.
+    % ----------------------------------
+    set( hZoomToggletool, 'state', 'off' );
+    set( hQCToggletool,   'state', 'off' );
+    set( hPanToggletool,  'state', 'off' );
+    
+    % Create a pointer to select the time limits
+    % ------------------------------------------
+    selTimePointer = ones(16)+1;
+    selTimePointer(1,:)       = 1; selTimePointer(16,:)      = 1;
+    selTimePointer(:,1)       = 1; selTimePointer(:,16)      = 1;
+    selTimePointer(1:4,8:9)   = 1; selTimePointer(13:16,8:9) = 1;
+    selTimePointer(8:9,1:4)   = 1; selTimePointer(8:9,13:16) = 1;
+    selTimePointer(5:12,5:12) = NaN; % Create a transparent region in the center
+    
+    % Activate clic mouse menu on second axes (salinity) for next rbbox
+    % ----------------------------------------------------------------
+    set(hMainFig,'WindowButtonDownFcn', @Time_SelectCallback);
+
+    % change cursor
+    % ---------------
+    set( hMainFig, 'Pointer', 'custom',...
+      'PointerShapeCData', selTimePointer, 'PointerShapeHotSpot',[9 9]);
+
+    % ----------------------------------------------------------------------
+    % nested function on mouse clic when Select Time toggle tool is selected
+    % ----------------------------------------------------------------------
+    function Time_SelectCallback(gcbo, eventdata)
+
+      % disable ButtonMotion on main fig during select
+      % prevent drawing to map
+      % ----------------------------------------------
+      set( hMainFig, 'WindowButtonMotionFcn', []);
+
+      % Retrieve named application data
+      % -------------------------------
+      tsg = getappdata( hMainFig, 'tsg_data');
+
+      % Selection of the data within the figure
+      % ---------------------------------------
+      point1    = get(gca,'CurrentPoint');    % button down detected
+      finalRect = rbbox;                      % return figure units
+      point2    = get(gca,'CurrentPoint');    % button up detected
+
+      point1 = point1(1,1:2);                 % extract x and y
+      point2 = point2(1,1:2);
+
+      p1 = min(point1,point2);
+      p2 = max(point1,point2);                % calculate locations
+
+      % get index on selected zone - Only on X axes (time)
+      % --------------------------------------------------
+      ind = find(tsg.DAYD >= p1(1,1) & tsg.DAYD <= p2(1,1));
+
+      % Write the date in the Editable uicontrol
+      % ----------------------------------------
+      set( hetDateMin, 'String', datestr(tsg.DAYD(ind(1)),   31));
+      set( hetDateMax, 'String', datestr(tsg.DAYD(ind(end)), 31));
+
+      % enable ButtonMotion on main fig after select QC area
+      % ----------------------------------------------------
+      set( hMainFig, 'WindowButtonMotionFcn', @MouseMotion);
+    end
+  end
+
+%% SelectTime_OffMenuCallback
+%----------------------------
+  function SelectTime_OffMenuCallback(hObject, eventdata)
+  % Callback function run when the ....
+  
+  % Desactivate time limit buttons
+  % ------------------------------
+  set( hTimelimitToggletool, 'state', 'off');
+
+  set( hMainFig, 'WindowButtonDownFcn', []);
+
+  set( hMainFig, 'Pointer', 'arrow');
+
+  end
+
 %% Clim_OffMenuCallback
   %------------------------------------------------------------------------
   % Callback function run when the Levitus climatology toolbar is unselected