pottemp.m 2.64 KiB
function THETA=pottemp(S,T0,P0,PR);
% POTTEMP Potential temperature from in-situ measurements.
% THETA=POTTEMP(S,TO,PO,PR) computes the potential temp.
% at a reference pressure PR (dbars) corresponding to the
% salinity S (ppt) and temperature TO (deg C) at pressure
% PO (dbars). The formula has been copied from the UNESCO
% algorithms (See comments for details).
%
% if PR is omitted, it is assumed to be zero.
%
% S, T0, P0 can be vectors or matrices. If S (T) is a scalar,
% and T (S) is a vector or matrix, we assume that the scalar
% value corresponds to all values in the vector.
%
% If P0 is a vector, and one of S,T0 is a matrix, we assume
% that elements of P0 correspond with all row entries in the
% matrices.
%
% REF: BRYDEN,H.,1973,DEEP-SEA RES.,20,401-408
% FOFONOFF,N.,1977,DEEP-SEA RES.,24,489-491
% Converted to 8 char length function names and for version 5.1 (AN)
if (nargin==3), PR=0; end;
[NS,MS]=size(S);
[NT,MT]=size(T0);
if ((NS==1 & MS==1) & (NT ~=1 | MT ~=1) ),
S=S*ones(NT,MT);
[NS,MS]=size(S);
elseif ((NT==1 & MT==1) & (NS ~=1 | MS ~=1) ),
T0=T0*ones(NS,MS);
[NT,MT]=size(T0);
end;
[NP,MP]=size(P0);
if (MP==1), P0=P0*ones(1,MT);
elseif (NP==1), P0=P0'*ones(1,MT); end;
% Notes: RP 29/Nov/91
%
% I have modified the FORTRAN code to make it Matlab compatible, but
% no numbers have been changed. In certain places "*" has been replaced
% with ".*" to allow vectorization.
%
% This routine calls adiabatt.m (renamed from ATG).
%C ***********************************
%C TO COMPUTE LOCAL POTENTIAL TEMPERATURE AT PR
%C USING BRYDEN 1973 POLYNOMIAL FOR ADIABATIC LAPSE RATE
%C AND RUNGE-KUTTA 4-TH ORDER INTEGRATION ALGORITHM.
%C UNITS:
%C PRESSURE P0 DECIBARS
%C TEMPERATURE T0 DEG CELSIUS (IPTS-68)
%C SALINITY S (IPSS-78)
%C REFERENCE PRS PR DECIBARS
%C POTENTIAL TMP. THETA DEG CELSIUS
%C CHECKVALUE: THETA= 36.89073 C,S=40 (IPSS-78),T0=40 DEG C,
%C P0=10000 DECIBARS,PR=0 DECIBARS
%C
%C SET-UP INTERMEDIATE TEMPERATURE AND PRESSURE VARIABLES
% IMPLICIT REAL*8 (A-H,O-Z)
P=P0;
T=T0;
%C**************
H = PR - P;
XK = H.*adiabatt(S,T,P) ;
T = T + 0.5*XK; Q = XK ;
P = P + 0.5*H ;
XK = H.*adiabatt(S,T,P) ;
T = T + 0.29289322*(XK-Q) ;
Q = 0.58578644*XK + 0.121320344*Q ;
XK = H.*adiabatt(S,T,P) ;
T = T + 1.707106781*(XK-Q);
Q = 3.414213562*XK - 4.121320344*Q;
P = P + 0.5*H ;
XK = H.*adiabatt(S,T,P) ;
THETA = T + (XK-2.0*Q)/6.0;