Changeset 3893

Timestamp:

2013-04-24T16:00:59+02:00 (11 years ago)

Author:

gm

Message:

dev_r3858_CNRS3_Ediag: #927 keep only Vallis and UNESCO eos introduce alpha beta and their primitive

Location:

branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM

Files:

• CONFIG/ORCA2_LIM/EXP00/namelist (modified) (1 diff)
• NEMO/OPA_SRC/DOM/istate.F90 (modified) (4 diffs)
• NEMO/OPA_SRC/DYN/dynspg_flt.F90 (modified) (1 diff)
• NEMO/OPA_SRC/LDF/ldfslp.F90 (modified) (9 diffs)
• NEMO/OPA_SRC/TRA/eosbn2.F90 (modified) (32 diffs)
• NEMO/OPA_SRC/TRD/trdini.F90 (modified) (1 diff)
• NEMO/OPA_SRC/TRD/trdken.F90 (modified) (1 diff)
• NEMO/OPA_SRC/TRD/trdpen.F90 (modified) (4 diffs)
• NEMO/OPA_SRC/oce.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/step.F90 (modified) (1 diff)

branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist

@@ -564 +564 @@
 &nameos        !   ocean physical parameters
 !-----------------------------------------------------------------------
-   nn_eos      =   0       !  type of equation of state and Brunt-Vaisala frequency
-                           !     = 0, UNESCO (formulation of Jackett and McDougall (1994) and of McDougall (1987) )
-                           !     = 1, linear: rho(T)   = rau0 * ( 1.028 - ralpha * T )
-                           !     = 2, linear: rho(T,S) = rau0 * ( rbeta * S - ralpha * T )
-   rn_alpha    =   2.0e-4  !  thermal expension coefficient (nn_eos= 1 or 2)
-   rn_beta     =   7.7e-4  !  saline  expension coefficient (nn_eos= 2)
+   nn_eos      =  0        !  type of equation of state and Brunt-Vaisala frequency
+                           !     =-1, UNESCO (formulation of Jackett and McDougall (1995) )
+                           !     = 0, UNESCO (modified formulation of Jackett and McDougall (1994) and of McDougall (1987) )
+                           !     = 1, simplified eos (based on Vallis, 2006): 
+                           !             prd(T,S,z)  =  - rn_alpha * (1+rn_cabbel*(T-rn_temp0)+rn_thermo*z) * (T-rn_temp0) 
+                           !                            + rn_beta  * (S-rn_salt0) 
+                           !                            + rn_gamma * z
+   rn_alph0    =   2.0e-4  !  thermal expension coefficient (nn_eos= 1)
+   rn_beta0    =   7.7e-4  !  saline  expension coefficient (nn_eos= 1)
+   rn_gamm0    =   4.39e-6 !  depth expansion coeff. (=0 for linear eos)
+   rn_cabbel   =   2.5e-2  !  cabbeling coeff. (=0 for linear eos)
+   rn_thermo   =   1.1e-4  !  thermobaric coeff. (=0 for linear eos)
 /
 !-----------------------------------------------------------------------

branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

@@ -32 +32 @@
    USE phycst          ! physical constants
    USE dtatsd          ! data temperature and salinity   (dta_tsd routine)
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O library
    USE zpshde          ! partial step: hor. derivative (zps_hde routine)
    USE eosbn2          ! equation of state            (eos bn2 routine)
@@ -42 +40 @@
    USE dynspg_ts       ! pressure gradient schemes
    USE sol_oce         ! ocean solver variables
+   !
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O library
    USE lib_mpp         ! MPP library
    USE restart         ! restart
@@ -56 +57 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.5 , NEMO Consortium (2013)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -81 +82 @@
       CALL dta_tsd_init                       ! Initialisation of T & S input data
 
-      rhd  (:,:,:  ) = 0.e0
-      rhop (:,:,:  ) = 0.e0
-      rn2  (:,:,:  ) = 0.e0 
-      tsa  (:,:,:,:) = 0.e0    
+      rhd  (:,:,:  ) = 0._wp
+      rhop (:,:,:  ) = 0._wp
+      rn2  (:,:,:  ) = 0._wp
+      tsa  (:,:,:,:) = 0._wp   
+      rab_b(:,:,:,:) = 0._wp
+      rab_n(:,:,:,:) = 0._wp
 
       IF( ln_rstart ) THEN                    ! Restart from a file

branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90

@@ -404 +404 @@
       IF( lrst_oce )   CALL flt_rst( kt, 'WRITE' )      ! write filtered free surface arrays in restart file
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt')
+      IF( nn_timing == 1 )   CALL timing_stop('dyn_spg_flt')
       !
    END SUBROUTINE dyn_spg_flt

branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

@@ -28 +28 @@
    USE zdfmxl         ! mixed layer depth
    USE eosbn2         ! equation of states
+   !
+   USE in_out_manager ! I/O manager
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE in_out_manager ! I/O manager
    USE prtctl         ! Print control
    USE wrk_nemo       ! work arrays
@@ -393 +394 @@
       REAL(wp) ::   zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim
       REAL(wp) ::   zdzrho_raw
-      REAL(wp) ::   zbeta0
       REAL(wp), POINTER, DIMENSION(:,:)     ::   z1_mlbw
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalbet
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zdxrho , zdyrho, zdzrho     ! Horizontal and vertical density gradients
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zti_mlb, ztj_mlb            ! for Griffies operator only
@@ -403 +402 @@
       !
       CALL wrk_alloc( jpi,jpj, z1_mlbw )
-      CALL wrk_alloc( jpi,jpj,jpk, zalbet )
       CALL wrk_alloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho,              klstart = 0  )
       CALL wrk_alloc( jpi,jpj,  2,2, zti_mlb, ztj_mlb,        kkstart = 0, klstart = 0  )
@@ -410 +408 @@
       !  Some preliminary calculation  !
       !--------------------------------!
-      !
-      CALL eos_alpbet( tsb, zalbet, zbeta0 )  !==  before local thermal/haline expension ratio at T-points  ==!
       !
       DO jl = 0, 1                            !==  unmasked before density i- j-, k-gradients  ==!
@@ -419 +415 @@
             DO jj = 1, jpjm1                  ! NB: not masked ==>  a minimum value is set
                DO ji = 1, fs_jpim1            ! vector opt.
-                  zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! i-gradient of T & S at u-point
-                  zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
-                  zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! j-gradient of T & S at v-point
-                  zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
-                  zdxrho_raw = ( - zalbet(ji+ip,jj   ,jk) * zdit + zbeta0*zdis ) / e1u(ji,jj)
-                  zdyrho_raw = ( - zalbet(ji   ,jj+jp,jk) * zdjt + zbeta0*zdjs ) / e2v(ji,jj)
-                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX(   repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign
+                  zdit = ( tsb(ji+1,jj  ,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! i-gradient of T & S at u-point
+                  zdis = ( tsb(ji+1,jj  ,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
+                  zdjt = ( tsb(ji  ,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! j-gradient of T & S at v-point
+                  zdjs = ( tsb(ji  ,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
+                  zdxrho_raw = ( - rab_b(ji+ip,jj   ,jk,jp_tem) * zdit + rab_b(ji+ip,jj   ,jk,jp_sal) * zdis ) / e1u(ji,jj)
+                  zdyrho_raw = ( - rab_b(ji   ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji   ,jj+jp,jk,jp_sal) * zdjs ) / e2v(ji,jj)
+                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign
                   zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
                END DO
@@ -431 +427 @@
          END DO
          !
-         IF( ln_zps.and.l_grad_zps ) THEN     ! partial steps: correction of i- & j-grad on bottom
+         IF( ln_zps .AND. l_grad_zps ) THEN     ! partial steps: correction of i- & j-grad on bottom
 # if defined key_vectopt_loop
             DO jj = 1, 1
@@ -442 +438 @@
                   zdit = gtsu(ji,jj,jp_tem)   ;   zdjt = gtsv(ji,jj,jp_tem)      ! i- & j-gradient of Temperature
                   zdis = gtsu(ji,jj,jp_sal)   ;   zdjs = gtsv(ji,jj,jp_sal)      ! i- & j-gradient of Salinity
-                  zdxrho_raw = ( - zalbet(ji+ip,jj   ,iku) * zdit + zbeta0*zdis ) / e1u(ji,jj)
-                  zdyrho_raw = ( - zalbet(ji   ,jj+jp,ikv) * zdjt + zbeta0*zdjs ) / e2v(ji,jj)
+                  zdxrho_raw = ( - rab_b(ji+ip,jj   ,iku,jp_tem) * zdit + rab_b(ji+ip,jj   ,iku,jp_sal) * zdis ) / e1u(ji,jj)
+                  zdyrho_raw = ( - rab_b(ji   ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji   ,jj+jp,ikv,jp_sal) * zdjs ) / e2v(ji,jj)
                   zdxrho(ji+ip,jj   ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign
                   zdyrho(ji   ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
@@ -463 +459 @@
                      zdks = 0._wp
                   ENDIF
-                  zdzrho_raw = ( - zalbet(ji   ,jj   ,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp)
-                  zdzrho(ji   ,jj   ,jk,  kp) =     - MIN( - repsln,      zdzrho_raw )    ! force zdzrho >= repsln
+                  zdzrho_raw = ( - rab_b(ji,jj,jk,jp_tem) * zdkt + rab_b(ji,jj,jk,jp_sal) * zdks ) / fse3w(ji,jj,jk+kp)
+                  zdzrho(ji,jj,jk,kp) = - MIN( - repsln, zdzrho_raw )    ! force zdzrho >= repsln
                  END DO
             END DO
@@ -608 +604 @@
       !
       CALL wrk_dealloc( jpi,jpj, z1_mlbw )
-      CALL wrk_dealloc( jpi,jpj,jpk, zalbet )
       CALL wrk_dealloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho,              klstart = 0  )
       CALL wrk_dealloc( jpi,jpj,  2,2, zti_mlb, ztj_mlb,        kkstart = 0, klstart = 0  )

branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

@@ -17 +17 @@
    !!            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
+   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_ab used in ldfslp
    !!            3.5  ! 2012-03  (F. Roquet, G. Madec)  add pen_ddt_dds used in trdpen
    !!             -   ! 2012-05  (F. Roquet)  add Vallis and original JM95 equation of state
-   !!             -   ! 2012-05  (F. Roquet)  add eos_alpha_beta
+   !!             -   ! 2012-05  (F. Roquet)  add eos_alpha_beta to be used in dynhpg
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   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/haline expansion ratio
-   !!   eos_expansion_ratio : calculates the partial derivative of in situ density with respect to T and S
-   !!   tfreez         : Compute the surface freezing temperature
+   !!   eos_insitu     : density anomaly (rho/rho0 - 1)
+   !!   eos_insitu_pot : density and surface referenced potential density anomalies
+   !!   eos_insitu_2d  : density anomaly for 2d fields
+   !!   eos_bn2        : Brunt-Vaisala frequency
+   !!   eos_rab        : partial derivative of density anomaly with respect to T and S
+   !!   tfreez         : surface freezing temperature
    !!   eos_init       : set eos parameters (namelist)
    !!----------------------------------------------------------------------
@@ -49 +47 @@
    !                   !! * Interface
    INTERFACE eos
-      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
+      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, &
+         &                eos_ab, eos_insitu_ab
    END INTERFACE
    INTERFACE bn2
-      MODULE PROCEDURE eos_bn2
+      MODULE PROCEDURE eos_bn2, eos_bn2_ab
    END INTERFACE
 
    PUBLIC   eos            ! called by step, istate, tranpc and zpsgrd modules
    PUBLIC   eos_init       ! called by istate module
-   PUBLIC   bn2            ! called by step module
-   PUBLIC   eos_alpbet     ! called by ldfslp module
-   PUBLIC   eos_alpha_beta ! used for density diagnostics in dyntra
+   PUBLIC   bn2            ! called by step module, update alpha/beta and compute bn2
    PUBLIC   eos_pen        ! used for pe diagnostics in trdpen
    PUBLIC   tfreez         ! called by sbcice_... modules
 
    !                                          !!* Namelist (nameos) *
-   INTEGER , PUBLIC ::   nn_eos   = 0         !: = -1/0/1/2/3 type of eq. of state and associated Brunt-Vaisala frequ.
-   REAL(wp), PUBLIC ::   rn_alpha = 2.0e-4_wp !: thermal expension coeff. (linear equation of state)
-   REAL(wp), PUBLIC ::   rn_beta  = 7.7e-4_wp !: saline  expension coeff. (linear equation of state)
-
-   REAL(wp) ::   vallis_ratio = 0   ! 1027/rau0
-   REAL(wp) ::   vallis_diff  = 0   ! (1027-rau0)/rau0
+   INTEGER , PUBLIC ::   nn_eos   = 0         !: = -1/0/1 type of eq. of state and Brunt-Vaisala frequ.
+
+   REAL(wp) ::   rn_alph0  = 1.67e-4 !: thermal expansion coeff. (simplified equation of state)
+   REAL(wp) ::   rn_beta0  = 7.8e-4 !: saline  expansion coeff. (simplified equation of state)
+   REAL(wp) ::   rn_gamm0  = 4.4e-6 !: depth expansion coeff. (simplified equation of state)
+   REAL(wp) ::   rn_cabbel = 3.0e-2 !: cabbeling coeff. (simplified equation of state)
+   REAL(wp) ::   rn_thermo = 1.1e-4 !: thermobaric coeff. (simplified equation of state)
+   REAL(wp) ::   offset    = 0      !: offset applied on Vallis density anomaly, so that rho(10,35,0)=1027
 
    !! * Substitutions
@@ -89 +88 @@
       !!       defined through the namelist parameter nn_eos.
       !!
-      !! ** Method  :   5 cases:
+      !! ** Method  :   3 cases:
       !!      nn_eos = -1 : Jackett and McDougall (1995) equation of state.
       !!         Check value: rho = 1041.83267 kg/m**3 for p=3000 dbar,
@@ -100 +99 @@
       !!         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
-      !!         with pressure                      p        decibars
+      !!              prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0
+      !!         with depth                         z        meters
       !!              potential temperature         t        deg celsius
       !!              salinity                      s        psu
@@ -107 +106 @@
       !!              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,
+      !!         Check value: rho = 1060.93298 kg/m**3 for z=10000 dbar,
       !!          t = 40 deg celcius, s=40 psu
-      !!      nn_eos = 1 : linear equation of state function of temperature only
-      !!              prd(t) = 0.0285 - rn_alpha * t
-      !!      nn_eos = 2 : linear equation of state function of temperature and
-      !!               salinity
-      !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
-      !!      nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006)
-      !!              prd(t,s,p) =  ( rho(t,s,p) - rau0 ) / rau0
-      !!              where the pressure p in decibars is approximated by the depth in meters.
+      !!      nn_eos = 1 : Vallis simplified equation of state (Vallis book, 2006)
+      !!              prd(t,s,z) =  ( rho(t,s,p) - rau0 ) / rau0
+      !!                         = -alpha*(1+cabbel*(T-T0)+thermo*z)*(T-T0) + beta*(S-S0) + gamma*z
+      !!              linear case function of T only: rn_alpha<>0, other coefficients = 0
+      !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0
+      !!              Vallis (2006) equation: use default values of coefficients 
       !!      Note that no boundary condition problem occurs in this routine
       !!      as pts are defined over the whole domain.
@@ -122 +119 @@
       !! ** Action  :   compute prd , the in situ density (no units)
       !!
-      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
-      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995
-      !!----------------------------------------------------------------------
+      !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1995
+      !!                Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
+      !!----------------------------------------------------------------------
+      !!
       REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
       !                                                      ! 2 : salinity               [psu]
       REAL(wp), DIMENSION(:,:,:)  , INTENT(  out) ::   prd   ! in situ density            [-]
-      !
+      !!
       INTEGER  ::   ji, jj, jk           ! dummy loop indices
       REAL(wp) ::   zt , zs , zh , zsr   ! local scalars
-      REAL(wp) ::   zr1, zr2, zr3, zr4   !   -      -
-      REAL(wp) ::   zrhop, ze, zbw, zb   !   -      -
-      REAL(wp) ::   zd , zc , zaw, za    !   -      -
-      REAL(wp) ::   zb1, za1, zkw, zk0   !   -      -
+      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      -
+      REAL(wp) ::   zn,  za1, za2, za    !   -      -
+      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      -
+      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
       !!----------------------------------------------------------------------
+
       !
       IF( nn_timing == 1 ) CALL timing_start('eos')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zws )
-      !
       SELECT CASE( nn_eos )
       !
       CASE( -1 )                !==  Jackett and McDougall (1995) formulation  ==!
+         !
+         CALL wrk_alloc( jpi, jpj, jpk, zws )
          !
 !CDIR NOVERRCHK
@@ -158 +157 @@
                   !
                   ! compute volumic mass pure water at atm pressure
-                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
                   ! seawater volumic mass atm pressure
-                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
                      &                      -4.0899e-3_wp ) *zt+0.824493_wp
-                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
-                  zr4= 4.8314e-4_wp
-                  !
-                  ! potential volumic mass (reference to the surface)
-                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
+                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+                  zn4= 4.8314e-4_wp
+                  zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1
                   !
                   ! add the compression terms
-                  ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
-                  zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
-                  zb = zbw + ze * zs
-                  !
-                  zd = -1.480266e-4_wp
-                  zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3
-                  zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp
-                  za = ( zd*zsr + zc ) *zs + zaw
-                  !
-                  zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp
-                  za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp
-                  zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp
-                  zk0= ( zb1*zsr + za1 )*zs + zkw
+                  za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
+                  za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
+                  za = za1 + za2 * zs
+                  !
+                  zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp
+                  zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3
+                  zb3 = 1.480266e-4_wp
+                  zb  = ( zb3*zsr + zb2 ) *zs + zb1
+                  !
+                  zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp
+                  zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp
+                  zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp
+                  zc  = ( zc3*zsr + zc2 )*zs + zc1
                   !
                   ! masked in situ density anomaly. Important: rau0=1035, should be 1027!
-                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
-                     &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk)
+                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   &
+                     &             - 1._wp  ) * tmask(ji,jj,jk)
                END DO
             END DO
          END DO
          !
+         CALL wrk_dealloc( jpi, jpj, jpk, zws )
+         !
       CASE( 0 )                !==  modified Jackett and McDougall (1995) formulation  ==!
+         !
+         CALL wrk_alloc( jpi, jpj, jpk, zws )
          !
 !CDIR NOVERRCHK
@@ -205 +206 @@
                   !
                   ! compute volumic mass pure water at atm pressure
-                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
                   ! seawater volumic mass atm pressure
-                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
                      &                      -4.0899e-3_wp ) *zt+0.824493_wp
-                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
-                  zr4= 4.8314e-4_wp
-                  !
-                  ! potential volumic mass (reference to the surface)
-                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
+                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+                  zn4= 4.8314e-4_wp
+                  zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1
                   !
                   ! add the compression terms
-                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
-                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
-                  zb = zbw + ze * zs
-                  !
-                  zd = -2.042967e-2_wp
-                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
-                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp
-                  za = ( zd*zsr + zc ) *zs + zaw
-                  !
-                  zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp
-                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp
-                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
-                  zk0= ( zb1*zsr + za1 )*zs + zkw
+                  za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
+                  za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
+                  za = za1 + za2 * zs
+                  !
+                  zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp
+                  zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp
+                  zb3= 2.042967e-2_wp
+                  zb = ( zb3*zsr + zb2 ) *zs + zb1
+                  !
+                  zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
+                  zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp
+                  zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp
+                  zc = ( zc3*zsr + zc2 )*zs + zc1
                   !
                   ! masked in situ density anomaly
-                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
-                     &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk)
+                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   &
+                     &             - 1._wp  ) * tmask(ji,jj,jk)
                END DO
             END DO
          END DO
          !
-      CASE( 1 )                !==  Linear formulation function of temperature only  ==!
-         DO jk = 1, jpkm1
-            prd(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
-         END DO
-         !
-      CASE( 2 )                !==  Linear formulation function of temperature and salinity  ==!
-         DO jk = 1, jpkm1
-            prd(:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
-         END DO
-         !
-      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==!
+         CALL wrk_dealloc( jpi, jpj, jpk, zws )
+         !
+      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==!
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -256 +247 @@
                   zh = fsdept(ji,jj,jk)        ! depth
                   ! masked in situ density anomaly
-                  prd(ji,jj,jk) = vallis_diff + vallis_ratio * (                                             &
-                     &   -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh  &
-                     &   ) * tmask(ji,jj,jk)
+                  prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt   &                                          &
+                     &   + rn_beta0 * zs + rn_gamm0 * zh  ) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -267 +257 @@
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos  : ', ovlap=1, kdim=jpk )
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zws )
-      !
       IF( nn_timing == 1 ) CALL timing_stop('eos')
       !
@@ -278 +266 @@
       !!                  ***  ROUTINE eos_insitu_pot  ***
       !!
-      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the
+      !! ** Purpose :   Compute the in situ density (ratio (rho-rau0)/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  :   5 cases:
-      !!      nn_eos = -1 : Jackett and McDougall (1995) equation of state.
-      !!      nn_eos = 0 : standard NEMO equation of state, based on 
-      !!         Jackett and McDougall (1995) 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)
-      !!      nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006)
-      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
-      !!              rhop(t,s)  = rho(t,s,0)
-      !!      Note that no boundary condition problem occurs in this routine
-      !!      as (tn,sn) or (ta,sa) are defined over the whole domain.
+      !! ** Method  :   Same as for eos_insitu
       !!
       !! ** Action  : - prd  , the in situ density (no units)
@@ -332 +286 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
       !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw   ! local scalars
-      REAL(wp) ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0               !   -      -
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars
+      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      -
+      REAL(wp) ::   zn,  za1, za2, za    !   -      -
+      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      -
+      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
       !!----------------------------------------------------------------------
@@ -357 +314 @@
                   !
                   ! compute volumic mass pure water at atm pressure
-                  zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   &
-                     &                          -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp
+                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
                   ! seawater volumic mass atm pressure
-                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   &
-                     &                                         -4.0899e-3_wp ) *zt+0.824493_wp
-                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
-                  zr4= 4.8314e-4_wp
-                  !
-                  ! potential volumic mass (reference to the surface)
-                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
-                  !
-                  ! save potential volumic mass
-                  prhop(ji,jj,jk) = zrhop * tmask(ji,jj,jk)
+                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                     &                      -4.0899e-3_wp ) *zt+0.824493_wp
+                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+                  zn4= 4.8314e-4_wp
+                  zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1
                   !
                   ! add the compression terms
-                  ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
-                  zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
-                  zb = zbw + ze * zs
-                  !
-                  zd = -1.480266e-4_wp
-                  zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3
-                  zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp
-                  za = ( zd*zsr + zc ) *zs + zaw
-                  !
-                  zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp
-                  za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp
-                  zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp
-                  zk0= ( zb1*zsr + za1 )*zs + zkw
-                  !
-                  ! masked in situ density anomaly
-                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
-                     &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk)
+                  za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
+                  za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
+                  za = za1 + za2 * zs
+                  !
+                  zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp
+                  zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3
+                  zb3 = 1.480266e-4_wp
+                  zb  = ( zb3*zsr + zb2 ) *zs + zb1
+                  !
+                  zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp
+                  zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp
+                  zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp
+                  zc  = ( zc3*zsr + zc2 )*zs + zc1
+                  !
+                  ! masked in situ density anomaly. Important: rau0=1035, should be 1027!
+                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   &
+                     &             - 1._wp  ) * tmask(ji,jj,jk)
+              prhop(ji,jj,jk) = zn * tmask(ji,jj,jk)
                END DO
             END DO
@@ -406 +359 @@
                   !
                   ! compute volumic mass pure water at atm pressure
-                  zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   &
-                     &                          -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp
+                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
                   ! seawater volumic mass atm pressure
-                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   &
-                     &                                         -4.0899e-3_wp ) *zt+0.824493_wp
-                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
-                  zr4= 4.8314e-4_wp
-                  !
-                  ! potential volumic mass (reference to the surface)
-                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
-                  !
-                  ! save potential volumic mass
-                  prhop(ji,jj,jk) = zrhop * tmask(ji,jj,jk)
+                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                     &                      -4.0899e-3_wp ) *zt+0.824493_wp
+                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+                  zn4= 4.8314e-4_wp
+                  zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1
                   !
                   ! add the compression terms
-                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
-                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
-                  zb = zbw + ze * zs
-                  !
-                  zd = -2.042967e-2_wp
-                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
-                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp
-                  za = ( zd*zsr + zc ) *zs + zaw
-                  !
-                  zb1=   (  -0.1909078_wp  *zt+7.390729_wp    ) *zt-55.87545_wp
-                  za1= ( (   2.326469e-3_wp*zt+1.553190_wp    ) *zt-65.00517_wp ) *zt + 1044.077_wp
-                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
-                  zk0= ( zb1*zsr + za1 )*zs + zkw
+                  za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
+                  za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
+                  za = za1 + za2 * zs
+                  !
+                  zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp
+                  zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp
+                  zb3= 2.042967e-2_wp
+                  zb = ( zb3*zsr + zb2 ) *zs + zb1
+                  !
+                  zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
+                  zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp
+                  zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp
+                  zc = ( zc3*zsr + zc2 )*zs + zc1
                   !
                   ! masked in situ density anomaly
-                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
-                     &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk)
+                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   &
+                     &             - 1._wp  ) * tmask(ji,jj,jk)
+              prhop(ji,jj,jk) = zn * tmask(ji,jj,jk)
                END DO
             END DO
          END DO
          !
-      CASE( 1 )                !==  Linear formulation = F( temperature )  ==!
-         DO jk = 1, jpkm1
-            prd  (:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk)
-            prhop(:,:,jk) = ( 1.e0_wp   +            prd (:,:,jk)       ) * rau0 * tmask(:,:,jk)
-         END DO
-         !
-      CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
-         DO jk = 1, jpkm1
-            prd  (:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk)
-            prhop(:,:,jk) = ( 1.e0_wp  + prd (:,:,jk) )                                       * rau0 * tmask(:,:,jk)
-         END DO
-         !
-      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==!
+      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==!
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -462 +399 @@
                   zh = fsdept(ji,jj,jk)        ! depth
                   ! masked in situ density anomaly
-                  prd(ji,jj,jk) = vallis_diff + vallis_ratio * (                                             &
-                     &   -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh  &
-                     &   ) * tmask(ji,jj,jk)
+                  prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt   &                                          &
+                     &   + rn_beta0 * zs + rn_gamm0 * zh  ) * tmask(ji,jj,jk)
                   ! masked in situ density anomaly
-                 prhop(ji,jj,jk) = ( 1.0_wp - ( 1.67e-4_wp + 0.5e-5_wp*zt ) *zt + 0.78e-3_wp *zs )   &
-                     &              * 1027._wp * tmask(ji,jj,jk)
+                 prhop(ji,jj,jk) = ( 1.0_wp + offset - rn_alph0 * ( 1._wp + rn_cabbel*zt ) *zt + rn_beta0 *zs )  &
+                     &              * rau0 * tmask(ji,jj,jk)
                  !
                END DO
@@ -492 +428 @@
       !!      defined through the namelist parameter nn_eos. * 2D field case
       !!
-      !! ** Method  :   5 cases:
-      !!      nn_eos = -1 : Jackett and McDougall (1995) equation of state.
-      !!      nn_eos = 0 : standard NEMO equation of state, based on 
-      !!         Jackett and McDougall (1995) 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
-      !!         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) = 0.0285 - rn_alpha * t
-      !!      nn_eos = 2 : linear equation of state function of temperature and
-      !!               salinity
-      !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
-      !!      nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006)
-      !!              prd(t,s,p) =  ( rho(t,s,p) - rau0 ) / rau0
-      !!              where the pressure p in decibars is approximated by the depth in meters.
-      !!      Note that no boundary condition problem occurs in this routine
-      !!      as pts are defined over the whole domain.
+      !! ** Method  :   Same as for eos_insitu
       !!
       !! ** Action  : - prd , the in situ density (no units)
@@ -532 +441 @@
       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density
       !!
-      INTEGER  ::   ji, jj                    ! dummy loop indices
-      REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw   ! temporary scalars
-      REAL(wp) ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zmask        !    -         -
+      INTEGER  ::   ji, jj              ! dummy loop indices
+      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars
+      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      -
+      REAL(wp) ::   zn,  za1, za2, za    !   -      -
+      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      -
+      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      -
       REAL(wp), POINTER, DIMENSION(:,:) :: zws
       !!----------------------------------------------------------------------
@@ -558 +470 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask
                zt    = pts  (ji,jj,jp_tem)     ! interpolated T
                zs    = pts  (ji,jj,jp_sal)     ! interpolated S
@@ -565 +476 @@
                !
                ! compute volumic mass pure water at atm pressure
-               zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   &
-                  &                        -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp
+               zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                  &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
                ! seawater volumic mass atm pressure
-               zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt   &
-                  &                                   -4.0899e-3_wp ) *zt+0.824493_wp
-               zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp
-               zr4 = 4.8314e-4_wp
-               !
-               ! potential volumic mass (reference to the surface)
-               zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
+               zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                  &                      -4.0899e-3_wp ) *zt+0.824493_wp
+               zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+               zn4= 4.8314e-4_wp
+               zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1
                !
                ! add the compression terms
-               ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
-               zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
-               zb = zbw + ze * zs
+               za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
+               za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
+               za = za1 + za2 * zs
                !
-               zd = -1.480266e-4_wp
-               zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3
-               zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp
-               za = ( zd*zsr + zc ) *zs + zaw
+               zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp
+               zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3
+               zb3 = 1.480266e-4_wp
+               zb  = ( zb3*zsr + zb2 ) *zs + zb1
                !
-               zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp
-               za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp
-               zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp
-               zk0= ( zb1*zsr + za1 )*zs + zkw
+               zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp
+               zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp
+               zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp
+               zc  = ( zc3*zsr + zc2 )*zs + zc1
                !
                ! masked in situ density anomaly
-               prd(ji,jj) = ( zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask
+               prd(ji,jj) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   &
+                  &             - 1._wp  ) * tmask(ji,jj,1)
             END DO
          END DO
@@ -607 +517 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask
                zt    = pts  (ji,jj,jp_tem)     ! interpolated T
                zs    = pts  (ji,jj,jp_sal)     ! interpolated S
@@ -614 +523 @@
                !
                ! compute volumic mass pure water at atm pressure
-               zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   &
-                  &                        -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp
+               zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                  &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
                ! seawater volumic mass atm pressure
-               zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt   &
-                  &                                   -4.0899e-3_wp ) *zt+0.824493_wp
-               zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp
-               zr4 = 4.8314e-4_wp
-               !
-               ! potential volumic mass (reference to the surface)
-               zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
+               zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                  &                      -4.0899e-3_wp ) *zt+0.824493_wp
+               zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+               zn4= 4.8314e-4_wp
+               zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1
                !
                ! add the compression terms
-               ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
-               zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
-               zb = zbw + ze * zs
+               za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
+               za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
+               za = za1 + za2 * zs
                !
-               zd =    -2.042967e-2_wp
-               zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
-               zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt -4.721788_wp
-               za = ( zd*zsr + zc ) *zs + zaw
+               zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp
+               zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp
+               zb3= 2.042967e-2_wp
+               zb = ( zb3*zsr + zb2 ) *zs + zb1
                !
-               zb1=     (-0.1909078_wp  *zt+7.390729_wp      ) *zt-55.87545_wp
-               za1=   ( ( 2.326469e-3_wp*zt+1.553190_wp      ) *zt-65.00517_wp ) *zt+1044.077_wp
-               zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp   ) *zt-30.41638_wp ) *zt   &
-                  &                             +2098.925_wp ) *zt+190925.6_wp
-               zk0= ( zb1*zsr + za1 )*zs + zkw
+               zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
+               zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp
+               zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp
+               zc = ( zc3*zsr + zc2 )*zs + zc1
                !
                ! masked in situ density anomaly
-               prd(ji,jj) = ( zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask
-            END DO
-         END DO
-         !
-      CASE( 1 )                !==  Linear formulation = F( temperature )  ==!
+               prd(ji,jj) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   &
+                  &             - 1._wp  ) * tmask(ji,jj,1)
+            END DO
+         END DO
+         !
+      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==!
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1)
-            END DO
-         END DO
-         !
-      CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1)
-            END DO
-         END DO
-         !
-      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==!
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask
                zt    = pts  (ji,jj,jp_tem) - 10._wp   ! interpolated T
                zs    = pts  (ji,jj,jp_sal) - 35._wp   ! interpolated S
                zh    = pdep (ji,jj)            ! depth at the partial step level
                ! masked in situ density anomaly
-               prd(ji,jj) = vallis_diff + vallis_ratio * (                                                &
-                  &   -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh  &
-                  &   ) * zmask
+               prd(ji,jj) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt  &                                           &
+                     &   + rn_beta0 * zs + rn_gamm0 * zh  ) * tmask(ji,jj,1)
                !
             END DO
@@ -686 +577 @@
 
 
-   SUBROUTINE eos_bn2( pts, pn2 )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE eos_bn2  ***
-      !!
-      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the time-
-      !!      step of the input arguments
-      !!
-      !! ** Method  :   5 cases:
-      !!       * nn_eos = -1  : The brunt-vaisala frequency is computed for
-      !!       the Jackett and McDougall (1995) equation of state:
-      !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
-      !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
-      !!      computed and used in zdfddm module :
-      !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
-      !!       * nn_eos = 0  : The brunt-vaisala frequency is based on the
-      !!            formulation of McDougall (1987).
-      !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
-      !!      computed and used in zdfddm module :
-      !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
-      !!       * nn_eos = 1  : linear equation of state (temperature only)
-      !!            N^2 = grav * rn_alpha * dk[ t ]/e3w
-      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
-      !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
-      !!      The use of potential density to compute N^2 introduces e r r o r
-      !!      in the sign of N^2 at great depths. We recommand the use of
-      !!      nn_eos = 0, except for academical studies.
-      !!        Macro-tasked on horizontal slab (jk-loop)
-      !!      N.B. N^2 is set to zero at the first level (JK=1) in inidtr
-      !!      and is never used at this level.
-      !!
-      !! ** Action  : - pn2 : the brunt-vaisala frequency
-      !!
-      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
-      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995
-      !!                McDougall, J. Phys. Oceano., 1987
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
-      !                                                               ! 2 : salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pn2   ! Brunt-Vaisala frequency    [s-1]
-      !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zgde3w, zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop      ! local scalars
-      REAL(wp) ::   ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM    !   -      -
-      REAL(wp) ::   zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom   !   -      -
-      REAL(wp) ::   zalpbet, zalpha, zbeta                                  !   -      -
-#if defined key_zdfddm
-      REAL(wp) ::   zds   ! local scalars
-#endif
+
+   SUBROUTINE eos_ab( pts, pab )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE eos_ab  ***
+      !!
+      !! ** Purpose :   Calculates the in situ thermal/haline expansion terms at T-points
+      !!
+      !! ** Method  :   calculates alpha and beta at T-points
+      !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state
+      !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state
+      !!       * nn_eos = 1  : Vallis equation of state (Vallis 2006)
+      !!                       alpha and beta ratios are returned as 3-D arrays.
+      !!
+      !! ** Action  : - pab : partial derivative of in situ density with respect to T & S at T-points
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! alpha and beta
+      !
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars
+      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      -
+      REAL(wp) ::   zn,  za1, za2, za    !   -      -
+      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      -
+      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      -
+      REAL(wp) ::   zm_inv, zdm, zdn     !   -      -
+      REAL(wp) ::   zdza, zdzb, zdzc     !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
       !!----------------------------------------------------------------------
-
-      !
-      IF( nn_timing == 1 ) CALL timing_start('bn2')
+      !
+      IF ( nn_timing == 1 )   CALL timing_start('eos_ab')
       !
       CALL wrk_alloc( jpi, jpj, jpk, zws )
       !
-      ! pn2 : interior points only (2=< jk =< jpkm1 )
-      ! --------------------------
-      !
-      SELECT CASE( nn_eos )
-      !
-      CASE( -1 )                !==  Jackett and McDougall (1995)  ==!
-         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
-         DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zgde3w = grav / fse3w(ji,jj,jk)
-                  zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )         ! potential temperature at w-pt
-                  zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) )         ! salinity at w-pt
-                  zsr = 0.5 * ( zws(ji,jj,jk) + zws(ji,jj,jk-1) )        ! square root of S at w-pt
-                  zh = fsdepw(ji,jj,jk)                                                ! depth in meters  at w-point
-                  !
-                  ! compute volumic mass pure water at atm pressure
-                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
-                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
-                  ! seawater volumic mass atm pressure
-                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
-                     &                      -4.0899e-3_wp ) *zt+0.824493_wp
-                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
-                  zr4= 4.8314e-4_wp
-                  !
-                  ! potential volumic mass (reference to the surface)
-                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
-                  !
-                  ! add the compression terms
-                  ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
-                  zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
-                  zb = zbw + ze * zs
-                  !
-                  zd = -1.480266e-4_wp
-                  zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3
-                  zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp
-                  za = ( zd*zsr + zc ) *zs + zaw
-                  !
-                  zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp
-                  za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp
-                  zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp
-                  zk0= ( zb1*zsr + za1 )*zs + zkw
-                  !
-                  zM= ( zb*zh - za )*zh + zk0
-                  zDenom= 1._wp - zh / zM
-                  ! compute in-situ density
-                  zrho = zrhop/zDenom
-              ! alpha
-                  zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs
-                  zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp)
-                  zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  &
-                     &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp)
-                  zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  &
-                     &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  &
-                     &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp)
-                  zdM = (zdzb*zh-zdza)*zh+zdzk0
-                  zdDenom = zh * zdM / (zM*zM)
-                  zalpha = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0
-                  ! beta
-                  zdzb  = ze
-                  zdza  = -3.0644505e-2_wp*zsr+zc
-                  zdzk0 = 1.5_wp*zb1*zsr+za1
-                  zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2
-                  zdM = (zdzb*zh-zdza)*zh+zdzk0
-                  zdDenom = zh * zdM / (zM*zM)
-                  zbeta =   ( zdrhop - zrho*zdDenom ) / zDenom / rau0
-                  ! alpha/beta
-                  zalpbet = zalpha / zbeta
-                  ! N2
-                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2
-                     &          * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   &
-                     &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) )
-#if defined key_zdfddm
-                  !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!!
-                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s])
-                  IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp
-                  rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
-#endif
-               END DO
-            END DO
-         END DO
-         !
-      CASE( 0 )                !==  McDougall (1987) formulation  ==!
-         DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zgde3w = grav / fse3w(ji,jj,jk)
-                  zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )         ! potential temperature at w-pt
-                  zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0  ! salinity anomaly (s-35) at w-pt
-                  zh = fsdepw(ji,jj,jk)                                                ! depth in meters  at w-point
-                  !
-                  zalpbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt  &   ! ratio alpha/beta
-                     &                                  - 0.203814e-03_wp ) * zt   &
-                     &                                  + 0.170907e-01_wp ) * zt   &
-                     &   +         0.665157e-01_wp                                 &
-                     &   +     ( - 0.678662e-05_wp * zs                            &
-                     &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   &
-                     &   +   ( ( - 0.302285e-13_wp * zh                            &
-                     &           - 0.251520e-11_wp * zs                            &
-                     &           + 0.512857e-12_wp * zt * zt              ) * zh   &
-                     &           - 0.164759e-06_wp * zs                            &
-                     &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   &
-                     &                                  + 0.380374e-04_wp ) * zh
-                     !
-                  zbeta  = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt      &   ! beta
-                     &                               - 0.301985e-05_wp ) * zt      &
-                     &   +       0.785567e-03_wp                                   &
-                     &   + (     0.515032e-08_wp * zs                              &
-                     &         + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs     &
-                     &   + ( (   0.121551e-17_wp * zh                              &
-                     &         - 0.602281e-15_wp * zs                              &
-                     &         - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh     &
-                     &                                + 0.408195e-10_wp   * zs     &
-                     &     + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt     &
-                     &                                - 0.121555e-07_wp ) * zh
-                     !
-                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2
-                     &          * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   &
-                     &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) )
-#if defined key_zdfddm
-                  !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!!
-                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s])
-                  IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp
-                  rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
-#endif
-               END DO
-            END DO
-         END DO
-         !
-      CASE( 1 )                !==  Linear formulation = F( temperature )  ==!
-         DO jk = 2, jpkm1
-            pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) )        &
-              &                / fse3w(:,:,jk) * tmask(:,:,jk)
-         END DO
-         !
-      CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
-         DO jk = 2, jpkm1
-            pn2(:,:,jk) = grav * (  rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) )      &
-               &                  - rn_beta  * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) )  )   &
-               &               / fse3w(:,:,jk) * tmask(:,:,jk)
-         END DO
-#if defined key_zdfddm
-         DO jk = 2, jpkm1                                 ! Rrau = (alpha / beta) (dk[t] / dk[s])
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
-                  IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp
-                  rrau(ji,jj,jk) = rn_alpha / rn_beta * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
-               END DO
-            END DO
-         END DO
-#endif
-         !
-      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==!
-         DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zgde3w = grav / fse3w(ji,jj,jk)
-                  zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) - 10._wp ! potential temperature at w-pt
-                  zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35._wp ! salinity anomaly (s-35) at w-pt
-                  zh = fsdepw(ji,jj,jk)                                                 ! depth in meters  at w-point
-                  !
-                  zalpha = 1027._wp * ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) / rau0 ! alpha
-                  zbeta  = 1027._wp * 0.78e-3_wp / rau0  ! beta
-                  zalpbet = zalpha / zbeta  ! ratio alpha/beta
-                  !
-                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2
-                     &          * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   &
-                     &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) )
-#if defined key_zdfddm
-                  !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!!
-                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s])
-                  IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp
-                  rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
-#endif
-               END DO
-            END DO
-         END DO
-         !
-      END SELECT
-
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk )
-#if defined key_zdfddm
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk )
-#endif
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, zws )
-      !
-      IF( nn_timing == 1 ) CALL timing_stop('bn2')
-      !
-   END SUBROUTINE eos_bn2
-
-
-   SUBROUTINE eos_alpbet( pts, palpbet, beta0 )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE eos_alpbet  ***
-      !!
-      !! ** Purpose :   Calculates the in situ thermal/haline expansion ratio at T-points
-      !!
-      !! ** Method  :   calculates alpha / beta ratio at T-points
-      !!       * nn_eos = -1  : alpha / beta ratio is computed for
-      !!       the Jackett and McDougall (1995) equation of state:
-      !!                       Scalar beta0 is returned = 1.
-      !!       * nn_eos = 0  : alpha / beta ratio using the formulation of McDougall (1987).
-      !!                       Scalar beta0 is returned = 1.
-      !!       * nn_eos = 1  : linear equation of state (temperature only)
-      !!                       The ratio is undefined, so we return alpha as palpbet
-      !!                       Scalar beta0 is returned = 0.
-      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
-      !!                       The alpha/beta ratio is returned as palpbet
-      !!                       Scalar beta0 is returned = 1.
-      !!       * nn_eos = 3  : Vallis equation of state (Vallis 2006)
-      !!                       Scalar beta0 is returned = 1.
-      !!
-      !! ** Action  : - palpbet : thermal/haline expansion ratio at T-points
-      !!            :   beta0   : 1. or 0.
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts       ! pot. temperature & salinity
-      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palpbet   ! thermal/haline expansion ratio
-      REAL(wp),                              INTENT(  out) ::   beta0     ! set = 1 except with case 1 eos, rho=rho(T)
-      !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop            ! temporary scalars
-      REAL(wp) ::   ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM  !    -         -
-      REAL(wp) ::   zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom !
-      REAL(wp) ::   zalpha, zbeta                                ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 ) CALL timing_start('eos_alpbet')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, zws )
-      !
       SELECT CASE ( nn_eos )
       !
       CASE ( -1 )               ! Jackett and McDougall (1995) formulation
+         !
          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
          DO jk = 1, jpkm1
@@ -988 +624 @@
                   !
                   ! compute volumic mass pure water at atm pressure
-                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
                   ! seawater volumic mass atm pressure
-                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
                      &                      -4.0899e-3_wp ) *zt+0.824493_wp
-                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
-                  zr4= 4.8314e-4_wp
-                  !
-                  ! potential volumic mass (reference to the surface)
-                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
+                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+                  zn4= 4.8314e-4_wp
+                  zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1
                   !
                   ! add the compression terms
-                  ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
-                  zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
-                  zb = zbw + ze * zs
-                  !
-                  zd = -1.480266e-4_wp
-                  zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3
-                  zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp
-                  za = ( zd*zsr + zc ) *zs + zaw
-                  !
-                  zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp
-                  za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp
-                  zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp
-                  zk0= ( zb1*zsr + za1 )*zs + zkw
-                  !
-                  zM= ( zb*zh - za )*zh + zk0
-                  zDenom= 1._wp - zh / zM
-                  ! compute in-situ density
-                  zrho = zrhop/zDenom
+                  za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
+                  za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
+                  za = za1 + za2 * zs
+                  !
+                  zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp
+                  zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3
+                  zb3 = 1.480266e-4_wp
+                  zb  = ( zb3*zsr + zb2 ) *zs + zb1
+                  !
+                  zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp
+                  zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp
+                  zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp
+                  zc  = ( zc3*zsr + zc2 )*zs + zc1
+                  !
+                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc )
               ! alpha
-                  zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs
-                  zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp)
-                  zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  &
-                     &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp)
-                  zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  &
-                     &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  &
-                     &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp)
-                  zdM = (zdzb*zh-zdza)*zh+zdzk0
-                  zdDenom = zh * zdM / (zM*zM)
-                  zalpha = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0
+                  zdza  = (2.78936e-8_wp*zt-1.202016e-6_wp) + (12.414646e-11_wp*zt+6.128773e-9_wp ) * zs
+                  zdzb  = (4.118662e-7_wp*zt-1.847318e-4_wp)*zs +   &
+                     &        ( 5.869245e-6_wp*zt-5.969284e-4_wp)*zt+2.212276e-2_wp
+                  zdzc  = ( (-9.239848e-3_wp*zt+9.085835e-2_wp)*zsr +   &
+                     &       (-15.252564e-4_wp*zt+12.566526e-2_wp)*zt-3.101089_wp ) * zs + &
+                     &       ((-16.761012e-4_wp*zt+28.946112e-2_wp)*zt-34.12206_wp )*zt+1444.304_wp
+                     
+                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  &
+                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  &
+                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt  &
+                     &                  -18.19058e-3_wp)*zt+6.793952e-2_wp)
+                  zdm = (zdza*zh+zdzb)*zh+zdzc
+                  !
+                  pab(ji,jj,jk,jp_tem) =  - ( &
+                     &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) &
+                     &   / rau0 * tmask(ji,jj,jk)
+                  !
                   ! beta
-                  zdzb  = ze
-                  zdza  = -3.0644505e-2_wp*zsr+zc
-                  zdzk0 = 1.5_wp*zb1*zsr+za1
-                  zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2
-                  zdM = (zdzb*zh-zdza)*zh+zdzk0
-                  zdDenom = zh * zdM / (zM*zM)
-                  zbeta =   ( zdrhop - zrho*zdDenom ) / zDenom / rau0
-                  ! alpha/beta
-                  palpbet(ji,jj,jk) = zalpha / zbeta * tmask(ji,jj,jk)
+                  zdza  = za2
+                  zdzb  = 2.220399e-4_wp*zsr+zb2
+                  zdzc  = 1.5_wp*zc3*zsr+zc2
+                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2
+                  zdm   = (zdza*zh+zdzb)*zh+zdzc
+                  pab(ji,jj,jk,jp_sal) =    ( &
+                     &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) &
+                     &   / rau0 * tmask(ji,jj,jk)
+                  !
                END DO
             END DO
          END DO
-         beta0 = 1._wp
-         !
-      CASE ( 0 )               ! McDougall (1987) formulation
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt = pts(ji,jj,jk,jp_tem)           ! potential temperature
-                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35)
-                  zh = fsdept(ji,jj,jk)               ! depth in meters
-                  !
-                  palpbet(ji,jj,jk) =                                              &
-                     &     ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &
-                     &                                  - 0.203814e-03_wp ) * zt   &
-                     &                                  + 0.170907e-01_wp ) * zt   &
-                     &   + 0.665157e-01_wp                                         &
-                     &   +     ( - 0.678662e-05_wp * zs                            &
-                     &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   &
-                     &   +   ( ( - 0.302285e-13_wp * zh                            &
-                     &           - 0.251520e-11_wp * zs                            &
-                     &           + 0.512857e-12_wp * zt * zt              ) * zh   &
-                     &           - 0.164759e-06_wp * zs                            &
-                     &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   &
-                     &                                  + 0.380374e-04_wp ) * zh
-               END DO
-            END DO
-         END DO
-         beta0 = 1._wp
-         !
-      CASE ( 1 )              !==  Linear formulation = F( temperature )  ==!
-         palpbet(:,:,:) = rn_alpha
-         beta0 = 0._wp
-         !
-      CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==!
-         palpbet(:,:,:) = rn_alpha / rn_beta
-         beta0 = 1._wp
-         !
-      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==!
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature
-                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35)
-                  zh = fsdepw(ji,jj,jk)                                                ! depth in meters  at w-point
-                  !
-                  zalpha = ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) ! alpha/vallis_ratio
-                  zbeta  = 0.78e-3_wp  ! beta/vallis_ratio
-                  !
-                  palpbet(ji,jj,jk) = zalpha / zbeta * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         beta0 = 1._wp
-         !
-      CASE DEFAULT
-         IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
-         nstop = nstop + 1
-         !
-      END SELECT
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, zws )
-      !
-      IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet')
-      !
-   END SUBROUTINE eos_alpbet
-
-
-   SUBROUTINE eos_alpha_beta( pts, palpha, pbeta )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE eos_alpha_beta  ***
-      !!
-      !! ** Purpose :   Calculates the in situ thermal/haline expansion terms at T-points
-      !!
-      !! ** Method  :   calculates alpha and beta at T-points
-      !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state
-      !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state
-      !!       * nn_eos = 1  : linear equation of state (temperature only)
-      !!                       palpha = rn_alpha
-      !!                       pbeta = 0
-      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
-      !!                       palpha = rn_alpha
-      !!                       pbeta = rn_beta
-      !!       * nn_eos = 3  : Vallis equation of state (Vallis 2006)
-      !!                       alpha and beta ratios are returned as 3-D arrays.
-      !!
-      !! ** Action  : - palpha : thermal expansion ratio at T-points
-      !!            : - pbeta  : haline expansion ratio at T-points
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts      ! pot. temperature & salinity
-      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palpha   ! alpha
-      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pbeta    ! beta
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop            ! temporary scalars
-      REAL(wp) ::   ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM  !    -         -
-      REAL(wp) ::   zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom !
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
-      !!----------------------------------------------------------------------
-      !
-      IF ( nn_timing == 1 )   CALL timing_start('eos_alpha_beta')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, zws )
-      !
-      SELECT CASE ( nn_eos )
-      !
-      CASE ( -1 )               ! Jackett and McDougall (1995) formulation
+         !
+      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation
          !
          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
@@ -1159 +693 @@
                   !
                   ! compute volumic mass pure water at atm pressure
-                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
                   ! seawater volumic mass atm pressure
-                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
                      &                      -4.0899e-3_wp ) *zt+0.824493_wp
-                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
-                  zr4= 4.8314e-4_wp
-                  !
-                  ! potential volumic mass (reference to the surface)
-                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
+                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+                  zn4= 4.8314e-4_wp
+                  zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1
                   !
                   ! add the compression terms
-                  ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
-                  zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
-                  zb = zbw + ze * zs
-                  !
-                  zd = -1.480266e-4_wp
-                  zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3
-                  zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp
-                  za = ( zd*zsr + zc ) *zs + zaw
-                  !
-                  zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp
-                  za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp
-                  zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp
-                  zk0= ( zb1*zsr + za1 )*zs + zkw
-                  !
-                  zM= ( zb*zh - za )*zh + zk0
-                  zDenom= 1._wp - zh / zM
-                  ! compute in-situ density
-                  zrho = zrhop/zDenom
+                  za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
+                  za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
+                  za = za1 + za2 * zs
+                  !
+                  zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp
+                  zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp
+                  zb3= 2.042967e-2_wp
+                  zb = ( zb3*zsr + zb2 ) *zs + zb1
+                  !
+                  zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
+                  zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp
+                  zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp
+                  zc = ( zc3*zsr + zc2 )*zs + zc1
+                  !
+                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc )
               ! alpha
-                  zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs
-                  zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp)
-                  zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  &
-                     &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp)
-                  zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  &
-                     &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  &
-                     &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp)
-                  zdM = (zdzb*zh-zdza)*zh+zdzk0
-                  zdDenom = zh * zdM / (zM*zM)
-                  palpha(ji,jj,jk) = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk)
+                  zdza  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs
+                  zdzb  = (14.535852e-5_wp*zt-2.598241e-3)*zs+(-17.81973e-6_wp*zt-5.025098e-3_wp)*zt+0.1028859_wp
+                  zdzc  = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  &
+                     &    +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp)
+                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  &
+                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  &
+                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp)
+                  zdm = (zdza*zh+zdzb)*zh+zdzc
+                  !
+                  pab(ji,jj,jk,jp_tem) =  - ( &
+                     &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) &
+                     &   / rau0 * tmask(ji,jj,jk)
                   ! beta
-                  zdzb  = ze
-                  zdza  = -3.0644505e-2_wp*zsr+zc
-                  zdzk0 = 1.5_wp*zb1*zsr+za1
-                  zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2
-                  zdM = (zdzb*zh-zdza)*zh+zdzk0
-                  zdDenom = zh * zdM / (zM*zM)
-                  pbeta(ji,jj,jk) =   ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk)
+                  zdza  = za2
+                  zdzb  = 3.0644505e-2_wp*zsr+zb2
+                  zdzc  = 1.5_wp*zc3*zsr+zc2
+                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2
+                  zdm   = (zdza*zh+zdzb)*zh+zdzc
+                  pab(ji,jj,jk,jp_sal) =    ( &
+                     &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) &
+                     &   / rau0 * tmask(ji,jj,jk)
+                  !
                END DO
             END DO
          END DO
          !
-      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation
+      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==!
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature anomaly (t-T0)
+                  zh = fsdept(ji,jj,jk)               ! depth in meters at t-point
+                  !
+                  pab(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*rn_cabbel*zt + rn_thermo*zh ) * tmask(ji,jj,jk)   ! alpha
+               END DO
+            END DO
+         END DO
+         pab(:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:)   ! beta
+         !
+      CASE DEFAULT
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
+         nstop = nstop + 1
+         !
+      END SELECT
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zws )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('eos_ab')
+      !
+   END SUBROUTINE eos_ab
+
+
+
+   SUBROUTINE eos_insitu_ab( pts, prd, pab )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE eos_insitu_ab  ***
+      !!
+      !! ** Purpose :   Calculates the in situ thermal/haline expansion terms
+      !!                  and the insitu density anomaly at T-points
+      !!
+      !! ** Method  :   calculates rhd, alpha, beta at T-points
+      !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state
+      !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state
+      !!       * nn_eos = 1  : Vallis equation of state (Vallis 2006)
+      !!                       alpha and beta ratios are returned as 3-D arrays.
+      !!
+      !! ** Action  : - pab : thermal and haline expansion ratios at T-points
+      !!              - prd , the density anomaly (no units)
+      !!
+      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
+      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity    [Celcius,psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   prd   ! density anomaly                [-]
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! partial derivative of density anomaly
+      !                                                               ! with respect to T and S        [1/Celcius,1/psu]
+      !
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars
+      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      -
+      REAL(wp) ::   zn,  za1, za2, za    !   -      -
+      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      -
+      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      -
+      REAL(wp) ::   zm_inv, zdm, zdn     !   -      -
+      REAL(wp) ::   zdza, zdzb, zdzc     !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
+      !!----------------------------------------------------------------------
+      !
+      IF ( nn_timing == 1 )   CALL timing_start('eos_insitu_ab')
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zws )
+      !
+      SELECT CASE ( nn_eos )
+      !
+      CASE ( -1 )               ! Jackett and McDougall (1995) formulation
          !
          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
@@ -1224 +826 @@
                   !
                   ! compute volumic mass pure water at atm pressure
-                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
                   ! seawater volumic mass atm pressure
-                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
                      &                      -4.0899e-3_wp ) *zt+0.824493_wp
-                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
-                  zr4= 4.8314e-4_wp
-                  !
-                  ! potential volumic mass (reference to the surface)
-                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
+                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+                  zn4= 4.8314e-4_wp
+                  zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1
                   !
                   ! add the compression terms
-                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
-                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
-                  zb = zbw + ze * zs
-                  !
-                  zd = -2.042967e-2_wp
-                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
-                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp
-                  za = ( zd*zsr + zc ) *zs + zaw
-                  !
-                  zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp
-                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp
-                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
-                  zk0= ( zb1*zsr + za1 )*zs + zkw
-                  !
-                  zM= ( zb*zh - za )*zh + zk0
-                  zDenom= 1._wp - zh / zM
-                  ! compute in-situ density
-                  zrho = zrhop/zDenom
+                  za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
+                  za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
+                  za = za1 + za2 * zs
+                  !
+                  zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp
+                  zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3
+                  zb3 = 1.480266e-4_wp
+                  zb  = ( zb3*zsr + zb2 ) *zs + zb1
+                  !
+                  zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp
+                  zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp
+                  zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp
+                  zc  = ( zc3*zsr + zc2 )*zs + zc1
+                  !
+                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc )
+                  !
+                  ! masked in situ density anomaly. Important: rau0=1035, should be 1027!
+                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh * zm_inv ) / rau0  - 1._wp  ) * tmask(ji,jj,jk)
+                  !
               ! alpha
-                  zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs
-                  zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp)
-                  zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  &
-                     &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp)
-                  zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  &
-                     &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  &
-                     &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp)
-                  zdM = (zdzb*zh-zdza)*zh+zdzk0
-                  zdDenom = zh * zdM / (zM*zM)
-                  palpha(ji,jj,jk) = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk)
+                  zdza  = (2.78936e-8_wp*zt-1.202016e-6_wp) + (12.414646e-11_wp*zt+6.128773e-9_wp ) * zs
+                  zdzb  = (4.118662e-7_wp*zt-1.847318e-4_wp)*zs +   &
+                     &        ( 5.869245e-6_wp*zt-5.969284e-4_wp)*zt+2.212276e-2_wp
+                  zdzc  = ( (-9.239848e-3_wp*zt+9.085835e-2_wp)*zsr +   &
+                     &       (-15.252564e-4_wp*zt+12.566526e-2_wp)*zt-3.101089_wp ) * zs + &
+                     &       ((-16.761012e-4_wp*zt+28.946112e-2_wp)*zt-34.12206_wp )*zt+1444.304_wp
+                     
+                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  &
+                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  &
+                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt  &
+                     &                  -18.19058e-3_wp)*zt+6.793952e-2_wp)
+                  zdm = (zdza*zh+zdzb)*zh+zdzc
+                  !
+                  pab(ji,jj,jk,jp_tem) =  - ( &
+                     &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) &
+                     &   / rau0 * tmask(ji,jj,jk)
+                  !
                   ! beta
-                  zdzb  = ze
-                  zdza  = -3.0644505e-2_wp*zsr+zc
-                  zdzk0 = 1.5_wp*zb1*zsr+za1
-                  zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2
-                  zdM = (zdzb*zh-zdza)*zh+zdzk0
-                  zdDenom = zh * zdM / (zM*zM)
-                  pbeta(ji,jj,jk) =   ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk)
+                  zdza  = za2
+                  zdzb  = 2.220399e-4_wp*zsr+zb2
+                  zdzc  = 1.5_wp*zc3*zsr+zc2
+                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2
+                  zdm   = (zdza*zh+zdzb)*zh+zdzc
+                  pab(ji,jj,jk,jp_sal) =    ( &
+                     &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) &
+                     &   / rau0 * tmask(ji,jj,jk)
+                  !
                END DO
             END DO
          END DO
          !
-      CASE ( 1 )              !==  Linear formulation = F( temperature )  ==!
-         palpha(:,:,:) = rn_alpha
-         pbeta (:,:,:) = 0._wp
-         !
-      CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==!
-         palpha(:,:,:) = rn_alpha
-         pbeta (:,:,:) = rn_beta
-         !
-      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==!
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature (t-10)
-                  zh = fsdept(ji,jj,jk)                                                ! depth in meters  at w-point
-                  palpha(ji,jj,jk) = vallis_ratio * &
-                      &  ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) * tmask(ji,jj,jk)   ! alpha
-                  !
-               END DO
-            END DO
-         END DO
-         pbeta (:,:,:) = vallis_ratio * 0.78e-3_wp * tmask(:,:,:)   ! beta
-         !
-      CASE DEFAULT
-         IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
-         nstop = nstop + 1
-         !
-      END SELECT
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, zws )
-      !
-      IF( nn_timing == 1 ) CALL timing_stop('eos_alpha_beta')
-      !
-   END SUBROUTINE eos_alpha_beta
-
-
-
-   SUBROUTINE eos_pen( pts, palphaPE, pbetaPE, pPEanom )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE eos_pen  ***
-      !!
-      !! ** Purpose :   Calculates alpha_PE, beta_PE and PE at T-points
-      !!
-      !! ** Method  :   PE is defined analytically as the vertical 
-      !!                   primitive of EOS times -g integrated between 0 and z>0.
-      !!                PEanom is the PE anomaly: PEanom = ( PE - rau0 gz ) / rau0 gz
-      !!                                           = -1/z * /int_z^0 rd , where rd density anomaly
-      !!                dalphaPE and dbetaPE are partial derivatives of PE anomaly with respect to T and S:
-      !!                    alphaPE = - 1/(rau0 gz) * dPE/dT = - dPEanom/dT
-      !!                    betaPE  =   1/(rau0 gz) * dPE/dS = - dPEanom/dS
-      !!
-      !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state.
-      !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state.
-      !!       * nn_eos = 1  : linear equation of state (temperature only)
-      !!                       palpha = rau0*g*rn_alpha*z
-      !!                       pbeta = 0
-      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
-      !!                       palpha = rau0*g*rn_alpha*z
-      !!                       pbeta = -rau0*g*rn_beta*z
-      !!       * nn_eos = 3  : Vallis equation of state (Vallis 2006)
-      !!
-      !! ** Action  : - pPEanom     : PE anomaly given at T-points
-      !!            : - palphaPE    : given at T-points
-      !!            : - pbetaPE     : given at T-points
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts       ! pot. temperature & salinity
-      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palphaPE  ! partial derivative of PE anom
-      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pbetaPE   ! with respect to T and S, resp.
-      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pPEanom   ! potential energy anomaly
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop            ! temporary scalars
-      REAL(wp) ::   ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM  !    -         -
-      REAL(wp) ::   zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom !
-      REAL(wp) ::   zap1, zdelta, zsdelta, zAp, zBp, zlogBp, zCp, zatanCp !
-      REAL(wp) ::   zddelta, zdAp, zdBp, zdlogBp, zdCp, zP, zdP, zEp      !
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
-      !!----------------------------------------------------------------------
-      !
-      IF ( nn_timing == 1 )   CALL timing_start('eos_pen')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, zws )
-      !
-      SELECT CASE ( nn_eos )
-      !
-      CASE ( -1 )               ! Jackett and McDougall (1995) formulation
+      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation
          !
          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
@@ -1374 +899 @@
                   !
                   ! compute volumic mass pure water at atm pressure
-                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
                   ! seawater volumic mass atm pressure
-                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
                      &                      -4.0899e-3_wp ) *zt+0.824493_wp
-                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
-                  zr4= 4.8314e-4_wp
-                  !
-                  ! potential volumic mass (reference to the surface)
-                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
+                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+                  zn4= 4.8314e-4_wp
+                  zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1
                   !
                   ! add the compression terms
-                  ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
-                  zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
-                  zb = zbw + ze * zs
-                  !
-                  zd = -1.480266e-4_wp
-                  zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3
-                  zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp
-                  za = ( zd*zsr + zc ) *zs + zaw
-                  !
-                  zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp
-                  za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp
-                  zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp
-                  zk0= ( zb1*zsr + za1 )*zs + zkw
-                  !
-                  zM= ( zb*zh - za )*zh + zk0
-                  zDenom= 1._wp - zh / zM
-                  ! compute in-situ density
-                  zrho = zrhop/zDenom
-                  !
-                  zap1    = 1._wp + za
-                  zdelta  = 4._wp*zb*zk0 - zap1*zap1
-                  zsdelta = sqrt( abs( zdelta ) )
-                  zAp     = zap1 / zb / zsdelta
-                  zBp     = ( zM - zh ) / zk0
-                  zlogBp  = log(abs(zBp))
-                  zCp     = zh * zsdelta / (2._wp*zk0-zh*zap1)
-                  zatanCp = atan(zCp)
-                  zP      = zh + zAp * zatanCp + 0.5_wp*zlogBp/zb
-                  ! compute potential energy
-                  pPEanom(ji,jj,jk)     = ( ( zrhop * zP ) / rau0 / zh - 1._wp ) * tmask(ji,jj,jk) !!wrong
-                  !
-              ! dPEdt
-                  zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs
-                  zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp)
-                  zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  &
-                     &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp)
-                  zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  &
-                     &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  &
-                     &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp)
-                  zdM = (zdzb*zh-zdza)*zh+zdzk0
-                  zdDenom = zh * zdM / (zM*zM)
-                  !
-                  zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza
-                  zdAp    = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta)
-                  zdlogBp = zdM/(zM-zh) - zdzk0/zk0
-                  zdCp    = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1))
-                  zdP     = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp) 
-                  !
-                  palphaPE(ji,jj,jk)    = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 
-                  !
-              ! dPEds
-                  zdzb  = ze
-                  zdza  = -3.0644505e-2_wp*zsr+zc
-                  zdzk0 = 1.5_wp*zb1*zsr+za1
-                  zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2
-                  zdM = (zdzb*zh-zdza)*zh+zdzk0
-                  zdDenom = zh * zdM / (zM*zM)
-                  !
-                  zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza
-                  zdAp    = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta)
-                  zdlogBp = zdM/(zM-zh) - zdzk0/zk0
-                  zdCp    = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1))
-                  zdP     = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp) 
-                  !
-                  pbetaPE(ji,jj,jk)    = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk)
-                  !
+                  za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
+                  za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
+                  za = za1 + za2 * zs
+                  !
+                  zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp
+                  zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp
+                  zb3= 2.042967e-2_wp
+                  zb = ( zb3*zsr + zb2 ) *zs + zb1
+                  !
+                  zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
+                  zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp
+                  zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp
+                  zc = ( zc3*zsr + zc2 )*zs + zc1
+                  !
+                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc )
+                  !
+                  ! masked in situ density anomaly (rho/rho0 - 1). Important: rau0=1035, should be 1027!
+                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh * zm_inv ) * r1_rau0  - 1._wp  ) * tmask(ji,jj,jk)
+                  !
+                  !                                 ! alpha
+                  zdza  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs
+                  zdzb  = (14.535852e-5_wp*zt-2.598241e-3)*zs+(-17.81973e-6_wp*zt-5.025098e-3_wp)*zt+0.1028859_wp
+                  zdzc  = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  &
+                     &    +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp)
+                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  &
+                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  &
+                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp)
+                  zdm = (zdza*zh+zdzb)*zh+zdzc
+                  !
+                  pab(ji,jj,jk,jp_tem) =  - (  zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv )   &
+                     &                    * r1_rau0 * tmask(ji,jj,jk)
+                  !
+                  !                                 ! beta
+                  zdza  = za2
+                  zdzb  = 3.0644505e-2_wp*zsr+zb2
+                  zdzc  = 1.5_wp*zc3*zsr+zc2
+                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2
+                  zdm   = (zdza*zh+zdzb)*zh+zdzc
+                  pab(ji,jj,jk,jp_sal) = (  zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv )   &
+                     &                 * r1_rau0 * tmask(ji,jj,jk)
                END DO
             END DO
          END DO
          !
-      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation
+      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==!
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature (t-T0)
+                  zh = fsdept(ji,jj,jk)                                                ! depth in meters  at t-point
+                  ! masked in situ density anomaly
+                  prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt   &                                          &
+                     &   + rn_beta0 * zs + rn_gamm0 * zh  ) * tmask(ji,jj,jk)
+                  ! alpha
+                  pab(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*rn_cabbel*zt + rn_thermo*zh )  &
+                      &                       * tmask(ji,jj,jk)   ! alpha
+                  !
+               END DO
+            END DO
+         END DO
+         ! beta
+         pab (:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:)   ! beta
+         !
+      CASE DEFAULT
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
+         nstop = nstop + 1
+         !
+      END SELECT
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zws )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('eos_insitu_ab')
+      !
+   END SUBROUTINE eos_insitu_ab
+
+
+
+   SUBROUTINE eos_pen( pts, pab_pe, pen )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE eos_pen  ***
+      !!
+      !! ** Purpose :   Calculates alpha_PE, beta_PE and PE at T-points
+      !!
+      !! ** Method  :   PE is defined analytically as the vertical 
+      !!                   primitive of EOS times -g integrated between 0 and z>0.
+      !!                pen is the PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz
+      !!                                           = -1/z * /int_z^0 rd , where rd density anomaly
+      !!                ab_pe are partial derivatives of PE anomaly with respect to T and S:
+      !!                    ab_pe(1) = - 1/(rau0 gz) * dPE/dT = - d(pen)/dT
+      !!                    ab_pe(2) =   1/(rau0 gz) * dPE/dS = - d(pen)/dS
+      !!
+      !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state.
+      !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state.
+      !!       * nn_eos = 1  : Vallis equation of state (Vallis 2006)
+      !!
+      !! ** Action  : - pen         : PE anomaly given at T-points
+      !!            : - pab_pe  : given at T-points
+      !!                    pab_pe(:,:,:,jp_tem) is alpha_pe
+      !!                    pab_pe(:,:,:,jp_sal) is beta_pe
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts     ! pot. temperature & salinity
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab_pe  ! alpha_pe and beta_pe
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pen     ! potential energy anomaly
+      !
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars
+      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      -
+      REAL(wp) ::   zn,  za1, za2, za    !   -      -
+      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      -
+      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      -
+      REAL(wp) ::   zm_inv, zdm, zdn     !   -      -
+      REAL(wp) ::   zdza, zdzb, zdzc     !   -      -
+      REAL(wp) ::   zdelta, zsdelta      !   -      -
+      REAL(wp) ::   zsgn, zAp, zlogBp    !   -      -
+      REAL(wp) ::   zCp, zatCp, zP       !   -      -
+      REAL(wp) ::   zddelta, zdAp        !   -      -
+      REAL(wp) ::   zdlogBp, zdCp, zdP   !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
+      !!----------------------------------------------------------------------
+      !
+      IF ( nn_timing == 1 )   CALL timing_start('eos_pen')
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zws )
+      !
+      SELECT CASE ( nn_eos )
+      !
+      CASE ( -1 )               ! Jackett and McDougall (1995) formulation
          !
          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
@@ -1464 +1044 @@
                   zt = pts   (ji,jj,jk,jp_tem)
                   zs = pts   (ji,jj,jk,jp_sal)
-                  zh = fsdept(ji,jj,jk)        ! depth
+                  zh = fsdept(ji,jj,jk)        ! depth (Important: zh must be different from 0)
                   zsr= zws   (ji,jj,jk)        ! square root salinity
                   !
                   ! compute volumic mass pure water at atm pressure
-                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
                   ! seawater volumic mass atm pressure
-                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
                      &                      -4.0899e-3_wp ) *zt+0.824493_wp
-                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
-                  zr4= 4.8314e-4_wp
-                  !
-                  ! potential volumic mass (reference to the surface)
-                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
+                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+                  zn4= 4.8314e-4_wp
+                  zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1
                   !
                   ! add the compression terms
-                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
-                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
-                  zb = zbw + ze * zs
-                  !
-                  zd = -2.042967e-2_wp
-                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
-                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp
-                  za = ( zd*zsr + zc ) *zs + zaw
-                  !
-                  zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp
-                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp
-                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
-                  zk0= ( zb1*zsr + za1 )*zs + zkw
-                  !
-                  zM= ( zb*zh - za )*zh + zk0
-                  zDenom= 1._wp - zh / zM
-                  ! compute in-situ density
-                  zrho = zrhop/zDenom
-                  !
-                  zap1    = 1._wp + za
-                  zdelta  = 4._wp*zb*zk0 - zap1*zap1
-                  zsdelta = sqrt( abs( zdelta ) )
-                  zAp     = zap1 / zb / zsdelta
-                  zBp     = ( zM - zh ) / zk0
-                  zlogBp  = log(abs(zBp))
-                  zCp     = zh * zsdelta / (2._wp*zk0-zh*zap1)
-                  zatanCp = atan(zCp)
-                  zP      = zh + zAp * zatanCp + 0.5_wp*zlogBp/zb
+                  za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp
+                  za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp
+                  za = za1 + za2 * zs
+                  !
+                  zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp
+                  zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3
+                  zb3 = 1.480266e-4_wp
+                  zb  = ( zb3*zsr + zb2 ) *zs + zb1
+                  !
+                  zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp
+                  zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp
+                  zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp
+                  zc  = ( zc3*zsr + zc2 )*zs + zc1
+                  !
+                  zdelta = 4._wp * za * zc - zb * zb
+                  ! ddt
+                  zdza  = (2.78936e-8_wp*zt-1.202016e-6_wp) + (12.414646e-11_wp*zt+6.128773e-9_wp ) * zs
+                  zdzb  = (4.118662e-7_wp*zt-1.847318e-4_wp)*zs +   &
+                     &        ( 5.869245e-6_wp*zt-5.969284e-4_wp)*zt+2.212276e-2_wp
+                  zdzc  = ( (-9.239848e-3_wp*zt+9.085835e-2_wp)*zsr +   &
+                     &       (-15.252564e-4_wp*zt+12.566526e-2_wp)*zt-3.101089_wp ) * zs + &
+                     &       ((-16.761012e-4_wp*zt+28.946112e-2_wp)*zt-34.12206_wp )*zt+1444.304_wp
+                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  &
+                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  &
+                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt  &
+                     &                  -18.19058e-3_wp)*zt+6.793952e-2_wp)
+                  zdm = (zdza*zh+zdzb)*zh+zdzc
+                  !
+                  ! stability condition for the calculation (add a small perturbation on T if needed)
+                  IF( ABS(zdelta) < 1.e-5 .OR. ABS(za)<1.e-9 ) THEN ;
+                     zt = zt + 1.e-3_wp
+                     !
+                     zn = zn + 1.e-3_wp * zdn
+                     za = za + 1.e-3_wp * zdza
+                     zb = zb + 1.e-3_wp * zdzb
+                     zc = zc + 1.e-3_wp * zdzc
+                     !
+                     zdelta = 4._wp * za * zc - zb * zb
+                     !
+                  ENDIF
+                  !
+                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc )
+                  zsgn   = SIGN( 1. , zdelta )
+                  zsdelta = SQRT( zsgn * zdelta )
+                  zAp     = - zb / za / zsdelta
+                  zlogBp  = - LOG( zc * zm_inv )
+                  zCp     = zh * zsdelta / (2._wp*zc+zh*zb)
+                  IF( zsgn > 0. ) THEN ;   zatCp = ATAN(zCp);
+                  ELSE  ;    zatCp = ATANH(zCp)
+                  ENDIF
+                  zP      = ( 0.5_wp*zlogBp/za + zAp * zatCp ) / zh
+                  !
                   ! compute potential energy
-                  pPEanom(ji,jj,jk)     = - grav * zrhop * zP * tmask(ji,jj,jk)
-                  !
-              ! dPEdt
-                  zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs
-                  zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp)
-                  zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  &
-                     &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp)
-                  zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  &
-                     &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  &
-                     &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp)
-                  zdM = (zdzb*zh-zdza)*zh+zdzk0
-                  zdDenom = zh * zdM / (zM*zM)
-                  !
-                  zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza
-                  zdAp    = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta)
-                  zdlogBp = zdM/(zM-zh) - zdzk0/zk0
-                  zdCp    = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1))
-                  zdP     = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp) 
-                  !
-                  palphaPE(ji,jj,jk)    = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk)
-                  !
-              ! dPEds
-                  zdzb  = ze
-                  zdza  = -3.0644505e-2_wp*zsr+zc
-                  zdzk0 = 1.5_wp*zb1*zsr+za1
-                  zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2
-                  zdM = (zdzb*zh-zdza)*zh+zdzk0
-                  zdDenom = zh * zdM / (zM*zM)
-                  !
-                  zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza
-                  zdAp    = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta)
-                  zdlogBp = zdM/(zM-zh) - zdzk0/zk0
-                  zdCp    = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1))
-                  zdP     = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp) 
-                  !
-                  pbetaPE(ji,jj,jk)    = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk)
+                  pen(ji,jj,jk)  = ( zn - rau0 + zn * zP ) / rau0 * tmask(ji,jj,jk)
+                  !
+              ! alphaPE
+                  zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta
+                  zdAp    = zAp * (zdzb/zb - zdza/za - zddelta)
+                  zdlogBp = zdm*zm_inv - zdzc/zc
+                  zdCp    = zCp * (  zddelta - (2._wp*zdzc+zh*zdzb) / (2._wp*zc+zh*zb)  )
+                  zdP     = (  0.5_wp * ( zdlogBp - zlogBp*zdza/za ) / za + zdAp*zatCp   &
+                     &       + zAp*zdCp / ( 1._wp + zsgn*zCp*zCp )  ) / zh
+                  !
+                  pab_pe(ji,jj,jk,jp_tem)    = - ( zdn * (1._wp + zP ) + zn * zdP ) / rau0 * tmask(ji,jj,jk)
+                  ! dds
+                  zdza  = za2
+                  zdzb  = 2.220399e-4_wp*zsr+zb2
+                  zdzc  = 1.5_wp*zc3*zsr+zc2
+                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2
+                  zdm   = (zdza*zh+zdzb)*zh+zdzc
+                  ! betaPE
+                  zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta
+                  zdAp    = zAp * (zdzb/zb - zdza/za - zddelta)
+                  zdlogBp = zdm*zm_inv - zdzc/zc
+                  zdCp    = zCp * (zddelta - (2._wp*zdzc+zh*zdzb)/(2._wp*zc+zh*zb))
+                  zdP     = ( 0.5_wp * ( zdlogBp - zlogBp*zdza/za ) / za + zdAp*zatCp + &
+                     &          zAp*zdCp/(1._wp+zsgn*zCp*zCp) ) / zh
+                  !
+                  pab_pe(ji,jj,jk,jp_sal)    = ( zdn * (1._wp + zP ) + zn * zdP )    &
+                                                    &           / rau0 * tmask(ji,jj,jk)
                   !
                END DO
@@ -1550 +1143 @@
          END DO
          !
-      CASE ( 1 )              !==  Linear formulation = F( temperature )  ==!
-         palphaPE(:,:,:) = rn_alpha * tmask(:,:,:)
-         pbetaPE (:,:,:) = 0._wp
-         DO jk = 1, jpkm1
-            pPEanom(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
-         END DO
-         !
-      CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==!
-         palphaPE(:,:,:) = rn_alpha * tmask(:,:,:)
-         pbetaPE (:,:,:) = rn_beta  * tmask(:,:,:)
-         DO jk = 1, jpkm1
-            pPEanom(:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
-         END DO
-         !
-      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==!
+      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation
+         !
+         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
          DO jk = 1, jpkm1
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature (t-10)
-                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35)
-                  zh = fsdept(ji,jj,jk)                                                ! depth in meters  at w-point
-                  !
-                  palphaPE(ji,jj,jk) = vallis_ratio * &
-                      &  ( 1.67e-4_wp * ( 1._wp + 0.55e-4_wp*zh ) + 1.e-5_wp*zt ) * tmask(ji,jj,jk)   ! alphaPE
-                  pPEanom(ji,jj,jk)     = vallis_diff + vallis_ratio *   &
-                      &   ( - ( 1.67e-4_wp*(1._wp+0.55e-4_wp*zh) + 0.5e-5_wp*zt )*zt + 0.78e-3_wp*zs + 2.195e-6_wp*zh ) &
-                      &   * tmask(ji,jj,jk)
+                  !
+                  zt = pts   (ji,jj,jk,jp_tem)
+                  zs = pts   (ji,jj,jk,jp_sal)
+                  zh = fsdept(ji,jj,jk)        ! depth (Important: zh must be different from 0)
+                  zsr= zws   (ji,jj,jk)        ! square root salinity
+                  !
+                  ! compute volumic mass pure water at atm pressure
+                  zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
+                  ! seawater volumic mass atm pressure
+                  zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                     &                      -4.0899e-3_wp ) *zt+0.824493_wp
+                  zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+                  zn4= 4.8314e-4_wp
+                  zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1
+                  !
+                  ! add the compression terms
+                  za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
+                  za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
+                  za = za1 + za2 * zs
+                  !
+                  zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp
+                  zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp
+                  zb3= 2.042967e-2_wp
+                  zb = ( zb3*zsr + zb2 ) *zs + zb1
+                  !
+                  zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
+                  zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp
+                  zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp
+                  zc = ( zc3*zsr + zc2 )*zs + zc1
+                  !
+                  zdelta = 4._wp * za * zc - zb * zb
+                  !
+                  ! dds
+                  zdza  = za2
+                  zdzb  = 3.0644505e-2_wp*zsr+zb2
+                  zdzc  = 1.5_wp*zc3*zsr+zc2
+                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2
+                  zdm   = (zdza*zh+zdzb)*zh+zdzc
+                  ! stability condition for the calculation (add a small perturbation on S if needed)
+                  IF( ABS(zdelta) < 1.e-5 .OR. ABS(za)<1.e-9 ) THEN ; !stability condition for the calculation
+                     zs = zs + 1.e-3_wp
+                     zsr= zsr + 0.5e-3_wp/zsr
+                     !
+                     zn = zn + 1.e-3_wp * zdn
+                     za = za + 1.e-3_wp * zdza
+                     zb = zb + 1.e-3_wp * zdzb
+                     zc = zc + 1.e-3_wp * zdzc
+                     !
+                     zdelta = 4._wp * za * zc - zb * zb
+                     !
+                  ENDIF
+                  !
+                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc )
+                  zsgn   = SIGN( 1. , zdelta )
+                  zsdelta = SQRT( zsgn * zdelta )
+                  zAp     = - zb / za / zsdelta
+                  zlogBp  = - LOG( zc * zm_inv )
+                  zCp     = zh * zsdelta / (2._wp*zc+zh*zb)
+                  IF( zsgn > 0. ) THEN ;   zatCp = ATAN(zCp);
+                  ELSE  ;    zatCp = ATANH(zCp)
+                  ENDIF
+                  zP      = ( 0.5_wp*zlogBp/za + zAp * zatCp ) / zh
+                  !
+                  ! compute potential energy
+                  pen(ji,jj,jk)  = ( zn - rau0 + zn * zP ) / rau0 * tmask(ji,jj,jk)
+                  !
+                  ! betaPE
+                  zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta
+                  zdAp    = zAp * (zdzb/zb - zdza/za - zddelta)
+                  zdlogBp = zdm*zm_inv - zdzc/zc
+                  zdCp    = zCp * (zddelta - (2._wp*zdzc+zh*zdzb)/(2._wp*zc+zh*zb))
+                  zdP     = ( 0.5_wp * ( zdlogBp - zlogBp*zdza/za ) / za + zdAp*zatCp + &
+                     &          zAp*zdCp/(1._wp+zsgn*zCp*zCp) ) / zh
+                  !
+                  pab_pe(ji,jj,jk,jp_sal)    = ( zdn * (1._wp + zP ) + zn * zdP )    &
+                                                    &           / rau0 * tmask(ji,jj,jk)
+                  !
+              ! ddt
+                  zdza  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs
+                  zdzb  = (14.535852e-5_wp*zt-2.598241e-3)*zs+(-17.81973e-6_wp*zt-5.025098e-3_wp)*zt+0.1028859_wp
+                  zdzc  = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  &
+                     &    +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp)
+                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  &
+                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  &
+                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp)
+                  zdm = (zdza*zh+zdzb)*zh+zdzc
+              ! alphaPE
+                  zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta
+                  zdAp    = zAp * (zdzb/zb - zdza/za - zddelta)
+                  zdlogBp = zdm*zm_inv - zdzc/zc
+                  zdCp    = zCp * (zddelta - (2._wp*zdzc+zh*zdzb)/(2._wp*zc+zh*zb))
+                  zdP     = ( 0.5_wp * ( zdlogBp - zlogBp*zdza/za ) / za + zdAp*zatCp + &
+                     &          zAp*zdCp/(1._wp+zsgn*zCp*zCp) ) / zh
+                  !
+                  pab_pe(ji,jj,jk,jp_tem) = - ( zdn * (1._wp + zP ) + zn * zdP ) / rau0 * tmask(ji,jj,jk)
                   !
                END DO
             END DO
          END DO
-         pbetaPE (:,:,:) = vallis_ratio * 0.78e-3_wp * tmask(:,:,:)   ! betaPE=beta
+         !
+      CASE( 1 )                !==  Vallis (2006) simplified EOS  ==!
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature (t-T0)
+                  zh = fsdept(ji,jj,jk)                                                ! depth in meters  at t-point
+                  pab_pe(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*(rn_cabbel*zt + rn_thermo*zh) )  &
+                      &          * tmask(ji,jj,jk)                                     ! alphaPE
+                  pen(ji,jj,jk) = ( - rn_alph0 * ( 1._wp + rn_cabbel*zt + 0.5_wp*rn_thermo*zh ) *zt   &                                          &
+                      &          + rn_beta0 * zs + 0.5_wp*rn_gamm0 * zh  ) * tmask(ji,jj,jk)
+                  !
+               END DO
+            END DO
+         END DO
+         pab_pe (:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:)   ! betaPE
          !
       CASE DEFAULT
@@ -1598 +1281 @@
 
 
+   SUBROUTINE eos_bn2_ab( pts, pn2 )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE eos_bn2_ab  ***
+      !!
+      !! ** Purpose :   Update alpha/beta then compute the local Brunt-Vaisala 
+      !!      frequency at the time-step of the input arguments
+      !!
+      !! ** Method  :   call eos_ab to obtain alpha and beta coefficients
+      !!                then interpolate them vertically, and compute pn2
+      !!        Macro-tasked on horizontal slab (jk-loop)
+      !!      N.B. N^2 is set to zero at the first level (JK=1) in inidtr
+      !!      and is never used at this level.
+      !!
+      !! ** Action  : - pn2 : the brunt-vaisala frequency
+      !!
+      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
+      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995
+      !!                McDougall, J. Phys. Oceano., 1987
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
+      !                                                               ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pn2   ! Brunt-Vaisala frequency    [s-1]
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zgde3w, za, zb, zr1          ! temporary scalars
+#if defined key_zdfddm
+      REAL(wp) ::   zds   ! local scalars
+#endif
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zab
+      !!----------------------------------------------------------------------
+      !
+      !
+      IF( nn_timing == 1 ) CALL timing_start('eos_bn2_ab')
+      !
+      ! pn2 : interior points only (2=< jk =< jpkm1 )
+      ! --------------------------
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, jpts, zab )
+      !
+      CALL eos_ab( pts, zab )
+      !
+      DO jk = 2, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               zgde3w = 1._wp / fse3w(ji,jj,jk)
+               zr1 = ( fsdepw(ji,jj,jk) - fsdept(ji,jj,jk) ) * zgde3w
+               !
+               za = ( zab(ji,jj,jk,jp_tem) +   &
+                  &   ( zab(ji,jj,jk-1,jp_tem) - zab(ji,jj,jk,jp_tem) ) * zr1 )   &
+                  &   * tmask(ji,jj,jk)
+               zb = ( zab(ji,jj,jk,jp_sal) +   &
+                  &   ( zab(ji,jj,jk-1,jp_sal) - zab(ji,jj,jk,jp_sal) ) * zr1 )   &
+                  &   * tmask(ji,jj,jk)
+               ! N2
+               pn2(ji,jj,jk) = grav * zgde3w                    &   ! N^2
+                  &          * (   za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   &
+                  &              - zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) )
+#if defined key_zdfddm
+               zds = zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )      ! Rrau = (alpha / beta) (dk[t] / dk[s])
+               IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp
+               rrau(ji,jj,jk) = za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
+#endif
+            END DO
+         END DO
+      END DO
+      !
+
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk )
+#if defined key_zdfddm
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk )
+#endif
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, jpts, zab )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('eos_bn2_ab')
+      !
+   END SUBROUTINE eos_bn2_ab
+
+
+   SUBROUTINE eos_bn2( pts, pab, pn2 )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE eos_bn2  ***
+      !!
+      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the 
+      !!                time-step of the input arguments
+      !!
+      !! ** Method  :   interpolate alpha/beta vertically, and compute pn2
+      !!        Macro-tasked on horizontal slab (jk-loop)
+      !!      N.B. N^2 is set to zero at the first level (JK=1) in inidtr
+      !!      and is never used at this level.
+      !!
+      !! ** Action  : - pn2 : the brunt-vaisala frequency
+      !!
+      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
+      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995
+      !!                McDougall, J. Phys. Oceano., 1987
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celcius,psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! partial derivative of in situ 
+      !                                                              ! density with respect to T and S [1/Celcius,1/psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency    [1/s^2]
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zgde3w, za, zb, zr1          ! temporary scalars
+#if defined key_zdfddm
+      REAL(wp) ::   zds   ! local scalars
+#endif
+      !!----------------------------------------------------------------------
+      !
+      !
+      IF( nn_timing == 1 ) CALL timing_start('eos_bn2')
+      !
+      ! pn2 : interior points only (2=< jk =< jpkm1 )
+      ! --------------------------
+      !
+      DO jk = 2, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zgde3w = 1._wp / fse3w(ji,jj,jk)
+               zr1 = ( fsdepw(ji,jj,jk) - fsdept(ji,jj,jk) ) * zgde3w
+               !
+               za = ( pab(ji,jj,jk,jp_tem) + ( pab(ji,jj,jk-1,jp_tem) - pab(ji,jj,jk,jp_tem) ) * zr1 ) * tmask(ji,jj,jk)
+               zb = ( pab(ji,jj,jk,jp_sal) + ( pab(ji,jj,jk-1,jp_sal) - pab(ji,jj,jk,jp_sal) ) * zr1 ) * tmask(ji,jj,jk)
+               !
+               pn2(ji,jj,jk) = grav * zgde3w                    &   ! N^2
+                  &          * (   za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   &
+                  &              - zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) )
+#if defined key_zdfddm
+               zds = zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )      ! Rrau = (alpha / beta) (dk[t] / dk[s])
+               IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp
+               rrau(ji,jj,jk) = za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
+#endif
+            END DO
+         END DO
+      END DO
+      !
+
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk )
+#if defined key_zdfddm
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk )
+#endif
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('eos_bn2')
+      !
+   END SUBROUTINE eos_bn2
+
+
+
    FUNCTION tfreez( psal ) RESULT( ptf )
       !!----------------------------------------------------------------------
@@ -1629 +1461 @@
       !! ** Method  :   Read the namelist nameos and control the parameters
       !!----------------------------------------------------------------------
-      NAMELIST/nameos/ nn_eos, rn_alpha, rn_beta
+      NAMELIST/nameos/ nn_eos, rn_alph0, rn_beta0, rn_gamm0, rn_cabbel, rn_thermo
       !!----------------------------------------------------------------------
       !
@@ -1640 +1472 @@
          WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '          Namelist nameos : set eos parameters'
-         WRITE(numout,*) '             flag for eq. of state and N^2  nn_eos   = ', nn_eos
-         WRITE(numout,*) '             thermal exp. coef. (linear)    rn_alpha = ', rn_alpha
-         WRITE(numout,*) '             saline  exp. coef. (linear)    rn_beta  = ', rn_beta
+         WRITE(numout,*) '             flag for eq. of state and N^2    nn_eos    = ', nn_eos
+         WRITE(numout,*) '             thermal exp. coef. (nn_eos=1)    rn_alph0  = ', rn_alph0
+         WRITE(numout,*) '             saline  exp. coef. (nn_eos=1)    rn_beta0  = ', rn_beta0
+         WRITE(numout,*) '             thermal exp. coef. (nn_eos=1)    rn_gamm0  = ', rn_gamm0
+         WRITE(numout,*) '             saline  exp. coef. (nn_eos=1)    rn_cabbel = ', rn_cabbel
+         WRITE(numout,*) '             thermal exp. coef. (nn_eos=1)    rn_thermo = ', rn_thermo
       ENDIF
       !
@@ -1656 +1491 @@
          IF(lwp) WRITE(numout,*) '                and McDougall (1987) Brunt-Vaisala frequency'
          !
-      CASE( 1 )                        !==  Linear formulation = F( temperature )  ==!
+      CASE( 1 )                        !==  Vallis (2006) formulation  ==!
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '          use of linear eos rho(T) = rau0 * ( 1.0285 - rn_alpha * T )'
-         IF( lk_zdfddm ) CALL ctl_stop( '          double diffusive mixing parameterization requires',   &
-              &                         ' that T and S are used as state variables' )
-         !
-      CASE( 2 )                        !==  Linear formulation = F( temperature , salinity )  ==!
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '          use of linear eos rho(T,S) = rau0 * ( rn_beta * S - rn_alpha * T )'
-         !
-      CASE( 3 )                        !==  Vallis (2006) formulation  ==!
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '          use of Vallis (2006) equation of state'
-         vallis_ratio = 1027._wp / rau0
-         vallis_diff = (1027._wp-rau0) / rau0
+         IF(lwp) WRITE(numout,*) '          use of simplified eos:    rhd(T,S,z) = '
+         IF(lwp) WRITE(numout,*) '             -alph0*(1+cabbel*(T-T0)+thermo*z)*(T-T0) + beta0*(S-S0) + gamm0*z'
+         IF( lk_zdfddm .AND. rn_beta0==0. )   &
+              &   CALL ctl_stop( '          double diffusive mixing parameterization requires that rn_beta<>0')
+         offset = 1027._wp / rau0 - 1._wp
          !
       CASE DEFAULT                     !==  ERROR in nn_eos  ==!

branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90

@@ -79 +79 @@
       !
 
+!!gm end
       IF( ln_tra_mxl .OR. ln_vor_trd )   CALL ctl_stop( 'ML tracer and Barotropic vorticity diags are still using old IOIPSL' )
 !!gm end

branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90

@@ -276 +276 @@
       ENDIF
       !
-      nkstp     = nit000 - 1
-!!gm check the value of nkstp  here nit000 should be OK
-      !
    END SUBROUTINE trd_ken_init
 

branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdpen.F90

@@ -33 +33 @@
    INTEGER ::   nkstp   ! current time step 
 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   alphaPE, betaPE   ! partial derivative of PEanom with respect to T and S
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   ab_pe   ! partial derivative of PE anomaly with respect to T and S
 
    !! * Substitutions
@@ -50 +50 @@
       !!                  ***  FUNCTION trd_tra_alloc  ***
       !!---------------------------------------------------------------------
-      ALLOCATE( alphaPE(jpi,jpj,jpk) , betaPE(jpi,jpj,jpk) , STAT= trd_pen_alloc )
+      ALLOCATE( ab_pe(jpi,jpj,jpk,jpts) , STAT= trd_pen_alloc )
       !
       IF( lk_mpp             )   CALL mpp_sum ( trd_pen_alloc )
-      IF( trd_pen_alloc /= 0 )   CALL ctl_warn('trd_pen_alloc: failed to allocate arrays')
+      IF( trd_pen_alloc /= 0 )   CALL ctl_warn( 'trd_pen_alloc: failed to allocate arrays' )
    END FUNCTION trd_pen_alloc
 
@@ -79 +79 @@
       IF ( kt /= nkstp ) THEN   ! full eos: set partial derivatives at the 1st call of kt time step
          nkstp = kt
-         CALL eos_pen( tsn, alphaPE, betaPE, zpe )
-         CALL iom_put( "alphaPE", alphaPE )
-         CALL iom_put( "betaPE", betaPE )
-         CALL iom_put( "PEanom", zpe )
+         CALL eos_pen( tsn, ab_PE, zpe )
+         CALL iom_put( "alphaPE", ab_pe(:,:,:,jp_tem) )
+         CALL iom_put( "betaPE" , ab_pe(:,:,:,jp_sal) )
+         CALL iom_put( "PEanom" , zpe )
       ENDIF
       !
       DO jk = 1, jpkm1
-         zpe(:,:,jk) = ( - alphaPE(:,:,jk) * ptrdx(:,:,jk)   &
-            &            + betaPE (:,:,jk) * ptrdy(:,:,jk)  )
+         zpe(:,:,jk) = ( - ab_pe(:,:,jk,jp_tem) * ptrdx(:,:,jk)   &
+            &            + ab_pe(:,:,jk,jp_sal) * ptrdy(:,:,jk)  )
       END DO
 
@@ -96 +96 @@
                                 IF( .NOT.lk_vvl ) THEN                   ! cst volume : adv flux through z=0 surface
                                    CALL wrk_alloc( jpi, jpj, z2d )
-                                   z2d(:,:) = wn(:,:,1) * ( - alphaPE(:,:,1) * tsn(:,:,1,jp_tem)    &
-                                      &                     + betaPE (:,:,1) * tsn(:,:,1,jp_sal)  ) / fse3t(:,:,1)
+                                   z2d(:,:) = wn(:,:,1) * ( - ab_pe(:,:,1,jp_tem) * tsn(:,:,1,jp_tem)    &
+                                      &                     + ab_pe(:,:,1,jp_sal) * tsn(:,:,1,jp_sal)  ) / fse3t(:,:,1)
                                    CALL iom_put( "petrd_sad" , z2d )
                                    CALL wrk_dealloc( jpi, jpj, z2d )

branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/oce.F90

@@ -20 +20 @@
    !! dynamics and tracer fields                            ! before ! now    ! after  ! the after trends becomes the fields
    !! --------------------------                            ! fields ! fields ! trends ! only after tra_zdf and dyn_spg
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   ub   ,  un    , ua     !: i-horizontal velocity        [m/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity        [m/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wn             !: vertical velocity            [m/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rotb ,  rotn           !: relative vorticity           [s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hdivb,  hdivn          !: horizontal divergence        [s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn            !: 4D T-S fields        [Celcius,psu] 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rn2b ,  rn2            !: brunt-vaisala frequency**2   [s-2]
-   !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  tsa             !: 4D T-S trends fields & work array 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   ub   ,  un    , ua     !: i-horizontal velocity          [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity          [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wn             !: vertical velocity              [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rotb ,  rotn           !: relative vorticity             [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hdivb,  hdivn          !: horizontal divergence          [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn   , tsa    !: 4D T-S fields                  [Celcius,psu] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   rab_b,  rab_n          !: thermal/haline expansion coef. [Celcius-1,psu-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rn2b ,  rn2            !: brunt-vaisala frequency**2     [s-2]
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0  [no units]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhop   !: potential volumic mass                           [kg/m3]
 
-   !! free surface                                      !  before  ! now    ! after  !
-   !! ------------                                      !  fields  ! fields ! trends !
+   !! free surface                                              !  before  ! now    ! after  !
+   !! ------------                                              !  fields  ! fields ! trends !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshb   , sshn   , ssha   !: sea surface height at t-point [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   sshu_b , sshu_n , sshu_a !: sea surface height at u-point [m]
@@ -75 +74 @@
          &      hdivb(jpi,jpj,jpk)      , hdivn(jpi,jpj,jpk)      ,                             &
          &      tsb  (jpi,jpj,jpk,jpts) , tsn  (jpi,jpj,jpk,jpts) , tsa(jpi,jpj,jpk,jpts) ,     &
+         &      rab_b(jpi,jpj,jpk,jpts) , rab_n(jpi,jpj,jpk,jpts) ,                             &
          &      rn2b (jpi,jpj,jpk)      , rn2  (jpi,jpj,jpk)                              , STAT=ierr(1) )
          !

branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -104 +104 @@
       ! Ocean physics update                (ua, va, tsa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL bn2( tsb, rn2b )        ! before Brunt-Vaisala frequency
-                         CALL bn2( tsn, rn2  )        ! now    Brunt-Vaisala frequency
+                         CALL eos( tsb, rab_b )              ! before thermal/haline expansion coef. 
+                         CALL bn2( tsb, rab_b, rn2b )        ! before Brunt-Vaisala frequency
+                         CALL eos( tsn, rab_n )              ! now thermal/haline expansion coef. 
+                         CALL bn2( tsn, rab_n, rn2  )        ! now    Brunt-Vaisala frequency
+                         
       !
       !  VERTICAL PHYSICS