;**************************************************************************
;**************************************************************************
undef("rho_alphabeta")
function rho_alphabeta(t2d[*][*]:numeric,s2d[*][*]:numeric,depth:numeric,iopt:integer) 

;-- based on Steve Yeager's IDL routine rhoalphabeta
;-- which in turn is based on POP state_mod.F (ccsm3_0_beta22) for 'mwjf'
;   iopt = 0  returns alpha
;   iopt \= 0  returns beta
;=========================================================================

local dims,nx,ny,c1,c1p5,c2,c3,c4,c5,c10,c1000,p001,mwjfnp0s0t0,mwjfnp0s0t1,\
      mwjfnp0s0t2,mwjfnp0s0t3,mwjfnp0s1t0,mwjfnp0s1t1,mwjfnp0s2t0,mwjfnp1s0t0,\
      mwjfnp1s0t2,mwjfnp1s1t0,mwjfnp2s0t0,mwjfnp2s0t2,mwjfdp0s0t0,mwjfdp0s0t1,\
      mwjfdp0s0t2,mwjfdp0s0t3,mwjfdp0s0t4,mwjfdp0s1t0,mwjfdp0s1t1,mwjfdp0s1t3,\
      mwjfdp0sqt0,mwjfdp0sqt2,mwjfdp1s0t0,mwjfdp2s0t3,mwjfdp3s0t1,pressure,p,\
      sqr,work1,work2,rhofull
begin

;========= define rho
  rhoout = new(dimsizes(t2d),typeof(t2d))

;========== define constants
  c1    = 1.0
  c1p5  = 1.5
  c2    = 2.0
  c3    = 3.0
  c4    = 4.0
  c5    = 5.0
  c10   = 10.0
  c1000 = 1000.0
      ;*** these constants will be used to construct the numerator
      ;*** factor unit change (kg/m^3 -> g/cm^3) into numerator terms
  p001 = 0.001
  mwjfnp0s0t0 =   9.99843699e+2 * p001
  mwjfnp0s0t1 =   7.35212840e+0 * p001
  mwjfnp0s0t2 =  -5.45928211e-2 * p001
  mwjfnp0s0t3 =   3.98476704e-4 * p001
  mwjfnp0s1t0 =   2.96938239e+0 * p001
  mwjfnp0s1t1 =  -7.23268813e-3 * p001
  mwjfnp0s2t0 =   2.12382341e-3 * p001
  mwjfnp1s0t0 =   1.04004591e-2 * p001
  mwjfnp1s0t2 =   1.03970529e-7 * p001
  mwjfnp1s1t0 =   5.18761880e-6 * p001
  mwjfnp2s0t0 =  -3.24041825e-8 * p001
  mwjfnp2s0t2 =  -1.23869360e-11 * p001
      ;*** these constants will be used to construct the denominator
  mwjfdp0s0t0 =   1.0e+0
  mwjfdp0s0t1 =   7.28606739e-3
  mwjfdp0s0t2 =  -4.60835542e-5
  mwjfdp0s0t3 =   3.68390573e-7
  mwjfdp0s0t4 =   1.80809186e-10
  mwjfdp0s1t0 =   2.14691708e-3
  mwjfdp0s1t1 =  -9.27062484e-6
  mwjfdp0s1t3 =  -1.78343643e-10
  mwjfdp0sqt0 =   4.76534122e-6
  mwjfdp0sqt2 =   1.63410736e-9
  mwjfdp1s0t0 =   5.30848875e-6
  mwjfdp2s0t3 =  -3.03175128e-16
  mwjfdp3s0t1 =  -1.27934137e-17

;=====pressure calculaton
; taken from gokhan's idl pressure.pro for references
;     this function computes pressure in bars from depth in meters
;     by using a mean density derived from depth-dependent global
;     average temperatures and salinities from Levitus_94, and
;     integrating using hydrostatic balance.
;
;     references:
;        Levitus, S., R. Burgett, and T.P. Boyer, World Ocean Atlas
;          1994, Volume 3: Salinity, NOAA Atlas NESDIS 3, US Dept. of
;          Commerce, 1994.
;        Levitus, S. and T.P. Boyer, World Ocean Atlas 1994,
;          Volume 4: Temperature, NOAA Atlas NESDIS 4, US Dept. of
;          Commerce, 1994.
;        Dukowicz, J. K., 2000: Reduction of Pressure and Pressure
;          Gradient Errors in Ocean Simulations, J. Phys. Oceanogr.,
;          submitted.

  if(depth.ne.0) then
    pressure = 0.059808*(exp(-0.025*depth) - 1.0)  + 0.100766*depth + \
               2.28405e-7*(depth^2)  
  else
    pressure = 0.
  end if
 
  p = pressure 
; CAS corrected this; verified w/Steve Yeager 9/8/09
; p = pressure only works for potential density; depth/pressure = 0. 
  p = pressure*c10

;========= compute the numerator of the MWFJ density  [P_1(S,T,p)]  

  mwjfnums0t0 = mwjfnp0s0t0 + p*(mwjfnp1s0t0 + p*mwjfnp2s0t0)
  mwjfnums0t1 = mwjfnp0s0t1
  mwjfnums0t2 = mwjfnp0s0t2 + p*(mwjfnp1s0t2 + p*mwjfnp2s0t2)
  mwjfnums0t3 = mwjfnp0s0t3
  mwjfnums1t0 = mwjfnp0s1t0 + p*mwjfnp1s1t0
  mwjfnums1t1 = mwjfnp0s1t1
  mwjfnums2t0 = mwjfnp0s2t0

  work1 = t2d
  work1 = mwjfnums0t0 + t2d * (mwjfnums0t1 + t2d * (mwjfnums0t2 + \
  mwjfnums0t3 * t2d )) + s2d  * (mwjfnums1t0 + \
  mwjfnums1t1 * t2d + mwjfnums2t0 * s2d)

;============= compute the denominator of  MWJF density [P_2(S,T,p)]

  sqr = sqrt(s2d)

  mwjfdens0t0 = mwjfdp0s0t0 + p*mwjfdp1s0t0
  mwjfdens0t1 = mwjfdp0s0t1 + p^3 * mwjfdp3s0t1
  mwjfdens0t2 = mwjfdp0s0t2
  mwjfdens0t3 = mwjfdp0s0t3 + p^2 * mwjfdp2s0t3
  mwjfdens0t4 = mwjfdp0s0t4
  mwjfdens1t0 = mwjfdp0s1t0
  mwjfdens1t1 = mwjfdp0s1t1
  mwjfdens1t3 = mwjfdp0s1t3
  mwjfdensqt0 = mwjfdp0sqt0
  mwjfdensqt2 = mwjfdp0sqt2

  work2 = t2d 
  work2 = mwjfdens0t0 + t2d * (mwjfdens0t1 + t2d  * (mwjfdens0t2 +   \
          t2d * (mwjfdens0t3 + mwjfdens0t4 * t2d ))) + \
          s2d * (mwjfdens1t0 + t2d * (mwjfdens1t1 + t2d*t2d*mwjfdens1t3) + \
          sqr * (mwjfdensqt0 + t2d*t2d*mwjfdensqt2))    

  denomk = work2
  denomk = c1/work2
        
  rhofull = work1
  rhofull = work1*denomk

  if (iopt.eq.0) then 
     work3 = t2d
     work3 =  mwjfnums0t1 + t2d * (c2*mwjfnums0t2 + \
             c3*mwjfnums0t3 * t2d) + mwjfnums1t1 * s2d
     work4 = s2d
     work4 = mwjfdens0t1 + s2d * mwjfdens1t1 + \
            t2d * (c2*(mwjfdens0t2 + s2d*sqr*mwjfdensqt2) + \
            t2d * (c3*(mwjfdens0t3 + s2d * mwjfdens1t3) + \
            t2d *  c4*mwjfdens0t4))
     drhodx = work4
     drhodx = (work3 - work1*denomk*work4)*denomk
     drhodx = - (c1/rhofull)*drhodx
     drhodx@long_name = "Alpha == -(1/rho)*(drho/dT)"
  else 
     work3 = t2d
     work3 =  mwjfnums1t0 + mwjfnums1t1 * t2d + c2*mwjfnums2t0 * s2d
     work4 = t2d
     work4 = mwjfdens1t0 + t2d * (mwjfdens1t1 +  \
          t2d*t2d*mwjfdens1t3) + c1p5*sqr*(mwjfdensqt0 + t2d*t2d*mwjfdensqt2)
     drhodx = work4
     drhodx =(work3-work1*denomk*work4)*denomk
     drhodx = (c1/rhofull)*drhodx
     drhodx@long_name = "Beta == (1/rho)*(drho/dS)"
   end if

;==== return either alpha or beta
  return (drhodx)   
end