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