Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/do_loop_substitute.h90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/ICE/do_loop_substitute.h90	(revision 15066)
+++ 	(revision )
@@ -1,77 +1,0 @@
-#if defined show_comments
-! These comments are not intended to be retained during preprocessing; i.e. do not define "show_comments"
-!!----------------------------------------------------------------------
-!! NEMO/OCE 4.x , NEMO Consortium (2020)
-!! Software governed by the CeCILL license (see ./LICENSE)
-!!----------------------------------------------------------------------
-! This header file contains preprocessor definitions and macros used in the do-loop substitutions introduced
-! between version 4.0 and 4.2. The primary aim of these macros is to assist in future applications of tiling
-! to improve performance. This is expected to be achieved by alternative versions of these macros in selected
-! locations. The initial introduction of these macros simply replaced all identifiable nested 2D- and 3D-loops
-! with single line statements (and adjusts indenting accordingly). Do loops were identifiable if they comformed
-! to either:
-!                                       DO jk = ....
-!   DO jj = ....                           DO jj = ...
-!      DO ji = ....                           DO ji = ...
-!         .                   OR                 .
-!         .                                      .
-!     END DO                                  END DO
-!   END DO                                 END DO
-!                                       END DO
-! and white-space variants thereof.
-!
-! Additionally, only loops with recognised jj and ji loops limits were treated; these were:
-! Lower limits of 1, 2 or fs_2
-! Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)
-!
-! The macro naming convention takes the form: DO_2D( L, R, B, T) where:
-!   L is the Left   offset from the PE's inner domain;
-!   R is the Right  offset from the PE's inner domain
-!   B is the Bottom offset from the PE's inner domain;
-!   T is the Top    offset from the PE's inner domain;
-!
-! So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:
-!
-!   DO jj = 2, jpj
-!      DO ji = 1, jpim1
-!         .
-!         .
-!      END DO
-!   END DO
-!
-! with:
-!
-!   DO_2D( 1, 0, 0, 1 )
-!      .
-!      .
-!   END_2D
-!
-! similar conventions apply to the 3D loops macros. jk loop limits are retained through macro arguments
-! and are not restricted. This includes the possibility of strides for which an extra set of DO_3DS
-! macros are defined.
-!
-! In the following definitions the inner PE domain is defined by start indices of (Nis0, Njs0) and end
-! indices of (Nie0, Nje0) where:
-!
-! Nis0 =   1 + nn_hls     Njs0 =   1 + nn_hls
-! Nie0 = jpi - nn_hls     Nje0 = jpj - nn_hls
-!
-#endif
-
-#define DO_2D(L, R, B, T) DO jj = ntsj-(B), ntej+(T)   ;   DO ji = ntsi-(L), ntei+(R)
-#define A1Di(H) ntsi-H:ntei+H
-#define A1Dj(H) ntsj-H:ntej+H
-#define A2D(H) A1Di(H),A1Dj(H)
-#define A1Di_T(T) (ntsi-nn_hls-1)*T+1:
-#define A1Dj_T(T) (ntsj-nn_hls-1)*T+1:
-#define A2D_T(T) A1Di_T(T),A1Dj_T(T)
-#define JPK  :
-#define JPTS  :
-#define KJPT  :
-
-#define DO_3D(L, R, B, T, ks, ke) DO jk = ks, ke   ;   DO_2D(L, R, B, T)
-
-#define DO_3DS(L, R, B, T, ks, ke, ki) DO jk = ks, ke, ki   ;   DO_2D(L, R, B, T)
-
-#define END_2D   END DO   ;   END DO
-#define END_3D   END DO   ;   END DO   ;   END DO
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcblk_algo_coare3p0.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 15066)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 15520)
@@ -91,5 +91,5 @@
    SUBROUTINE turb_coare3p0( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, &
       &                      Cd, Ch, Ce, t_zu, q_zu, Ubzu,                               &
-      &                      nb_iter, Cdn, Chn, Cen,                                     & ! optional output
+      &                      nb_iter, CdN, ChN, CeN,                                     & ! optional output
       &                      Qsw, rad_lw, slp, pdT_cs,                                   & ! optionals for cool-skin (and warm-layer)
       &                      pdT_wl, pHz_wl )                                              ! optionals for warm-layer only
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcblk_algo_coare3p6.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 15066)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 15520)
@@ -81,5 +81,5 @@
    SUBROUTINE turb_coare3p6( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, &
       &                      Cd, Ch, Ce, t_zu, q_zu, Ubzu,                               &
-      &                      nb_iter, Cdn, Chn, Cen,                                     & ! optional output
+      &                      nb_iter, CdN, ChN, CeN,                                     & ! optional output
       &                      Qsw, rad_lw, slp, pdT_cs,                                   & ! optionals for cool-skin (and warm-layer)
       &                      pdT_wl, pHz_wl )                                              ! optionals for warm-layer only
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 15066)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 15520)
@@ -85,5 +85,5 @@
    SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, &
       &                      Cd, Ch, Ce, t_zu, q_zu, Ubzu,                            &
-      &                      nb_iter, Cdn, Chn, Cen,                                           & ! optional output
+      &                      nb_iter, CdN, ChN, CeN,                                           & ! optional output
       &                      Qsw, rad_lw, slp, pdT_cs,                                & ! optionals for cool-skin (and warm-layer)
       &                      pdT_wl, pHz_wl )                                           ! optionals for warm-layer only
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 15066)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/eosbn2.F90	(revision 15520)
@@ -56,11 +56,9 @@
    !                  !! * Interface
    INTERFACE eos
-      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
+      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d
    END INTERFACE
    !
    INTERFACE eos_rab
-      MODULE PROCEDURE rab_3d_wp, rab_2d, rab_0d
-      MODULE PROCEDURE rab_3d_mixed
+      MODULE PROCEDURE rab_3d, rab_2d, rab_0d
    END INTERFACE
    !
@@ -192,5 +190,5 @@
 CONTAINS
 
-   SUBROUTINE eos_insitu_wp( pts, prd, pdep )
+   SUBROUTINE eos_insitu( pts, prd, pdep )
       !!
       REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
@@ -199,18 +197,9 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pdep  ! depth                      [m]
       !!
-      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 )
+      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 
@@ -318,115 +307,7 @@
       IF( ln_timing )   CALL timing_stop('eos-insitu')
       !
-   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 )
+   END SUBROUTINE eos_insitu_t
+  
+   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
       !!
       REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
@@ -436,20 +317,8 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pdep   ! depth                      [m]
       !!
-      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 )
+      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 )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_insitu_pot  ***
@@ -605,163 +474,7 @@
       IF( ln_timing )   CALL timing_stop('eos-pot')
       !
-   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
-
-
+   END SUBROUTINE eos_insitu_pot_t
+
+  
    SUBROUTINE eos_insitu_2d( pts, pdep, prd )
       !!
@@ -947,5 +660,5 @@
 
 
-   SUBROUTINE rab_3d_wp( pts, pab, Kmm )
+   SUBROUTINE rab_3d( pts, pab, Kmm )
       !!
       INTEGER                     , INTENT(in   ) ::   Kmm   ! time level index
@@ -953,18 +666,8 @@
       REAL(wp), DIMENSION(:,:,:,:), INTENT(  out) ::   pab   ! thermal/haline expansion ratio
       !!
-      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 )
+      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 )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE rab_3d  ***
@@ -1070,112 +773,5 @@
       IF( ln_timing )   CALL timing_stop('rab_3d')
       !
-   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, 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 , ztm        ! local scalars
-      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('rab_3d')
-      !
-      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  = 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
-            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,jk,jp_tem) = zn * r1_rho0 * ztm
-            !
-            ! 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,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
-
+   END SUBROUTINE rab_3d_t
 
    SUBROUTINE rab_2d( pts, pdep, pab, Kmm )
