!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

	 module state_ccsm_mod

!BOP
! !MODULE: state_mod
!
! !DESCRIPTION:
!  This module contains routines necessary for computing the density 
!  from model temperature and salinity using an equation of state.
!
!  The module supports four forms of EOS:
!  \begin{enumerate}
!     \item The UNESCO equation of state computed using the 
!           potential-temperature-based bulk modulus from Jackett and 
!           McDougall, JTECH, Vol.12, pp 381-389, April, 1995.
!     \item The faster and more accurate alternative to the UNESCO eos
!           of McDougall, Wright, Jackett and Feistel (hereafter 
!           MWJF, 2001 submission to JTECH).
!     \item a simple linear EOS based on constant expansion coeffs
!  \end{enumerate}
!
! !REVISION HISTORY:
!  CVS:$Id: state_mod.F90,v 1.14 2003/02/11 23:33:15 pwjones Exp $
!  CVS:$Name:  $

! !USES:

	   use kinds_mod
!	   use constants
	   implicit none
	   private
	   save

! !PUBLIC MEMBER FUNCTIONS:

	   public :: state_ccsm,init_state

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  valid ranges and pressure as function of depth
!
!-----------------------------------------------------------------------

	real (r8) :: & 
      	tmin, tmax,        &! valid temperature range
      	smin, smax          ! valid salinity    range

	   real (r8), dimension(:), allocatable :: & 
      	pressz              ! ref pressure (bars) at each level

!-----------------------------------------------------------------------
!
!  choices for eos type and valid range checks
!
!-----------------------------------------------------------------------

	   integer (int_kind), parameter :: &
      state_type_jmcd       = 1,    &! integer ids for state choice
      state_type_mwjf       = 2,    &
      state_type_linear     = 4

	   integer (int_kind) ::    &
      state_itype           ! input state type chosen

!-----------------------------------------------------------------------
!
!  UNESCO EOS constants and JMcD bulk modulus constants
!
!-----------------------------------------------------------------------

   !*** for density of fresh water (standard UNESCO)

	   real (r8), parameter ::           &
      unt0 =   999.842594d0,         &
      unt1 =  6.793952d-2,           &
      unt2 = -9.095290d-3,           &
      unt3 =  1.001685d-4,           &
      unt4 = -1.120083d-6,           &
      unt5 =  6.536332d-9

   !*** for dependence of surface density on salinity (UNESCO)

	   real (r8), parameter ::           &
      uns1t0 =  0.824493d0 ,         &
      uns1t3 = -8.2467d-7,           &
      uns1t1 = -4.0899d-3,           &
      uns1t2 =  7.6438d-5,           &
      uns1t4 =  5.3875d-9,           &
      unsqt0 = -5.72466d-3,          &
      unsqt1 =  1.0227d-4,           &
      unsqt2 = -1.6546d-6,           &
      uns2t0 =  4.8314d-4

   !*** from Table A1 of Jackett and McDougall

	   real (r8), parameter ::           &
      bup0s0t0 =  1.965933d+4,       &
      bup0s0t1 =  1.444304d+2,       &
      bup0s0t2 = -1.706103d0   ,     &
      bup0s0t3 =  9.648704d-3,       &
      bup0s0t4 = -4.190253d-5

	   real (r8), parameter ::           &
      bup0s1t0 =  5.284855d+1,       &
      bup0s1t1 = -3.101089d-1,       &
      bup0s1t2 =  6.283263d-3,       &
      bup0s1t3 = -5.084188d-5

	   real (r8), parameter ::           &
      bup0sqt0 =  3.886640d-1,       &
      bup0sqt1 =  9.085835d-3,       &
      bup0sqt2 = -4.619924d-4

	   real (r8), parameter ::           &
      bup1s0t0 =  3.186519d0   ,     &
      bup1s0t1 =  2.212276d-2,       &
      bup1s0t2 = -2.984642d-4,       &
      bup1s0t3 =  1.956415d-6 

	   real (r8), parameter ::           &
      bup1s1t0 =  6.704388d-3,       &
      bup1s1t1 = -1.847318d-4,       &
      bup1s1t2 =  2.059331d-7,       &
      bup1sqt0 =  1.480266d-4 

	   real (r8), parameter ::           &
      bup2s0t0 =  2.102898d-4,       &
      bup2s0t1 = -1.202016d-5,       &
      bup2s0t2 =  1.394680d-7,       &
      bup2s1t0 = -2.040237d-6,       &
      bup2s1t1 =  6.128773d-8,       &
      bup2s1t2 =  6.207323d-10

!-----------------------------------------------------------------------
!
!  MWJF EOS coefficients
!
!-----------------------------------------------------------------------

   !*** these constants will be used to construct the numerator
   !*** factor unit change (kg/m^3 -> g/cm^3) into numerator terms

	   real (r8), parameter ::                   &
      mwjfnp0s0t0 =   9.99843699d+2 * 1.d-3, &
      mwjfnp0s0t1 =   7.35212840d+0 * 1.d-3, &
      mwjfnp0s0t2 =  -5.45928211d-2 * 1.d-3, &
      mwjfnp0s0t3 =   3.98476704d-4 * 1.d-3, &
      mwjfnp0s1t0 =   2.96938239d+0 * 1.d-3, &
      mwjfnp0s1t1 =  -7.23268813d-3 * 1.d-3, &
      mwjfnp0s2t0 =   2.12382341d-3 * 1.d-3, &
      mwjfnp1s0t0 =   1.04004591d-2 * 1.d-3, &
      mwjfnp1s0t2 =   1.03970529d-7 * 1.d-3, &
      mwjfnp1s1t0 =   5.18761880d-6 * 1.d-3, &
      mwjfnp2s0t0 =  -3.24041825d-8 * 1.d-3, &
      mwjfnp2s0t2 =  -1.23869360d-11* 1.d-3

   !*** these constants will be used to construct the denominator

	   real (kind=r8), parameter ::       &
      mwjfdp0s0t0 =   1.0d0,          &
      mwjfdp0s0t1 =   7.28606739d-3,  &
      mwjfdp0s0t2 =  -4.60835542d-5,  &
      mwjfdp0s0t3 =   3.68390573d-7,  &
      mwjfdp0s0t4 =   1.80809186d-10, &
      mwjfdp0s1t0 =   2.14691708d-3,  &
      mwjfdp0s1t1 =  -9.27062484d-6,  &
      mwjfdp0s1t3 =  -1.78343643d-10, &
      mwjfdp0sqt0 =   4.76534122d-6,  &
      mwjfdp0sqt2 =   1.63410736d-9,  &
      mwjfdp1s0t0 =   5.30848875d-6,  &
      mwjfdp2s0t3 =  -3.03175128d-16, &
      mwjfdp3s0t1 =  -1.27934137d-17

   !*** MWJF numerator coefficients including pressure

	   real (r8) ::                                                        &
      mwjfnums0t0, mwjfnums0t1, mwjfnums0t2, mwjfnums0t3,              &
      mwjfnums1t0, mwjfnums1t1, mwjfnums2t0,                           &
      mwjfdens0t0, mwjfdens0t1, mwjfdens0t2, mwjfdens0t3, mwjfdens0t4, &
      mwjfdens1t0, mwjfdens1t1, mwjfdens1t3,                           &
      mwjfdensqt0, mwjfdensqt2

!-----------------------------------------------------------------------
!
!  parameters for linear eos
!
!-----------------------------------------------------------------------

	   real (r8), parameter ::       & 
      T_leos_ref = 19.0d0,       &! reference T for linear eos (deg C)
      S_leos_ref = 0.035d0,      &! reference S for linear eos (msu)
      rho_leos_ref = 1.025022d0, &! ref dens (g/cm3) at ref T,S and 0 bar
      alf = 2.55d-4,             &! expansion coeff -(drho/dT) (gr/cm^3/K)
      bet = 7.64d-1               ! expansion coeff (drho/dS) (gr/cm^3/msu)

!EOC
!***********************************************************************

	 contains

!***********************************************************************
!BOP
! !IROUTINE: state
! !INTERFACE:

	 subroutine state_ccsm(k, kk, TEMPK, SALTK, &
                  RHOOUT, MASK, fillValue, RHOFULL, DRHODT, DRHODS)

! !DESCRIPTION:
!  Returns the density of water at level k from equation of state
!  $\rho = \rho(d,\theta,S)$ where $d$ is depth, $\theta$ is
!  potential temperature, and $S$ is salinity. the density can be
!  returned as a perturbation (RHOOUT) or as the full density
!  (RHOFULL). Note that only the polynomial EOS choice will return
!  a perturbation density; in other cases the full density is returned
!  regardless of which argument is requested.
!
!  This routine also computes derivatives of density with respect
!  to temperature and salinity at level k from equation of state
!  if requested (ie the optional arguments are present).
!
!  If $k = kk$ are equal the density for level k is returned.
!  If $k \neq kk$ the density returned is that for a parcel
!  adiabatically displaced from level k to level kk.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

	   integer (int_kind), intent(in) :: &
      k,                    &! depth level index
      kk                     ! level to which water is adiabatically 
                             ! displaced

	   real (r4), dimension(:,:), intent(in) :: & 
      TEMPK,             &! temperature (deg C) at level k
      SALTK               ! salinity    (msu)   at level k

	   logical, dimension(:,:), intent(in), optional :: &
      MASK                   ! for masking land points

	   real (r4), intent(in), optional :: &
      fillValue              ! optional value for filling land points
                             !   zero is used if no fill value supplied

! !OUTPUT PARAMETERS:

	   real (r8), dimension(:,:), optional, intent(out) :: & 
      RHOOUT,  &! perturbation density of water
      RHOFULL, &! full density of water
      DRHODT,  &! derivative of density with respect to temperature
      DRHODS    ! derivative of density with respect to salinity

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables:
!
!-----------------------------------------------------------------------

	   real (r8), dimension(size(TEMPK,dim=1),size(TEMPK,dim=2)) :: &
      TQ,SQ,             &! adjusted T,S
      BULK_MOD,          &! Bulk modulus
      RHO_S,             &! density at the surface
      DRDT0,             &! d(density)/d(temperature), for surface
      DRDS0,             &! d(density)/d(salinity   ), for surface
      DKDT,              &! d(bulk modulus)/d(pot. temp.)
      DKDS,              &! d(bulk modulus)/d(salinity  )
      SQR,DENOMK,        &! work arrays
      WORK1, WORK2, WORK3, WORK4, T2

	   real (r8) :: p,p2,fill !&
!      p, p2     &! temporary pressure scalars
!      fill       ! value to use for filling land points

!-----------------------------------------------------------------------
!
!  force T,S to be in valid range
!
!-----------------------------------------------------------------------

       TQ = min(DBLE(TEMPK),tmax)
       TQ = max(TQ,   tmin)

       SQ = min(DBLE(SALTK),smax)
       SQ = max(SQ,   smin)


       if (present(fillValue)) then
          fill = fillValue
       else
          fill = 0.d0
       endif

!-----------------------------------------------------------------------
!
!  now compute density or expansion coefficients
!
!-----------------------------------------------------------------------

       select case (state_itype)

!-----------------------------------------------------------------------
!
!  McDougall, Wright, Jackett, and Feistel EOS
!  test value : rho = 1.033213242 for
!  S = 35.0 PSU, theta = 20.0, pressz = 200.0
!
!-----------------------------------------------------------------------

       case (state_type_mwjf)


!	 TQ(230,150) = 20.
!	 SQ(230,150) = 35.
          p   = 10.d0*pressz(kk)
!	print*,'pressz(kk)=',pressz(kk)

!	  p = 200.
          SQ  = 1000.d0*SQ
	 SQ(230,150) = 35.
          SQR = sqrt(SQ)

      !***
      !*** first calculate numerator of MWJF 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 = mwjfnums0t0 + TQ * (mwjfnums0t1 + TQ * (mwjfnums0t2 + &
              mwjfnums0t3 * TQ)) + SQ * (mwjfnums1t0 +              &
              mwjfnums1t1 * TQ + mwjfnums2t0 * SQ)

      !***
      !*** now calculate denominator of MWJF density [P_2(S,T,p)]
      !***

          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 = mwjfdens0t0 + TQ * (mwjfdens0t1 + TQ * (mwjfdens0t2 +    &
           TQ * (mwjfdens0t3 + mwjfdens0t4 * TQ))) +                   &
           SQ * (mwjfdens1t0 + TQ * (mwjfdens1t1 + TQ*TQ*mwjfdens1t3)+ &
           SQR * (mwjfdensqt0 + TQ*TQ*mwjfdensqt2))

          DENOMK = 1.d0/WORK2

      if (present(RHOOUT)) then
         RHOOUT  = merge(WORK1*DENOMK,fill, MASK)
!         RHOOUT  = WORK1*DENOMK
!	 print*,RHOOUT(230,150),TQ(230,150),SQ(230,150)

      endif

      if (present(RHOFULL)) then
         RHOFULL = merge(WORK1*DENOMK,fill,MASK)
      endif

      if (present(DRHODT)) then
         WORK3 = &! dP_1/dT
                 mwjfnums0t1 + TQ * (2.d0*mwjfnums0t2 +    &
                 3.d0*mwjfnums0t3 * TQ) + mwjfnums1t1 * SQ

         WORK4 = &! dP_2/dT
                 mwjfdens0t1 + SQ * mwjfdens1t1 +               &
                 TQ * (2.d0*(mwjfdens0t2 + SQ*SQR*mwjfdensqt2) +  &
                 TQ * (3.d0*(mwjfdens0t3 + SQ * mwjfdens1t3) +    &
                 TQ *  4.d0*mwjfdens0t4))

         DRHODT = merge((WORK3 - WORK1*DENOMK*WORK4)*DENOMK,fill,MASK)
      endif

      if (present(DRHODS)) then
         WORK3 = &! dP_1/dS
                 mwjfnums1t0 + mwjfnums1t1 * TQ + 2.d0*mwjfnums2t0 * SQ

         WORK4 = mwjfdens1t0 +   &! dP_2/dS
                 TQ * (mwjfdens1t1 + TQ*TQ*mwjfdens1t3) +   &
                 1.5d0*SQR*(mwjfdensqt0 + TQ*TQ*mwjfdensqt2)

         DRHODS = merge((WORK3 - WORK1*DENOMK*WORK4)*DENOMK * 1000.d0,&
		  fill,MASK)
      endif

!-----------------------------------------------------------------------
!
!  Jackett and McDougall EOS
!
!-----------------------------------------------------------------------

   case (state_type_jmcd)

      p   = pressz(kk)
      p2  = p*p

      SQ  = 1000.d0*SQ
      SQR = sqrt(SQ)
      T2  = TQ*TQ


      !***
      !*** first calculate surface (p=0) values from UNESCO eqns.
      !***

      WORK1 = uns1t0 + uns1t1*TQ + & 
             (uns1t2 + uns1t3*TQ + uns1t4*T2)*T2
      WORK2 = SQR*(unsqt0 + unsqt1*TQ + unsqt2*T2)

      RHO_S = unt1*TQ + (unt2 + unt3*TQ + (unt4 + unt5*TQ)*T2)*T2 &
                      + (uns2t0*SQ + WORK1 + WORK2)*SQ


      !***
      !*** now calculate bulk modulus at pressure p from 
      !*** Jackett and McDougall formula
      !***

      WORK3 = bup0s1t0 + bup0s1t1*TQ +                    &
             (bup0s1t2 + bup0s1t3*TQ)*T2 +                &
              p *(bup1s1t0 + bup1s1t1*TQ + bup1s1t2*T2) + &
              p2*(bup2s1t0 + bup2s1t1*TQ + bup2s1t2*T2)
      WORK4 = SQR*(bup0sqt0 + bup0sqt1*TQ + bup0sqt2*T2 + &
                   bup1sqt0*p)

      BULK_MOD  = bup0s0t0 + bup0s0t1*TQ +                    &
                  (bup0s0t2 + bup0s0t3*TQ + bup0s0t4*T2)*T2 + &
                  p *(bup1s0t0 + bup1s0t1*TQ +                &
                     (bup1s0t2 + bup1s0t3*TQ)*T2) +           &
                  p2*(bup2s0t0 + bup2s0t1*TQ + bup2s0t2*T2) + &
                  SQ*(WORK3 + WORK4)

      DENOMK = 1.d0/(BULK_MOD - p)

      !***
      !*** now calculate required fields
      !***

      if (present(RHOOUT)) then
         RHOOUT  = merge(((unt0 + RHO_S)*BULK_MOD*DENOMK)*1.d-3, &
                         fill, MASK)
      endif

      if (present(RHOFULL)) then
         RHOFULL = merge(((unt0 + RHO_S)*BULK_MOD*DENOMK)*1.d-3, &
                         fill, MASK)
      endif

      if (present(DRHODT)) then
         DRDT0 =  unt1 + 2.d0*unt2*TQ +                      &
                  (3.d0*unt3 + 4.d0*unt4*TQ + 5.d0*unt5*T2)*T2 + &
                  (uns1t1 + 2.d0*uns1t2*TQ +                 &
                   (3.d0*uns1t3 + 4.d0*uns1t4*TQ)*T2 +         &
                   (unsqt1 + 2.d0*unsqt2*TQ)*SQR )*SQ

         DKDT  = bup0s0t1 + 2.d0*bup0s0t2*TQ +                       &
                 (3.d0*bup0s0t3 + 4.d0*bup0s0t4*TQ)*T2 +               &
                 p *(bup1s0t1 + 2.d0*bup1s0t2*TQ + 3.d0*bup1s0t3*T2) + &
                 p2*(bup2s0t1 + 2.d0*bup2s0t2*TQ) +                  &
                 SQ*(bup0s1t1 + 2.d0*bup0s1t2*TQ + 3.d0*bup0s1t3*T2 +  &
                     p  *(bup1s1t1 + 2.d0*bup1s1t2*TQ) +             &
                     p2 *(bup2s1t1 + 2.d0*bup2s1t2*TQ) +             &
                     SQR*(bup0sqt1 + 2.d0*bup0sqt2*TQ))

         DRHODT = merge ((DENOMK*(DRDT0*BULK_MOD -                    &
                           pressz(kk)*(unt0+RHO_S)*DKDT*DENOMK))*1.d-3, &
                         fill, MASK)
      endif

      if (present(DRHODS)) then
         DRDS0  = 2.d0*uns2t0*SQ + WORK1 + 1.5d0*WORK2
         DKDS = WORK3 + 1.5d0*WORK4

         DRHODS = merge (DENOMK*(DRDS0*BULK_MOD -                    &
                          pressz(kk)*(unt0+RHO_S)*DKDS*DENOMK),       &
                         fill, MASK)

      endif
 
!-----------------------------------------------------------------------
!
!  linear EOS
!
!-----------------------------------------------------------------------

   case (state_type_linear)

      if (present(RHOOUT )) RHOOUT  = bet*(SALTK - S_leos_ref) - &
                                      alf*(TEMPK - T_leos_ref)
      if (present(RHOFULL)) RHOFULL = rho_leos_ref + & 
                                      bet*(SALTK - S_leos_ref) - &
                                      alf*(TEMPK - T_leos_ref)
      if (present(DRHODT )) DRHODT  = -alf
      if (present(DRHODS )) DRHODS  =  bet

!-----------------------------------------------------------------------

   end select

!-----------------------------------------------------------------------
!EOC

 end subroutine state_ccsm

!***********************************************************************
!BOP
! !IROUTINE: init_state
! !INTERFACE:

 subroutine init_state(state_choice, zt)

! !DESCRIPTION:
!  Initializes eos choice and initializes coefficients for
!  equation of state.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   character (*), intent(in) :: &
      state_choice          ! choice of which EOS to use
                            ! valid choices are: mwjf, jmcd, linear

   real (r4), dimension(:), intent(in) :: &
      zt                    ! depth (cm) of each k level

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) ::  & 
      k                    ! vertical loop index

!-----------------------------------------------------------------------
!
!  set EOS type from input choice
!
!-----------------------------------------------------------------------

   select case (state_choice(1:4))
   case ('jmcd')
      state_itype = state_type_jmcd
   case ('mwjf')
      state_itype = state_type_mwjf
   case ('line')
      state_itype = state_type_linear
   case default
      stop 'unknown choice for EOS'
   end select

!-----------------------------------------------------------------------
!
!  calculate pressure as function of depth
!
!-----------------------------------------------------------------------

   allocate(pressz(size(zt)))

   do k=1,size(zt)
      pressz(k) = pressure(zt(k)*0.01d0)  ! convert depth to meters
   end do

!-----------------------------------------------------------------------
!
!  initialize various arrays, constants for each state type
!
!-----------------------------------------------------------------------

   select case (state_itype)

!-----------------------------------------------------------------------
!
!  McDougall, Wright, Jackett, and Feistel
!
!-----------------------------------------------------------------------

   case (state_type_mwjf)

      tmin =  -2.0d0  ! limited   on the low  end
      tmax = 999.0d0  ! unlimited on the high end
      smin =   0.0d0  ! limited   on the low  end
      smax = 0.999d0  ! unlimited on the high end

!-----------------------------------------------------------------------
!
!  Jackett and McDougall
!
!-----------------------------------------------------------------------

   case (state_type_jmcd)

      tmin = -2.0d0  ! valid pot. temp. range for level k
      tmax = 40.0d0 
      smin =  0.0d0  ! valid salinity   range for level k
      smax = 0.042d0 

!-----------------------------------------------------------------------
!
!  linear EOS - intialize valid range
!
!-----------------------------------------------------------------------

   case (state_type_linear)

      tmin = -2.0d0  ! valid pot. temp. range for level k
      tmax = 40.0d0 
      smin =  0.0d0  ! valid salinity   range for level k
      smax = 0.042d0 

   end select

!-----------------------------------------------------------------------
!EOC

 end subroutine init_state

!***********************************************************************
!BOP
! !IROUTINE: pressure
! !INTERFACE:

 function pressure(depth)

! !DESCRIPTION:
!  This function computes pressure in bars from depth in meters
!  using a mean density derived from depth-dependent global 
!  average temperatures and salinities from Levitus 1994, 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.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   real (r8), intent(in) :: depth    ! depth in meters

! !OUTPUT PARAMETERS:

   real (r8) :: pressure   ! pressure [bars]

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  convert depth in meters to pressure in bars
!
!-----------------------------------------------------------------------

   pressure = 0.059808d0*(exp(-0.025d0*depth) - 1.d0)     &
            + 0.100766d0*depth + 2.28405d-7*depth**2

!-----------------------------------------------------------------------
!EOC

 end function pressure

!***********************************************************************

 end module state_ccsm_mod

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

