Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/dia25h.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/dia25h.F90	(revision 14986)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/dia25h.F90	(revision 15540)
@@ -25,10 +25,10 @@
    ! variables for calculating 25-hourly means
    INTEGER , SAVE ::   cnt_25h           ! Counter for 25 hour means
-   REAL(wp), SAVE ::   r1_25 = 0.04_wp   ! =1/25 
-   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   tn_25h  , sn_25h
-   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   sshn_25h 
-   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   un_25h  , vn_25h  , wn_25h
-   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avt_25h , avm_25h
-   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   en_25h  , rmxln_25h
+   REAL(dp), SAVE ::   r1_25 = 0.04_wp   ! =1/25 
+   REAL(dp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   tn_25h  , sn_25h
+   REAL(dp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   sshn_25h 
+   REAL(dp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   un_25h  , vn_25h  , wn_25h
+   REAL(dp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avt_25h , avm_25h
+   REAL(dp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   en_25h  , rmxln_25h
 
    !!----------------------------------------------------------------------
@@ -130,9 +130,9 @@
       INTEGER                          ::   iyear0, nimonth0,iday0            ! start year,imonth,day
       LOGICAL ::   ll_print = .FALSE.    ! =T print and flush numout
-      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi   ! local scalars
+      REAL(dp)                         ::   zsto, zout, zmax, zjulian, zmdi   ! local scalars
       INTEGER                          ::   i_steps                           ! no of timesteps per hour
-      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                ! workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                              ! workspace
-      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                             ! workspace
+      REAL(dp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                ! workspace
+      REAL(dp), DIMENSION(jpi,jpj,jpk) ::   zw3d                              ! workspace
+      REAL(dp), DIMENSION(jpi,jpj,3)   ::   zwtmb                             ! workspace
       !!----------------------------------------------------------------------
 
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diaar5.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diaar5.F90	(revision 14986)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diaar5.F90	(revision 15540)
@@ -30,6 +30,6 @@
    PUBLIC   dia_ar5_hst    ! heat/salt transport
 
-   REAL(wp)                         ::   vol0         ! ocean volume (interior domain)
-   REAL(wp)                         ::   area_tot     ! total ocean surface (interior domain)
+   REAL(dp)                         ::   vol0         ! ocean volume (interior domain)
+   REAL(dp)                         ::   area_tot     ! total ocean surface (interior domain)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   thick0       ! ocean thickness (interior domain)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sn0          ! initial salinity
@@ -78,6 +78,8 @@
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z2d, zpe                   ! 2D workspace
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: z3d, zrhd, ztpot, zgdept   ! 3D workspace (zgdept: needed to use the substitute)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)      :: z2d! 2D workspace
+      REAL(dp), ALLOCATABLE, DIMENSION(:,:)      :: zpe! 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)    :: zgdept! 3D workspace (zgdept: needed to use the substitute)
+      REAL(dp), ALLOCATABLE, DIMENSION(:,:,:)    :: z3d, zrhd, ztpot! 3D workspace (zgdept: needed to use the substitute)
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace
 
@@ -137,5 +139,5 @@
             zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
          END DO
-         CALL eos( ztsn, zrhd, zgdept)                       ! now in situ density using initial salinity
+         CALL eos( CASTDP(ztsn), zrhd, CASTDP(zgdept))                       ! now in situ density using initial salinity
          !
          zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
@@ -232,5 +234,5 @@
             ztpot(:,:,jpk) = 0._wp
             DO jk = 1, jpkm1
-               ztpot(:,:,jk) = eos_pt_from_ct( CASTWP(ts(:,:,jk,jp_tem,Kmm)), CASTWP(ts(:,:,jk,jp_sal,Kmm)) )
+               ztpot(:,:,jk) = eos_pt_from_ct( ts(:,:,jk,jp_tem,Kmm), ts(:,:,jk,jp_sal,Kmm) )
             END DO
             !
@@ -315,9 +317,9 @@
       INTEGER                         , INTENT(in )  :: ktra  ! tracer index
       CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in)   :: puflx  ! u-flux of advection/diffusion
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in)   :: pvflx  ! v-flux of advection/diffusion
+      REAL(dp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in)   :: puflx  ! u-flux of advection/diffusion
+      REAL(dp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in)   :: pvflx  ! v-flux of advection/diffusion
       !
       INTEGER    ::  ji, jj, jk
-      REAL(wp), DIMENSION(A2D(nn_hls))  :: z2d
+      REAL(dp), DIMENSION(A2D(nn_hls))  :: z2d
 
       z2d(:,:) = puflx(:,:,1)
@@ -360,5 +362,5 @@
       INTEGER  ::   ji, jj, jk  ! dummy loop indices
       REAL(wp) ::   zztmp
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity
+      REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zvol0
       !
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diacfl.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diacfl.F90	(revision 14986)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diacfl.F90	(revision 15540)
@@ -27,5 +27,5 @@
    !
    INTEGER, DIMENSION(3) ::   nCu_loc, nCv_loc, nCw_loc   ! U, V, and W run max locations in the global domain
-   REAL(wp)              ::   rCu_max, rCv_max, rCw_max   ! associated run max Courant number 
+   REAL(dp)              ::   rCu_max, rCv_max, rCw_max   ! associated run max Courant number 
 
    PUBLIC   dia_cfl       ! routine called by step.F90
@@ -53,7 +53,7 @@
       !
       INTEGER                          ::   ji, jj, jk                       ! dummy loop indices
-      REAL(wp)                         ::   zCu_max, zCv_max, zCw_max        ! local scalars
+      REAL(dp)                         ::   zCu_max, zCv_max, zCw_max        ! local scalars
       INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace
+      REAL(dp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace
       LOGICAL , DIMENSION(jpi,jpj,jpk) ::   llmsk
       !!----------------------------------------------------------------------
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diadct.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diadct.F90	(revision 14986)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diadct.F90	(revision 15540)
@@ -61,5 +61,5 @@
 
    TYPE COORD_SECTION
-      REAL(wp) :: lon,lat
+      REAL(dp) :: lon,lat
    END TYPE COORD_SECTION
 
@@ -73,11 +73,11 @@
       INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section
       CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class
-      REAL(wp), DIMENSION(nb_class_max)            :: zsigi             ! in-situ   density classes    (99 if you don't want)
-      REAL(wp), DIMENSION(nb_class_max)            :: zsigp             ! potential density classes    (99 if you don't want)
-      REAL(wp), DIMENSION(nb_class_max)            :: zsal              ! salinity classes   (99 if you don't want)
-      REAL(wp), DIMENSION(nb_class_max)            :: ztem              ! temperature classes(99 if you don't want)
-      REAL(wp), DIMENSION(nb_class_max)            :: zlay              ! level classes      (99 if you don't want)
-      REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output
-      REAL(wp)                                         :: slopeSection  ! slope of the section
+      REAL(dp), DIMENSION(nb_class_max)            :: zsigi             ! in-situ   density classes    (99 if you don't want)
+      REAL(dp), DIMENSION(nb_class_max)            :: zsigp             ! potential density classes    (99 if you don't want)
+      REAL(dp), DIMENSION(nb_class_max)            :: zsal              ! salinity classes   (99 if you don't want)
+      REAL(dp), DIMENSION(nb_class_max)            :: ztem              ! temperature classes(99 if you don't want)
+      REAL(dp), DIMENSION(nb_class_max)            :: zlay              ! level classes      (99 if you don't want)
+      REAL(dp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output
+      REAL(dp)                                         :: slopeSection  ! slope of the section
       INTEGER                                          :: nb_point      ! number of points in the section
       TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! list of points in the sections
@@ -86,6 +86,6 @@
    TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
  
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d 
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d  
+   REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d 
+   REAL(dp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d  
 
 
@@ -202,6 +202,6 @@
      INTEGER              , DIMENSION(1)    ::   ish     ! work array for mpp_sum
      INTEGER              , DIMENSION(3)    ::   ish2    !   "
-     REAL(wp), ALLOCATABLE, DIMENSION(:)    ::   zwork   !   "  
-     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)::   zsum    !   "
+     REAL(dp), ALLOCATABLE, DIMENSION(:)    ::   zwork   !   "  
+     REAL(dp), ALLOCATABLE, DIMENSION(:,:,:)::   zsum    !   "
      !!---------------------------------------------------------------------    
      !
@@ -587,7 +587,7 @@
      !
      INTEGER ::   jk, jseg, jclass,jl, isgnu, isgnv    ! loop on level/segment/classes/ice categories
-     REAL(wp)::   zumid, zvmid, zumid_ice, zvmid_ice   ! U/V ocean & ice velocity on a cell segment 
-     REAL(wp)::   zTnorm                               ! transport of velocity through one cell's sides 
-     REAL(wp)::   ztn, zsn, zrhoi, zrhop, zsshn, zdep  ! temperature/salinity/potential density/ssh/depth at u/v point
+     REAL(dp)::   zumid, zvmid, zumid_ice, zvmid_ice   ! U/V ocean & ice velocity on a cell segment 
+     REAL(dp)::   zTnorm                               ! transport of velocity through one cell's sides 
+     REAL(dp)::   ztn, zsn, zrhoi, zrhop, zsshn, zdep  ! temperature/salinity/potential density/ssh/depth at u/v point
      TYPE(POINT_SECTION) ::   k
       !!--------------------------------------------------------
@@ -680,5 +680,5 @@
                   zsn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) ) 
                   zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop) 
-                  zrhoi = interp(Kmm,k%I,k%J,jk,'V',CASTDP(rhd*rho0+rho0)) 
+                  zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rho0+rho0) 
                   zsshn =  0.5*( ssh(k%I,k%J,Kmm) + ssh(k%I,k%J+1,Kmm)    ) * vmask(k%I,k%J,1) 
                CASE(2,3) 
@@ -686,5 +686,5 @@
                   zsn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) ) 
                   zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop) 
-                  zrhoi = interp(Kmm,k%I,k%J,jk,'U',CASTDP(rhd*rho0+rho0)) 
+                  zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rho0+rho0) 
                   zsshn =  0.5*( ssh(k%I,k%J,Kmm) + ssh(k%I+1,k%J,Kmm)    ) * umask(k%I,k%J,1)  
                END SELECT 
@@ -794,5 +794,5 @@
      TYPE(POINT_SECTION) :: k 
      INTEGER  :: jk,jseg,jclass                        ! dummy variables for looping on level/segment/classes  
-     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zdep ! temperature/salinity/ssh/potential density /depth at u/v point 
+     REAL(dp) :: ztn, zsn, zrhoi, zrhop, zsshn, zdep ! temperature/salinity/ssh/potential density /depth at u/v point 
      !!------------------------------------------------------------- 
  
@@ -853,5 +853,5 @@
                  zsn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) ) 
                  zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop) 
-                 zrhoi = interp(Kmm,k%I,k%J,jk,'V',CASTDP(rhd*rho0+rho0))
+                 zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rho0+rho0)
 
               CASE(2,3) 
@@ -859,5 +859,5 @@
                  zsn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) ) 
                  zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop) 
-                 zrhoi = interp(Kmm,k%I,k%J,jk,'U',CASTDP(rhd*rho0+rho0)) 
+                 zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rho0+rho0) 
                  zsshn =  0.5*( ssh(k%I,k%J,Kmm)    + ssh(k%I+1,k%J,Kmm)    ) * umask(k%I,k%J,1)  
               END SELECT 
@@ -987,8 +987,8 @@
      INTEGER               :: jclass             ! Dummy loop
      CHARACTER(len=2)      :: classe             ! Classname 
-     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
-     REAL(wp)              :: zslope             ! section's slope coeff
+     REAL(dp)              :: zbnd1,zbnd2        ! Class bounds
+     REAL(dp)              :: zslope             ! section's slope coeff
      !
-     REAL(wp), DIMENSION(nb_type_class)::   zsumclasses   ! 1D workspace 
+     REAL(dp), DIMENSION(nb_type_class)::   zsumclasses   ! 1D workspace 
      !!------------------------------------------------------------- 
 
@@ -1171,12 +1171,12 @@
   CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
   REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab         ! variable to compute at (ki, kj, kk )
-  REAL(wp)                                     :: interp       ! interpolated variable 
+  REAL(dp)                                     :: interp       ! interpolated variable 
 
   !*local declations
   INTEGER :: ii1, ij1, ii2, ij2                                ! local integer
-  REAL(wp):: ze3t, ze3, zwgt1, zwgt2, zbis, zdepu            ! local real
-  REAL(wp):: zet1, zet2                                        ! weight for interpolation 
-  REAL(wp):: zdep1,zdep2                                       ! differences of depth
-  REAL(wp):: zmsk                                              ! mask value
+  REAL(dp):: ze3t, ze3, zwgt1, zwgt2, zbis, zdepu            ! local real
+  REAL(dp):: zet1, zet2                                        ! weight for interpolation 
+  REAL(dp):: zdep1,zdep2                                       ! differences of depth
+  REAL(dp):: zmsk                                              ! mask value
   !!----------------------------------------------------------------------
 
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diadetide.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diadetide.F90	(revision 14986)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diadetide.F90	(revision 15540)
@@ -22,5 +22,5 @@
    LOGICAL, PUBLIC                               ::   lk_diadetide
    INTEGER                                       ::   ndiadetide
-   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::   tdiadetide
+   REAL(dp), SAVE, ALLOCATABLE, DIMENSION(:)     ::   tdiadetide
 
    PUBLIC ::   dia_detide_init, dia_detide
@@ -92,6 +92,6 @@
 
       INTEGER, INTENT(in)          ::   kt
-      REAL(wp), DIMENSION(jpi,jpj) ::   zwght_2D
-      REAL(wp)                     ::   zwght, ztmp
+      REAL(dp), DIMENSION(jpi,jpj) ::   zwght_2D
+      REAL(dp)                     ::   zwght, ztmp
       INTEGER                      ::   jn
 
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diahsb.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diahsb.F90	(revision 14986)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diahsb.F90	(revision 15540)
@@ -40,15 +40,15 @@
    LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets
 
-   REAL(wp) ::   surf_tot              ! ocean surface
-   REAL(wp) ::   frc_t, frc_s, frc_v   ! global forcing trends
-   REAL(wp) ::   frc_wn_t, frc_wn_s    ! global forcing trends
+   REAL(dp) ::   surf_tot              ! ocean surface
+   REAL(dp) ::   frc_t, frc_s, frc_v   ! global forcing trends
+   REAL(dp) ::   frc_wn_t, frc_wn_s    ! global forcing trends
    !
    REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   surf 
    REAL(dp), DIMENSION(:,:)  , ALLOCATABLE  :: surf_ini
-   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE :: ssh_ini                            !
-   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   !
+   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE :: ssh_ini                            !
+   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   !
    REAL(dp), DIMENSION(:,:,:), ALLOCATABLE  :: sc_loc_ini, e3t_ini               !
-   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini                       !
-   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   tmask_ini
+   REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini                       !
+   REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::   tmask_ini
 
    !! * Substitutions
@@ -78,5 +78,5 @@
       REAL(wp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
       REAL(wp)   ::   zdiff_hc1   , zdiff_sc1     !  -         -     -        - 
-      REAL(wp)   ::   zdiff_v1          ! volume variation
+      REAL(dp)   ::   zdiff_v1          ! volume variation
       REAL(dp)   ::   zdiff_v2
       REAL(wp)   ::   zerr_hc1    , zerr_sc1      ! heat and salt content misfit
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diahth.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diahth.F90	(revision 14986)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diahth.F90	(revision 15540)
@@ -31,11 +31,11 @@
    
    ! note: following variables should move to local variables once iom_put is always used 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hth    !: depth of the max vertical temperature gradient [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd20   !: depth of 20 C isotherm                         [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd26   !: depth of 26 C isotherm                         [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd28   !: depth of 28 C isotherm                         [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc3   !: heat content of first 300 m                    [W]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc7   !: heat content of first 700 m                    [W]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc20  !: heat content of first 2000 m                   [W]
+   REAL(dp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hth    !: depth of the max vertical temperature gradient [m]
+   REAL(dp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd20   !: depth of 20 C isotherm                         [m]
+   REAL(dp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd26   !: depth of 26 C isotherm                         [m]
+   REAL(dp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd28   !: depth of 28 C isotherm                         [m]
+   REAL(dp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc3   !: heat content of first 300 m                    [W]
+   REAL(dp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc7   !: heat content of first 700 m                    [W]
+   REAL(dp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc20  !: heat content of first 2000 m                   [W]
 
 
@@ -262,5 +262,5 @@
          IF( iom_use ('hc300') ) THEN  
             zzdep = 300.
-            CALL  dia_hth_htc( Kmm, zzdep, CASTWP(ts(:,:,:,jp_tem,Kmm)), htc3 )
+            CALL  dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc3 )
             CALL iom_put( 'hc300', rho0_rcp * htc3 )  ! vertically integrated heat content (J/m2)
          ENDIF
@@ -271,5 +271,5 @@
          IF( iom_use ('hc700') ) THEN  
             zzdep = 700.
-            CALL  dia_hth_htc( Kmm, zzdep, CASTWP(ts(:,:,:,jp_tem,Kmm)), htc7 )
+            CALL  dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc7 )
             CALL iom_put( 'hc700', rho0_rcp * htc7 )  ! vertically integrated heat content (J/m2)
   
@@ -281,5 +281,5 @@
          IF( iom_use ('hc2000') ) THEN  
             zzdep = 2000.
-            CALL  dia_hth_htc( Kmm, zzdep, CASTWP(ts(:,:,:,jp_tem,Kmm)), htc20 )
+            CALL  dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc20 )
             CALL iom_put( 'hc2000', rho0_rcp * htc20 )  ! vertically integrated heat content (J/m2)  
          ENDIF
@@ -296,8 +296,8 @@
       INTEGER , INTENT(in) :: Kmm      ! ocean time level index
       REAL(wp), INTENT(in) :: ptem
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdept     
+      REAL(dp), DIMENSION(jpi,jpj), INTENT(out) :: pdept     
       !
       INTEGER  :: ji, jj, jk, iid
-      REAL(wp) :: zztmp, zzdep
+      REAL(dp) :: zztmp, zzdep
       INTEGER, DIMENSION(jpi,jpj) :: iktem
       
@@ -337,9 +337,9 @@
       INTEGER , INTENT(in) :: Kmm      ! ocean time level index
       REAL(wp), INTENT(in) :: pdep     ! depth over the heat content
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pt   
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phtc  
+      REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pt   
+      REAL(dp), DIMENSION(jpi,jpj), INTENT(inout) :: phtc  
       !
       INTEGER  :: ji, jj, jk, ik
-      REAL(wp), DIMENSION(jpi,jpj) :: zthick
+      REAL(dp), DIMENSION(jpi,jpj) :: zthick
       INTEGER , DIMENSION(jpi,jpj) :: ilevel
 
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diamlr.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diamlr.F90	(revision 14986)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diamlr.F90	(revision 15540)
@@ -409,5 +409,5 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj) ::   zadatrj2d
+      REAL(dp), DIMENSION(jpi,jpj) ::   zadatrj2d
       !!----------------------------------------------------------------------
 
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/dianam.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/dianam.F90	(revision 14986)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/dianam.F90	(revision 15540)
@@ -63,5 +63,5 @@
       REAL(dp)           ::   zsec1, zsec2                     ! not used
       REAL(dp)           ::          zjul                      ! temporary scalars
-      REAL(wp)           ::   zdrun                      ! temporary scalars
+      REAL(dp)           ::   zdrun                      ! temporary scalars
       !!----------------------------------------------------------------------
 
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diaptr.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diaptr.F90	(revision 14986)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diaptr.F90	(revision 15540)
@@ -44,7 +44,7 @@
    PUBLIC   dia_ptr_hst    ! called from tra_ldf/tra_adv routines
 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hstr_adv, hstr_ldf, hstr_eiv   !: Heat/Salt TRansports(adv, diff, Bolus.)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hstr_ove, hstr_btr, hstr_vtr   !: heat Salt TRansports(overturn, baro, merional)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   pvtr_int, pzon_int             !: Other zonal integrals
+   REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hstr_adv, hstr_ldf, hstr_eiv   !: Heat/Salt TRansports(adv, diff, Bolus.)
+   REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hstr_ove, hstr_btr, hstr_vtr   !: heat Salt TRansports(overturn, baro, merional)
+   REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   pvtr_int, pzon_int             !: Other zonal integrals
 
    LOGICAL, PUBLIC    ::   l_diaptr       !: tracers  trend flag
@@ -52,10 +52,10 @@
    INTEGER, PARAMETER ::   jp_vtr = 4
 
-   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
-   REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rho0 x Cp)
-   REAL(wp) ::   rc_ggram = 1.e-9_wp   ! conversion from g    to Gg  (further x rho0)
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk34 ! mask out Southern Ocean (=0 south of 34°S)
+   REAL(dp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
+   REAL(dp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rho0 x Cp)
+   REAL(dp) ::   rc_ggram = 1.e-9_wp   ! conversion from g    to Gg  (further x rho0)
+
+   REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks
+   REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk34 ! mask out Southern Ocean (=0 south of 34°S)
 
    LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag
@@ -63,4 +63,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "single_precision_substitute.h90"
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
@@ -113,13 +114,13 @@
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
-      REAL(wp), DIMENSION(jpj)      ::  zvsum, ztsum, zssum   ! 1D workspace
+      REAL(dp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
+      REAL(dp), DIMENSION(jpj)      ::  zvsum, ztsum, zssum   ! 1D workspace
       !
       !overturning calculation
-      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   sjk, r1_sjk, v_msf  ! i-mean i-k-surface and its inverse
-      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   zt_jk, zs_jk        ! i-mean T and S, j-Stream-Function
-
-      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   z4d1, z4d2
-      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   z3dtr
+      REAL(dp), DIMENSION(:,:,:  ), ALLOCATABLE ::   sjk, r1_sjk, v_msf  ! i-mean i-k-surface and its inverse
+      REAL(dp), DIMENSION(:,:,:  ), ALLOCATABLE ::   zt_jk, zs_jk        ! i-mean T and S, j-Stream-Function
+
+      REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE ::   z4d1, z4d2
+      REAL(dp), DIMENSION(:,:,:  ), ALLOCATABLE ::   z3dtr
       !!----------------------------------------------------------------------
       !
@@ -349,9 +350,9 @@
       INTEGER                     , INTENT(in)           :: Kmm          ! time level index
       REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in), OPTIONAL :: pvtr         ! j-effective transport
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE            :: zmask        ! 3D workspace
-      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE          :: zts          ! 4D workspace
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE            :: sjk, v_msf   ! Zonal sum: i-k surface area, j-effective transport
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE            :: zt_jk, zs_jk ! Zonal sum: i-k surface area * (T, S)
-      REAL(wp)                                           :: zsfc, zvfc   ! i-k surface area
+      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE            :: zmask        ! 3D workspace
+      REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE          :: zts          ! 4D workspace
+      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE            :: sjk, v_msf   ! Zonal sum: i-k surface area, j-effective transport
+      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE            :: zt_jk, zs_jk ! Zonal sum: i-k surface area * (T, S)
+      REAL(dp)                                           :: zsfc, zvfc   ! i-k surface area
       INTEGER  ::   ji, jj, jk, jn                                       ! dummy loop indices
       !!----------------------------------------------------------------------
@@ -363,5 +364,5 @@
 
             DO jn = 1, nbasin
-               v_msf(:,:,jn) = ptr_sjk( pvtr(:,:,:), btmsk34(:,:,jn) )
+               v_msf(:,:,jn) =ptr_sjk( CASTDP(pvtr(:,:,:)), btmsk34(:,:,jn) )
             ENDDO
 
@@ -538,6 +539,6 @@
       INTEGER                         , INTENT(in )  :: ktra  ! tracer index
       CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv'
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in)   :: pvflx ! 3D input array of advection/diffusion
-      REAL(wp), DIMENSION(A1Dj(nn_hls),nbasin)                 :: zsj   !
+      REAL(dp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in)   :: pvflx ! 3D input array of advection/diffusion
+      REAL(dp), DIMENSION(A1Dj(nn_hls),nbasin)                 :: zsj   !
       INTEGER                                        :: jn    !
 
@@ -574,11 +575,11 @@
       !! ** Action  : phstr
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpj,nbasin) , INTENT(inout)         ::  phstr  !
-      REAL(wp), DIMENSION(A1Dj(nn_hls),nbasin), INTENT(in)            ::  pva    !
+      REAL(dp), DIMENSION(jpj,nbasin) , INTENT(inout)         ::  phstr  !
+      REAL(dp), DIMENSION(A1Dj(nn_hls),nbasin), INTENT(in)            ::  pva    !
       INTEGER                                               ::  jj
 #if ! defined key_mpi_off
       INTEGER, DIMENSION(1)           ::  ish1d
       INTEGER, DIMENSION(2)           ::  ish2d
-      REAL(wp), DIMENSION(jpj*nbasin) ::  zwork
+      REAL(dp), DIMENSION(jpj*nbasin) ::  zwork
 #endif
 
@@ -610,11 +611,11 @@
       !! ** Action  : phstr
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpj,jpk,nbasin) , INTENT(inout)     ::  phstr  !
-      REAL(wp), DIMENSION(A1Dj(nn_hls),jpk,nbasin), INTENT(in)        ::  pva    !
+      REAL(dp), DIMENSION(jpj,jpk,nbasin) , INTENT(inout)     ::  phstr  !
+      REAL(dp), DIMENSION(A1Dj(nn_hls),jpk,nbasin), INTENT(in)        ::  pva    !
       INTEGER                                               ::  jj, jk
 #if ! defined key_mpi_off
       INTEGER, DIMENSION(1)              ::  ish1d
       INTEGER, DIMENSION(3)              ::  ish3d
-      REAL(wp), DIMENSION(jpj*jpk*nbasin)  ::  zwork
+      REAL(dp), DIMENSION(jpj*jpk*nbasin)  ::  zwork
 #endif
 
@@ -675,9 +676,9 @@
       !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(A2D(nn_hls),jpk)  ::   pvflx  ! mask flux array at V-point
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj)  ::   pmsk   ! Optional 2D basin mask
+      REAL(dp), INTENT(in), DIMENSION(A2D(nn_hls),jpk)  ::   pvflx  ! mask flux array at V-point
+      REAL(dp), INTENT(in), DIMENSION(jpi,jpj)  ::   pmsk   ! Optional 2D basin mask
       !
       INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
-      REAL(wp), DIMENSION(A1Dj(nn_hls)) :: p_fval  ! function value
+      REAL(dp), DIMENSION(A1Dj(nn_hls)) :: p_fval  ! function value
       !!--------------------------------------------------------------------
       !
@@ -704,5 +705,5 @@
       !
       INTEGER                  ::   ji,jj       ! dummy loop arguments
-      REAL(wp), DIMENSION(A1Dj(nn_hls)) :: p_fval ! function value
+      REAL(dp), DIMENSION(A1Dj(nn_hls)) :: p_fval ! function value
       !!--------------------------------------------------------------------
       !
@@ -723,9 +724,9 @@
       !! ** Action  : - p_fval: j-cumulated sum of pva
       !!----------------------------------------------------------------------
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)  ::   pva   ! mask flux array at V-point
+      REAL(dp) , INTENT(in), DIMENSION(jpi,jpj)  ::   pva   ! mask flux array at V-point
       !
       INTEGER                  ::   ji,jj,jc       ! dummy loop arguments
       INTEGER                  ::   ijpj        ! ???
-      REAL(wp), DIMENSION(jpi,jpj) :: p_fval ! function value
+      REAL(dp), DIMENSION(jpi,jpj) :: p_fval ! function value
       !!--------------------------------------------------------------------
       !
@@ -754,9 +755,9 @@
       !!
       IMPLICIT none
-      REAL(wp) , INTENT(in), DIMENSION(A2D(nn_hls),jpk) ::   pta    ! mask flux array at V-point
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
+      REAL(dp) , INTENT(in), DIMENSION(A2D(nn_hls),jpk) ::   pta    ! mask flux array at V-point
+      REAL(dp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
       !!
       INTEGER                           :: ji, jj, jk ! dummy loop arguments
-      REAL(wp), DIMENSION(A1Dj(nn_hls),jpk) :: p_fval     ! return function value
+      REAL(dp), DIMENSION(A1Dj(nn_hls),jpk) :: p_fval     ! return function value
       !!--------------------------------------------------------------------
       !
Index: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diawri.F90	(revision 14986)
+++ NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diawri.F90	(revision 15540)
@@ -123,5 +123,5 @@
       REAL(wp)::   zztmp2, zztmpy   !   -      -
       REAL(wp)::   ze3
-      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace
+      REAL(dp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
       !!----------------------------------------------------------------------
@@ -1092,5 +1092,5 @@
       !!
       INTEGER :: inum, jk
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept       ! 3D workspace for qco substitution
+      REAL(dp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept       ! 3D workspace for qco substitution
       !!----------------------------------------------------------------------
       ! 
