/*      REAL FUNCTION THETA (S,T,P,PR)                                         *
C--------------------------------------------------------------------------
C
C                            Function = THETA
C
C     Function to compute local potential temperature at PR using Bryden's
C     1973 polynomial for adiabatic lapse rate and Runge-Kutta 4th order
C     integration algorithm.
C
C     Input parameters:
C            S      -  Salinity (DSS-78)
C            T0     -  Temperature in degrees C
C            P0     -  Pressure in decibars
C            PR     -  Reference pressure in decibars
C     Output value:
C            THETA  -  Potential temperature in degrees C
C
C     Reference:  Bryden, H., 1973, Deep-Sea Research,20,401-408
C                 Fofonoff,N., 1977, Deep-Sea research,24,489-491
C
C     Checkvalue: THETA= 36.89073 C,S=40 (PSS-78),T0=40 DEG C,
C                 P0=10000 DECIBARS,PR=0 DECIBARS
C
C--------------------------------------------------------------------------
*/
#include <stdio.h>
#include <math.h>
#include "ocean.h"

#if PROTOTYPE_ALLOWED
double THETA (double S, double T0, double P0, double PR)
#else
double THETA (S, T0, P0, PR)
double S, T0, P0, PR;
#endif
{
     double P, T, H, XK, Q, theta;
/*    =====>  Set up intermediate temperature and pressure values.... */

      P=P0;
      T=T0;
      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);
      theta = T + (XK-2.0*Q)/6.0;

     return(theta);
}
#ifdef EXE
#if PROTOTYPE_ALLOWED
void main(int argc, char *argv[])
#else
void main(argc, argv)
int argc;
char *argv[];
#endif
{
   double P = 10000.0, T = 40.0, S = 40.0, PR = 0.0;

   if (argc != 5)
   {
      printf("\n USAGE: theta <salinity> <temperature> <pressure> <reference-pressure>");
      printf("\n EXAMPLE: theta 40.0 40.0 10000.0 0.0\n");
   }
   else
   {
      S = atof(argv[1]);
      T = atof(argv[2]);
      P = atof(argv[3]);
      PR = atof(argv[4]);
   }
   printf("\n potential temperature = %f deg C", THETA(S, T, P, PR));
   printf("\n for S = %g (PSS-78), T = %g deg C, P = %g dbar, ref P = %g dbar\n\n", S, T, P, PR);
}
#endif