function THETA=potentmp(S,T0,P0,PR); % POTENTMP Potential temperature from in-situ measurements. % THETA=POTENTMP(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 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 adiabat.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.*adiabat(S,T,P) ; T = T + 0.5*XK; Q = XK ; P = P + 0.5*H ; XK = H.*adiabat(S,T,P) ; T = T + 0.29289322*(XK-Q) ; Q = 0.58578644*XK + 0.121320344*Q ; XK = H.*adiabat(S,T,P) ; T = T + 1.707106781*(XK-Q); Q = 3.414213562*XK - 4.121320344*Q; P = P + 0.5*H ; XK = H.*adiabat(S,T,P) ; THETA = T + (XK-2.0*Q)/6.0;