function th=theta(S,T,P,PR)
%function th=theta(S,T,P,PR)
%
% ******************************************************
% TO COMPUTE LOCAL POTENTIAL TEMPERATURE AT PR
% USING BRYDEN 1973 POLYNOMIAL FOR ADIABATIC LAPSE RATE
% AND RUNGE-KUTTA 4-TH ORDER INTEGRATION ALGORITHM.
% REF: BRYDEN,H.,1973,DEEP-SEA RES.,20,401-408
% FOFONOFF,N.,1977,DEEP-SEA RES.,24,489-491
% ***
% UNITS:      
%       PRESSURE        P0       DECIBARS
%       TEMPERATURE     T0       DEG CELSIUS (IPTS-68)
%       SALINITY        S        (IPSS-78)
%       REFERENCE PRS   PR       DECIBARS
%       POTENTIAL TMP.  THETA    DEG CELSIUS 
% ***
% CHECKVALUE: THETA= 36.89073 C,S=40 (IPSS-78),T0=40 DEG C,
% P0=10000 DECIBARS,PR=0 DECIBARS
% ***             

      H = PR - P;
      XK = H.*atg(S,T,P);
      T = T + 0.5*XK;
      Q = XK ; 
      P = P + 0.5*H; 
      XK = H.*atg(S,T,P);
      T = T + 0.29289322*(XK-Q); 
      Q = 0.58578644*XK + 0.121320344*Q; 
      XK = H.*atg(S,T,P);
      T = T + 1.707106781*(XK-Q);
      Q = 3.414213562*XK - 4.121320344*Q;
      P = P + 0.5*H; 
      XK = H.*atg(S,T,P);
      th = T + (XK-2.0*Q)/6.0;
