function[] = ctd_plot(varargin); % function[] = ctd_plot(varargin); % % This script makes plots for CTD casts % % ctd_plot('file',OPTIONS) % % CTD_PLOT(INPUTFILE) will create plot T,S,Sigma as a function of Z. % This is the default setup. % % Possible Options: % % 'out' -- "'out',OUTPUTFILE" will print the figure to a colour % postscript % % 'cast' -- "'cast','tscr'" or any combination of those letters will % plot what are selected of (t)emperature, (s)alinity, % (c)onductivity, and (r)ho, i.e. sigma. Default is TSR. % % 'vert' -- "'vert','z/r/p' will plot the profiles against Z (depth) % Sigma (r), or Pressure(p). Default is Z. % %---------------------------------------- % % e.g. ctd_plot('BED0301.bdn','out','stn1.ps','vert','p','cast','ts') % % This command will open "BED0301.bdn" and plot T and S vs. P, % printing the output to the file stn1.ps % %---------------------------------------- % % This function assumes the following columnated structure for the input file: % % column1: depth % column2: pressure % column3: temperature % column4: salinity % % Currently limits on Temperature, Salinity, Conductivity, Depth, etc., % are set to match the Bedford Basin Data set. See the beginning of the % code to edit them. %--------------------------------------------------------------------- % % AUTHOR: Ramzi Mirshak, Department of Oceanography, Dalhousie University % CONTACT: ramzi.mirshak@phys.ocean.dal.ca % DATE: 23 October, 2003 %--------------------------------------------------------------------- % Limits for variables to be plotted salimit = [28 32]; % Salinity (PSU) -- original values [28 32] delimit = [-1 65]; % Depth (m) -- original values [-1 65] colimit = [0.6 0.9]; % Conductivity (R_p) -- original values [0.6 0.9] rhlimit = [20 26]; % Density (Sigma) -- original values [20 26] telimit = [0 15]; % Temperature -- original values [0 15] % Number of tickmarks in across the axis for different variables samarks = 5; % Salinity comarks = 4; % Conductivity rhmarks = 4; % Density temarks = 4; % Temperature %===================================================================== %================== DO NOT EDIT BELOW THIS LINE!! ================== %===================================================================== error(nargchk(1,7,nargin)); if((nargin + 1)/2 ~= ceil((nargin+1)/2)) error('Bad choice of inputs'); end % Flags set default values... Tflag = 1; % things to be plotted Sflag = 1; Cflag = 0; Rflag = 1; Pflag = 0; % things to plot against R2flag =0; outflag = 0; % print as postscript? f_in = varargin{1}; if (nargin > 1) switch(varargin{2}) case 'out' outflag = 1; output = varargin{3}; case 'vert'; vertical = 1; if(varargin{3} == 'p') Pflag = 1; elseif(varargin{3} == 'r') R2flag = 1; end case 'cast' cast = 1; Tflag = 0; Sflag = 0; Rflag = 0; look = findstr(varargin{3},'t'); if(~isempty(look)), Tflag = 1; end; look = findstr(varargin{3},'s'); if(~isempty(look)), Sflag = 1; end; look = findstr(varargin{3},'c'); if(~isempty(look)), Cflag = 1; end; look = findstr(varargin{3},'r'); if(~isempty(look)), Rflag = 1; end; otherwise error(['Unknown Option: ',varargin{2}]); end end if (nargin > 3) switch(varargin{4}) case 'out' outflag = 1; output = varargin{5}; case 'vert'; if(varargin{5} == 'p') Pflag = 1; elseif(varargin{5} == 'r') R2flag = 1; end case 'cast' Tflag = 0; Sflag = 0; Rflag = 0; look = findstr(varargin{5},'t'); if(~isempty(look)), Tflag = 1; end; look = findstr(varargin{5},'s'); if(~isempty(look)), Sflag = 1; end; look = findstr(varargin{5},'c'); if(~isempty(look)), Cflag = 1; end; look = findstr(varargin{5},'r'); if(~isempty(look)), Rflag = 1; end; otherwise error(['Unknown Option: ',varargin{4}]); end end if (nargin > 5) switch(varargin{6}) case 'out' outflag = 1; output = varargin{7}; case 'vert'; if(varargin{7} == 'p') Pflag = 1; elseif(varargin{7} == 'r') R2flag = 1; end case 'cast' Tflag = 0; Sflag = 0; Rflag = 0; look = findstr(varargin{7},'t'); if(~isempty(look)), Tflag = 1; end; look = findstr(varargin{7},'s'); if(~isempty(look)), Sflag = 1; end; look = findstr(varargin{7},'c'); if(~isempty(look)), Cflag = 1; end; look = findstr(varargin{7},'r'); if(~isempty(look)), Rflag = 1; end; otherwise error(['Unknown Option: ',varargin{4}]); end end %addpath /home/mirshak/matlab/oceans %addpath /home/mirshak/matlab/oceans_dk % ticklines above plots for SAlinity, TEmperature, COnductivity, and RHo satick = [salimit(1):(diff(salimit)/(samarks-1)):salimit(2)]; tetick = [telimit(1):(diff(telimit)/(temarks-1)):telimit(2)]; cotick = [colimit(1):(diff(colimit)/(comarks-1)):colimit(2)]; rhtick = [rhlimit(1):(diff(rhlimit)/(rhmarks-1)):rhlimit(2)]; cast_in = load(f_in); Z = cast_in(:,1); P = cast_in(:,2); T = cast_in(:,3); S = cast_in(:,4); C = sw_cndr(S,T,P); R = swstate(S,T,P); if(R2flag) depth_var = R; delimit = rhlimit; depth_msg = 'Density, \sigma_T (kg m^{-3})'; elseif(Pflag); depth_var = P; depth_msg = 'Pressure (dbar)'; else depth_var = Z; depth_msg = 'Depth (m)'; end clf %===================================================================== % Temperature if(Tflag) set(axes,'position',[0.1 0.1 0.8 0.68]); plot(T,depth_var,'r','linewidth',2); axis([telimit delimit]); zticks = get(gca,'ytick'); set(gca,'visible','off'); axis ij; set(axes,'position',[0.1 0.73 0.8 0.1]); set(gca,'visible','off'); axis([telimit 0 1]); axis ij; line([telimit],[0.35 0.35],'color','r'); for i = 1:length(tetick) eval(['line ([',num2str(tetick(i)),' ',num2str(tetick(i)),'], ', ... '[0.35 0.44],''color'',''r'');']); eval(['text (',num2str(tetick(i)),',0.2,''',num2str(tetick(i)), ... ''',''horiz'',''center'',''fontsize'',12,''color'',''r'');']) end text(telimit(1) + diff(telimit)/2,0.05,'Temperature (^oC)' ... ,'horiz','center','fontsize',14,'color','r'); end %--------------------------------------------------------------------- %===================================================================== % Salinity if(Sflag) set(axes,'position',[0.1 0.1 0.8 0.68]); plot(S,depth_var,'g','linewidth',2); axis([salimit delimit]); zticks = get(gca,'ytick'); set(gca,'visible','off'); axis ij; set(axes,'position',[0.1 0.78 0.8 0.1]); set(gca,'visible','off'); axis([salimit 0 1]); axis ij; line([salimit],[0.35 0.35],'color','g'); for i = 1:length(satick) eval(['line ([',num2str(satick(i)),' ',num2str(satick(i)),'], ', ... '[0.35 0.44],''color'',''g'');']); eval(['text (',num2str(satick(i)),',0.2,''',num2str(satick(i)), ... ''',''horiz'',''center'',''fontsize'',12,''color'',''g'');']) end text(salimit(1) + diff(salimit)/2,0.05,'Salinity' ... ,'horiz','center','fontsize',14,'color','g'); end %--------------------------------------------------------------------- %===================================================================== % Density if(Rflag) set(axes,'position',[0.1 0.1 0.8 0.68]); plot(R,depth_var,'k','linewidth',2); axis([rhlimit delimit]); zticks = get(gca,'ytick'); set(gca,'visible','off'); axis ij; set(axes,'position',[0.1 0.83 0.8 0.1]); set(gca,'visible','off'); axis([rhlimit 0 1]); axis ij; line([rhlimit],[0.35 0.35],'color','k'); for i = 1:length(rhtick) eval(['line ([',num2str(rhtick(i)),' ',num2str(rhtick(i)),'], ', ... '[0.35 0.44],''color'',''k'');']); eval(['text (',num2str(rhtick(i)),',0.2,''',num2str(rhtick(i)), ... ''',''horiz'',''center'',''fontsize'',12,''color'',''k'');']) end text(rhlimit(1) + diff(rhlimit)/2,0.1,'Density, \sigma_T (kg m^{-3})', ... 'horiz','center','fontsize',14,'color','k'); end %--------------------------------------------------------------------- %===================================================================== % Conductivity if(Cflag) set(axes,'position',[0.1 0.1 0.8 0.68]); plot(C,depth_var,'b','linewidth',2); axis([colimit delimit]); zticks = get(gca,'ytick'); set(gca,'visible','off'); axis ij; set(axes,'position',[0.1 0.88 0.8 0.1]); set(gca,'visible','off'); axis([colimit 0 1]); axis ij; line([colimit],[0.35 0.35],'color','b'); for i = 1:length(cotick) eval(['line ([',num2str(cotick(i)),' ',num2str(cotick(i)),'], ', ... '[0.35 0.44],''color'',''b'');']); eval(['text (',num2str(cotick(i)),',0.2,''',num2str(cotick(i)), ... ''',''horiz'',''center'',''fontsize'',12,''color'',''b'');']) end text(colimit(1) + diff(colimit)/2,0.05,'Conductivity Ratio, R_p', ... 'horiz','center','fontsize',14,'color','b'); end %--------------------------------------------------------------------- set(axes,'position',[0.04 0.1 0.1 0.68]); axis([0 1 delimit]); line([0.5 0.5],[zticks(1) zticks(end)],'color','k'); for i = 1:length(zticks) line([0.5 0.6],[zticks(i) zticks(i)],'color','k'); text(0.4, zticks(i), num2str(zticks(i)),'vert','middle','fontsize',12, ... 'color','k','horiz','right'); end text(0.05, mean(zticks([1 end])),depth_msg,'fontsize',14,'color','k', ... 'horiz','center','rotation',90); axis ij set(gca,'visible','off'); orient tall wysiwyg if(outflag) eval(['print -dpsc ',output]); end %--------------------------------------------------------------------- %--------------------------------------------------------------------- %--- BELOW THIS LINE ARE FUNCTIONS THAT ARE REQUIRED BY THE ABOVE ---- %---- SCRIPT. DO NOT ALTER OR EDIT THESE IN ANY WAY WHATSOEVER!! ---- %--------------------------------------------------------------------- %--------------------------------------------------------------------- function [SVAN,SIGMA]=swstate(S,T,P0,ftn); % SWSTATE State equation for seawater % % SIGMA=SWSTATE(S,T,P) returns the density anomaly SIGMA (kg/m^3) % given the salinity S (ppt), temperature T (deg C) and pressure % P (dbars). % % [SVAN,SIGMA]=SWSTATE(S,T,P) returns the specific volume % anomaly SVAN (m^3/kg*1e-8) and the density anomaly SIGMA (kg/m^3). % % [dVdT,dRdT]=SWSTATE(S,T,P,'dT') returns derivatives w.r.t. % temperature of the volume and density. % % [dVdS,dRdS]=SWSTATE(S,T,P,'dS') returns derivatives w.r.t. % salinity. % % [dVdP,dRdP]=SWSTATE(S,T,P,'dP') returns derivatives w.r.t. % pressure. % % In all cases, if only one output variable is specified, the % relevant variable will be SIGMA rather than SVAN. % % All elements can be scalars, vectors, or matrices but should be % the same size. %Notes: RP (WHOI 2/dec/91) % % This stuff is directly copied from the UNESCO algorithms, with some % minor changes to make it Matlab compatible (like adding ";" and changing % "*" to ".*" when necessary. % % RP (WHOI 3/dec/91) % % Added first derivative calculations. % % RP (IOS) 4/Apr/96 % 8-char names *and* made the default thing to return SIGMA rather than % SVAN (which no-one uses). derivT=0; derivS=0; derivP=0; if (nargin==4), if (ftn=='dT'), derivT=1; elseif (ftn=='dS'), derivS=1; elseif (ftn=='dP'), derivP=1; else error('swstate: Unrecognized option!'); end; end; % ****************************************************** % SPECIFIC VOLUME ANOMALY (STERIC ANOMALY) BASED ON 1980 EQUATION % OF STATE FOR SEAWATER AND 1978 PRACTICAL SALINITY SCALE. % REFERENCES % MILLERO, ET AL (1980) DEEP-SEA RES.,27A,255-264 % MILLERO AND POISSON 1981,DEEP-SEA RES.,28A PP 625-629. % BOTH ABOVE REFERENCES ARE ALSO FOUND IN UNESCO REPORT 38 (1981) % UNITS: % PRESSURE P0 DECIBARS % TEMPERATURE T DEG CELSIUS (IPTS-68) % SALINITY S (IPSS-78) % SPEC. VOL. ANA. SVAN M**3/KG *1.0E-8 % DENSITY ANA. SIGMA KG/M**3 % ****************************************************************** % CHECK VALUE: SVAN=981.3021 E-8 M**3/KG. FOR S = 40 (IPSS-78) , % T = 40 DEG C, P0= 10000 DECIBARS. % CHECK VALUE: SIGMA = 59.82037 KG/M**3 FOR S = 40 (IPSS-78) , % T = 40 DEG C, P0= 10000 DECIBARS. %HECK VALUE: FOR S = 40 (IPSS-78) , T = 40 DEG C, P0= 10000 DECIBARS. % DR/DP DR/DT DR/DS % DRV(1,7) DRV(2,3) DRV(1,8) % % FINITE DIFFERENCE WITH 3RD ORDER CORRECTION DONE IN DOUBLE PRECSION % % 3.46969238E-3 -.43311722 .705110777 % % EXPLICIT DIFFERENTIATION SINGLE PRECISION FORMULATION EOS80 % % 3.4696929E-3 -.4331173 .7051107 % % (RP...I think this ---------^^^^^^ should be -.4431173!); % ******************************************************* % DATA R3500=1028.1063; R4=4.8314E-4; DR350=28.106331; % CONVERT PRESSURE TO BARS AND TAKE SQUARE ROOT SALINITY. P=P0/10.; SAL=S; SR = sqrt(abs(S)); % ********************************************************* % PURE WATER DENSITY AT ATMOSPHERIC PRESSURE % BIGG P.H.,(1967) BR. J. APPLIED PHYSICS 8 PP 521-537. % R1 = ((((6.536332E-9*T-1.120083E-6).*T+1.001685E-4).*T ... -9.095290E-3).*T+6.793952E-2).*T-28.263737; % SEAWATER DENSITY ATM PRESS. % COEFFICIENTS INVOLVING SALINITY R2 = (((5.3875E-9*T-8.2467E-7).*T+7.6438E-5).*T-4.0899E-3).*T+8.24493E-1; R3 = (-1.6546E-6*T+1.0227E-4).*T-5.72466E-3; % INTERNATIONAL ONE-ATMOSPHERE EQUATION OF STATE OF SEAWATER SIG = (R4*S + R3.*SR + R2).*S + R1; % SPECIFIC VOLUME AT ATMOSPHERIC PRESSURE V350P = 1.0/R3500; SVA = -SIG*V350P./(R3500+SIG); SIGMA=SIG+DR350; V0 = 1.0./(1000.0 + SIGMA); % SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS SVAN=SVA*1.0E+8; if (derivS), % These are derivatives for (S,T,0). R4S=9.6628E-4; RHO1 = 1000.0 + SIGMA; RHOS=R4S*SAL+1.5.*R3.*SR+R2; V0S=-RHOS./(RHO1.*RHO1); elseif (derivT), R1 =(((3.268166E-8*T-4.480332E-6).*T+3.005055E-4).*T... -1.819058E-2).*T+6.793952E-2; R2 = ((2.155E-8*T-2.47401E-6).*T+1.52876E-4).*T-4.0899E-3; R3 = -3.3092E-6*T+1.0227E-4; RHO1 = 1000.0 + SIGMA; RHOT = (R3.*SR + R2).*SAL + R1; V0T = -RHOT./(RHO1.*RHO1); end; % ****************************************************************** % ****** NEW HIGH PRESSURE EQUATION OF STATE FOR SEAWATER ******** % ****************************************************************** % MILLERO, ET AL , 1980 DSR 27A, PP 255-264 % CONSTANT NOTATION FOLLOWS ARTICLE %******************************************************** % COMPUTE COMPRESSION TERMS E = (9.1697E-10*T+2.0816E-8).*T-9.9348E-7; BW = (5.2787E-8*T-6.12293E-6).*T+3.47718E-5; B = BW + E.*S; % Bulk Modulus (almost) % CORRECT B FOR ANAMOLY BIAS CHANGE Bout = B + 5.03217E-5; if (derivS), DBDS=E; elseif (derivT), BW = 1.05574E-7*T-6.12293E-6; E = 1.83394E-9*T +2.0816E-8; BT = BW + E.*SAL; end; % D = 1.91075E-4; C = (-1.6078E-6*T-1.0981E-5).*T+2.2838E-3; AW = ((-5.77905E-7*T+1.16092E-4).*T+1.43713E-3).*T-0.1194975; A = (D*SR + C).*S + AW; % CORRECT A FOR ANAMOLY BIAS CHANGE Aout = A + 3.3594055; if (derivS), DADS=2.866125E-4*SR+C; elseif (derivT), C = -3.2156E-6*T -1.0981E-5; AW = (-1.733715E-6*T+2.32184E-4).*T+1.43713E-3; AT = C.*SAL + AW; end; B1 = (-5.3009E-4*T+1.6483E-2).*T+7.944E-2; A1 = ((-6.1670E-5*T+1.09987E-2).*T-0.603459).*T+54.6746; KW = (((-5.155288E-5*T+1.360477E-2).*T-2.327105).*T+148.4206).*T-1930.06; K0 = (B1.*SR + A1).*S + KW; if (derivS), K0S=1.5*B1.*SR+A1; KS=(DBDS.*P+DADS).*P+K0S; elseif (derivT), B1 = -1.06018E-3*T+1.6483E-2; % APRIL 9 1984 CORRECT A1 BIAS FROM -.603457 !!! A1 = (-1.8501E-4*T+2.19974E-2).*T-0.603459; KW = ((-2.0621152E-4*T+4.081431E-2).*T-4.65421).*T+148.4206; K0T = (B1.*SR+A1).*SAL + KW; KT = (BT.*P + AT).*P + K0T; end; % EVALUATE PRESSURE POLYNOMIAL % *********************************************** % K EQUALS THE SECANT BULK MODULUS OF SEAWATER % DK=K(S,T,P)-K(35,0,P) % K35=K(35,0,P) % *********************************************** DK = (B.*P + A).*P + K0; K35 = (5.03217E-5*P+3.359406).*P+21582.27; GAM=P./K35; PK = 1.0 - GAM; SVA = SVA.*PK + (V350P+SVA).*P.*DK./(K35.*(K35+DK)); % SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS SVAN=SVA*1.0E+8; % Volume anomaly V350P = V350P.*PK; % **************************************************** % COMPUTE DENSITY ANAMOLY WITH RESPECT TO 1000.0 KG/M**3 % 1) DR350: DENSITY ANAMOLY AT 35 (IPSS-78), 0 DEG. C AND 0 DECIBARS % 2) DR35P: DENSITY ANAMOLY 35 (IPSS-78), 0 DEG. C , PRES. VARIATION % 3) DVAN : DENSITY ANAMOLY VARIATIONS INVOLVING SPECFIC VOL. ANAMOLY % ******************************************************************** % CHECK VALUE: SIGMA = 59.82037 KG/M**3 FOR S = 40 (IPSS-78), % T = 40 DEG C, P0= 10000 DECIBARS. % ******************************************************* DR35P=GAM./V350P; DVAN=SVA./(V350P.*(V350P+SVA)); SIGMA=DR350+DR35P-DVAN; % Density anomaly K=K35+DK; VP=1.0-P./K; V = (1.) ./(SIGMA+1000.0); if (derivS), VS=V0S.*VP+V0.*P.*KS./(K.*K); SVAN=VS; % dVdS SIGMA=-VS./(V.*V); % dRdS elseif (derivT), VT = V0T.*VP + V0.*P.*KT./(K.*K); SVAN=VT; % dVdT SIGMA=-VT./(V.*V); % dRdT elseif (derivP), DKDP = 2.0*Bout.*P + Aout; % CORRECT DVDP TO PER DECIBAR BY MULTIPLE *.1 DVDP = -.1*V0.*(1.0 - P.*DKDP./K)./K; SVAN=DVDP; % dVdP SIGMA=-DVDP./(V.*V); % dRdP end; if (nargout==1), SVAN=SIGMA; end; 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 %----------------------------------------------------------------------- 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 %---------------------------------------------------------------------- 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 %----------------------------------------------------------------------- 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; rt = c0 + (c1 + (c2 + (c3 + c4.*T).*T).*T).*T; return %-------------------------------------------------------------------- function wysiwyg %WYSIWYG -- this function is called with no args and merely % changes the size of the figure on the screen to equal % the size of the figure that would be printed, % according to the papersize attribute. Use this function % to give a more accurate picture of what will be % printed. % Dan(K) Braithwaite, Dept. of Hydrology U.of.A 11/93 unis = get(gcf,'units'); ppos = get(gcf,'paperposition'); set(gcf,'units',get(gcf,'paperunits')); pos = get(gcf,'position'); pos(3:4) = ppos(3:4); set(gcf,'position',pos); set(gcf,'units',unis);