#include"cppdefs.h" !----------------------------------------------------------------------- !BOP ! ! !MODULE: eqstate --- the equation of state \label{sec:eqstate} ! ! !INTERFACE: MODULE eqstate ! ! !DESCRIPTION: ! Computes the density, $\mean{\rho}$, and buoyancy from the ! salinity, $S$, the temperature, $\Theta$, and the thermodynamic ! pressure, $P$, according to an \emph{equation of state}, ! \begin{equation} ! \label{DefEOS} ! \mean{\rho} = \hat{\rho} (S,\Theta,P) ! \point ! \end{equation} ! ! The following remark on the thermodynamic interpretation of ! density, temperature, and pressure is useful here. If $\Theta$ is ! identified with the in-situ temperature, and $P$ with the in-situ ! pressure, then $\mean{\rho}$ will be the in-situ density. On the ! other hand, if $P$ is identified with the surface pressure, and ! $\Theta$ with the potential temperature, the same equation of ! state, \eq{DefEOS}, will yield $\mean{\rho}$ as the potential ! density. Note that the quantity {\tt sigma\_t} found in the GOTM output ! is simply computed from $\mean{\rho}$ - 1000 kg m$^{-3}$, and may ! therefore adopt different meanings. ! ! ! At present, two different models for the equation of state ("modes"), ! and four different "methods" how to evalute the equation of state ! are implemented. ! ! Modes: ! \begin{enumerate} ! \item The UNESCO equation of state according to \cite{FofonoffMillard83} ! \item The \cite{JACKETTea05} equation of state ! \end{enumerate} ! Methods: ! \begin{enumerate} ! \item the full equation of state --- including pressure effects ! \item the full equation of state --- without pressure effects ! \item the linearised equation of state ! \item a general linear form of the equation of state ! \end{enumerate} !USES: IMPLICIT NONE ! default: all is private. private ! ! !PUBLIC MEMBER FUNCTIONS: public init_eqstate,eqstate1,eos_alpha,eos_beta,unesco,rho_feistel ! ! !REVISION HISTORY: ! Original author(s): Hans Burchard & Karsten Bolding ! ! $Log: eqstate.F90,v $ ! Revision 20.0 2013/12/14 00:14:01 fms ! Merged revision 1.1.2.1 onto trunk ! ! Revision 1.1.2.1 2012/05/15 16:01:18 smg ! initial cvs ci to mom5 ! AUTHOR:Griffies ! REVIEWERS: ! TEST STATUS: ! CHANGES PUBLIC INTERFACES? ! CHANGES ANSWERS? ! ! Revision 1.1.2.1.390.1 2012/04/23 20:33:41 smg ! updated to gotm-2012.03.09 CVS tagged code. ! AUTHOR:Martin Schmidt ! REVIEWERS:Griffies ! TEST STATUS: ! CHANGES PUBLIC INTERFACES? ! CHANGES ANSWERS? ! ! Revision 1.9 2010-09-17 12:53:53 jorn ! extensive code clean-up to ensure proper initialization and clean-up of all variables ! ! Revision 1.8 2008-08-01 07:33:14 lars ! modified documentation ! ! Revision 1.7 2007-01-06 11:49:13 kbk ! namelist file extension changed .inp --> .nml ! ! Revision 1.6 2005/06/27 13:44:07 kbk ! modified + removed traling blanks ! ! Revision 1.5 2003/03/28 09:20:36 kbk ! added new copyright to files ! ! Revision 1.4 2003/03/28 08:06:33 kbk ! removed tabs ! ! Revision 1.3 2003/03/10 08:54:16 gotm ! Improved documentation and cleaned up code ! ! Revision 1.2 2001/11/27 19:44:32 gotm ! Fixed an initialisation bug ! ! Revision 1.1.1.1 2001/02/12 15:55:58 gotm ! initial import into CVS ! !EOP ! ! private data memebers integer :: eq_state_method, eq_state_mode REALTYPE :: T0,S0,p0,dtr0,dsr0 logical :: init_linearised ! !----------------------------------------------------------------------- contains !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: Read the namelist {\tt eqstate} ! ! !INTERFACE: subroutine init_eqstate(namlst) ! ! !DESCRIPTION: ! Here, the namelist {\tt eqstate} in the namelist file {\tt gotmrun.nml} ! is read. ! ! !USES: IMPLICIT NONE ! ! !INPUT PARAMETERS: integer, optional, intent(in) :: namlst ! ! !REVISION HISTORY: ! Original author(s): Hans Burchard & Karsten Bolding ! !EOP ! ! !LOCAL VARIABLES: namelist /eqstate/ eq_state_mode,eq_state_method,T0,S0,p0,dtr0,dsr0 ! !----------------------------------------------------------------------- !BOC LEVEL1 'init_eqstate' init_linearised = .true. if(present(namlst)) then read(namlst,nml=eqstate,err=80) end if return 80 FATAL 'I could not read "eqstate" namelist' stop 'init_eqstate' end subroutine init_eqstate !EOC !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: Select an equation of state ! ! !INTERFACE: REALTYPE function eqstate1(S,T,p,g,rho_0) ! ! !DESCRIPTION: ! Calculates the in-situ buoyancy according to the selected method. ! {\tt S} is salinity $S$ in psu, {\tt T} is ! potential temperature $\theta$ in $^{\circ}$C (ITS-90), {\tt p} is ! gauge pressure (absolute pressure - 10.1325 bar), {\tt g} is the ! gravitational acceleration in m\,s$^{-2}$ and {\tt rho\_0} the reference ! density in kg\,m$^{-3}$. {\tt eqstate1} is the in-situ-density ! in kg\,m$^{-3}$. ! For {\tt eq\_state\_method}=1, the UNESCO equation of state is used, ! for {\tt eq\_state\_method}=2, the \cite{JACKETTea05} equation ! of state is used. Here, some care is needed, since the UNESCO equation ! used bar for pressure and the \cite{JACKETTea05} uses dbar for pressure. ! For values of ! {\tt eq\_state\_method} ranging from 1 to 4, one of the following methods ! will be used. ! ! \begin{enumerate} ! \item the full equation of state for sea water ! including pressure dependence. ! \item the equation of state for sea water ! with the pressure evaluated at the sea surface as ! reference level. This is the choice ! for computations based on potential temperature and density. ! \item a linearised equation of state. ! The parameters {\tt T0}, ! {\tt S0} and {\tt p0} have to be specified in the namelist. ! \item a linear equation of state with prescribed {\tt rho0}, {\tt T0}, ! {\tt S0}, {\tt dtr0}, {\tt dsr0} according to ! \begin{equation} ! \label{eosLinear} ! \rho = \rho_0 + \text{\tt dtr0} (T - T_0) ! + \text{\tt dsr0} (S - S_0) ! \point ! \end{equation} ! \end{enumerate} ! ! !USES: IMPLICIT NONE ! ! !INPUT PARAMETERS: REALTYPE,intent(in) :: S,T,p REALTYPE,optional,intent(in) :: g,rho_0 ! ! !REVISION HISTORY: ! Original author(s): Hans Burchard & Karsten Bolding ! !EOP ! ! !LOCAL VARIABLES: REALTYPE :: x REALTYPE, save :: rh0,dtr,dsr REALTYPE :: dTT,dSS logical :: press ! !----------------------------------------------------------------------- !BOC select case (eq_state_mode) case(1) select case (eq_state_method) case (1) press=.true. x=unesco(S,T,p,press) case (2) press=.false. x=unesco(S,T,p,press) case (3) if (init_linearised) then press=.true. ! This allows for choosing potentials other than p=0 dTT=0.001 dSS=0.001 rh0= unesco(S0,T0,p0,press) dtr=(unesco(S0,T0+0.5*dTT,p0,press) & -unesco(S0,T0-0.5*dTT,p0,press))/dTT dsr=(unesco(S0+0.5*dSS,T0,p0,press) & -unesco(S0-0.5*dSS,T0,p0,press))/dSS init_linearised=.false. end if x=rh0+dtr*(T-T0)+dsr*(S-S0) case (4) x=rho_0+dtr0*(T-T0)+dsr0*(S-S0) case default end select case(2) select case (eq_state_method) case (1) press=.true. x=rho_feistel(S,T,p*10.,press) case (2) press=.false. x=rho_feistel(S,T,p*10.,press) case (3) if (init_linearised) then press=.true. ! This allows for choosing potentials other than p=0 dTT=0.001 dSS=0.001 rh0= rho_feistel(S0,T0,p0*10.,press) dtr=(rho_feistel(S0,T0+0.5*dTT,p0*10.,press) & -rho_feistel(S0,T0-0.5*dTT,p0*10.,press))/dTT dsr=(rho_feistel(S0+0.5*dSS,T0,p0*10.,press) & -rho_feistel(S0-0.5*dSS,T0,p0*10.,press))/dSS init_linearised=.false. end if x=rh0+dtr*(T-T0)+dsr*(S-S0) case (4) x=rho_0+dtr0*(T-T0)+dsr0*(S-S0) case default end select case default end select eqstate1=-g*(x-rho_0)/rho_0 return end function eqstate1 !EOC !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: Compute thermal expansion coefficient ! ! !INTERFACE: REALTYPE function eos_alpha(S,T,p,g,rho_0) ! ! !DESCRIPTION: ! Computes the thermal expansion coefficient defined by ! \begin{equation} ! \label{eosAlpha} ! \alpha = ! - \dfrac{1}{\rho_0} ! \left( \partder{\rho_{is}}{T} \right)_S ! = ! \dfrac{1}{g} ! \left( \partder{B_{is}}{T} \right)_S ! \comma ! \end{equation} ! where $B_{is}$ denotes the in-situ buoyancy. The computation is carried ! out by a finite difference approximation of \eq{eosAlpha}, ! requiring two evaluations of the equation of state. ! Note, that comparing \eq{eosAlpha} with \eq{eosLinear} it follows that ! $\alpha = - \text{\tt dtr0}/\rho_0$. ! ! !USES: IMPLICIT NONE ! ! !INPUT PARAMETERS: REALTYPE,intent(in) :: S,T,p REALTYPE,optional,intent(in) :: g,rho_0 ! ! !REVISION HISTORY: ! Original author(s): Lars Umlauf ! !EOP ! ! !LOCAL VARIABLES: ! REALTYPE,parameter :: delta = 0.01 REALTYPE :: buoy_a,buoy_b !----------------------------------------------------------------------- !BOC buoy_a = eqstate1(S,T+0.5*delta,p,g,rho_0) buoy_b = eqstate1(S,T-0.5*delta,p,g,rho_0) eos_alpha = (buoy_a - buoy_b) / (g*delta) end function eos_alpha !EOC !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: Compute saline contraction coefficient ! ! !INTERFACE: REALTYPE function eos_beta(S,T,p,g,rho_0) ! ! !DESCRIPTION: ! Computes the saline contractioncoefficient defined by ! \begin{equation} ! \label{eosBeta} ! \beta = ! \dfrac{1}{\rho_0} ! \left( \partder{\rho_{is}}{S} \right)_T ! = ! - \dfrac{1}{g} ! \left( \partder{B_{is}}{S} \right)_T ! \comma ! \end{equation} ! where $B_{is}$ denotes the in-situ buoyancy. The computation is carried ! out by a finite difference approximation of \eq{eosBeta}, ! requiring two evaluations of the equation of state. ! Note, that comparing \eq{eosBeta} with \eq{eosLinear} it follows that ! $\beta = \text{\tt dsr0}/\rho_0$. ! ! ! !USES: IMPLICIT NONE ! ! !INPUT PARAMETERS: REALTYPE,intent(in) :: S,T,p REALTYPE,optional,intent(in) :: g,rho_0 ! ! !REVISION HISTORY: ! Original author(s): Lars Umlauf ! !EOP ! ! !LOCAL VARIABLES: ! REALTYPE,parameter :: delta = 0.01 REALTYPE :: buoy_a,buoy_b !----------------------------------------------------------------------- !BOC buoy_a = eqstate1(S+0.5*delta,T,p,g,rho_0) buoy_b = eqstate1(S-0.5*delta,T,p,g,rho_0) eos_beta = -(buoy_a - buoy_b) / (g*delta) end function eos_beta !EOC !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: The UNESCO equation of state ! ! !INTERFACE: REALTYPE function unesco(S,T,p,UNPress) ! ! !DESCRIPTION: ! Computes the in-situ density in \eq{DefEOS} according to the ! UNESCO equation of state for sea water (see \cite{FofonoffMillard83}). ! The pressure ! dependence can be switched on ({\tt UNPress=.true.}) or off ! ({\tt UNPress=.false.}). {\tt S} is salinity $S$ in psu, {\tt T} is ! potential temperature $\theta$ in $^{\circ}$C (ITS-90), {\tt p} is ! gauge pressure (absolute pressure - 10.1325 bar) and ! {\tt unesco} is the in-situ density in kg\,m$^{-3}$. ! The check value is {\tt unesco(35,25,1000) = 1062.53817} . ! ! !USES: IMPLICIT NONE ! ! !INPUT PARAMETERS: REALTYPE, intent(in) :: S,T,p LOGICAL, intent(in) :: UNPress ! ! !REVISION HISTORY: ! Original author(s): Hans Burchard & Karsten Bolding ! !EOP ! ! !LOCAL VARIABLES: REALTYPE :: x,K REALTYPE :: T2,T3,T4,T5,S15,S2,S3,p2 ! !----------------------------------------------------------------------- !BOC T2 = T*T T3 = T*T2 T4 = T2*T2 T5 = T*T4 S15= S**1.5 S2 = S*S S3 = S*S2 x=999.842594+6.793952e-02*T-9.09529e-03*T2+1.001685e-04*T3 x=x-1.120083e-06*T4+6.536332e-09*T5 x=x+S*(0.824493-4.0899e-03*T+7.6438e-05*T2-8.2467e-07*T3) x=x+S*5.3875e-09*T4 x=x+sqrt(S3)*(-5.72466e-03+1.0227e-04*T-1.6546e-06*T2) x=x+4.8314e-04*S2 if ((UNPress).and.(p.gt.0)) then p2=p*p K= 19652.21 & +148.4206 *T -2.327105 *T2 & + 1.360477E-2*T3 -5.155288E-5 *T4 & + 3.239908 *p +1.43713E-3 *T *p & + 1.16092E-4 *T2*p -5.77905E-7 *T3*p & + 8.50935E-5 *p2 -6.12293E-6 *T *p2 & + 5.2787E-8 *T2*p2 & + 54.6746 *S -0.603459 *T *S & + 1.09987E-2 *T2 *S -6.1670E-5 *T3 *S & + 7.944E-2 *S15+1.6483E-2 *T *S15 & - 5.3009E-4 *T2 *S15+2.2838E-3 *p *S & - 1.0981E-5 *T *p *S -1.6078E-6 *T2*p *S & + 1.91075E-4 *p *S15-9.9348E-7 *p2*S & + 2.0816E-8 *T *p2*S +9.1697E-10 *T2*p2*S x=x/(1.-p/K) end if unesco=x return end function unesco !EOC !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: The \cite{JACKETTea05} equation of state ! ! !INTERFACE: REALTYPE function rho_feistel(s,th,p,UNPress) ! ! !DESCRIPTION: ! Computes the in-situ density in \eq{DefEOS} according to the ! \cite{JACKETTea05} equation of state for sea water, which is based ! on the Gibbs potential developed by \cite{FEISTEL03}. The pressure ! dependence can be switched on ({\tt UNPress=.true.}) or off ! ({\tt UNPress=.false.}). {\tt s} is salinity $S$ in psu, {\tt th} is ! potential temperature $\theta$ in $^{\circ}$C (ITS-90), {\tt p} is ! gauge pressure (absolute pressure - 10.1325 dbar) and ! {\tt rho\_feistel} is the in-situ density in kg\,m$^{-3}$. ! The check value is {\tt rho\_feistel(20,20,1000) = 1017.728868019642} . ! ! !USES: IMPLICIT NONE ! ! !INPUT PARAMETERS: REALTYPE, intent(in) :: s,th,p LOGICAL, intent(in) :: UNPress ! ! !REVISION HISTORY: ! Original author(s): Hans Burchard & Karsten Bolding ! !EOP ! ! !LOCAL VARIABLES: REALTYPE :: th2,sqrts,pth,anum,aden ! !----------------------------------------------------------------------- !BOC th2 = th*th sqrts = sqrt(s) anum = 9.9984085444849347d+02 + & th*( 7.3471625860981584d+00 + & th*(-5.3211231792841769d-02 + & th* 3.6492439109814549d-04)) + & s*( 2.5880571023991390d+00 - & th* 6.7168282786692355d-03 + & s* 1.9203202055760151d-03) aden = 1.0000000000000000d+00 + & th*( 7.2815210113327091d-03 + & th*(-4.4787265461983921d-05 + & th*( 3.3851002965802430d-07 + & th* 1.3651202389758572d-10))) + & s*( 1.7632126669040377d-03 - & th*( 8.8066583251206474d-06 + & th2* 1.8832689434804897d-10) + & sqrts*( 5.7463776745432097d-06 + & th2* 1.4716275472242334d-09)) if((UNPress).and.(p.gt.0.d0)) then pth = p*th anum = anum + p*( 1.1798263740430364d-02 + & th2* 9.8920219266399117d-08 + & s* 4.6996642771754730d-06 - & p*( 2.5862187075154352d-08 + & th2* 3.2921414007960662d-12)) aden = aden + p*( 6.7103246285651894d-06 - & pth*(th2* 2.4461698007024582d-17 + & p* 9.1534417604289062d-18)) end if rho_feistel = anum/aden ! Note: this function should always be run in double precision ! (since rho is returned rather than sigma = rho-1.0d3) return end function rho_feistel !EOC end module eqstate !----------------------------------------------------------------------- ! Copyright by the GOTM-team under the GNU Public License - www.gnu.org !-----------------------------------------------------------------------