      subroutine state(tx,sx,depth,rho)
c
c=======================================================================
c                                                                   ====
c  Nonlinear or linear equation of state for seawater.              ====
c                                                                   ====
c  On Input:                                                        ====
c                                                                   ====
c     T       temperature (real array; degress Celcius)             ====
c     S       salinity anomaly from SMEAN (real array; PSS)         ====
c             (SMEAN is added locally in order to compute density)  ====
c     DEPTH   depths (real array; cm)                               ====
c                                                                   ====
c  On Output:                                                       ====
c                                                                   ====
c     RHO     density anomaly (real array; nondimensional)          ====
c                                                                   ====
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <rhomean.h>
c
c-----------------------------------------------------------------------
c  Define local data
c-----------------------------------------------------------------------
c
      integer i,k
      FLOAT
     *      depth(imt,km),tx(imt,km),sx(imt,km),rho(imt,km)
#ifdef unesco
      FLOAT
     *      a,aw,a1,b,bw,b1,c,d,d350,dk,dr35p,dr350,dvan,e,gam,kw,k0,
     *      k35,pk,r1,r2,r3,r4,r3500,t,s,sig,sigma,sr,sva,v350p,z
c
      save dr350,r3500,r4
c
      data dr350,r3500,r4 /28.106331,1028.1063,4.8314e-4/
# endif
c
c-----------------------------------------------------------------------
c  Begin executable code.
c-----------------------------------------------------------------------
c
#ifndef leqstate
# ifdef unesco
c  International equation of state for seawater (Millero et al., 1980:
c  A new high-pressure equation of state for seawater. Deep-Sea Res.,
c  27a, 255-264).
c
      do 20 k=1,km
        do 10 i=1,imt
          z=depth(i,k)*0.001
          t=tx(i,k)
          s=sx(i,k)+smean
          sr=sqrt(abs(s))
          r1=((((6.536332e-9*t-1.120083e-6)*t+1.001685e-4)*t
     *       -9.095290e-3)*t+6.793952e-2)*t-28.263737
          r2=(((5.3875e-9*t-8.2467e-7)*t+7.6438e-5)*t-4.0899e-3)*t
     *       +8.24493e-1
          r3=(-1.6546e-6*t+1.0227e-4)*t-5.72466e-3
          sig=(r4*s+r3*sr+r2)*s+r1
          v350p=1.0/r3500
          sva=-sig*v350p/(r3500+sig)
          sigma=dr350-(sva/(v350p*(v350p+sva)))
          e=(9.1697e-10*t+2.0816e-8)*t-9.9348e-7
          bw=(5.2787e-8*t-6.12293e-6)*t+3.47718e-5
          b=bw+e*s
          d=1.91075e-4
          c=(-1.6078e-6*t-1.0981e-5)*t+2.2838e-3
          aw=((-5.77905e-7*t+1.16092e-4)*t+1.43713e-3)*t-0.1194975
          a=(d*sr+c)*s+aw
          b1=(-5.3009e-4*t+1.6483e-2)*t+7.944e-2
          a1=((-6.1670e-5*t+1.09987e-2)*t-0.603459)*t+54.6746
          kw=(((-5.155288e-5*t+1.360477e-2)*t-2.327105)*t+148.4206)*t
     *       -1930.06
          k0=(b1*sr+a1)*s+kw
          dk=(b*z+a)*z+k0
          k35=(5.03217e-5*z+3.359406)*z+21582.27
          gam=z/k35
          pk=1.0-gam
          sva=sva*pk+(v350p+sva)*z*dk/(k35*(k35+dk))
          v350p=v350p*pk
          d350=gam/pk
          dr35p=r3500*d350
          dvan=sva/(v350p*(v350p+sva))
          sigma=dr350+dr35p-dvan
          rho(i,k)=sigma/rho0
  10    continue
  20  continue
# else
c  Simplified nonlinear equation of state.
c
      do 20 k=1,km
        do 10 i=1,imt
          rho(i,k)=5891.+depth(i,k)/1013.+3.*(sx(i,k)+smean)+
     *             (38.-0.375*tx(i,k))*tx(i,k)
  10    continue
  20  continue
      do 40 k=1,km
        do 30 i=1,imt
          rho(i,k)=(1779.5+(11.25-0.0745*tx(i,k))*tx(i,k)
     *             -(3.8+0.01*tx(i,k))*(sx(i,k)+smean))/rho(i,k)
  30    continue
  40  continue
      do 60 k=1,km
        do 50 i=1,imt
          rho(i,k)=1.0/(0.698+rho(i,k))-1.02
  50    continue
  60  continue
# endif
#else
c  Linear equation of state. The user needs to select the appropriate
c  values for the mean temperature TBAR, thermal expansion coefficient
c  TCOEF, mean salinity SBAR, and saline contraction coefficeint SCOEF.
c  These values are specified in common block RHOMEAN.
c
      do 10 k=1,km
      do 10 i=1,imt
        rho(i,k)=-tcoef*(tx(i,k)-tbar)+scoef*((sx(i,k)+smean)-sbar)
  10  continue
#endif
      return
      end
