      real function rho_wright(depth_cm,theta,S_moma)
c ----------------------------------------------------------------
c     from Wright, 1997. JOAT, pp735-740. Uses Eckart formula, gives
c     rho in terms of pot'l temperature
c   here we use 'extended' formula, valid over reduced range, [col 4, p736]
c   p < 5.10^7 Pa, -2 < theta < 30., 28 <S <38
c
c     inputs: depth_cm    depth in cm 
c             theta   pot'l temp in deg
c             S_moma  salinity in moma units [(S-35)*10^-3 ]
c     returns wright  density in g/cm^3
c
c     working of function is in MKS units ..
c ----------------------------------------------------------------
      implicit none
      real depth_cm, theta, S_moma
Cf2py intent (in) depth_cm
Cf2py intent (in) theta
Cf2py intent (in) S_moma


#ifdef linear_rho
#include "lin_dncoef.h"
      rho_wright = dRhodT*(theta - T0)
#ifndef no_salt 
     &          + dRhods*(S_moma - s0)
#endif /* .not. no_salt */



#else /* .not. linear_rho */
      real a(0:2), b(0:5), c(0:5)
      data a/7.057924d-4, 3.480336d-7, -1.112733d-7/
      data b/5.790749d8, 3.516535d6, -4.002714d4
     &      , 2.084372d2, 5.944068d5, -9.643486d3/
      data c/1.704853d5, 7.904722d2, -7.984422d0
     &      ,5.140652d-2, -2.302158d2, -3.079464 /
      save a,b,c

      real gravRho_d100,rho00
      parameter (gravRho_d100 = 100.0,rho00=1.0)

      real p,S,alpha0,p0,lambda

      p = gravRho_d100*depth_cm
      S = S_moma*1000. + 35.

      alpha0 = a(0) + a(1)*theta + a(2)*S
      p0 = b(0) +b(4)*S +
     &      theta*(b(1) + b(5)*S + theta*(b(2) + theta*b(3)))
      lambda = c(0) + c(4)*S +
     &      theta*(c(1) + c(5)*S + theta*(c(2) + theta*c(3)))

      rho_wright = 1.d-3*(p+p0)/(lambda + alpha0*(p+p0)) - rho00
#endif /* .not. linear_rho */
      return
      end

      subroutine drho_wright(depth_cm,theta,S_moma,drho_W_dT,ndyn)
c ----------------------------------------------------------------
c     from Wright, 1997. JOAT, pp735-740. Uses Eckart formula, gives
c     rho in terms of pot'l temperature
c   here we use 'extended' formula, valid over reduced range, [col 4, p736]
c   p < 5.10^7 Pa, -2 < theta < 30., 28 <S <38
c
c     inputs: depth_cm    depth in cm 
c             theta   pot'l temp in deg
c             S_moma  salinity in moma units [(S-35)*10^-3 ]
c             ndyn  = 2 unless rho function of theta only
c     returns drho_W_dT(ndyn)  drhodt in g/(cm^3 deg), [drhods in g/(cm^3 frac salinity)]
c
c     working of function is in MKS units ..
c ----------------------------------------------------------------
      implicit none
      integer ndyn
      real depth_cm, theta, S_moma,drho_W_dT(ndyn)
Cf2py intent (in) ndyn
Cf2py intent (in) depth_cm
Cf2py intent (in) theta
Cf2py intent (in) S_moma
Cf2py intent (out) drho_W_dT


#ifdef linear_rho
# include "lin_dncoef.h"
      drho_W_dT(1) = dRhodT
# ifndef no_salt 
      drho_W_dT(2) = dRhods
# endif /* .not. no_salt */



#else /* .not. linear_rho */
      real a(0:2), b(0:5), c(0:5)
      real scale(2)
      data a/7.057924d-4, 3.480336d-7, -1.112733d-7/
      data b/5.790749d8, 3.516535d6, -4.002714d4
     &      , 2.084372d2, 5.944068d5, -9.643486d3/
      data c/1.704853d5, 7.904722d2, -7.984422d0
     &      ,5.140652d-2, -2.302158d2, -3.079464d0 /
      data scale/1.0d0, 1.0e3/
      save a,b,c,scale

      real gravRho_d100,rho00,c3,c2
      parameter (gravRho_d100 = 100.0,rho00=1.0,c3=3.0,c2=2.0)

      integer m
      real p,S,alpha0,p0,lambda
      real rdenom, dp0dT(2), dldT(2)

      p = gravRho_d100*depth_cm
      S = S_moma*1000. + 35.

      alpha0 = a(0) + a(1)*theta + a(2)*S
      p0 = b(0) + b(4)*S +
     &      theta*(b(1) + b(5)*S + theta*(b(2) + theta*b(3)))
      lambda = c(0) + c(4)*S +
     &      theta*(c(1) + c(5)*S + theta*(c(2) + theta*c(3)))
      
      dp0dT(1) = b(1) + theta*(c2*b(2) + c3*b(3)*theta) + S*b(5)
      dp0dT(2) = b(4) + b(5)*theta

      dldT(1) = c(1) + theta*(c2*c(2) + c3*c(3)*theta) + S*c(5)
      dldT(2) = c(4) + c(5)*theta

      rdenom = 1./( 1.e3*(lambda + alpha0*(p+p0))**2 )

c   scale multiplies salinity gradient by 1000. so we have /frac S
c    rather than /mil

      do m=1,ndyn
         drho_W_dT(m) = rdenom*scale(m)*
     &         (lambda*dp0dT(m) - (p+p0)*(dldT(m) + (p+p0)*a(m)) )
         end do

#endif /* .not. linear_rho */
      return
      end


      real function dP_wright(depth_cm_top,depth_cm_below,theta,S_moma)
c ----------------------------------------------------------------
c     finds pressure diff dP_wright, [actually int rho_dz] in cgs [g.cm^-2]
c     based on formula for rho from Wright, 1997. JOAT, pp735-740.
c
c   here we use 'extended' formula, valid over reduced range, [col 4, p736]
c   p < 5.10^7 Pa, -2 < theta < 30., 28 <S <38
c
c     inputs: depth_cm_top     depth in cm of upper level 
c             depth_cm_below   depth in cm of lower level 
c             theta            pot'l temp in deg
c             S_moma           salinity in moma units [(S-35)*10^-3 ]
c     returns dP_wright        int density dz in g/cm^2
c
c     working of function is in MKS units ..
c ----------------------------------------------------------------
      implicit none
      real depth_cm_top, depth_cm_below, theta, S_moma
Cf2py intent (in) depth_cm_top
Cf2py intent (in) depth_cm_below
Cf2py intent (in) theta
Cf2py intent (in) S_moma


#ifdef linear_rho
#include "lin_dncoef.h"
      dP_wright = (depth_cm_below - depth_cm_top)*(
     &              dRhodT*(theta - T0)
#ifndef no_salt 
     &             + dRhods*(S_moma - s0)
#endif /* .not. no_salt */
     &         )



#else /* .not. linear_rho */
      real a(0:2), b(0:5), c(0:5)
      data a/7.057924d-4, 3.480336d-7, -1.112733d-7/
      data b/5.790749d8, 3.516535d6, -4.002714d4
     &      , 2.084372d2, 5.944068d5, -9.643486d3/
      data c/1.704853d5, 7.904722d2, -7.984422d0
     &      ,5.140652d-2, -2.302158d2, -3.079464 /
      save a,b,c

      real gravRho_d100,rho00
      parameter (gravRho_d100 = 100.0,rho00=1.0)

      real p_top,p_below,S,alpha0,p0,lambda,l

      p_top = gravRho_d100*depth_cm_top
      p_below = gravRho_d100*depth_cm_below
      S = S_moma*1000. + 35.

      alpha0 = a(0) + a(1)*theta + a(2)*S
      p0 = b(0) +b(4)*S +
     &      theta*(b(1) + b(5)*S + theta*(b(2) + theta*b(3)))
      lambda = c(0) + c(4)*S +
     &      theta*(c(1) + c(5)*S + theta*(c(2) + theta*c(3)))

      l = lambda/alpha0
      dP_wright = (1.d-3/(alpha0*gravRho_d100))*
     &        (p_below - p_top - l*log((l+p0+p_below)/(l+p0+p_top)))
     &  -  rho00*(depth_cm_below - depth_cm_top)        

#endif /* .not. linear_rho */
      return
      end

      subroutine drho_wright_n(depth_cm,theta,S_moma
     &      ,drho_W_dT,drho_W_ds,klength)
c ----------------------------------------------------------------
c     from Wright, 1997. JOAT, pp735-740. Uses Eckart formula, gives
c     rho in terms of pot'l temperature
c   here we use 'extended' formula, valid over reduced range, [col 4, p736]
c   p < 5.10^7 Pa, -2 < theta < 30., 28 <S <38
c
c     inputs: depth_cm    depth in cm 
c             theta   pot'l temp in deg
c             S_moma  salinity in moma units [(S-35)*10^-3 ]
c     returns drho/dt, drhods in g/cm^3 per degree & per *frac* salinity
c
c     working of function is in MKS units ..
c ----------------------------------------------------------------
      implicit none
      integer k,klength
      real depth_cm(klength), theta(klength), S_moma(klength)
     &      ,drho_W_dT(klength),drho_W_ds(klength)
Cf2py intent (in) klength
Cf2py intent (in) depth_cm
Cf2py intent (in) theta
Cf2py intent (in) S_moma
Cf2py intent (out) drho_W_dT,drho_W_ds


#ifdef linear_rho
#include "lin_dncoef.h"
      do k=1,klength
         drho_W_dT(k) = dRhodT
#ifndef no_salt 
         drho_W_dT(k) = dRhods
#endif /* .not. no_salt */
         end do


#else /* .not. linear_rho */
      real a(0:2), b(0:5), c(0:5)
      data a/7.057924d-4, 3.480336d-7, -1.112733d-7/
      data b/5.790749d8, 3.516535d6, -4.002714d4
     &      , 2.084372d2, 5.944068d5, -9.643486d3/
      data c/1.704853d5, 7.904722d2, -7.984422d0
     &      ,5.140652d-2, -2.302158d2, -3.079464d0 /
      save a,b,c

      real gravRho_d100,rho00,c3,c2
      parameter (gravRho_d100 = 100.0,rho00=1.0,c3=3.0,c2=2.0)

      real p,S,T,alpha0,p0,lambda
      real rdenom, dp0dT(2), dldT(2)

      do k=1,klength
         p = gravRho_d100*depth_cm(k)
         S = S_moma(k)*1000. + 35.
         T = theta(k)

       	 alpha0 = a(0) + a(1)*T + a(2)*S
       	 p0 = b(0) + b(4)*S +
     & 	      T*(b(1) + b(5)*S + T*(b(2) + T*b(3)))
       	 lambda = c(0) + c(4)*S +
     & 	      T*(c(1) + c(5)*S + T*(c(2) + T*c(3)))
       	
       	 dp0dT(1) = b(1) + T*(c2*b(2) + c3*b(3)*T) + S*b(5)
       	 dp0dT(2) = b(4) + b(5)*T
	 
       	 dldT(1) = c(1) + T*(c2*c(2) + c3*c(3)*T) + S*c(5)
       	 dldT(2) = c(4) + c(5)*T
	 
c       factor 1.e3 converts from MKS to cgs .
       	 rdenom = 1./( 1.e3*(lambda + alpha0*(p+p0))**2 )
	 
	 
       	 drho_W_dT(k) = rdenom*
     & 	         (lambda*dp0dT(1) - (p+p0)*(dldT(1) + (p+p0)*a(1)) )
#ifndef no_salt 
c       factor 1.e3 converts from per mil to per unit of fractional S
       	 drho_W_ds(k) = 1.e3*rdenom*
     & 	         (lambda*dp0dT(2) - (p+p0)*(dldT(2) + (p+p0)*a(2)) )
#endif /* .not. no_salt */
         end do

#endif /* .not. linear_rho */
      return
      end

      subroutine wright_sub(depth_cm,theta,S_moma, rho, n)
c ----------------------------------------------------------------
c     from Wright, 1997. JOAT, pp735-740. Uses Eckart formula, gives
c     rho in terms of pot'l temperature
c   here we use 'extended' formula, valid over reduced range, [col 4, p736]
c   p < 5.10^7 Pa, -2 < theta < 30., 28 <S <38
c
c     inputs: depth_cm    depth in cm 
c             theta   pot'l temp in deg
c             S_moma  salinity in moma units [(S-35)*10^-3 ]
c     returns wright  density in g/cm^3
c
c     working of function is in MKS units ..
c ----------------------------------------------------------------
      implicit none
      integer n
      real depth_cm(n), theta(n), S_moma(n), rho(n)
Cf2py intent (in) n
Cf2py intent (in) depth_cm
Cf2py intent (in) theta
Cf2py intent (in) S_moma
Cf2py intent (out) rho


#ifdef linear_rho
#include "lin_dncoef.h"
      integer i
      do i=1,n
         rho(i) = dRhodT*(theta(i) - T0)
#ifndef no_salt 
     &          + dRhods*(S_moma(i) - s0)
#endif /* .not. no_salt */
         end do


#else /* .not. linear_rho */
      real a(0:2), b(0:5), c(0:5)
      data a/7.057924d-4, 3.480336d-7, -1.112733d-7/
      data b/5.790749d8, 3.516535d6, -4.002714d4
     &      , 2.084372d2, 5.944068d5, -9.643486d3/
      data c/1.704853d5, 7.904722d2, -7.984422d0
     &      ,5.140652d-2, -2.302158d2, -3.079464 /
      save a,b,c

      real gravRho_d100,rho00
      parameter (gravRho_d100 = 100.0,rho00=1.0)
      integer nmax,i
      parameter (nmax = 1000)

      real p,S,alpha0,p0,lambda

      do i=1,n
         p = gravRho_d100*depth_cm(i)
         S = S_moma(i)*1000. + 35.
   
         alpha0 = a(0) + a(1)*theta(i) + a(2)*S
         p0 = b(0) +b(4)*S +
     &    theta(i)*(b(1) + b(5)*S + theta(i)*(b(2) + theta(i)*b(3)))
         lambda = c(0) + c(4)*S +
     &    theta(i)*(c(1) + c(5)*S + theta(i)*(c(2) + theta(i)*c(3)))
   
         rho(i) =1.d-3*(p+p0)/(lambda + alpha0*(p+p0)) - rho00
         end do
#endif /* .not. linear_rho */
      return
      end

