;library of oceanography functions for NCL
;many functions taken from the CSIRO SEAWATER matlab library
;see  Fofonoff, P. & Millard, R.C. Unesco 1983. Algorithms for computation of fundamental properties of seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44.



 ;constants
 pi = 3.14159
 Re = 6371220.0  ;mean radius of earth (m)
 g  = 9.80665    ;standard gravity



;list of functions

;fm2meter(z) - convert dept in fathoms to meters
;pres2depth(P,lat) - convert pressure (db) to depth (m)
;depth2pres(z,lat) - convert depth (m) to pressure (db)
;dT_adiab_sw - calculate adiabatic temperature gradient
;theta_sw(T,S,P,Pr) - calculate potential temp. for sea water
;cp_sw - calculate constant pressure specific heat
;sw_smow - density of Standard Mean Ocean Water
;sw_dens0 - calculate seawater density at atmos. surface pressure
;sw_seck - calculate Secant Bulk Modulus (K) of seawater
;sw_dens(T,S,P) - calculate seawater density
;sw_svan - calc specific volume anomaly (only use if you don't already have density)





;======================
;fm2meter
;=======================

;convert dept in fathoms to meters
 undef("fm2m")

 procedure fm2m (Z:numeric)

 begin
   
   Z = Z * 1.822
 end

;====================================
;pres2depth
;====================================

;convert pressure (db) to depth (m)
;(see Saunders, P.M., 1981. Practical conversion of pressure to depth. J. Phys. Ocean, 11, 573-574.

  undef("pres2depth")
  
  function pres2depth (P:numeric, lat:numeric)

;p   - depth (m)
;lat - latitude (deg. North); must be same size as p

  local dimd, diml,z, c1, c2, c3, c4, d2rad, gam, y, rad, gy, z, bline, tline

  begin

    ;check array sizes match
    
    dimd = dimsizes(P)
    diml = dimsizes(lat)


    if (dimsizes(dimd).ne.dimsizes(diml).or.any(dimd.ne.diml) ) then
      
      print("Error. ocean_funcs.ncl: pres2depth")
      print("P and lat arrays must be same size")
      exit
    end if

    d2rad = pi/180.
    c1    = 9.72659
    c2    = -2.2512e-5
    c3    = 2.279e-10
    c4    = -1.82e-15
    gam   = 2.184e-6

    
    y     = abs(lat)
    rad   = sin(d2rad*y)
    rad   = rad^2

    gy    = 9.780318 * (1. + (rad*5.2788e-3) + (rad^2*2.36e-5) )
   

    ;calc. depth
    bline =gy + (gam*.5*P)
    tline = (c1*P) + (c2*P^2) + (c3*P^3) + (c4*P^4)

    z = tline/bline
   

    return(z)


  end

;====================================
;depth2pres
;====================================

;convert depth (m) to pressure (db)
;(see Saunders, P.M., 1981. Practical conversion of pressure to depth. J. Phys. Ocean, 11, 573-574.

  undef("depth2pres")
  
  function depth2pres (d:numeric, lat:numeric)

;d   - depth (m)
;lat - latitude (deg. North); must be same size as d

  local dimd, diml, c1, c2, Y, deg2rad, p

  begin

    ;check array sizes match
    
    dimd = dimsizes(d)
    diml = dimsizes(lat)


    if (dimsizes(dimd).ne.dimsizes(diml).or.any(dimd.ne.diml) ) then
      
      print("Error. ocean_funcs.ncl: depth2pres")
      print("d and lat arrays must be same size")
      exit
    end if

    
    c2      = 2.21e-6

    ;calc constant c1
    deg2rad = pi/180.
    Y       = sin(abs(lat)*deg2rad)
    c1      = (5.92 + (5.25 * Y^2.)) * 1.e-3

    ;calc. pressure


    p       = ( (1-c1) - sqrt( (1-c1)^2 - (4*c2*d) )  )/(2 * c2)


    return(p)


  end


;=========================================
;dT_adiab_sw
;==========================================

;calculate adiabatic temperature gradient (viz Fofonoff & Millard, 1983)

  undef("dT_adiab_sw")

  function dT_adiab_sw(T:numeric, S:numeric, P:numeric)

;T - temperature (C)
;S - salinity (psu)
;P - pressure (db)
;all three arrays must have the same dimensions


  local dimp, dimt, dims, T68, a0, a1, a2, a3, b0, b1, c0, c1, c2, c3, e0, e1, e2, out

  begin

    ;check array sizes
    dimp = dimsizes(P)
    dimt = dimsizes(T)
    dims = dimsizes(S)
    

    if (dimsizes(dimp).ne.dimsizes(dimt).or.any(dimp.ne.dimt).or.dimsizes(dims).ne.dimsizes(dimt).or.any(dims.ne.dimt) ) then
      
      print("Error. ocean_funcs.ncl: dT_adiab_sw")
      print("T, S and P arrays must be same size")
      exit
    end if


    
    ;constants
    a0 =  3.5803E-5
    a1 =  8.5258E-6
    a2 = -6.836E-8
    a3 =  6.6228E-10

    b0 =  1.8932E-6
    b1 = -4.2393E-8

    c0 =  1.8741E-8
    c1 = -6.7795E-10
    c2 =  8.733E-12
    c3 = -5.4481E-14

    d0 = -1.1351E-10
    d1 =  2.7759E-12

    e0 = -4.6206E-13
    e1 =  1.8676E-14
    e2 = -2.1687E-16

    T68 = T * 1.00024  ;convert to 1968 temperature scale


    out = a0 + (a1 + (a2 + a3*T68)*T68)*T68 + (b0 + b1*T68)*(S-35) + ( (c0 + (c1 + (c2 + c3*T68)*T68)*T68) + (d0 + d1*T68)*(S-35) )*P + (  e0 + (e1 + e2*T68)*T68 )*P*P

    return(out)
         

  end


  

;==============================================
;theta_sw
;===============================================

;calculate potential temp. for sea water, from salinity/pressure/temp.



  undef("theta_sw")

  function theta_sw( T:numeric, S:numeric, P:numeric, Pr:numeric)

;T  - temperature (Celsius)
;S  - salinity (psu): must be same size as T and P
;P  - pressure (db) : must be same size as T and S
;Pr - reference pressure (db), scalar OR same size as T and P

  local c68, dP, dth, th, q, dimpr, dimt


  begin

    c68  = 1.00024   ;conversion constant to 1968 T scale

    ;array size check is done by function dT_adiab_sw
    dimpr = dimsizes(Pr)
    
    ;reference pressure dimensions
    if (dimsizes(dimpr).ne.1.and.dimpr(0).ne.1) then
      dimt = dimsizes(T)
      if (dimsizes(dimpr).ne.dimsizes(dimt).or.any(dimpr.ne.dimt) ) then
      
        print("Error. ocean_funcs.ncl: theta_sw")
        print("Pr array must scalar or same dimensions as T,S,P")
        exit
      end if
    end if
    

    dP = Pr - P  ;pressure difference

    
    ;1st iteration
    dth = dP * dT_adiab_sw(T, S, P)
    th  = (T * c68) + (0.5 * dth)       
    q   = dth


    ; 2nd interation
    dth = dP * dT_adiab_sw(th/c68, S, (P+ (0.5*dP)) )
    
    th  = th + ( (1- (1/sqrt(2)) )*(dth-q) )
    
    q   = ( (2-sqrt(2))*dth) + ( ( (3/sqrt(2)) - 2) * q )
    


    ;3rd iteration
    dth = dP * dT_adiab_sw(th/c68, S, (P + (0.5*dP)) )

    th  = th + ( (1 + (1/sqrt(2)) )*(dth-q) )
    
    q   = ( (2+ sqrt(2))*dth) + ( ((-3/sqrt(2)) - 2) * q )


    ;4th interation
    dth = dP * dT_adiab_sw(th/c68, S, (P + dP) )

    th  = ( th + (dth - (2*q))/6 )/ c68
    
    return(th)

  end

;********************************
;cp_sw
;*******************************


 ;calculate constant pressure specific heat (cp) for seawater, from T, P and S 

  undef("cp_sw")

  function cp_sw( T:numeric, S:numeric, P:numeric)

;T  - temperature (Celsius)
;S  - salinity (psu): must be same size as T and P
;P  - pressure (db) : must be scalar or same size as T and S


  local dimp, dimt, dims, T1,T2, T3, T4, c0,c1,c2,c3,c4,a0,a1,a2,a3,a4,b0,b1,b2,b3,b4,d0,d1,d2,d3,d4,e0,e1,e2,f0,f1,f2,f3,g0,h0,h1,h2,j1, A, B, cp_st0, cp_0t0, d1_cp, d2_cp, cp, Pbar


  begin

 ;check array sizes

    dimp = dimsizes(P)
    dimt = dimsizes(T)
    dims = dimsizes(S)

    if (product(dimp).ne.1.and.( dimsizes(dimp).ne.dimsizes(dimt).or.any(dimp.ne.dimt))) then
      print("Error. ocean_funcs.ncl: cp_sw")
      print("P must be scalar or same size as T, S arrays")
      exit
    end if


    if (dimsizes(dims).ne.dimsizes(dimt).or.any(dims.ne.dimt)) then
      print("Error. ocean_funcs.ncl: cp_sw")
      print("T, S arrays must be same size")
      exit
    end if


;check valid ranges

    if (any(S.gt.40..or.S.lt.0.)) then
      print("Warning: ocean_funcs.ncl: cp_sw")
      print("S is outside valid range of 0-40")
    end if

    if (any(T.gt.35..or.T.lt.0.)) then
      print("Warning: ocean_funcs.ncl: cp_sw")
      print("T is outside valid range of 0-35C")
    end if

;convert P from dB to B, and convert T to 1968 T-scale

    Pbar = P/10.
    T1   = T * 1.00024


;specific heat at P = 0


    ;temperature powers
    T2 = T1^2
    T3 = T1^3
    T4 = T1^4


    ;empirical constants
    c0 = 4217.4
    c1 = -3.720283
    c2 = 0.1412855
    c3 = -2.654387e-3
    c4 = 2.093236e-5

    a0 = -7.643575
    a1 = 0.1072763
    a2 = -1.38385e-3

    b0 = 0.1770383
    b1 = -4.07718e-3
    b2 = 5.148e-5


    
    cp_0t0 = c0 +(c1*T1) + (c2*T2) + (c3*T3) + (c4*T4)
    
    A      = a0 + (a1*T1) + (a2*T2)

    B      = b0 + (b1*T1) + (b2*T2)

    cp_st0 = cp_0t0 + (A*S) + (B*S^1.5)




;pressure dependance

    a0 = -4.9592e-1
    a1 =  1.45747e-2
    a2 = -3.13885e-4
    a3 =  2.0357e-6
    a4 =  1.7168e-8

    b0 =  2.4931e-4
    b1 = -1.08645e-5
    b2 =  2.87533e-7
    b3 = -4.0027e-9
    b4 =  2.2956e-11

    c0 = -5.422e-8
    c1 =  2.6380e-9
    c2 = -6.5637e-11
    c3 =  6.136e-13


    d1_cp = (Pbar * (a0 + (a1*T1) + (a2*T2) + (a3*T3) + (a4*T4) ) ) + ( Pbar^2 * (b0 + (b1*T1) + (b2*T2) + (b3*T3) + (b4*T4) )) + (Pbar^3 * (c0 + (c1*T1) + (c2*T2) + (c3*T3) ))


    d0 =  4.9247e-3
    d1 = -1.28315e-4
    d2 =  9.802e-7
    d3 =  2.5941e-8
    d4 = -2.9179e-10

    e0 = -1.2331e-4
    e1 = -1.517e-6
    e2 =  3.122e-8

    f0 = -2.9558e-6
    f1 =  1.17054e-7
    f2 = -2.3905e-9
    f3 =  1.8448e-11

    g0 =  9.971e-8

    h0 =  5.540e-10
    h1 = -1.7682e-11
    h2 =  3.513e-13

    j1 = -1.4300e-12
    

    d2_cp = Pbar *( (S *(d0 +(d1*T1)+(d2*T2)+(d3*T3)+(d4*T4)))  + (S^1.5 *(e0 +(e1*T1) + (e2*T2))))  + \
 (Pbar^2 *( (S *(f0+(f1*T1)+(f2*T2)+(f3*T3))) + (g0*S^1.5))  ) + \
(Pbar^3 * ( (S *(h0+(h1*T1)+(h2*T2))) +(j1*T1*S^1.5)) )



    cp = cp_st0 + d1_cp + d2_cp


    return(cp)

   end

;*************************
;sw_smow
;**************************

 ;density of Standard Mean Ocean Water (pure water)

  undef("sw_smow")

  function sw_smow( T:numeric )

;T  - temperature (Celsius)


  local a0, a1, a2, a3, a4, a5, T68, dens

  begin
  
  ;coefficients
    a0 = 999.842594
    a1 = 6.793952e-2
    a2 = -9.095290e-3
    a3 = 1.001685e-4
    a4 = -1.120083e-6
    a5 = 6.536332e-9

    T68 = T*1.00024

    dens = a0 +(a1*T68) + (a2*T68^2) + (a3*T68^3) + (a4*T68^4) + (a5*T68^5)

    return(dens)

  end

;**************************
;sw_dens0
;**********************

 ;calculate seawater density at atmos. surface pressure


  undef("sw_dens0")

  function sw_dens0( T:numeric, S:numeric)

;T  - temperature (Celsius)
;S  - salinity (psu): must be same size as T 


  local dims,dimt, T68, b0, b1, b2, b3, b4, c0, c1, c2, d0, dens


;check dimension sizes

  begin

    dims = dimsizes(S)
    dimt = dimsizes(T)

    if (dimsizes(dims).ne.dimsizes(dimt).or.any(dims.ne.dimt)) then
      print("Error. ocean_funcs.ncl: sw_dens0")
      print("T, S arrays must be same size")
      exit
    end if


;constants

  b0 = 8.24493e-1
  b1 = -4.0899e-3
  b2 = 7.6438e-5
  b3 = -8.2467e-7
  b4 = 5.3875e-9

  c0 = -5.72466e-3
  c1 = 1.0227e-4
  c2 = -1.6546e-6

  d0 = 4.8314e-4


  T68 = T * 1.00024
  
  dens = S*(b0 + (b1*T68)+(b2*T68^2)+(b3*T68^3)+(b4*T68^4) ) + S^1.5*(c0 + (c1*T68)+(c2*T68^2) ) + (d0 * S^2)

  dens = dens + sw_smow(T68)

  return(dens)

  end


;**************************
;sw_seck
;************************

 ;calculate Secant Bulk Modulus (K) of seawater

  undef("sw_seck")

  function sw_seck( T:numeric, S:numeric, P:numeric)


;T  - temperature (Celsius)
;S  - salinity (psu): must be same size as T 
;P  - Pressure (db): must be same size as T

  local dims, dimt, dimp, h0, h1, h2, h3, T68, AW, k0, k1, k2, BW, e0, e1, e2, e3, e4, KW, j0, i0, i1, i2, A, m0, m1, m2, f0, f1, f2, f3, g0, g1, g2, B, K0, K, Patm

  begin


    dims = dimsizes(S)
    dimt = dimsizes(T)


    ;check dimension sizes

    if (dimsizes(dims).ne.dimsizes(dimt).or.any(dims.ne.dimt)) then
      print("Error. ocean_funcs.ncl: sw_seck")
      print("T, S arrays must be same size")
      exit
    end if


    dimp = dimsizes(P)

    if (dimsizes(dimp).ne.dimsizes(dimt).or.any(dimp.ne.dimt)) then
      print("Error. ocean_funcs.ncl: sw_seck")
      print("T, P arrays must be same size")
      exit
    end if


    ;compression terms

    T68  = T * 1.00024
    Patm = P/10. ;vonverty to mb

    h3 = -5.77905E-7
    h2 = 1.16092E-4
    h1 = 1.43713E-3
    h0 = 3.239908


    AW  = h0 + (h1*T68) + (h2*T68^2) + (h3*T68^3)


    k2 =  5.2787E-8
    k1 = -6.12293E-6
    k0 =  8.50935E-5

    BW  = k0 + (k1 + k2*T68)*T68

    e4 = -5.155288E-5
    e3 = 1.360477E-2
    e2 = -2.327105
    e1 = 148.4206
    e0 = 19652.21

    KW  = e0 + (e1 + (e2 + (e3 + e4*T68)*T68)*T68) *T68 


    ;K at stmos. pressure


    j0 = 1.91075E-4

    i2 = -1.6078E-6
    i1 = -1.0981E-5
    i0 =  2.2838E-3


    A  = AW + S*(i0 + (i1*T68) + (i2*T68^2) ) + (j0*S^1.5)


    m2 =  9.1697E-10
    m1 = 2.0816E-8
    m0 = -9.9348E-7

    B = BW + (m0 + (m1*T68) + (m2*T68^2) )*S  ;   eqn 18

    f3 =  -6.1670E-5
    f2 =  1.09987E-2
    f1 =  -0.603459
    f0 = 54.6746

    g2 = -5.3009E-4
    g1 = 1.6483E-2
    g0 = 7.944E-2

    K0 = KW + S*(f0 + (f1*T68)+ (f2*T68^2) + (f3*T68^3) )  + S^1.5*(g0 + (g1*T68) + (g2*T68^2) )   ;eqn 16



  ; K at T, S, P
    K = K0 + (A*Patm) + (B*Patm^2)   ; eqn 15
    return(K)

  end

;************************
;sw_dens
;**********************


 ;calculate seawater density from T, S, P

  undef("sw_dens")

  function sw_dens( T:numeric, S:numeric, P:numeric)


;T  - temperature (Celsius)
;S  - salinity (psu): must be same size as T 
;P  - Pressure (db): must be same size as T

  local dims, dimt, dimp, dens0, K, Patm, dens

  begin


    dims = dimsizes(S)
    dimt = dimsizes(T)


    ;check dimension sizes

    if (dimsizes(dims).ne.dimsizes(dimt).or.any(dims.ne.dimt)) then
      print("Error. ocean_funcs.ncl: sw_dens")
      print("T, S arrays must be same size")
      exit
    end if


    dimp = dimsizes(P)

    if (dimsizes(dimp).ne.dimsizes(dimt).or.any(dimp.ne.dimt)) then
      print("Error. ocean_funcs.ncl: sw_dens")
      print("T, P arrays must be same size")
      exit
    end if

;check valid ranges

    if (any(S.gt.42..or.S.lt.0.)) then
      print("Warning: ocean_funcs.ncl: sw_dens")
      print("S is outside valid range of 0-42")
    end if

    if (any(T.gt.40..or.T.lt.-2.)) then
      print("Warning: ocean_funcs.ncl: sw_dens")
      print("T is outside valid range of -2 to 40C")
    end if

    if (any(P.gt.10000..or.P.lt.0.)) then
      print("Warning: ocean_funcs.ncl: sw_dens")
      print("P is outside valid range of 0 to 10000db")
    end if



    dens0    = sw_dens0(T,S)
    
    K         = sw_seck(T,S,P)

    Patm     = P/10.  ;convert db to mb


    dens       = dens0/(1-Patm/K)

    return(dens)


    end 
;*******************************
;sw_svan
;****************************

    undef("sw_svan")
    function sw_svan(T:numeric,S:numeric,P:numeric)


;calc specific volume anomaly

;T  - temperature (Celsius)
;S  - salinity (psu): must be same size as T 
;P  - Pressure (db): must be same size as T

  local dims, dimt, dimp, dens0, dens, svan

  begin


    dims = dimsizes(S)
    dimt = dimsizes(T)


    ;check dimension sizes

    if (dimsizes(dims).ne.dimsizes(dimt).or.any(dims.ne.dimt)) then
      print("Error. ocean_funcs.ncl: sw_svan")
      print("T, S arrays must be same size")
      exit
    end if


    dimp = dimsizes(P)

    if (dimsizes(dimp).ne.dimsizes(dimt).or.any(dimp.ne.dimt)) then
      print("Error. ocean_funcs.ncl: sw_svan")
      print("T, P arrays must be same size")
      exit
    end if



    rho     = sw_dens(T, S, P)
    rho0   = sw_dens(conform(P,0.,-1), conform(P,35.,-1), P)


    svan   = (1/rho) - (1/rho0)


    return(svan)

  end