function [the]=theta(p0,t0,s,pr)
% theta fcn ****** april 6 1981 *****
% ***********************************
% to compute local potential temperature at p0 referenced to 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.,?977,deep-sea res.?24,489-491
%%checkvalue: theta= 36.89072 deg c for s=40 nsu,t=40 deg c,
%%p=10000 dbars,pr=0 dbars
%%units of inputs p0 and pr are decibars
%% units changed to decibars for use at mit
% modified by N.L. Bindoff for use in MATLAB
%
% input
% p0  m*n matrix of pressures (db)
% t0  m*n matrix of insitu temperatures (degrees C)
% s   m*n vector of salinity (psu)
% pr  scalar
%
% Output
% the m*n matrix potential temperature with respect to the pr.
%
%
      [m,n]=size(p0);
      xk=zeros(m,n);
      h = pr*ones(m,n) - p0;
      xk = h.*adtg(p0,t0,s)/1000.0;
      t0 = t0 + 0.50*xk;
      q = xk;
      p0 = p0 + 0.50*h;
      p=p0;
      t=t0;
      xk = h.*adtg(p,t,s)/1000.0;
      t0 = t0 + 0.292893220*(xk-q);
      q = 0.58578644*xk + 0.121320344*q;
      t=t0;
      xk = h.*adtg(p,t,s)/1000.0;
      t0 = t0 + 1.707106781*(xk-q);
      q = 3.414213562*xk - 4.1213203440*q;
      p0 = p0 + 0.50*h;
      p=p0;
      t=t0;
      xk = h.*adtg(p,t,s)/1000.0;
      the = t0 + (xk-2.00*q)/6.0;





