Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/C1D/step_c1d.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/C1D/step_c1d.F90	(revision 15065)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/C1D/step_c1d.F90	(revision 15066)
@@ -72,6 +72,6 @@
                          CALL eos_rab( CASTWP(ts(:,:,:,:,Nbb)), rab_b, Nnn )  ! before local thermal/haline expension ratio at T-points
                          CALL eos_rab( CASTWP(ts(:,:,:,:,Nnn)), rab_n, Nnn )  ! now    local thermal/haline expension ratio at T-points
-                         CALL bn2( CASTWP(ts(:,:,:,:,Nbb)), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
-                         CALL bn2( CASTWP(ts(:,:,:,:,Nnn)), rab_n, rn2 , Nnn ) ! now    Brunt-Vaisala frequency
+                         CALL bn2( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
+                         CALL bn2( ts(:,:,:,:,Nnn), rab_n, rn2 , Nnn ) ! now    Brunt-Vaisala frequency
 
       !  VERTICAL PHYSICS
@@ -108,5 +108,5 @@
       IF( ln_zdfosm  )  CALL tra_osm( kstp, Nnn     , ts, Nrhs  )  ! OSMOSIS non-local tracer fluxes
                         CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa   )         ! vertical mixing
-                        CALL eos( CASTEWP(ts(:,:,:,:,Nnn)), rhd, rhop, gdept_0(:,:,:) )  ! now potential density for zdfmxl
+                        CALL eos( CASTWP(ts(:,:,:,:,Nnn)), rhd, rhop, gdept_0(:,:,:) )  ! now potential density for zdfmxl
       IF( ln_zdfnpc )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa   )         ! applied non penetrative convective adjustment on (t,s)
                         CALL tra_atf( kstp, Nbb, Nnn, Naa, ts )                 ! time filtering of "now" tracer arrays
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DOM/domqco.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DOM/domqco.F90	(revision 15065)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DOM/domqco.F90	(revision 15066)
@@ -118,9 +118,9 @@
       !                                ! Horizontal interpolation of e3t
 #if defined key_RK3
-      CALL dom_qco_r3c( CASTWP(ssh(:,:,Kbb)), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb), r3f(:,:) )
-      CALL dom_qco_r3c( CASTWP(ssh(:,:,Kmm)), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm)           )
+      CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb), r3f(:,:) )
+      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm)           )
 #else
-      CALL dom_qco_r3c( CASTWP(ssh(:,:,Kbb)), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb)           )
-      CALL dom_qco_r3c( CASTWP(ssh(:,:,Kmm)), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )
+      CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb)           )
+      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )
 #endif
       ! dom_qco_r3c defines over [nn_hls, nn_hls-1, nn_hls, nn_hls-1]
@@ -142,5 +142,5 @@
       !!              - compute the ratio ssh/h_0 at t-,u-,v-pts, (f-pt optional)
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:)          , INTENT(in   )  ::   pssh               ! sea surface height   [m]
+      REAL(dp), DIMENSION(:,:)          , INTENT(in   )  ::   pssh               ! sea surface height   [m]
       REAL(wp), DIMENSION(:,:)          , INTENT(  out)  ::   pr3t, pr3u, pr3v   ! ssh/h0 ratio at t-, u-, v-,points  [-]
       REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT(  out)  ::   pr3f               ! ssh/h0 ratio at f-point   [-]
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DYN/dynspg_ts.F90	(revision 15065)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DYN/dynspg_ts.F90	(revision 15066)
@@ -87,8 +87,4 @@
 
 
-   INTERFACE dyn_drg_init
-      MODULE PROCEDURE dyn_drg_init_wp, dyn_drg_init_mixed
-   END INTERFACE
-
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -288,5 +284,5 @@
          END_2D
       ELSE				     !* remove baroclinic drag AND provide the barotropic drag coefficients
-         CALL dyn_drg_init( Kbb, Kmm, CASTWP(puu), CASTWP(pvv), puu_b, pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v )
+         CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b, pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v )
       ENDIF
       !
@@ -1347,5 +1343,5 @@
      
 
-   SUBROUTINE dyn_drg_init_wp( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
+   SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_drg_init  ***
@@ -1358,5 +1354,5 @@
       !!----------------------------------------------------------------------
       INTEGER                             , INTENT(in   ) ::  Kbb, Kmm           ! ocean time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) ::  puu, pvv           ! ocean velocities and RHS of momentum equation
+      REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) ::  puu, pvv           ! ocean velocities and RHS of momentum equation
       REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(in   ) ::  puu_b, pvv_b       ! barotropic velocities at main time levels
       REAL(wp), DIMENSION(jpi,jpj)        , INTENT(inout) ::  pu_RHSi, pv_RHSi   ! baroclinic part of the barotropic RHS
@@ -1451,111 +1447,5 @@
       ENDIF
       !
-   END SUBROUTINE dyn_drg_init_wp
-
-   SUBROUTINE dyn_drg_init_mixed( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dyn_drg_init  ***
-      !!                    
-      !! ** Purpose : - add the baroclinic top/bottom drag contribution to 
-      !!              the baroclinic part of the barotropic RHS
-      !!              - compute the barotropic drag coefficients
-      !!
-      !! ** Method  :   computation done over the INNER domain only 
-      !!----------------------------------------------------------------------
-      INTEGER                             , INTENT(in   ) ::  Kbb, Kmm           ! ocean time level indices
-      REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) ::  puu, pvv           ! ocean velocities and RHS of momentum equation
-      REAL(sp), DIMENSION(jpi,jpj,jpt)    , INTENT(in   ) ::  puu_b, pvv_b       ! barotropic velocities at main time levels
-      REAL(sp), DIMENSION(jpi,jpj)        , INTENT(inout) ::  pu_RHSi, pv_RHSi   ! baroclinic part of the barotropic RHS
-      REAL(sp), DIMENSION(jpi,jpj)        , INTENT(  out) ::  pCdU_u , pCdU_v    ! barotropic drag coefficients
-      !
-      INTEGER  ::   ji, jj   ! dummy loop indices
-      INTEGER  ::   ikbu, ikbv, iktu, iktv
-      REAL(sp) ::   zztmp
-      REAL(sp), DIMENSION(jpi,jpj) ::   zu_i, zv_i
-      !!----------------------------------------------------------------------
-      !
-      !                    !==  Set the barotropic drag coef.  ==!
-      !
-      IF( ln_isfcav.OR.ln_drgice_imp ) THEN          ! top+bottom friction (ocean cavities)
-         
-         DO_2D( 0, 0, 0, 0 )
-            pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
-            pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
-         END_2D
-      ELSE                          ! bottom friction only
-         DO_2D( 0, 0, 0, 0 )
-            pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
-            pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
-         END_2D
-      ENDIF
-      !
-      !                    !==  BOTTOM stress contribution from baroclinic velocities  ==!
-      !
-      IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
-         
-         DO_2D( 0, 0, 0, 0 )
-            ikbu = mbku(ji,jj)       
-            ikbv = mbkv(ji,jj)    
-            zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
-            zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
-         END_2D
-      ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
-         
-         DO_2D( 0, 0, 0, 0 )
-            ikbu = mbku(ji,jj)       
-            ikbv = mbkv(ji,jj)    
-            zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
-            zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
-         END_2D
-      ENDIF
-      !
-      IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
-         zztmp = -1._wp / rDt_e
-         DO_2D( 0, 0, 0, 0 )
-            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
-                 &                              r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
-            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 & 
-                 &                              r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp  )
-         END_2D
-      ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
-         
-         DO_2D( 0, 0, 0, 0 )
-            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj)
-            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj)
-         END_2D
-      END IF
-      !
-      !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case)
-      !
-      IF( ln_isfcav.OR.ln_drgice_imp ) THEN
-         !
-         IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
-            
-            DO_2D( 0, 0, 0, 0 )
-               iktu = miku(ji,jj)
-               iktv = mikv(ji,jj)
-               zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
-               zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
-            END_2D
-         ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
-            
-            DO_2D( 0, 0, 0, 0 )
-               iktu = miku(ji,jj)
-               iktv = mikv(ji,jj)
-               zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
-               zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
-            END_2D
-         ENDIF
-         !
-         !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
-         
-         DO_2D( 0, 0, 0, 0 )
-            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj)
-            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj)
-         END_2D
-         !
-      ENDIF
-      !
-   END SUBROUTINE dyn_drg_init_mixed
+   END SUBROUTINE dyn_drg_init
 
    SUBROUTINE ts_bck_interp( jn, ll_init,       &   ! <== in
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/eosbn2.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/eosbn2.F90	(revision 15065)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/eosbn2.F90	(revision 15066)
@@ -56,9 +56,11 @@
    !                  !! * Interface
    INTERFACE eos
-      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d
+      MODULE PROCEDURE eos_insitu_wp, eos_insitu_pot_wp, eos_insitu_2d, eos_insitu_pot_2d
+      MODULE PROCEDURE eos_insitu_mixed, eos_insitu_pot_mixed
    END INTERFACE
    !
    INTERFACE eos_rab
-      MODULE PROCEDURE rab_3d, rab_2d, rab_0d
+      MODULE PROCEDURE rab_3d_wp, rab_2d, rab_0d
+      MODULE PROCEDURE rab_3d_mixed
    END INTERFACE
    !
@@ -190,5 +192,5 @@
 CONTAINS
 
-   SUBROUTINE eos_insitu( pts, prd, pdep )
+   SUBROUTINE eos_insitu_wp( pts, prd, pdep )
       !!
       REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
@@ -197,10 +199,20 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pdep  ! depth                      [m]
       !!
-      CALL eos_insitu_t( pts, is_tile(pts), prd, is_tile(prd), pdep, is_tile(pdep) )
-   END SUBROUTINE eos_insitu
-
-   SUBROUTINE eos_insitu_t( pts, ktts, prd, ktrd, pdep, ktdep )
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE eos_insitu  ***
+      CALL eos_insitu_t_wp( pts, is_tile(pts), prd, is_tile(prd), pdep, is_tile(pdep) )
+   END SUBROUTINE eos_insitu_wp
+
+   SUBROUTINE eos_insitu_mixed( pts, prd, pdep )
+      !!
+      REAL(dp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
+      !                                                      ! 2 : salinity               [psu]
+      REAL(sp), DIMENSION(:,:,:)  , INTENT(  out) ::   prd   ! in situ density            [-]
+      REAL(sp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pdep  ! depth                      [m]
+      !!
+      CALL eos_insitu_t_mixed( pts, is_tile(pts), prd, is_tile(prd), pdep, is_tile(pdep) )
+   END SUBROUTINE 
+
+   SUBROUTINE eos_insitu_t_wp( pts, ktts, prd, ktrd, pdep, ktdep )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE 
       !!
       !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from
@@ -306,8 +318,115 @@
       IF( ln_timing )   CALL timing_stop('eos-insitu')
       !
-   END SUBROUTINE eos_insitu_t
-
-
-   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
+   END SUBROUTINE eos_insitu_t_wp
+   SUBROUTINE eos_insitu_t_mixed( pts, ktts, prd, ktrd, pdep, ktdep )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE eos_insitu  ***
+      !!
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from
+      !!       potential temperature and salinity using an equation of state
+      !!       selected in the nameos namelist
+      !!
+      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0
+      !!         with   prd    in situ density anomaly      no units
+      !!                t      TEOS10: CT or EOS80: PT      Celsius
+      !!                s      TEOS10: SA or EOS80: SP      TEOS10: g/kg or EOS80: psu
+      !!                z      depth                        meters
+      !!                rho    in situ density              kg/m^3
+      !!                rho0   reference density            kg/m^3
+      !!
+      !!     ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
+      !!         Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg
+      !!
+      !!     ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z).
+      !!         Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu
+      !!
+      !!     ln_seos : simplified equation of state
+      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0
+      !!              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 like equation: use default values of coefficients
+      !!
+      !! ** Action  :   compute prd , the in situ density (no units)
+      !!
+      !! References :   Roquet et al, Ocean Modelling, in preparation (2014)
+      !!                Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
+      !!                TEOS-10 Manual, 2010
+      !!----------------------------------------------------------------------
+      INTEGER                                 , INTENT(in   ) ::   ktts, ktrd, ktdep
+      REAL(dp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
+      !                                                                  ! 2 : salinity               [psu]
+      REAL(sp), DIMENSION(A2D_T(ktrd) ,JPK     ), INTENT(  out) ::   prd   ! in situ density            [-]
+      REAL(sp), DIMENSION(A2D_T(ktdep),JPK     ), INTENT(in   ) ::   pdep  ! depth                      [m]
+      !
+      INTEGER  ::   ji, jj, jk                ! dummy loop indices
+      REAL(sp) ::   zt , zh , zs , ztm        ! local scalars
+      REAL(sp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('eos-insitu')
+      !
+      SELECT CASE( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
+         !
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+            !
+            zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
+            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            ztm = tmask(ji,jj,jk)                                         ! tmask
+            !
+            zn3 = EOS013*zt   &
+               &   + EOS103*zs+EOS003
+               !
+            zn2 = (EOS022*zt   &
+               &   + EOS112*zs+EOS012)*zt   &
+               &   + (EOS202*zs+EOS102)*zs+EOS002
+               !
+            zn1 = (((EOS041*zt   &
+               &   + EOS131*zs+EOS031)*zt   &
+               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
+               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
+               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+               !
+            zn0 = (((((EOS060*zt   &
+               &   + EOS150*zs+EOS050)*zt   &
+               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked)
+            !
+         END_3D
+         !
+      CASE( np_seos )                !==  simplified EOS  ==!
+         !
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
+            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
+            zh  = pdep (ji,jj,jk)
+            ztm = tmask(ji,jj,jk)
+            !
+            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
+               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
+               &  - rn_nu * zt * zs
+               !
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
+         END_3D
+         !
+      END SELECT
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=REAL(prd, wp), clinfo1=' eos-insitu  : ', kdim=jpk )
+      !
+      IF( ln_timing )   CALL timing_stop('eos-insitu')
+      !
+   END SUBROUTINE eos_insitu_t_mixed
+
+
+   SUBROUTINE eos_insitu_pot_wp( pts, prd, prhop, pdep )
       !!
       REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
@@ -317,9 +436,20 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pdep   ! depth                      [m]
       !!
-      CALL eos_insitu_pot_t( pts, is_tile(pts), prd, is_tile(prd), prhop, is_tile(prhop), pdep, is_tile(pdep) )
-   END SUBROUTINE eos_insitu_pot
-
-
-   SUBROUTINE eos_insitu_pot_t( pts, ktts, prd, ktrd, prhop, ktrhop, pdep, ktdep )
+      CALL eos_insitu_pot_t_wp( pts, is_tile(pts), prd, is_tile(prd), prhop, is_tile(prhop), pdep, is_tile(pdep) )
+   END SUBROUTINE eos_insitu_pot_wp
+
+   SUBROUTINE eos_insitu_pot_mixed( pts, prd, prhop, pdep )
+      !!
+      REAL(dp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
+      !                                                       ! 2 : salinity               [psu]
+      REAL(sp), DIMENSION(:,:,:)  , INTENT(  out) ::   prd    ! in situ density            [-]
+      REAL(dp), DIMENSION(:,:,:)  , INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      REAL(sp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pdep   ! depth                      [m]
+      !!
+      CALL eos_insitu_pot_t_mixed( pts, is_tile(pts), prd, is_tile(prd), prhop, is_tile(prhop), pdep, is_tile(pdep) )
+   END SUBROUTINE eos_insitu_pot_mixed
+
+
+   SUBROUTINE eos_insitu_pot_t_wp( pts, ktts, prd, ktrd, prhop, ktrhop, pdep, ktdep )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_insitu_pot  ***
@@ -475,5 +605,161 @@
       IF( ln_timing )   CALL timing_stop('eos-pot')
       !
-   END SUBROUTINE eos_insitu_pot_t
+   END SUBROUTINE eos_insitu_pot_t_wp
+
+   SUBROUTINE eos_insitu_pot_t_mixed( pts, ktts, prd, ktrd, prhop, ktrhop, pdep, ktdep )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE eos_insitu_pot  ***
+      !!
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the
+      !!      potential volumic mass (Kg/m3) from potential temperature and
+      !!      salinity fields using an equation of state selected in the
+      !!     namelist.
+      !!
+      !! ** Action  : - prd  , the in situ density (no units)
+      !!              - prhop, the potential volumic mass (Kg/m3)
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) ::   ktts, ktrd, ktrhop, ktdep
+      REAL(dp), DIMENSION(A2D_T(ktts)  ,JPK,JPTS), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
+      !                                                                    ! 2 : salinity               [psu]
+      REAL(sp), DIMENSION(A2D_T(ktrd)  ,JPK     ), INTENT(  out) ::   prd    ! in situ density            [-]
+      REAL(dp), DIMENSION(A2D_T(ktrhop),JPK     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      REAL(sp), DIMENSION(A2D_T(ktdep) ,JPK     ), INTENT(in   ) ::   pdep   ! depth                      [m]
+      !
+      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices
+      INTEGER  ::   jdof
+      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars
+      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      -
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('eos-pot')
+      !
+      SELECT CASE ( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
+         !
+         ! Stochastic equation of state
+         IF ( ln_sto_eos ) THEN
+            ALLOCATE(zn0_sto(1:2*nn_sto_eos))
+            ALLOCATE(zn_sto(1:2*nn_sto_eos))
+            ALLOCATE(zsign(1:2*nn_sto_eos))
+            DO jsmp = 1, 2*nn_sto_eos, 2
+              zsign(jsmp)   = 1._wp
+              zsign(jsmp+1) = -1._wp
+            END DO
+            !
+            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+               !
+               ! compute density (2*nn_sto_eos) times:
+               ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)
+               ! (2) for t-dt, s-ds (with the opposite fluctuation)
+               DO jsmp = 1, nn_sto_eos*2
+                  jdof   = (jsmp + 1) / 2
+                  zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth
+                  zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature
+                  zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp)
+                  zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity
+                  ztm    = tmask(ji,jj,jk)                                         ! tmask
+                  !
+                  zn3 = EOS013*zt   &
+                     &   + EOS103*zs+EOS003
+                     !
+                  zn2 = (EOS022*zt   &
+                     &   + EOS112*zs+EOS012)*zt   &
+                     &   + (EOS202*zs+EOS102)*zs+EOS002
+                     !
+                  zn1 = (((EOS041*zt   &
+                     &   + EOS131*zs+EOS031)*zt   &
+                     &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
+                     &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
+                     &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+                     !
+                  zn0_sto(jsmp) = (((((EOS060*zt   &
+                     &   + EOS150*zs+EOS050)*zt   &
+                     &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+                     &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+                     &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+                     &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+                     &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+                     !
+                  zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp)
+               END DO
+               !
+               ! compute stochastic density as the mean of the (2*nn_sto_eos) densities
+               prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp
+               DO jsmp = 1, nn_sto_eos*2
+                  prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface
+                  !
+                  prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked)
+               END DO
+               prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
+               prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos
+            END_3D
+            DEALLOCATE(zn0_sto,zn_sto,zsign)
+         ! Non-stochastic equation of state
+         ELSE
+            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+               !
+               zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
+               zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
+               zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+               ztm = tmask(ji,jj,jk)                                         ! tmask
+               !
+               zn3 = EOS013*zt   &
+                  &   + EOS103*zs+EOS003
+                  !
+               zn2 = (EOS022*zt   &
+                  &   + EOS112*zs+EOS012)*zt   &
+                  &   + (EOS202*zs+EOS102)*zs+EOS002
+                  !
+               zn1 = (((EOS041*zt   &
+                  &   + EOS131*zs+EOS031)*zt   &
+                  &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
+                  &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
+                  &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+                  !
+               zn0 = (((((EOS060*zt   &
+                  &   + EOS150*zs+EOS050)*zt   &
+                  &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+                  &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+                  &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+                  &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+                  &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+                  !
+               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+               !
+               prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
+               !
+               prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked)
+            END_3D
+         ENDIF
+
+      CASE( np_seos )                !==  simplified EOS  ==!
+         !
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
+            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
+            zh  = pdep (ji,jj,jk)
+            ztm = tmask(ji,jj,jk)
+            !                                                     ! potential density referenced at the surface
+            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   &
+               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   &
+               &  - rn_nu * zt * zs
+            prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
+            !                                                     ! density anomaly (masked)
+            zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm
+            !
+         END_3D
+         !
+      END SELECT
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=REAL(prd, wp), clinfo1=' eos-pot: ', &
+         &                                  tab3d_2=REAL(prhop, wp), clinfo2=' pot : ', kdim=jpk )
+      !
+      IF( ln_timing )   CALL timing_stop('eos-pot')
+      !
+   END SUBROUTINE eos_insitu_pot_t_mixed
 
 
@@ -661,5 +947,5 @@
 
 
-   SUBROUTINE rab_3d( pts, pab, Kmm )
+   SUBROUTINE rab_3d_wp( pts, pab, Kmm )
       !!
       INTEGER                     , INTENT(in   ) ::   Kmm   ! time level index
@@ -667,9 +953,18 @@
       REAL(wp), DIMENSION(:,:,:,:), INTENT(  out) ::   pab   ! thermal/haline expansion ratio
       !!
-      CALL rab_3d_t( pts, is_tile(pts), pab, is_tile(pab), Kmm )
-   END SUBROUTINE rab_3d
-
-
-   SUBROUTINE rab_3d_t( pts, ktts, pab, ktab, Kmm )
+      CALL rab_3d_t_wp( pts, is_tile(pts), pab, is_tile(pab), Kmm )
+   END SUBROUTINE rab_3d_wp
+
+   SUBROUTINE rab_3d_mixed( pts, pab, Kmm )
+      !!
+      INTEGER                     , INTENT(in   ) ::   Kmm   ! time level index
+      REAL(dp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! pot. temperature & salinity
+      REAL(sp), DIMENSION(:,:,:,:), INTENT(  out) ::   pab   ! thermal/haline expansion ratio
+      !!
+      CALL rab_3d_t_mixed( pts, is_tile(pts), pab, is_tile(pab), Kmm )
+   END SUBROUTINE rab_3d_mixed
+
+
+   SUBROUTINE rab_3d_t_wp( pts, ktts, pab, ktab, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE rab_3d  ***
@@ -775,40 +1070,27 @@
       IF( ln_timing )   CALL timing_stop('rab_3d')
       !
-   END SUBROUTINE rab_3d_t
-
-
-   SUBROUTINE rab_2d( pts, pdep, pab, Kmm )
-      !!
-      INTEGER                   , INTENT(in   ) ::   Kmm   ! time level index
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pts    ! pot. temperature & salinity
-      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pdep   ! depth                  [m]
-      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pab    ! thermal/haline expansion ratio
-      !!
-      CALL rab_2d_t(pts, is_tile(pts), pdep, is_tile(pdep), pab, is_tile(pab), Kmm)
-   END SUBROUTINE rab_2d
-
-
-   SUBROUTINE rab_2d_t( pts, ktts, pdep, ktdep, pab, ktab, Kmm )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE rab_2d  ***
-      !!
-      !! ** Purpose :   Calculates thermal/haline expansion ratio for a 2d field (unmasked)
+   END SUBROUTINE rab_3d_t_wp
+
+   SUBROUTINE rab_3d_t_mixed( pts, ktts, pab, ktab, Kmm )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE rab_3d  ***
+      !!
+      !! ** Purpose :   Calculates thermal/haline expansion ratio at T-points
+      !!
+      !! ** Method  :   calculates alpha / beta at T-points
       !!
       !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
       !!----------------------------------------------------------------------
-      INTEGER                            , INTENT(in   ) ::   Kmm   ! time level index
-      INTEGER                            , INTENT(in   ) ::   ktts, ktdep, ktab
-      REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in   ) ::   pts    ! pot. temperature & salinity
-      REAL(wp), DIMENSION(A2D_T(ktdep)    ), INTENT(in   ) ::   pdep   ! depth                  [m]
-      REAL(wp), DIMENSION(A2D_T(ktab),JPTS), INTENT(  out) ::   pab    ! thermal/haline expansion ratio
+      INTEGER                                , INTENT(in   ) ::   Kmm   ! time level index
+      INTEGER                                , INTENT(in   ) ::   ktts, ktab
+      REAL(dp), DIMENSION(A2D_T(ktts),JPK,JPTS), INTENT(in   ) ::   pts   ! pot. temperature & salinity
+      REAL(sp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT(  out) ::   pab   ! thermal/haline expansion ratio
       !
       INTEGER  ::   ji, jj, jk                ! dummy loop indices
-      REAL(wp) ::   zt , zh , zs              ! local scalars
+      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars
       REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
       !!----------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('rab_2d')
-      !
-      pab(:,:,:) = 0._wp
+      IF( ln_timing )   CALL timing_start('rab_3d')
       !
       SELECT CASE ( neos )
@@ -816,9 +1098,10 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-            !
-            zh  = pdep(ji,jj) * r1_Z0                                  ! depth
-            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
-            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+            !
+            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
+            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            ztm = tmask(ji,jj,jk)                                         ! tmask
             !
             ! alpha
@@ -841,5 +1124,5 @@
             zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
             !
-            pab(ji,jj,jp_tem) = zn * r1_rho0
+            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm
             !
             ! beta
@@ -862,4 +1145,122 @@
             zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
             !
+            pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm
+            !
+         END_3D
+         !
+      CASE( np_seos )                  !==  simplified EOS  ==!
+         !
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+            zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
+            zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
+            zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point
+            ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask
+            !
+            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
+            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha
+            !
+            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
+            pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta
+            !
+         END_3D
+         !
+      CASE DEFAULT
+         WRITE(ctmp1,*) '          bad flag value for neos = ', neos
+         CALL ctl_stop( 'rab_3d:', ctmp1 )
+         !
+      END SELECT
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=REAL(pab(:,:,:,jp_tem), wp), clinfo1=' rab_3d_t: ', &
+         &                                  tab3d_2=REAL(pab(:,:,:,jp_sal),  wp), clinfo2=' rab_3d_s : ', kdim=jpk )
+      !
+      IF( ln_timing )   CALL timing_stop('rab_3d')
+      !
+   END SUBROUTINE rab_3d_t_mixed
+
+
+   SUBROUTINE rab_2d( pts, pdep, pab, Kmm )
+      !!
+      INTEGER                   , INTENT(in   ) ::   Kmm   ! time level index
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pts    ! pot. temperature & salinity
+      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pdep   ! depth                  [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pab    ! thermal/haline expansion ratio
+      !!
+      CALL rab_2d_t(pts, is_tile(pts), pdep, is_tile(pdep), pab, is_tile(pab), Kmm)
+   END SUBROUTINE rab_2d
+
+
+   SUBROUTINE rab_2d_t( pts, ktts, pdep, ktdep, pab, ktab, Kmm )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE rab_2d  ***
+      !!
+      !! ** Purpose :   Calculates thermal/haline expansion ratio for a 2d field (unmasked)
+      !!
+      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
+      !!----------------------------------------------------------------------
+      INTEGER                            , INTENT(in   ) ::   Kmm   ! time level index
+      INTEGER                            , INTENT(in   ) ::   ktts, ktdep, ktab
+      REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in   ) ::   pts    ! pot. temperature & salinity
+      REAL(wp), DIMENSION(A2D_T(ktdep)    ), INTENT(in   ) ::   pdep   ! depth                  [m]
+      REAL(wp), DIMENSION(A2D_T(ktab),JPTS), INTENT(  out) ::   pab    ! thermal/haline expansion ratio
+      !
+      INTEGER  ::   ji, jj, jk                ! dummy loop indices
+      REAL(wp) ::   zt , zh , zs              ! local scalars
+      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('rab_2d')
+      !
+      pab(:,:,:) = 0._wp
+      !
+      SELECT CASE ( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
+         !
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            !
+            zh  = pdep(ji,jj) * r1_Z0                                  ! depth
+            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            !
+            ! alpha
+            zn3 = ALP003
+            !
+            zn2 = ALP012*zt + ALP102*zs+ALP002
+            !
+            zn1 = ((ALP031*zt   &
+               &   + ALP121*zs+ALP021)*zt   &
+               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
+               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
+               !
+            zn0 = ((((ALP050*zt   &
+               &   + ALP140*zs+ALP040)*zt   &
+               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
+               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
+               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
+               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            pab(ji,jj,jp_tem) = zn * r1_rho0
+            !
+            ! beta
+            zn3 = BET003
+            !
+            zn2 = BET012*zt + BET102*zs+BET002
+            !
+            zn1 = ((BET031*zt   &
+               &   + BET121*zs+BET021)*zt   &
+               &   + (BET211*zs+BET111)*zs+BET011)*zt   &
+               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
+               !
+            zn0 = ((((BET050*zt   &
+               &   + BET140*zs+BET040)*zt   &
+               &   + (BET230*zs+BET130)*zs+BET030)*zt   &
+               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
+               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
+               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
             pab(ji,jj,jp_sal) = zn / zs * r1_rho0
             !
@@ -997,5 +1398,5 @@
       !!
       INTEGER                              , INTENT(in   ) ::  Kmm   ! time level index
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu]
+      REAL(dp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu]
       REAL(wp), DIMENSION(:,:,:,:)         , INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celsius-1,psu-1]
       REAL(wp), DIMENSION(:,:,:)           , INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2]
@@ -1021,5 +1422,5 @@
       INTEGER                                , INTENT(in   ) ::  Kmm   ! time level index
       INTEGER                                , INTENT(in   ) ::  ktab, ktn2
-      REAL(wp), DIMENSION(jpi,jpj,  jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu]
+      REAL(dp), DIMENSION(jpi,jpj,  jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu]
       REAL(wp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celsius-1,psu-1]
       REAL(wp), DIMENSION(A2D_T(ktn2),JPK     ), INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2]
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/tranpc.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/tranpc.F90	(revision 15065)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/tranpc.F90	(revision 15066)
@@ -101,5 +101,5 @@
          !
          CALL eos_rab( CASTWP(pts(:,:,:,:,Kaa)), zab, Kmm )         ! after alpha and beta (given on T-points)
-         CALL bn2    ( CASTWP(pts(:,:,:,:,Kaa)), zab, zn2, Kmm )    ! after Brunt-Vaisala  (given on W-points)
+         CALL bn2    ( pts(:,:,:,:,Kaa), zab, zn2, Kmm )    ! after Brunt-Vaisala  (given on W-points)
          !
          IF( .NOT. l_istiled .OR. ntile == 1 ) nnpcc = 0         ! Do only on the first tile
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/ZDF/zdfosm.F90	(revision 15065)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/ZDF/zdfosm.F90	(revision 15066)
@@ -3295,5 +3295,5 @@
       ! w-level of the mixing and mixed layers
       CALL eos_rab( CASTWP(ts(:,:,:,:,Kmm)), rab_n, Kmm )
-      CALL bn2( CASTWP(ts(:,:,:,:,Kmm)), rab_n, rn2, Kmm )
+      CALL bn2( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm )
       imld_rst(:,:) = nlb10            ! Initialization to the number of w ocean point
       hbl(:,:) = 0.0_wp                ! Here hbl used as a dummy variable, integrating vertically N^2
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/step.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/step.F90	(revision 15065)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/step.F90	(revision 15066)
@@ -171,6 +171,6 @@
                          CALL eos_rab( CASTWP(ts(:,:,:,:,Nbb)), rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points
                          CALL eos_rab( CASTWP(ts(:,:,:,:,Nnn)), rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points
-                         CALL bn2    ( CASTWP(ts(:,:,:,:,Nbb)), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
-                         CALL bn2    ( CASTWP(ts(:,:,:,:,Nnn)), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency
+                         CALL bn2    ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
+                         CALL bn2    ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency
 
       !  VERTICAL PHYSICS
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/stpmlf.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/stpmlf.F90	(revision 15065)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/stpmlf.F90	(revision 15066)
@@ -176,8 +176,8 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !  THERMODYNAMICS
-                         CALL eos_rab( CASTWP(ts(:,:,:,:,Nbb)), rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points
-                         CALL eos_rab( CASTWP(ts(:,:,:,:,Nnn)), rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points
-                         CALL bn2    ( CASTWP(ts(:,:,:,:,Nbb)), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
-                         CALL bn2    ( CASTWP(ts(:,:,:,:,Nnn)), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency
+                         CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points
+                         CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points
+                         CALL bn2    ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
+                         CALL bn2    ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency
 
       !  VERTICAL PHYSICS
@@ -217,7 +217,7 @@
                          CALL ssh_nxt    ( kstp, Nbb, Nnn, ssh,  Naa )   ! after ssh (includes call to div_hor)
       IF( .NOT.lk_linssh ) THEN
-                         CALL dom_qco_r3c( CASTWP(ssh(:,:,Naa)), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa)           )   ! "after" ssh/h_0 ratio at t,u,v pts
+                         CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa)           )   ! "after" ssh/h_0 ratio at t,u,v pts
          IF( ln_dynspg_exp )   &
-            &            CALL dom_qco_r3c( CASTWP(ssh(:,:,Nnn)), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn), r3f(:,:) )   ! spg_exp : needed only for "now" ssh/h_0 ratio at f point
+            &            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn), r3f(:,:) )   ! spg_exp : needed only for "now" ssh/h_0 ratio at f point
       ENDIF
                          CALL wzv        ( kstp, Nbb, Nnn, Naa, ww  )    ! Nnn cross-level velocity
@@ -227,5 +227,5 @@
                             zgdept(:,:,jk) = gdept(:,:,jk,Nnn)
                          END DO
-                         CALL eos        ( CASTWP(ts(:,:,:,:,Nnn)), rhd, rhop, zgdept ) ! now in situ density for hpg computation
+                         CALL eos        ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept ) ! now in situ density for hpg computation
                          DEALLOCATE( zgdept )
 
@@ -268,5 +268,5 @@
                                       ! as well as vertical scale factors and vertical velocity need to be updated
                             CALL div_hor    ( kstp, Nbb, Nnn )                  ! Horizontal divergence  (2nd call in time-split case)
-            IF(.NOT.lk_linssh) CALL dom_qco_r3c( CASTWP(ssh(:,:,Naa)), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! update ssh/h_0 ratio at t,u,v,f pts
+            IF(.NOT.lk_linssh) CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! update ssh/h_0 ratio at t,u,v,f pts
          ENDIF
                             CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion
@@ -303,5 +303,5 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                          CALL ssh_atf    ( kstp, Nbb, Nnn, Naa, ssh )            ! time filtering of "now" sea surface height
-      IF(.NOT.lk_linssh) CALL dom_qco_r3c( CASTWP(ssh(:,:,Nnn)), r3t_f, r3u_f, r3v_f )   ! "now" ssh/h_0 ratio from filtrered ssh
+      IF(.NOT.lk_linssh) CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f )   ! "now" ssh/h_0 ratio from filtrered ssh
 #if defined key_top
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
