Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 7277)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 7278)
@@ -75,6 +75,15 @@
 
    !                               !!** Namelist nameos **
-   INTEGER , PUBLIC ::   nn_eos     ! = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
-   LOGICAL , PUBLIC ::   ln_useCT   ! determine if eos_pt_from_ct is used to compute sst_m
+   LOGICAL , PUBLIC ::   ln_TEOS10   ! determine if eos_pt_from_ct is used to compute sst_m
+   LOGICAL , PUBLIC ::   ln_EOS80   ! determine if eos_pt_from_ct is used to compute sst_m
+   LOGICAL , PUBLIC ::   ln_SEOS   ! determine if eos_pt_from_ct is used to compute sst_m
+
+   ! Parameters
+   LOGICAL , PUBLIC    ::   l_useCT         ! =T in ln_TEOS10=T (i.e. use eos_pt_from_ct to compute sst_m), =F otherwise
+   INTEGER , PUBLIC    ::   neos            ! Identifier for equation of state used
+
+   INTEGER , PARAMETER ::   np_teos10 = -1  ! parameter for using TEOS10
+   INTEGER , PARAMETER ::   np_eos80  =  0  ! parameter for using EOS80
+   INTEGER , PARAMETER ::   np_seos   = 1   ! parameter for using Simplified Equation of state
 
    !                               !!!  simplified eos coefficients (default value: Vallis 2006)
@@ -184,5 +193,5 @@
       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
       !!       potential temperature and salinity using an equation of state
-      !!       defined through the namelist parameter nn_eos.
+      !!       selected in the nameos namelist
       !!
       !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0
@@ -194,11 +203,11 @@
       !!                rau0   reference density            kg/m^3
       !!
-      !!     nn_eos = -1 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
+      !!     ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
       !!         Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg
       !!
-      !!     nn_eos =  0 : polynomial EOS-80 equation of state is used for rho(t,s,z).
+      !!     ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z).
       !!         Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu
       !!
-      !!     nn_eos =  1 : simplified equation of state
+      !!     ln_seos : simplified equation of state
       !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rau0
       !!              linear case function of T only: rn_alpha<>0, other coefficients = 0
@@ -224,7 +233,7 @@
       IF( nn_timing == 1 )   CALL timing_start('eos-insitu')
       !
-      SELECT CASE( nn_eos )
-      !
-      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
+      SELECT CASE( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
          DO jk = 1, jpkm1
@@ -266,5 +275,5 @@
          END DO
          !
-      CASE( 1 )                !==  simplified EOS  ==!
+      CASE( np_seos )                !==  simplified EOS  ==!
          !
          DO jk = 1, jpkm1
@@ -300,6 +309,6 @@
       !! ** 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.
+      !!      salinity fields using an equation of state selected in the
+      !!     namelist.
       !!
       !! ** Action  : - prd  , the in situ density (no units)
@@ -322,7 +331,7 @@
       IF( nn_timing == 1 )   CALL timing_start('eos-pot')
       !
-      SELECT CASE ( nn_eos )
-      !
-      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
+      SELECT CASE ( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
          ! Stochastic equation of state
@@ -430,5 +439,5 @@
          ENDIF
          
-      CASE( 1 )                !==  simplified EOS  ==!
+      CASE( np_seos )                !==  simplified EOS  ==!
          !
          DO jk = 1, jpkm1
@@ -467,5 +476,5 @@
       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
       !!      potential temperature and salinity using an equation of state
-      !!      defined through the namelist parameter nn_eos. * 2D field case
+      !!      selected in the nameos namelist. * 2D field case
       !!
       !! ** Action  : - prd , the in situ density (no units) (unmasked)
@@ -486,7 +495,7 @@
       prd(:,:) = 0._wp
       !
-      SELECT CASE( nn_eos )
-      !
-      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
+      SELECT CASE( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
          DO jj = 1, jpjm1
@@ -527,5 +536,5 @@
          CALL lbc_lnk( prd, 'T', 1. )                    ! Lateral boundary conditions
          !
-      CASE( 1 )                !==  simplified EOS  ==!
+      CASE( np_seos )                !==  simplified EOS  ==!
          !
          DO jj = 1, jpjm1
@@ -576,7 +585,7 @@
       IF( nn_timing == 1 )   CALL timing_start('rab_3d')
       !
-      SELECT CASE ( nn_eos )
-      !
-      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
+      SELECT CASE ( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
          DO jk = 1, jpkm1
@@ -635,5 +644,5 @@
          END DO
          !
-      CASE( 1 )                  !==  simplified EOS  ==!
+      CASE( np_seos )                  !==  simplified EOS  ==!
          !
          DO jk = 1, jpkm1
@@ -657,5 +666,5 @@
       CASE DEFAULT
          IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
+         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos
          nstop = nstop + 1
          !
@@ -668,4 +677,5 @@
       !
    END SUBROUTINE rab_3d
+
 
    SUBROUTINE rab_2d( pts, pdep, pab )
@@ -690,7 +700,7 @@
       pab(:,:,:) = 0._wp
       !
-      SELECT CASE ( nn_eos )
-      !
-      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
+      SELECT CASE ( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
          DO jj = 1, jpjm1
@@ -750,5 +760,5 @@
          CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. )                    
          !
-      CASE( 1 )                  !==  simplified EOS  ==!
+      CASE( np_seos )                  !==  simplified EOS  ==!
          !
          DO jj = 1, jpjm1
@@ -773,5 +783,5 @@
       CASE DEFAULT
          IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
+         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos
          nstop = nstop + 1
          !
@@ -806,7 +816,7 @@
       pab(:) = 0._wp
       !
-      SELECT CASE ( nn_eos )
-      !
-      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
+      SELECT CASE ( neos )
+      !
+      CASE( np_teos10, np_eos80 )      !==  polynomial TEOS-10 / EOS-80 ==!
          !
          !
@@ -859,9 +869,9 @@
          !
          !
-      CASE( 1 )                  !==  simplified EOS  ==!
+      CASE( np_seos )                  !==  simplified EOS  ==!
          !
          zt    = pts(jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
          zs    = pts(jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
-         zh    = pdep                    ! depth at the partial step level
+         zh    = pdep                   ! depth at the partial step level
          !
          zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
@@ -873,5 +883,5 @@
       CASE DEFAULT
          IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
+         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos
          nstop = nstop + 1
          !
@@ -1005,15 +1015,16 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(out  )           ::   ptf    ! freezing temperature [Celsius]
       !
-      INTEGER  ::   ji, jj   ! dummy loop indices
-      REAL(wp) ::   zt, zs   ! local scalars
-      !!----------------------------------------------------------------------
-      !
-      SELECT CASE ( nn_eos )
-      !
-      CASE ( -1, 1 )                !==  CT,SA (TEOS-10 formulation) ==!
-         !
+      INTEGER  ::   ji, jj          ! dummy loop indices
+      REAL(wp) ::   zt, zs, z1_S0   ! local scalars
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE ( neos )
+      !
+      CASE ( np_teos10, np_seos )      !==  CT,SA (TEOS-10 and S-EOS formulations) ==!
+         !
+         z1_S0 = 1._wp / 35.16504_wp
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zs= SQRT( ABS( psal(ji,jj) ) * r1_S0 )           ! square root salinity
+               zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity
                ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
                   &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
@@ -1024,5 +1035,5 @@
          IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
          !
-      CASE ( 0 )                     !==  PT,SP (UNESCO formulation)  ==!
+      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==!
          !
          ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
@@ -1033,5 +1044,5 @@
       CASE DEFAULT
          IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
+         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos
          nstop = nstop + 1
          !
@@ -1039,4 +1050,5 @@
       !
   END SUBROUTINE eos_fzp_2d
+
 
   SUBROUTINE eos_fzp_0d( psal, ptf, pdep )
@@ -1059,9 +1071,9 @@
       !!----------------------------------------------------------------------
       !
-      SELECT CASE ( nn_eos )
-      !
-      CASE ( -1, 1 )                !==  CT,SA (TEOS-10 formulation) ==!
-         !
-         zs  = SQRT( ABS( psal ) * r1_S0 )           ! square root salinity
+      SELECT CASE ( neos )
+      !
+      CASE ( np_teos10, np_seos )      !==  CT,SA (TEOS-10 and S-EOS formulations) ==!
+         !
+         zs  = SQRT( ABS( psal ) / 35.16504_wp )           ! square root salinity
          ptf = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
                   &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
@@ -1070,5 +1082,5 @@
          IF( PRESENT( pdep ) )   ptf = ptf - 7.53e-4 * pdep
          !
-      CASE ( 0 )                     !==  PT,SP (UNESCO formulation)  ==!
+      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==!
          !
          ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal )   &
@@ -1079,5 +1091,5 @@
       CASE DEFAULT
          IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
+         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos
          nstop = nstop + 1
          !
@@ -1118,7 +1130,7 @@
       IF( nn_timing == 1 )   CALL timing_start('eos_pen')
       !
-      SELECT CASE ( nn_eos )
-      !
-      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
+      SELECT CASE ( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
          DO jk = 1, jpkm1
@@ -1183,5 +1195,5 @@
          END DO
          !
-      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==!
+      CASE( np_seos )                !==  Vallis (2006) simplified EOS  ==!
          !
          DO jk = 1, jpkm1
@@ -1205,5 +1217,5 @@
       CASE DEFAULT
          IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
+         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos
          nstop = nstop + 1
          !
@@ -1224,6 +1236,7 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ios   ! local integer
-      !!
-      NAMELIST/nameos/ nn_eos, ln_useCT, rn_a0, rn_b0, rn_lambda1, rn_mu1,   &
+      INTEGER  ::   ioptio   ! local integer
+      !!
+      NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS, rn_a0, rn_b0, rn_lambda1, rn_mu1,   &
          &                                             rn_lambda2, rn_mu2, rn_nu
       !!----------------------------------------------------------------------
@@ -1245,19 +1258,24 @@
          WRITE(numout,*) 'eos_init : equation of state'
          WRITE(numout,*) '~~~~~~~~'
-         WRITE(numout,*) '          Namelist nameos : set eos parameters'
-         WRITE(numout,*) '             flag for eq. of state and N^2  nn_eos   = ', nn_eos
-         IF( ln_useCT )   THEN
-            WRITE(numout,*) '             model uses Conservative Temperature'
-            WRITE(numout,*) '             Important: model must be initialized with CT and SA fields'
-         ELSE
-            WRITE(numout,*) '             model does not use Conservative Temperature'
-         ENDIF
+         WRITE(numout,*) '   Namelist nameos : Chosen the Equation Of Seawater (EOS)'
+         WRITE(numout,*) '      TEOS-10 : rho=F(Conservative Temperature, Absolute  Salinity, depth)   ln_TEOS10 = ', ln_TEOS10
+         WRITE(numout,*) '      EOS-80  : rho=F(Potential    Temperature, Practical Salinity, depth)   ln_EOS80  = ', ln_EOS80
+         WRITE(numout,*) '      S-EOS   : rho=F(Conservative Temperature, Absolute  Salinity, depth)   ln_SEOS   = ', ln_SEOS
       ENDIF
-      !
-      SELECT CASE( nn_eos )         ! check option
-      !
-      CASE( -1 )                       !==  polynomial TEOS-10  ==!
+
+      ! Check options for equation of state & set neos based on logical flags
+      ioptio = 0
+      IF( ln_TEOS10 ) THEN   ;   ioptio = ioptio+1   ;   neos = np_teos10   ;   ENDIF
+      IF( ln_EOS80  ) THEN   ;   ioptio = ioptio+1   ;   neos = np_eos80    ;   ENDIF
+      IF( ln_SEOS   ) THEN   ;   ioptio = ioptio+1   ;   neos = np_seos     ;   ENDIF
+      IF( ioptio /= 1 )   CALL ctl_stop("Exactly one equation of state option must be selected")
+      !
+      SELECT CASE( neos )         ! check option
+      !
+      CASE( np_teos10 )                       !==  polynomial TEOS-10  ==!
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '          use of TEOS-10 equation of state (cons. temp. and abs. salinity)'
+         !
+         l_useCT = .TRUE.                          ! model temperature is Conservative temperature 
          !
          rdeltaS = 32._wp
@@ -1446,9 +1464,10 @@
          BPE002 = 1.7269476440e-04_wp
          !
-      CASE( 0 )                        !==  polynomial EOS-80 formulation  ==!
+      CASE( np_eos80 )                        !==  polynomial EOS-80 formulation  ==!
          !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '          use of EOS-80 equation of state (pot. temp. and pract. salinity)'
          !
+         l_useCT = .FALSE.                         ! model temperature is Potential temperature
          rdeltaS = 20._wp
          r1_S0  = 1._wp/40._wp
@@ -1636,5 +1655,5 @@
          BPE002 = 5.3661089288e-04_wp
          !
-      CASE( 1 )                        !==  Simplified EOS     ==!
+      CASE( np_seos )                        !==  Simplified EOS     ==!
          IF(lwp) THEN
             WRITE(numout,*)
@@ -1651,7 +1670,8 @@
             WRITE(numout,*) '               Caution: rn_beta0=0 incompatible with ddm parameterization '
          ENDIF
-         !
-      CASE DEFAULT                     !==  ERROR in nn_eos  ==!
-         WRITE(ctmp1,*) '          bad flag value for nn_eos = ', nn_eos
+         l_useCT = .TRUE.          ! Use conservative temperature
+         !
+      CASE DEFAULT                     !==  ERROR in neos  ==!
+         WRITE(ctmp1,*) '          bad flag value for neos = ', neos, '. You should never see this error'
          CALL ctl_stop( ctmp1 )
          !
@@ -1662,4 +1682,13 @@
       r1_rcp      = 1._wp / rcp
       r1_rau0_rcp = 1._wp / rau0_rcp 
+      !
+      IF(lwp) THEN
+         IF( l_useCT )   THEN
+            WRITE(numout,*) '             model uses Conservative Temperature'
+            WRITE(numout,*) '             Important: model must be initialized with CT and SA fields'
+         ELSE
+            WRITE(numout,*) '             model does not use Conservative Temperature'
+         ENDIF
+      ENDIF
       !
       IF(lwp) WRITE(numout,*)
