forcing-perturbation/prep_restarts/comp_rhop/eosbn2.f90

MODULE eosbn2
   !!==============================================================================
   !!                       ***  MODULE  eosbn2  ***
   !! Ocean diagnostic variable : equation of state - in situ and potential density
   !!                                               - Brunt-Vaisala frequency 
   !!==============================================================================
   !! History :  OPA  ! 1989-03  (O. Marti)  Original code
   !!            6.0  ! 1994-07  (G. Madec, M. Imbard)  add bn2
   !!            6.0  ! 1994-08  (G. Madec)  Add Jackett & McDougall eos
   !!            7.0  ! 1996-01  (G. Madec)  statement function for e3
   !!            8.1  ! 1997-07  (G. Madec)  density instead of volumic mass
   !!             -   ! 1999-02  (G. Madec, N. Grima) semi-implicit pressure gradient
   !!            8.2  ! 2001-09  (M. Ben Jelloul)  bugfix on linear eos
   !!   NEMO     1.0  ! 2002-10  (G. Madec)  add eos_init
   !!             -   ! 2002-11  (G. Madec, A. Bozec)  partial step, eos_insitu_2d
   !!             -   ! 2003-08  (G. Madec)  F90, free form
   !!            3.0  ! 2006-08  (G. Madec)  add tfreez function
   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_alpbet used in ldfslp
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   eos            : generic interface of the equation of state
   !!   eos_insitu     : Compute the in situ density
   !!   eos_insitu_pot : Compute the insitu and surface referenced potential
   !!                    volumic mass
   !!   eos_insitu_2d  : Compute the in situ density for 2d fields
   !!   eos_bn2        : Compute the Brunt-Vaisala frequency
   !!   eos_alpbet     : calculates the in situ thermal and haline expansion coeff.
   !!   tfreez         : Compute the surface freezing temperature
   !!   eos_init       : set eos parameters (namelist)
   !!----------------------------------------------------------------------
!   USE dom_oce         ! ocean space and time domain
!   USE phycst          ! physical constants
!   USE zdfddm          ! vertical physics: double diffusion
!   USE in_out_manager  ! I/O manager
!   USE lib_mpp         ! MPP library
!   USE prtctl          ! Print control

   IMPLICIT NONE

   !                   !! * Interface 
   !                                          !!* Namelist (nameos) *
   INTEGER ::   nn_eos   = 0         !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
   DOUBLE PRECISION ::   rn_alpha = 2.0e-4    !: thermal expension coeff. (linear equation of state)
   DOUBLE PRECISION ::   rn_beta  = 7.7e-4    !: saline  expension coeff. (linear equation of state)

   DOUBLE PRECISION ::   ralpbet              !: alpha / beta ratio
   INTEGER, PARAMETER :: jpi=1442, jpj=1021, jpk=46, jpts=2
   DOUBLE PRECISION :: rau0 = 1035 
   
CONTAINS

   SUBROUTINE eos_insitu_pot( pt, ps, pd, pmask, prhop )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE eos_insitu_pot  ***
      !!           
      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the
      !!      potential volumic mass (Kg/m3) from potential temperature and
      !!      salinity fields using an equation of state defined through the 
      !!     namelist parameter nn_eos.
      !!
      !! ** Method  :
      !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state.
      !!         the in situ density is computed directly as a function of
      !!         potential temperature relative to the surface (the opa t
      !!         variable), salt and pressure (assuming no pressure variation
      !!         along geopotential surfaces, i.e. the pressure p in decibars
      !!         is approximated by the depth in meters.
      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
      !!              rhop(t,s)  = rho(t,s,0)
      !!         with pressure                      p        decibars
      !!              potential temperature         t        deg celsius
      !!              salinity                      s        psu
      !!              reference volumic mass        rau0     kg/m**3
      !!              in situ volumic mass          rho      kg/m**3
      !!              in situ density anomalie      prd      no units
      !!
      !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
      !!          t = 40 deg celcius, s=40 psu
      !!
      !!      nn_eos = 1 : linear equation of state function of temperature only
      !!              prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - rn_alpha * t
      !!              rhop(t,s)  = rho(t,s)
      !!
      !!      nn_eos = 2 : linear equation of state function of temperature and
      !!               salinity
      !!              prd(t,s) = ( rho(t,s) - rau0 ) / rau0 
      !!                       = rn_beta * s - rn_alpha * tn - 1.
      !!              rhop(t,s)  = rho(t,s)
      !!      Note that no boundary condition problem occurs in this routine
      !!      as (tn,sn) or (ta,sa) are defined over the whole domain.
      !!
      !! ** Action  : - prd  , the in situ density (no units)
      !!              - prhop, the potential volumic mass (Kg/m3)
      !!
      !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
      !!                Brown and Campana, Mon. Weather Rev., 1978
      !!----------------------------------------------------------------------
      !!
      DOUBLE PRECISION, DIMENSION(jpi,jpj,jpk     ), INTENT(   in) ::   pt    !  potential temperature  [Celcius]
      DOUBLE PRECISION, DIMENSION(jpi,jpj,jpk     ), INTENT(   in) ::   ps    !  salinity               [psu]
      DOUBLE PRECISION, DIMENSION(jpi,jpj,jpk     ), INTENT(   in) ::   pmask
      DOUBLE PRECISION, DIMENSION(jpk             ), INTENT(   in) ::   pd   
      DOUBLE PRECISION, DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
      !
      INTEGER  ::   ji, jj, jk   ! dummy loop indices
      DOUBLE PRECISION, DIMENSION(jpi,jpj,jpk     ) ::   zrd    ! in situ density            [-]
      DOUBLE PRECISION, DIMENSION(jpi,jpj,jpk     ) :: zws
      DOUBLE PRECISION ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw   ! local scalars
      DOUBLE PRECISION ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zrau0r       !   -      -
      !!----------------------------------------------------------------------

      SELECT CASE ( nn_eos )
      !
      CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==!
         zrau0r = 1.e0 / rau0
!CDIR NOVERRCHK
         zws(:,:,:) = SQRT( ABS( ps(:,:,:) ) )
         !  
         DO jk = 1, (jpk-1)
            DO jj = 1, jpj
               DO ji = 1, jpi
                  zt = pt    (ji,jj,jk)
                  zs = ps    (ji,jj,jk)
                  zh = pd    (jk)        ! depth
                  zsr= zws   (ji,jj,jk)        ! square root salinity
                  !
                  ! compute volumic mass pure water at atm pressure
                  zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4 )*zt   &
                     &                          -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594
                  ! seawater volumic mass atm pressure
                  zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt   &
                     &                                         -4.0899e-3 ) *zt+0.824493
                  zr3= ( -1.6546e-6*zt+1.0227e-4 )    *zt-5.72466e-3
                  zr4= 4.8314e-4
                  !
                  ! potential volumic mass (reference to the surface)
                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
                  !
                  ! save potential volumic mass
                  prhop(ji,jj,jk) = zrhop * pmask(ji,jj,jk)
                  !
                  ! add the compression terms
                  ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6
                  zbw= (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4
                  zb = zbw + ze * zs
                  !
                  zd = -2.042967e-2
                  zc =   (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896
                  zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788
                  za = ( zd*zsr + zc ) *zs + zaw
                  !
                  zb1=   (  -0.1909078  *zt+7.390729    ) *zt-55.87545
                  za1= ( (   2.326469e-3*zt+1.553190    ) *zt-65.00517 ) *zt + 1044.077
                  zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6
                  zk0= ( zb1*zsr + za1 )*zs + zkw
                  !
                  ! masked in situ density anomaly
                  zrd(ji,jj,jk) = (  zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
                     &             - rau0  ) * zrau0r * pmask(ji,jj,jk)
               END DO
            END DO
         END DO
         !
      CASE( 1 )                !==  Linear formulation = F( temperature )  ==!
         DO jk = 1, (jpk-1)
            zrd  (:,:,jk) = ( 0.0285 - rn_alpha * pt (:,:,jk) )        * pmask(:,:,jk)
            prhop(:,:,jk) = ( 1.e0   +            zrd (:,:,jk)       ) * rau0 * pmask(:,:,jk)
         END DO
         !
      CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
         DO jk = 1, (jpk-1)
            zrd  (:,:,jk) = ( rn_beta  * ps(:,:,jk) - rn_alpha * pt(:,:,jk) )        * pmask(:,:,jk)
            prhop(:,:,jk) = ( 1.e0  + zrd (:,:,jk) )                                       * rau0 * pmask(:,:,jk)
         END DO
         !
      END SELECT
      !
      !
   END SUBROUTINE eos_insitu_pot


   !!======================================================================
END MODULE eosbn2