diff --git a/ext/AGRIF/AGRIF_FILES/modbcfunction.F90 b/ext/AGRIF/AGRIF_FILES/modbcfunction.F90
index e6add0329e709fbd64311d119d638024a7ad141d..ae7de45aa4f67ad6544725c894575005c8d1f88c 100644
--- a/ext/AGRIF/AGRIF_FILES/modbcfunction.F90
+++ b/ext/AGRIF/AGRIF_FILES/modbcfunction.F90
@@ -447,14 +447,14 @@ subroutine Agrif_Find_Nearest ( tabvarsindic, fineloc, parentloc )
     type(Agrif_Variable), pointer :: parent_var
     type(Agrif_Variable), pointer :: child_var
     integer :: i
-    integer, dimension(6)           :: nb_child     !< Number of cells on the child grid
-    integer, dimension(6)           :: ub_child     !< Upper bound on the child grid
-    integer, dimension(6)           :: lb_child     !< Lower bound on the child grid
-    integer, dimension(6)           :: lb_parent    !< Lower bound on the parent grid
-    real, dimension(6)              :: s_child      !< Child  grid position (s_root = 0)
-    real, dimension(6)              :: s_parent     !< Parent grid position (s_root = 0)
-    real, dimension(6)              :: ds_child     !< Child  grid dx (ds_root = 1)
-    real, dimension(6)              :: ds_parent    !< Parent grid dx (ds_root = 1)
+    integer,      dimension(6)      :: nb_child     !< Number of cells on the child grid
+    integer,      dimension(6)      :: ub_child     !< Upper bound on the child grid
+    integer,      dimension(6)      :: lb_child     !< Lower bound on the child grid
+    integer,      dimension(6)      :: lb_parent    !< Lower bound on the parent grid
+    real(kind=8), dimension(6)      :: s_child      !< Child  grid position (s_root = 0)
+    real(kind=8), dimension(6)      :: s_parent     !< Parent grid position (s_root = 0)
+    real(kind=8), dimension(6)      :: ds_child     !< Child  grid dx (ds_root = 1)
+    real(kind=8), dimension(6)      :: ds_parent    !< Parent grid dx (ds_root = 1)
     integer                         :: nbdim        !< Number of dimensions
     real, dimension(6) :: xfineloc
 !
diff --git a/ext/AGRIF/AGRIF_FILES/modupdate.F90 b/ext/AGRIF/AGRIF_FILES/modupdate.F90
index 3bdb5b1b32c32a04011776637e240d6645e4d5c5..c647fcb7930460be33cc07e7cb3b51fdb5f3e63d 100644
--- a/ext/AGRIF/AGRIF_FILES/modupdate.F90
+++ b/ext/AGRIF/AGRIF_FILES/modupdate.F90
@@ -1904,8 +1904,8 @@ subroutine Agrif_Update_1D_Recursive ( type_update,         &
     integer,                            intent(in)  :: type_update            !< Type of update (copy or average)
     integer,                            intent(in)  :: indmin, indmax
     integer,                            intent(in)  :: lb_child, ub_child
-    real(kind=8),                               intent(in)  ::  s_child,  s_parent
-    real(kind=8),                               intent(in)  :: ds_child, ds_parent
+    real(kind=8),                       intent(in)  ::  s_child,  s_parent
+    real(kind=8),                       intent(in)  :: ds_child, ds_parent
     real, dimension(indmin:indmax),     intent(out) :: tempP
     real, dimension(lb_child:ub_child), intent(in)  :: tempC
 !---------------------------------------------------------------------------------------------------
@@ -1933,11 +1933,11 @@ subroutine Agrif_Update_2D_Recursive ( type_update,         &
                                         s_child,  s_parent, &
                                        ds_child, ds_parent )
 !---------------------------------------------------------------------------------------------------
-    integer, dimension(2),          intent(in)  :: type_update            !< Type of update (copy or average)
-    integer, dimension(2),          intent(in)  :: indmin, indmax
-    integer, dimension(2),          intent(in)  :: lb_child, ub_child
-    real(kind=8),    dimension(2),          intent(in)  ::  s_child,  s_parent
-    real(kind=8),    dimension(2),          intent(in)  :: ds_child, ds_parent
+    integer,     dimension(2),      intent(in)  :: type_update            !< Type of update (copy or average)
+    integer,     dimension(2),      intent(in)  :: indmin, indmax
+    integer,     dimension(2),      intent(in)  :: lb_child, ub_child
+    real(kind=8),dimension(2),      intent(in)  ::  s_child,  s_parent
+    real(kind=8),dimension(2),      intent(in)  :: ds_child, ds_parent
     real,    dimension(          &
         indmin(1):indmax(1),     &
         indmin(2):indmax(2)),       intent(out) :: tempP
@@ -2047,11 +2047,11 @@ subroutine Agrif_Update_2D_Recursive_ok ( type_update, &
                                        lb_child, ub_child,                    &
                                        s_child, s_parent, ds_child, ds_parent )
 !---------------------------------------------------------------------------------------------------
-    INTEGER, DIMENSION(2), intent(in)   :: type_update            !< Type of update (copy or average)
-    INTEGER, DIMENSION(2), intent(in)   :: indmin, indmax
-    INTEGER, DIMENSION(2), intent(in)   :: lb_child, ub_child
-    REAL,    DIMENSION(2), intent(in)   :: s_child,  s_parent
-    REAL,    DIMENSION(2), intent(in)   :: ds_child, ds_parent
+    INTEGER,     DIMENSION(2), intent(in)   :: type_update            !< Type of update (copy or average)
+    INTEGER,     DIMENSION(2), intent(in)   :: indmin, indmax
+    INTEGER,     DIMENSION(2), intent(in)   :: lb_child, ub_child
+    REAL(kind=8),DIMENSION(2), intent(in)   :: s_child,  s_parent
+    REAL(kind=8),DIMENSION(2), intent(in)   :: ds_child, ds_parent
     REAL,    DIMENSION(                 &
                 indmin(1):indmax(1),    &
                 indmin(2):indmax(2)),           intent(out) :: tempP
@@ -2101,11 +2101,11 @@ subroutine Agrif_Update_3D_Recursive ( type_update,         &
                                         s_child,  s_parent, &
                                        ds_child, ds_parent )
 !---------------------------------------------------------------------------------------------------
-    integer, dimension(3),          intent(in)  :: type_update            !< Type of update (copy or average)
-    integer, dimension(3),          intent(in)  :: indmin, indmax
-    integer, dimension(3),          intent(in)  :: lb_child, ub_child
-    real,    dimension(3),          intent(in)  ::  s_child,  s_parent
-    real,    dimension(3),          intent(in)  :: ds_child, ds_parent
+    integer,     dimension(3),      intent(in)  :: type_update            !< Type of update (copy or average)
+    integer,     dimension(3),      intent(in)  :: indmin, indmax
+    integer,     dimension(3),      intent(in)  :: lb_child, ub_child
+    real(kind=8),dimension(3),      intent(in)  ::  s_child,  s_parent
+    real(kind=8),dimension(3),      intent(in)  :: ds_child, ds_parent
     real,    dimension(          &
         indmin(1):indmax(1),     &
         indmin(2):indmax(2),     &
@@ -2213,11 +2213,11 @@ subroutine Agrif_Update_4D_Recursive ( type_update,         &
                                         s_child,  s_parent, &
                                        ds_child, ds_parent )
 !---------------------------------------------------------------------------------------------------
-    integer, dimension(4),          intent(in)  :: type_update            !< Type of update (copy or average)
-    integer, dimension(4),          intent(in)  :: indmin, indmax
-    integer, dimension(4),          intent(in)  :: lb_child, ub_child
-    real,    dimension(4),          intent(in)  ::  s_child,  s_parent
-    real,    dimension(4),          intent(in)  :: ds_child, ds_parent
+    integer,     dimension(4),      intent(in)  :: type_update            !< Type of update (copy or average)
+    integer,     dimension(4),      intent(in)  :: indmin, indmax
+    integer,     dimension(4),      intent(in)  :: lb_child, ub_child
+    real(kind=8),dimension(4),      intent(in)  ::  s_child,  s_parent
+    real(kind=8),dimension(4),      intent(in)  :: ds_child, ds_parent
     real,    dimension(          &
         indmin(1):indmax(1),     &
         indmin(2):indmax(2),     &
@@ -2285,11 +2285,11 @@ subroutine Agrif_Update_5D_Recursive ( type_update,         &
                                         s_child,  s_parent, &
                                        ds_child, ds_parent )
 !---------------------------------------------------------------------------------------------------
-    integer, dimension(5),          intent(in)  :: type_update            !< Type of update (copy or average)
-    integer, dimension(5),          intent(in)  :: indmin, indmax
-    integer, dimension(5),          intent(in)  :: lb_child, ub_child
-    real,    dimension(5),          intent(in)  ::  s_child,  s_parent
-    real,    dimension(5),          intent(in)  :: ds_child, ds_parent
+    integer,     dimension(5),      intent(in)  :: type_update            !< Type of update (copy or average)
+    integer,     dimension(5),      intent(in)  :: indmin, indmax
+    integer,     dimension(5),      intent(in)  :: lb_child, ub_child
+    real(kind=8),dimension(5),      intent(in)  ::  s_child,  s_parent
+    real(kind=8),dimension(5),      intent(in)  :: ds_child, ds_parent
     real,    dimension(          &
         indmin(1):indmax(1),     &
         indmin(2):indmax(2),     &
@@ -2364,11 +2364,11 @@ subroutine Agrif_Update_6D_Recursive ( type_update,         &
                                         s_child,  s_parent, &
                                        ds_child, ds_parent )
 !---------------------------------------------------------------------------------------------------
-    integer, dimension(6),          intent(in)  :: type_update            !< Type of update (copy or average)
-    integer, dimension(6),          intent(in)  :: indmin, indmax
-    integer, dimension(6),          intent(in)  :: lb_child, ub_child
-    real,    dimension(6),          intent(in)  ::  s_child,  s_parent
-    real,    dimension(6),          intent(in)  :: ds_child, ds_parent
+    integer,     dimension(6),      intent(in)  :: type_update            !< Type of update (copy or average)
+    integer,     dimension(6),      intent(in)  :: indmin, indmax
+    integer,     dimension(6),      intent(in)  :: lb_child, ub_child
+    real(kind=8),dimension(6),      intent(in)  ::  s_child,  s_parent
+    real(kind=8),dimension(6),      intent(in)  :: ds_child, ds_parent
     real,    dimension(          &
         indmin(1):indmax(1),     &
         indmin(2):indmax(2),     &
@@ -2454,8 +2454,8 @@ subroutine Agrif_UpdateBase ( type_update,              &
     integer,                            intent(in) :: lb_child, ub_child
     real, dimension(indmin:indmax),     intent(out):: parent_tab
     real, dimension(lb_child:ub_child), intent(in) :: child_tab
-    real,                               intent(in) :: s_parent,  s_child
-    real,                               intent(in) :: ds_parent, ds_child
+    real(kind=8),                       intent(in) :: s_parent,  s_child
+    real(kind=8),                       intent(in) :: ds_parent, ds_child
 !---------------------------------------------------------------------------------------------------
     integer :: np       ! Length of parent array
     integer :: nc       ! Length of child  array
diff --git a/ext/AGRIF/AGRIF_FILES/modupdatebasic.F90 b/ext/AGRIF/AGRIF_FILES/modupdatebasic.F90
index 900f0386f97c1dc62b574a4172085baf9759eb9d..76367bcfe79fbba0a5831ad285fd57b26f351332 100644
--- a/ext/AGRIF/AGRIF_FILES/modupdatebasic.F90
+++ b/ext/AGRIF/AGRIF_FILES/modupdatebasic.F90
@@ -231,16 +231,15 @@ subroutine Agrif_basicupdate_max1d ( x, y, np, nc, s_parent, s_child, ds_parent,
     REAL, DIMENSION(np), intent(out)    :: x
     REAL, DIMENSION(nc), intent(in)     :: y
     INTEGER,             intent(in)     :: np,nc
-    REAL,                intent(in)     :: s_parent,  s_child
-    REAL,                intent(in)     :: ds_parent, ds_child
+    REAL(kind=8),        intent(in)     :: s_parent,  s_child
+    REAL(kind=8),        intent(in)     :: ds_parent, ds_child
 !
     INTEGER :: i, ii, locind_child_left, coeffraf
-    REAL    :: xpos, invcoeffraf
+    REAL(kind=8) :: xpos
     INTEGER :: nbnonnuls
     INTEGER :: diffmod
 !
     coeffraf = nint(ds_parent/ds_child)
-    invcoeffraf = 1./coeffraf
 !
     if (coeffraf == 1) then
         locind_child_left = 1 + nint((s_parent - s_child)/ds_child)
diff --git a/ext/PPR/src/bfun1d.h90 b/ext/PPR/src/bfun1d.h90
index 279634678bb4805d2e674e5479a3e09bbaa0b257..1a40a981c6d3a4d69449580957363d1507761767 100644
--- a/ext/PPR/src/bfun1d.h90
+++ b/ext/PPR/src/bfun1d.h90
@@ -59,46 +59,46 @@
     !------------------------------------ -1th-order basis !
             select case (ndof)
             case (+1)
-                bfun(1) = sval**1 / 1.e0
+                bfun(1) = sval**1 / 1.e0_8
                 
             case (+2)
-                bfun(1) = sval**1 / 1.e0
-                bfun(2) = sval**2 / 2.e0
+                bfun(1) = sval**1 / 1.e0_8
+                bfun(2) = sval**2 / 2.e0_8
                 
             case (+3)
-                bfun(1) = sval**1 / 1.e0
-                bfun(2) = sval**2 / 2.e0
-                bfun(3) = sval**3 / 3.e0
+                bfun(1) = sval**1 / 1.e0_8
+                bfun(2) = sval**2 / 2.e0_8
+                bfun(3) = sval**3 / 3.e0_8
                 
             case (+4)
-                bfun(1) = sval**1 / 1.e0
-                bfun(2) = sval**2 / 2.e0
-                bfun(3) = sval**3 / 3.e0
-                bfun(4) = sval**4 / 4.e0
+                bfun(1) = sval**1 / 1.e0_8
+                bfun(2) = sval**2 / 2.e0_8
+                bfun(3) = sval**3 / 3.e0_8
+                bfun(4) = sval**4 / 4.e0_8
                 
             case (+5)
-                bfun(1) = sval**1 / 1.e0
-                bfun(2) = sval**2 / 2.e0
-                bfun(3) = sval**3 / 3.e0
-                bfun(4) = sval**4 / 4.e0
-                bfun(5) = sval**5 / 5.e0
+                bfun(1) = sval**1 / 1.e0_8
+                bfun(2) = sval**2 / 2.e0_8
+                bfun(3) = sval**3 / 3.e0_8
+                bfun(4) = sval**4 / 4.e0_8
+                bfun(5) = sval**5 / 5.e0_8
                 
             case (+6)
-                bfun(1) = sval**1 / 1.e0
-                bfun(2) = sval**2 / 2.e0
-                bfun(3) = sval**3 / 3.e0
-                bfun(4) = sval**4 / 4.e0
-                bfun(5) = sval**5 / 5.e0
-                bfun(6) = sval**6 / 6.e0
+                bfun(1) = sval**1 / 1.e0_8
+                bfun(2) = sval**2 / 2.e0_8
+                bfun(3) = sval**3 / 3.e0_8
+                bfun(4) = sval**4 / 4.e0_8
+                bfun(5) = sval**5 / 5.e0_8
+                bfun(6) = sval**6 / 6.e0_8
                 
             case (+7)
-                bfun(1) = sval**1 / 1.e0
-                bfun(2) = sval**2 / 2.e0
-                bfun(3) = sval**3 / 3.e0
-                bfun(4) = sval**4 / 4.e0
-                bfun(5) = sval**5 / 5.e0
-                bfun(6) = sval**6 / 6.e0
-                bfun(7) = sval**7 / 7.e0
+                bfun(1) = sval**1 / 1.e0_8
+                bfun(2) = sval**2 / 2.e0_8
+                bfun(3) = sval**3 / 3.e0_8
+                bfun(4) = sval**4 / 4.e0_8
+                bfun(5) = sval**5 / 5.e0_8
+                bfun(6) = sval**6 / 6.e0_8
+                bfun(7) = sval**7 / 7.e0_8
             
             end select
 
@@ -106,46 +106,46 @@
     !------------------------------------ +0th-order basis !
             select case (ndof)
             case (+1)
-                bfun(1) =           1.e0
+                bfun(1) =           1.e0_8
                 
             case (+2)
-                bfun(1) =           1.e0
-                bfun(2) = sval**1 * 1.e0
+                bfun(1) =           1.e0_8
+                bfun(2) = sval**1 * 1.e0_8
                 
             case (+3)
-                bfun(1) =           1.e0
-                bfun(2) = sval**1 * 1.e0
-                bfun(3) = sval**2 * 1.e0
+                bfun(1) =           1.e0_8
+                bfun(2) = sval**1 * 1.e0_8
+                bfun(3) = sval**2 * 1.e0_8
                 
             case (+4)
-                bfun(1) =           1.e0
-                bfun(2) = sval**1 * 1.e0
-                bfun(3) = sval**2 * 1.e0
-                bfun(4) = sval**3 * 1.e0
+                bfun(1) =           1.e0_8
+                bfun(2) = sval**1 * 1.e0_8
+                bfun(3) = sval**2 * 1.e0_8
+                bfun(4) = sval**3 * 1.e0_8
                 
             case (+5)
-                bfun(1) =           1.e0
-                bfun(2) = sval**1 * 1.e0
-                bfun(3) = sval**2 * 1.e0
-                bfun(4) = sval**3 * 1.e0
-                bfun(5) = sval**4 * 1.e0
+                bfun(1) =           1.e0_8
+                bfun(2) = sval**1 * 1.e0_8
+                bfun(3) = sval**2 * 1.e0_8
+                bfun(4) = sval**3 * 1.e0_8
+                bfun(5) = sval**4 * 1.e0_8
                 
             case (+6)
-                bfun(1) =           1.e0
-                bfun(2) = sval**1 * 1.e0
-                bfun(3) = sval**2 * 1.e0
-                bfun(4) = sval**3 * 1.e0
-                bfun(5) = sval**4 * 1.e0
-                bfun(6) = sval**5 * 1.e0
+                bfun(1) =           1.e0_8
+                bfun(2) = sval**1 * 1.e0_8
+                bfun(3) = sval**2 * 1.e0_8
+                bfun(4) = sval**3 * 1.e0_8
+                bfun(5) = sval**4 * 1.e0_8
+                bfun(6) = sval**5 * 1.e0_8
                 
             case (+7)
-                bfun(1) =           1.e0
-                bfun(2) = sval**1 * 1.e0
-                bfun(3) = sval**2 * 1.e0
-                bfun(4) = sval**3 * 1.e0
-                bfun(5) = sval**4 * 1.e0
-                bfun(6) = sval**5 * 1.e0
-                bfun(7) = sval**6 * 1.e0
+                bfun(1) =           1.e0_8
+                bfun(2) = sval**1 * 1.e0_8
+                bfun(3) = sval**2 * 1.e0_8
+                bfun(4) = sval**3 * 1.e0_8
+                bfun(5) = sval**4 * 1.e0_8
+                bfun(6) = sval**5 * 1.e0_8
+                bfun(7) = sval**6 * 1.e0_8
             
             end select
 
@@ -153,46 +153,46 @@
     !------------------------------------ +1st-order basis !
             select case (ndof)
             case (+1)
-                bfun(1) =           0.e0
+                bfun(1) =           0.e0_8
                 
             case (+2)
-                bfun(1) =           0.e0
-                bfun(2) =           1.e0
+                bfun(1) =           0.e0_8
+                bfun(2) =           1.e0_8
                 
             case (+3)
-                bfun(1) =           0.e0
-                bfun(2) =           1.e0
-                bfun(3) = sval**1 * 2.e0
+                bfun(1) =           0.e0_8
+                bfun(2) =           1.e0_8
+                bfun(3) = sval**1 * 2.e0_8
                 
             case (+4)
-                bfun(1) =           0.e0
-                bfun(2) =           1.e0
-                bfun(3) = sval**1 * 2.e0
-                bfun(4) = sval**2 * 3.e0
+                bfun(1) =           0.e0_8
+                bfun(2) =           1.e0_8
+                bfun(3) = sval**1 * 2.e0_8
+                bfun(4) = sval**2 * 3.e0_8
                 
             case (+5)
-                bfun(1) =           0.e0
-                bfun(2) =           1.e0
-                bfun(3) = sval**1 * 2.e0
-                bfun(4) = sval**2 * 3.e0
-                bfun(5) = sval**3 * 4.e0
+                bfun(1) =           0.e0_8
+                bfun(2) =           1.e0_8
+                bfun(3) = sval**1 * 2.e0_8
+                bfun(4) = sval**2 * 3.e0_8
+                bfun(5) = sval**3 * 4.e0_8
                 
             case (+6)
-                bfun(1) =           0.e0
-                bfun(2) =           1.e0
-                bfun(3) = sval**1 * 2.e0
-                bfun(4) = sval**2 * 3.e0
-                bfun(5) = sval**3 * 4.e0
-                bfun(6) = sval**4 * 5.e0
+                bfun(1) =           0.e0_8
+                bfun(2) =           1.e0_8
+                bfun(3) = sval**1 * 2.e0_8
+                bfun(4) = sval**2 * 3.e0_8
+                bfun(5) = sval**3 * 4.e0_8
+                bfun(6) = sval**4 * 5.e0_8
                 
             case (+7)
-                bfun(1) =           0.e0
-                bfun(2) =           1.e0
-                bfun(3) = sval**1 * 2.e0
-                bfun(4) = sval**2 * 3.e0
-                bfun(5) = sval**3 * 4.e0
-                bfun(6) = sval**4 * 5.e0
-                bfun(7) = sval**5 * 6.e0
+                bfun(1) =           0.e0_8
+                bfun(2) =           1.e0_8
+                bfun(3) = sval**1 * 2.e0_8
+                bfun(4) = sval**2 * 3.e0_8
+                bfun(5) = sval**3 * 4.e0_8
+                bfun(6) = sval**4 * 5.e0_8
+                bfun(7) = sval**5 * 6.e0_8
             
             end select
 
@@ -200,46 +200,46 @@
     !------------------------------------ +2nd-order basis !
             select case (ndof)
             case (+1)
-                bfun(1) =           0.e0
+                bfun(1) =           0.e0_8
                 
             case (+2)
-                bfun(1) =           0.e0
-                bfun(2) =           0.e0
+                bfun(1) =           0.e0_8
+                bfun(2) =           0.e0_8
                 
             case (+3)
-                bfun(1) =           0.e0
-                bfun(2) =           0.e0
-                bfun(3) =           2.e0
+                bfun(1) =           0.e0_8
+                bfun(2) =           0.e0_8
+                bfun(3) =           2.e0_8
                 
             case (+4)
-                bfun(1) =           0.e0
-                bfun(2) =           0.e0
-                bfun(3) =           2.e0
-                bfun(4) = sval**1 * 6.e0
+                bfun(1) =           0.e0_8
+                bfun(2) =           0.e0_8
+                bfun(3) =           2.e0_8
+                bfun(4) = sval**1 * 6.e0_8
                 
             case (+5)
-                bfun(1) =           0.e0
-                bfun(2) =           0.e0
-                bfun(3) =           2.e0
-                bfun(4) = sval**1 * 6.e0
-                bfun(5) = sval**2 *12.e0
+                bfun(1) =           0.e0_8
+                bfun(2) =           0.e0_8
+                bfun(3) =           2.e0_8
+                bfun(4) = sval**1 * 6.e0_8
+                bfun(5) = sval**2 *12.e0_8
                 
             case (+6)
-                bfun(1) =           0.e0
-                bfun(2) =           0.e0
-                bfun(3) =           2.e0
-                bfun(4) = sval**1 * 6.e0
-                bfun(5) = sval**2 *12.e0
-                bfun(6) = sval**3 *20.e0
+                bfun(1) =           0.e0_8
+                bfun(2) =           0.e0_8
+                bfun(3) =           2.e0_8
+                bfun(4) = sval**1 * 6.e0_8
+                bfun(5) = sval**2 *12.e0_8
+                bfun(6) = sval**3 *20.e0_8
                 
             case (+7)
-                bfun(1) =           0.e0
-                bfun(2) =           0.e0
-                bfun(3) =           2.e0
-                bfun(4) = sval**1 * 6.e0
-                bfun(5) = sval**2 *12.e0
-                bfun(6) = sval**3 *20.e0
-                bfun(7) = sval**4 *30.e0
+                bfun(1) =           0.e0_8
+                bfun(2) =           0.e0_8
+                bfun(3) =           2.e0_8
+                bfun(4) = sval**1 * 6.e0_8
+                bfun(5) = sval**2 *12.e0_8
+                bfun(6) = sval**3 *20.e0_8
+                bfun(7) = sval**4 *30.e0_8
             
             end select
 
diff --git a/ext/PPR/src/ffsl1d.h90 b/ext/PPR/src/ffsl1d.h90
index c2a3d9fa438ce31622a1c4cf17391e5f4ee350b7..45b140740af61bcfc6180e43774b324d121e93b0 100644
--- a/ext/PPR/src/ffsl1d.h90
+++ b/ext/PPR/src/ffsl1d.h90
@@ -250,7 +250,7 @@
         
         do  ipos = +2 , npos - 1
         
-            if (uvel(ipos) .gt. +0.e0) then
+            if (uvel(ipos) .gt. +0.e0_8) then
             
     !----------- integrate profile over upwind cell IPOS-1 !
  
@@ -264,8 +264,8 @@
         &        * tDEL / SPAC(    +1)
             end if
  
-            ss11 = +1.e0 - 2.e0 * uCFL
-            ss22 = +1.e0
+            ss11 = +1.e0_8 - 2.e0_8 * uCFL
+            ss22 = +1.e0_8
 
             call bfun1d(-1,mdof,ss11,vv11)
             call bfun1d(-1,mdof,ss22,vv22)
@@ -285,7 +285,7 @@
             end do
                       
             else &
-        &   if (uvel(ipos) .lt. -0.e0) then
+        &   if (uvel(ipos) .lt. -0.e0_8) then
             
     !----------- integrate profile over upwind cell IPOS+0 !
             
@@ -299,8 +299,8 @@
         &        * tDEL / SPAC(    +1)
             end if
 
-            ss11 = -1.e0 - 2.e0 * uCFL
-            ss22 = -1.e0
+            ss11 = -1.e0_8 - 2.e0_8 * uCFL
+            ss22 = -1.e0_8
         
             call bfun1d(-1,mdof,ss11,vv11)
             call bfun1d(-1,mdof,ss22,vv22)
diff --git a/ext/PPR/src/inv.h90 b/ext/PPR/src/inv.h90
index daceabb9a808cc6348e29c44604b21daedda93e4..491a1192322ec33479ce35ec93e9f73aa8f7f2f8 100644
--- a/ext/PPR/src/inv.h90
+++ b/ext/PPR/src/inv.h90
@@ -139,7 +139,7 @@
 
     !-------------------------------- C = C + scal * A * B !
 
-        if (scal .eq. +1.e0) then
+        if (scal .eq. +1.e0_8) then
 
         cmat(1,1) = cmat(1,1) & 
              + ( amat(1,1) * bmat(1,1) &
@@ -156,7 +156,7 @@
                + amat(2,2) * bmat(2,2) )
 
         else &
-        if (scal .eq. -1.e0) then
+        if (scal .eq. -1.e0_8) then
 
         cmat(1,1) = cmat(1,1) & 
              - ( amat(1,1) * bmat(1,1) &
@@ -210,7 +210,7 @@
 
     !-------------------------------- C = C + scal * A * B !
 
-        if (scal .eq. +1.e0) then
+        if (scal .eq. +1.e0_8) then
 
         cmat(1,1) = cmat(1,1) &
              + ( amat(1,1) * bmat(1,1) &
@@ -252,7 +252,7 @@
                + amat(3,3) * bmat(3,3) )
 
         else &
-        if (scal .eq. -1.e0) then
+        if (scal .eq. -1.e0_8) then
 
         cmat(1,1) = cmat(1,1) &
              - ( amat(1,1) * bmat(1,1) &
@@ -499,23 +499,23 @@
 
     !---------------------------------------- L = C * A^-1 !
 
-        lmat(1,1) = +0.e0
-        lmat(1,2) = +0.e0
-        lmat(2,1) = +0.e0
-        lmat(2,2) = +0.e0
+        lmat(1,1) = +0.e0_8
+        lmat(1,2) = +0.e0_8
+        lmat(2,1) = +0.e0_8
+        lmat(2,2) = +0.e0_8
 
         call mul_2x2(amat(3,1),adim,ainv,LDIM, &
-                    +1.e0,lmat,LDIM)
+                    +1.e0_8,lmat,LDIM)
         
     !---------------------------------------- U = A^-1 * B !       
 
-        umat(1,1) = +0.e0
-        umat(1,2) = +0.e0
-        umat(2,1) = +0.e0
-        umat(2,2) = +0.e0
+        umat(1,1) = +0.e0_8
+        umat(1,2) = +0.e0_8
+        umat(2,1) = +0.e0_8
+        umat(2,2) = +0.e0_8
  
         call mul_2x2(ainv,LDIM,amat(1,3),adim, &
-                    +1.e0,umat,LDIM)
+                    +1.e0_8,umat,LDIM)
         
     !-------------------------------- S = D - C * A^-1 * B !
 
@@ -525,7 +525,7 @@
         smat(2,2) = amat(4,4)
 
         call mul_2x2(lmat,LDIM,amat(1,3),adim, &
-                    -1.e0/adet,smat,LDIM)
+                    -1.e0_8/adet,smat,LDIM)
 
         call inv_2x2(smat,LDIM,sinv,LDIM,sdet)
 
@@ -640,33 +640,33 @@
 
     !---------------------------------------- L = C * A^-1 !
 
-        lmat(1,1) = +0.e0
-        lmat(1,2) = +0.e0
-        lmat(1,3) = +0.e0
-        lmat(2,1) = +0.e0
-        lmat(2,2) = +0.e0
-        lmat(2,3) = +0.e0
-        lmat(3,1) = +0.e0
-        lmat(3,2) = +0.e0
-        lmat(3,3) = +0.e0
+        lmat(1,1) = +0.e0_8
+        lmat(1,2) = +0.e0_8
+        lmat(1,3) = +0.e0_8
+        lmat(2,1) = +0.e0_8
+        lmat(2,2) = +0.e0_8
+        lmat(2,3) = +0.e0_8
+        lmat(3,1) = +0.e0_8
+        lmat(3,2) = +0.e0_8
+        lmat(3,3) = +0.e0_8
 
         call mul_3x3(amat(4,1),adim,ainv,LDIM, &
-                    +1.e0,lmat,LDIM)
+                    +1.e0_8,lmat,LDIM)
         
     !---------------------------------------- U = A^-1 * B !       
 
-        umat(1,1) = +0.e0
-        umat(1,2) = +0.e0
-        umat(1,3) = +0.e0
-        umat(2,1) = +0.e0
-        umat(2,2) = +0.e0
-        umat(2,3) = +0.e0
-        umat(3,1) = +0.e0
-        umat(3,2) = +0.e0
-        umat(3,3) = +0.e0
+        umat(1,1) = +0.e0_8
+        umat(1,2) = +0.e0_8
+        umat(1,3) = +0.e0_8
+        umat(2,1) = +0.e0_8
+        umat(2,2) = +0.e0_8
+        umat(2,3) = +0.e0_8
+        umat(3,1) = +0.e0_8
+        umat(3,2) = +0.e0_8
+        umat(3,3) = +0.e0_8
  
         call mul_3x3(ainv,LDIM,amat(1,4),adim, &
-                    +1.e0,umat,LDIM)
+                    +1.e0_8,umat,LDIM)
         
     !-------------------------------- S = D - C * A^-1 * B !
 
@@ -681,7 +681,7 @@
         smat(3,3) = amat(6,6)
 
         call mul_3x3(lmat,LDIM,amat(1,4),adim, &
-                    -1.e0/adet,smat,LDIM)
+                    -1.e0_8/adet,smat,LDIM)
 
         call inv_3x3(smat,LDIM,sinv,LDIM,sdet)
 
diff --git a/ext/PPR/src/oscl1d.h90 b/ext/PPR/src/oscl1d.h90
index 30f50982f55c0bb862f73334e0219a6102f4621a..178365709b2fe1245ff5ee2055fff8f3dd79802b 100644
--- a/ext/PPR/src/oscl1d.h90
+++ b/ext/PPR/src/oscl1d.h90
@@ -66,8 +66,8 @@
     !------------------------------- at least 3 grid-cells !
         do  ipos = +1, npos-1
         do  ivar = +1, nvar-0
-            oscl(1,ivar,ipos) = +0.e0
-            oscl(2,ivar,ipos) = +0.e0
+            oscl(1,ivar,ipos) = +0.e0_8
+            oscl(2,ivar,ipos) = +0.e0_8
         end do
         end do
         end if
@@ -143,25 +143,25 @@
         hhlc = hhll + hhcc
         hhmm = hhll + hhcc + hhrr
 
-        cmat(1,1) = -(hhcc+2.e0*hhrr)/(hhlc*hhmm)
+        cmat(1,1) = -(hhcc+2.e0_8*hhrr)/(hhlc*hhmm)
         cmat(1,2) = -(hhll-hhrr)* &
-        &          (3.e0*hhcc+2.e0*(hhll+hhrr))/&
+        &          (3.e0_8*hhcc+2.e0_8*(hhll+hhrr))/&
         &            (hhlc*hhrc*hhmm)
-        cmat(1,3) = +(hhcc+2.e0*hhll)/(hhrc*hhmm)
+        cmat(1,3) = +(hhcc+2.e0_8*hhll)/(hhrc*hhmm)
 
-        cmat(2,1) = +3.e0/(hhlc*hhmm)
-        cmat(2,2) = -3.e0*(2.e0*hhcc+hhll+hhrr)/&
+        cmat(2,1) = +3.e0_8/(hhlc*hhmm)
+        cmat(2,2) = -3.e0_8*(2.e0_8*hhcc+hhll+hhrr)/&
         &            (hhlc*hhrc*hhmm)
-        cmat(2,3) = +3.e0/(hhrc*hhmm)
+        cmat(2,3) = +3.e0_8/(hhrc*hhmm)
 
         do  ivar = 1, nvar
 
-            oscl(1,ivar,ipos) = +1.e0 * ( &
+            oscl(1,ivar,ipos) = +1.e0_8 * ( &
         & + cmat(1,1)*fdat(1,ivar,ipos-1) &
         & + cmat(1,2)*fdat(1,ivar,ipos+0) &
         & + cmat(1,3)*fdat(1,ivar,ipos+1) )
 
-            oscl(2,ivar,ipos) = +2.e0 * ( &
+            oscl(2,ivar,ipos) = +2.e0_8 * ( &
         & + cmat(2,1)*fdat(1,ivar,ipos-1) &
         & + cmat(2,2)*fdat(1,ivar,ipos+0) &
         & + cmat(2,3)*fdat(1,ivar,ipos+1) )
@@ -176,8 +176,8 @@
         hhcc = max(delx(head+1),dmin)
         hhrr = max(delx(head+2),dmin)
 
-        cmat(1,1) = -2.e0 / (hhll+hhcc)
-        cmat(1,2) = +2.e0 / (hhll+hhcc)
+        cmat(1,1) = -2.e0_8 / (hhll+hhcc)
+        cmat(1,2) = +2.e0_8 / (hhll+hhcc)
 
         do  ivar = 1, nvar
 
@@ -185,7 +185,7 @@
         & + cmat(1,1)*fdat(1,ivar,head+0) &
         & + cmat(1,2)*fdat(1,ivar,head+1)
              
-            oscl(2,ivar,head) = +0.e0
+            oscl(2,ivar,head) = +0.e0_8
 
         end do
 
@@ -195,8 +195,8 @@
         hhcc = max(delx(tail-1),dmin)
         hhrr = max(delx(tail-0),dmin)
 
-        cmat(1,2) = -2.e0 / (hhrr+hhcc)
-        cmat(1,3) = +2.e0 / (hhrr+hhcc)
+        cmat(1,2) = -2.e0_8 / (hhrr+hhcc)
+        cmat(1,3) = +2.e0_8 / (hhrr+hhcc)
 
         do  ivar = 1, nvar
 
@@ -204,7 +204,7 @@
         & + cmat(1,2)*fdat(1,ivar,tail-1) &
         & + cmat(1,3)*fdat(1,ivar,tail+0)
 
-            oscl(2,ivar,tail) = +0.e0
+            oscl(2,ivar,tail) = +0.e0_8
 
         end do
         
@@ -268,7 +268,7 @@
         & + .50d+0 * fdat(1,ivar,head+1) &
         & - .50d+0 * fdat(1,ivar,head+0)
         
-            oscl(2,ivar,head) = +0.e0
+            oscl(2,ivar,head) = +0.e0_8
 
         end do
 
@@ -280,7 +280,7 @@
         & + .50d+0 * fdat(1,ivar,tail+0) &
         & - .50d+0 * fdat(1,ivar,tail-1)
             
-            oscl(2,ivar,tail) = +0.e0
+            oscl(2,ivar,tail) = +0.e0_8
 
         end do
         
diff --git a/ext/PPR/src/p1e.h90 b/ext/PPR/src/p1e.h90
index 39b639357f49aceba928005a09220803a9a00830..467c3507f633ecdcadbe178ac532a353a7fad249 100644
--- a/ext/PPR/src/p1e.h90
+++ b/ext/PPR/src/p1e.h90
@@ -83,10 +83,10 @@
         do  ivar = 1,nvar
         
             edge(ivar,1) = fdat(1,ivar,1)
-            dfdx(ivar,1) = 0.e0
+            dfdx(ivar,1) = 0.e0_8
             
             edge(ivar,2) = fdat(1,ivar,1)
-            dfdx(ivar,2) = 0.e0
+            dfdx(ivar,2) = 0.e0_8
             
         end do
         end if
@@ -103,7 +103,7 @@
     
     !--------------- reconstruction: constant grid-spacing !
             
-            dd10 = delx(+1) * 2.e0
+            dd10 = delx(+1) * 2.e0_8
             
             do  ivar = +1, nvar
 
@@ -175,8 +175,8 @@
             edge(ivar,tail+1) = &
         &       fdat(+1,ivar,tail+0)
         
-            dfdx(ivar,head-1) = 0.e0
-            dfdx(ivar,tail+1) = 0.e0
+            dfdx(ivar,head-1) = 0.e0_8
+            dfdx(ivar,tail+1) = 0.e0_8
 
         end do
     
diff --git a/ext/PPR/src/p3e.h90 b/ext/PPR/src/p3e.h90
index 6becd63358c1a4c9d1538386b65f59e91d164aea..a651873433e6cfbbc90e9242501c086b4af7bb57 100644
--- a/ext/PPR/src/p3e.h90
+++ b/ext/PPR/src/p3e.h90
@@ -120,24 +120,24 @@
             do  ivar = 1, nvar
 
                 edge(ivar,ipos) = ( &
-        &     -      1.e0 * &
+        &     -      1.e0_8 * &
         &       fdat(1,ivar,ipos-2) &
-        &     +      7.e0 * &
+        &     +      7.e0_8 * &
         &       fdat(1,ivar,ipos-1) &
-        &     +      7.e0 * &
+        &     +      7.e0_8 * &
         &       fdat(1,ivar,ipos+0) &
-        &     -      1.e0 * &
-        &       fdat(1,ivar,ipos+1) ) / 12.e0
+        &     -      1.e0_8 * &
+        &       fdat(1,ivar,ipos+1) ) / 12.e0_8
             
                 dfdx(ivar,ipos) = ( &
-        &     +      1.e0 * &
+        &     +      1.e0_8 * &
         &       fdat(1,ivar,ipos-2) &
-        &     -     15.e0 * &
+        &     -     15.e0_8 * &
         &       fdat(1,ivar,ipos-1) &
-        &     +     15.e0 * &
+        &     +     15.e0_8 * &
         &       fdat(1,ivar,ipos+0) &
-        &     -      1.e0 * &
-        &       fdat(1,ivar,ipos+1) ) / 12.e0
+        &     -      1.e0_8 * &
+        &       fdat(1,ivar,ipos+1) ) / 12.e0_8
 
                 dfdx(ivar,ipos) = &
         &       dfdx(ivar,ipos) / delx(+1)
@@ -165,7 +165,7 @@
             xmap(-2) = -( delh(-2) &
         &              +  delh(-1) ) / xhat
             xmap(-1) = -  delh(-1)   / xhat
-            xmap(+0) = +  0.e0
+            xmap(+0) = +  0.e0_8
             xmap(+1) = +  delh(+0)   / xhat
             xmap(+2) = +( delh(+0) &
         &              +  delh(+1) ) / xhat
diff --git a/ext/PPR/src/p5e.h90 b/ext/PPR/src/p5e.h90
index fe4b446f7c1beec2e8a120c8902510a805d3c421..1e3f65fe93b366ec1d75b1948232e19474bc66a0 100644
--- a/ext/PPR/src/p5e.h90
+++ b/ext/PPR/src/p5e.h90
@@ -120,31 +120,31 @@
             do  ivar = 1, nvar
 
                 edge(ivar,ipos) = &
-        &     + ( 1.e0 / 60.e0) * & 
+        &     + ( 1.e0_8 / 60.e0_8) * & 
         &       fdat(1,ivar,ipos-3) &
-        &     - ( 8.e0 / 60.e0) * &
+        &     - ( 8.e0_8 / 60.e0_8) * &
         &       fdat(1,ivar,ipos-2) &
-        &     + (37.e0 / 60.e0) * &
+        &     + (37.e0_8 / 60.e0_8) * &
         &       fdat(1,ivar,ipos-1) &
-        &     + (37.e0 / 60.e0) * &
+        &     + (37.e0_8 / 60.e0_8) * &
         &       fdat(1,ivar,ipos+0) &
-        &     - ( 8.e0 / 60.e0) * &
+        &     - ( 8.e0_8 / 60.e0_8) * &
         &       fdat(1,ivar,ipos+1) &
-        &     + ( 1.e0 / 60.e0) * &
+        &     + ( 1.e0_8 / 60.e0_8) * &
         &       fdat(1,ivar,ipos+2)
             
                 dfdx(ivar,ipos) = &
-        &     - ( 1.e0 / 90.e0) * & 
+        &     - ( 1.e0_8 / 90.e0_8) * & 
         &       fdat(1,ivar,ipos-3) &
-        &     + ( 5.e0 / 36.e0) * &
+        &     + ( 5.e0_8 / 36.e0_8) * &
         &       fdat(1,ivar,ipos-2) &
-        &     - (49.e0 / 36.e0) * &
+        &     - (49.e0_8 / 36.e0_8) * &
         &       fdat(1,ivar,ipos-1) &
-        &     + (49.e0 / 36.e0) * &
+        &     + (49.e0_8 / 36.e0_8) * &
         &       fdat(1,ivar,ipos+0) &
-        &     - ( 5.e0 / 36.e0) * &
+        &     - ( 5.e0_8 / 36.e0_8) * &
         &       fdat(1,ivar,ipos+1) &
-        &     + ( 1.e0 / 90.e0) * &
+        &     + ( 1.e0_8 / 90.e0_8) * &
         &       fdat(1,ivar,ipos+2)
 
                 dfdx(ivar,ipos) = &
@@ -184,7 +184,7 @@
             xmap(-2) = -( delh(-2) &
         &              +  delh(-1) ) / xhat
             xmap(-1) = -  delh(-1)   / xhat
-            xmap(+0) = +  0.e0
+            xmap(+0) = +  0.e0_8
             xmap(+1) = +  delh(+0)   / xhat
             xmap(+2) = +( delh(+0) &
         &              +  delh(+1) ) / xhat
diff --git a/ext/PPR/src/pbc.h90 b/ext/PPR/src/pbc.h90
index 1d09e3a9bc810bec8824ca176c579f7117e94339..13debd46f9e7adc0162e57b36b5c8e920f473370 100644
--- a/ext/PPR/src/pbc.h90
+++ b/ext/PPR/src/pbc.h90
@@ -244,8 +244,8 @@
 
         xmap(-1) =-(delh(-1) + &
         &           delh(+0)*0.5d0)/xhat
-        xmap(+0) = -1.e0
-        xmap(+1) = +1.e0
+        xmap(+0) = -1.e0_8
+        xmap(+1) = +1.e0_8
         xmap(+2) = (delh(+1) + &
         &           delh(+0)*0.5d0)/xhat
 
@@ -453,7 +453,7 @@
         if (bcon(ivar)%bcopt.eq.bopt)  then
 
             eval(-1) = &
-        &   fdat(1,ivar,head-1) * 1.e0
+        &   fdat(1,ivar,head-1) * 1.e0_8
             eval(+0) = &
         &   fdat(1,ivar,head-1) * .5d0 + &
         &   fdat(1,ivar,head+0) * .5d0
@@ -575,8 +575,8 @@
 
         xmap(-1) =-(delh(-1) + &
         &           delh(+0)*0.5d0)/xhat
-        xmap(+0) = -1.e0
-        xmap(+1) = +1.e0
+        xmap(+0) = -1.e0_8
+        xmap(+1) = +1.e0_8
         xmap(+2) = (delh(+1) + &
         &           delh(+0)*0.5d0)/xhat
 
@@ -790,7 +790,7 @@
         &   fdat(1,ivar,tail+0) * .5d0 + &
         &   fdat(1,ivar,tail+1) * .5d0
             eval(+2) = &
-        &   fdat(1,ivar,tail+1) * 1.e0
+        &   fdat(1,ivar,tail+1) * 1.e0_8
    
             gval(+0) = &
         &   fdat(1,ivar,tail+0) * .5d0 - &
diff --git a/ext/PPR/src/ppm.h90 b/ext/PPR/src/ppm.h90
index f8158111cfd1c1e9b6a285e4ed357fe433e51124..923c80f956f0e3f296a11efca8380d84397127eb 100644
--- a/ext/PPR/src/ppm.h90
+++ b/ext/PPR/src/ppm.h90
@@ -94,8 +94,8 @@
         do  ivar = +1, nvar
             fhat(1,ivar,+1) = &
         &   fdat(1,ivar,+1)
-            fhat(2,ivar,+1) = 0.e0
-            fhat(3,ivar,+1) = 0.e0
+            fhat(2,ivar,+1) = 0.e0_8
+            fhat(3,ivar,+1) = 0.e0_8
         end do
         end if
 
@@ -275,8 +275,8 @@
             mono = +1
 
             lhat(1) = ff00
-            lhat(2) = 0.e0
-            lhat(3) = 0.e0
+            lhat(2) = 0.e0_8
+            lhat(3) = 0.e0_8
              
             return
               
diff --git a/ext/PPR/src/pqm.h90 b/ext/PPR/src/pqm.h90
index c834c3daa5886871c672841bd99239ce62c09ce2..19db6503c361f0e99950fb726fae12fa237cec60 100644
--- a/ext/PPR/src/pqm.h90
+++ b/ext/PPR/src/pqm.h90
@@ -99,10 +99,10 @@
     !----- default to reduced order if insufficient points !
         do  ivar = +1, nvar
             fhat(1,ivar,+1) = fdat(1,ivar,+1)
-            fhat(2,ivar,+1) = 0.e0
-            fhat(3,ivar,+1) = 0.e0
-            fhat(4,ivar,+1) = 0.e0
-            fhat(5,ivar,+1) = 0.e0
+            fhat(2,ivar,+1) = 0.e0_8
+            fhat(3,ivar,+1) = 0.e0_8
+            fhat(4,ivar,+1) = 0.e0_8
+            fhat(5,ivar,+1) = 0.e0_8
         end do
         end if
 
@@ -313,10 +313,10 @@
             mono = +1
 
             lhat(1) = ff00
-            lhat(2) = 0.e0
-            lhat(3) = 0.e0
-            lhat(4) = 0.e0
-            lhat(5) = 0.e0
+            lhat(2) = 0.e0_8
+            lhat(3) = 0.e0_8
+            lhat(4) = 0.e0_8
+            lhat(5) = 0.e0_8
              
             return
               
diff --git a/ext/PPR/src/rmap1d.h90 b/ext/PPR/src/rmap1d.h90
index fa3486e8b022199df66f9d71fb53aa52e3515c41..dc10cd74ee1530d390157a00135a916d2076ef52 100644
--- a/ext/PPR/src/rmap1d.h90
+++ b/ext/PPR/src/rmap1d.h90
@@ -244,8 +244,8 @@
 
     !------------------------------------- initializations !
 
-        vvlo(+1:+5) = 0.e0
-        vvhi(+1:+5) = 0.e0
+        vvlo(+1:+5) = 0.e0_8
+        vvhi(+1:+5) = 0.e0_8
 
     !------------- remap FDAT from XPOS to XNEW using FHAT !
 
@@ -285,7 +285,7 @@
             do  idof = +1,ndof
             do  ivar = +1,nvar
                 
-                fnew(idof,ivar,kpos) = 0.e0
+                fnew(idof,ivar,kpos) = 0.e0_8
             
             end do
             end do
@@ -374,7 +374,7 @@
     !   wrt. column-wise conservation. Use KBN approach to
     !   account for FP roundoff.
 
-        sdat = 0.e0; serr = 0.e0
+        sdat = 0.e0_8; serr = 0.e0_8
         do  ipos = +1, npos-1
         do  ivar = +1, nvar-0
         
@@ -407,7 +407,7 @@
 
         sdat =  sdat + serr
 
-        snew = 0.e0; serr = 0.e0
+        snew = 0.e0_8; serr = 0.e0_8
         do  ipos = +1, nnew-1
         do  ivar = +1, nvar-0
         
@@ -450,7 +450,7 @@
     
         do  ivar = +1, nvar-0
 
-            if (serr(ivar) .gt. 0.e0) then
+            if (serr(ivar) .gt. 0.e0_8) then
 
                 vmin = kmin(ivar)
 
@@ -459,7 +459,7 @@
         &  serr(ivar)/(xnew(vmin+1)-xnew(vmin+0))
 
             else &
-        &   if (serr(ivar) .lt. 0.e0) then
+        &   if (serr(ivar) .lt. 0.e0_8) then
 
                 vmax = kmax(ivar)
 
diff --git a/ext/PPR/src/util1d.h90 b/ext/PPR/src/util1d.h90
index 6b4353ba76ee5c76aaf0541e8f191770ba0ecfd1..10cd4c58c573e501145b647580300ee00ff85160 100644
--- a/ext/PPR/src/util1d.h90
+++ b/ext/PPR/src/util1d.h90
@@ -110,7 +110,7 @@
 
             call random_number (rand)
 
-            rand = 2.e0 * (rand-.5d0)
+            rand = 2.e0_8 * (rand-.5d0)
             
             move = rand *  move
 
diff --git a/ext/PPR/src/weno1d.h90 b/ext/PPR/src/weno1d.h90
index d74c76faa23d142eba27ced7929cb5dc4e017c46..e65272ee7e9fd89709fc53cf1d0cb05fa94e1ca7 100644
--- a/ext/PPR/src/weno1d.h90
+++ b/ext/PPR/src/weno1d.h90
@@ -111,7 +111,7 @@
         wsum = wval(1) + wval(2) + ZERO
         wval(1) = wval(1) / wsum
     !   wval(2) = wval(2) / wsum
-        wval(2) =-wval(1) + 1.e0 ! wval(2)/wsum but robust !
+        wval(2) =-wval(1) + 1.e0_8 ! wval(2)/wsum but robust !
 
         return
 
@@ -168,15 +168,15 @@
         
         if (ipos-halo.lt.head) then
         
-            omax = 1.e0 
-            omin = 0.e0 ; return
+            omax = 1.e0_8 
+            omin = 0.e0_8 ; return
         
         end if
         
         if (ipos+halo.gt.tail) then
             
-            omax = 1.e0 
-            omin = 0.e0 ; return
+            omax = 1.e0_8 
+            omin = 0.e0_8 ; return
         
         end if
 
@@ -203,7 +203,7 @@
       
     !---------------------------------------- "lower" part !
       
-        delh = 0.e0
+        delh = 0.e0_8
 
         do  hpos = ipos-1, imin, -1
         
@@ -235,7 +235,7 @@
       
     !---------------------------------------- "upper" part !     
       
-        delh = 0.e0
+        delh = 0.e0_8
         
         do  hpos = ipos+1, imax, +1
         
@@ -319,15 +319,15 @@
         
         if (ipos-halo.lt.head) then
         
-            omax = 1.e0 
-            omin = 0.e0 ; return
+            omax = 1.e0_8 
+            omin = 0.e0_8 ; return
         
         end if
         
         if (ipos+halo.gt.tail) then
             
-            omax = 1.e0 
-            omin = 0.e0 ; return
+            omax = 1.e0_8 
+            omin = 0.e0_8 ; return
         
         end if
 
@@ -343,21 +343,21 @@
         dfx1 = oscl(1,ivar,ipos)
         dfx2 = oscl(2,ivar,ipos)
       
-        oval = (2.e0**1*dfx1)**2 &
-        &    + (2.e0**2*dfx2)**2
+        oval = (2.e0_8**1*dfx1)**2 &
+        &    + (2.e0_8**2*dfx2)**2
       
         omin = oval
         omax = oval
       
     !---------------------------------------- "lower" part !
       
-        delh = 0.e0
+        delh = 0.e0_8
 
         do  hpos = ipos-1, imin, -1
         
     !------------------ calc. derivatives centred on IPOS. !     
 
-        delh = delh + 2.e0
+        delh = delh + 2.e0_8
         
         dfx1 = oscl(1,ivar,hpos)
         dfx2 = oscl(2,ivar,hpos)
@@ -366,8 +366,8 @@
 
     !------------------ indicator: NORM(H^N * D^N/DX^N(F)) !
 
-        oval = (2.e0**1*dfx1)**2 &
-    &        + (2.e0**2*dfx2)**2
+        oval = (2.e0_8**1*dfx1)**2 &
+    &        + (2.e0_8**2*dfx2)**2
 
         if (oval .lt. omin) then
             omin = oval 
@@ -380,13 +380,13 @@
       
     !---------------------------------------- "upper" part !     
       
-        delh = 0.e0
+        delh = 0.e0_8
         
         do  hpos = ipos+1, imax, +1
         
     !------------------ calc. derivatives centred on IPOS. !     
 
-        delh = delh - 2.e0
+        delh = delh - 2.e0_8
         
         dfx1 = oscl(1,ivar,hpos)
         dfx2 = oscl(2,ivar,hpos)
@@ -395,8 +395,8 @@
     
     !------------------ indicator: NORM(H^N * D^N/DX^N(F)) !
 
-        oval = (2.e0**1*dfx1)**2 &
-    &        + (2.e0**2*dfx2)**2
+        oval = (2.e0_8**1*dfx1)**2 &
+    &        + (2.e0_8**2*dfx2)**2
   
         if (oval .lt. omin) then
             omin = oval 
diff --git a/src/ICE/icedyn_adv_umx.F90 b/src/ICE/icedyn_adv_umx.F90
index 29f1db15e33ae2ffae9d7c8644ff2c47d8837745..8e79fb8c8650820bc3a86a96d98d067e9b12115c 100644
--- a/src/ICE/icedyn_adv_umx.F90
+++ b/src/ICE/icedyn_adv_umx.F90
@@ -1157,7 +1157,7 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj    ) :: zbup, zbdo
       REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_ups, ztj_ups
       !!----------------------------------------------------------------------
-      zbig = 1.e+40_wp
+      zbig = 1.e+20_wp   ! works ok with simple/double precison
 
       ! antidiffusive flux : high order minus low order
       ! --------------------------------------------------
diff --git a/src/ICE/icedyn_rhg_vp.F90 b/src/ICE/icedyn_rhg_vp.F90
index c087f97a8c2cad4cd03c234b646906e34eeaf515..d59afb869fc488c29ce88b6ce92633a44a07674d 100644
--- a/src/ICE/icedyn_rhg_vp.F90
+++ b/src/ICE/icedyn_rhg_vp.F90
@@ -1115,7 +1115,7 @@ CONTAINS
 
          END_2D
 
-!!$         CALL lbc_lnk( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. )
+!!$         CALL lbc_lnk( 'icedyn_rhg_vp', zs1, 'T', 1._wp, zs2, 'T', 1._wp, zs12, 'T', 1._wp )
       
       ENDIF
       
@@ -1132,7 +1132,7 @@ CONTAINS
                
          END_2D
 
-         CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. )
+         CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1._wp )
          
       ENDIF
       
@@ -1169,8 +1169,8 @@ CONTAINS
          END_2D
          
          !
-         CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, &
-!            &                          ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
+         CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1._wp, ztauy_oi, 'V', -1._wp, ztaux_ai, 'U', -1._wp, ztauy_ai, 'V', -1._wp ) !, &
+!            &                          ztaux_bi, 'U', -1._wp, ztauy_bi, 'V', -1._wp )
          !
          CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
          CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
@@ -1200,7 +1200,7 @@ CONTAINS
             zsig_II(ji,jj)   =   0.5_wp * SQRT ( zs2(ji,jj) * zs2(ji,jj) + 4. * zs12(ji,jj) * zs12(ji,jj) )
          END_2D
 
-!!$         CALL lbc_lnk( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.)
+!!$         CALL lbc_lnk( 'icedyn_rhg_vp', zsig_I, 'T', 1._wp, zsig_II, 'T', 1._wp)
          
          IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,   zsig_I(:,:)  * zmsk00(:,:) ) ! Normal stress
          IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' ,   zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
@@ -1241,7 +1241,7 @@ CONTAINS
                
          END_2D
          !
-         ! CALL lbc_lnk( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.)
+         ! CALL lbc_lnk( 'icedyn_rhg_vp', zsig1_p, 'T', 1._wp, zsig2_p, 'T', 1._wp)
          !
          CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
          CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
@@ -1265,8 +1265,8 @@ CONTAINS
                            &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
          END_2D
          !
-         CALL lbc_lnk( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., &
-            &                           zCorU, 'U', -1., zCorV, 'V', -1. )
+         CALL lbc_lnk( 'icedyn_rhg_vp', zspgU, 'U', -1._wp, zspgV, 'V', -1._wp, &
+            &                           zCorU, 'U', -1._wp, zCorV, 'V', -1._wp )
          !
          CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
          CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
@@ -1296,7 +1296,7 @@ CONTAINS
 
          END_2D
             
-         CALL lbc_lnk( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. )
+         CALL lbc_lnk( 'icedyn_rhg_vp', zfU, 'U', -1._wp, zfV, 'V', -1._wp )
          
          CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
          CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
@@ -1326,9 +1326,9 @@ CONTAINS
 
          END_2D
          
-         CALL lbc_lnk( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
-            &                           zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
-            &                           zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
+         CALL lbc_lnk( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1._wp, zdiag_ymtrp_ice, 'V', -1._wp, &
+            &                           zdiag_xmtrp_snw, 'U', -1._wp, zdiag_ymtrp_snw, 'V', -1._wp, &
+            &                           zdiag_xatrp    , 'U', -1._wp, zdiag_yatrp    , 'V', -1._wp )
 
          CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
          CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport 
@@ -1612,14 +1612,14 @@ CONTAINS
 
          END_2D
          
-         IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_res,  'U',  1., zv_res , 'V',  1. )
+         IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_res,  'U',  1._wp, zv_res , 'V',  1._wp )
 
          DO_2D( 0, 0, 0, 0 ) !clem check bounds
         
                pvel_res(ji,jj) = 0.25_wp * ( zu_res(ji-1,jj) + zu_res(ji,jj) + zv_res(ji,jj-1) + zv_res(ji,jj) )
          
          END_2D
-         CALL lbc_lnk( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1. )
+         CALL lbc_lnk( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1._wp )
 
       ELSE
 
diff --git a/src/ICE/icethd.F90 b/src/ICE/icethd.F90
index 3971c4415240950ed6423023821b958247370d3f..fbd87f717158afabb8540444462742d147a4925e 100644
--- a/src/ICE/icethd.F90
+++ b/src/ICE/icethd.F90
@@ -342,11 +342,22 @@ CONTAINS
          CALL tab_2d_1d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
          !
          ! --- Change units of e_i, e_s from J/m2 to J/m3 --- !
+         ! Here we make sure that we don't divide by very small, but physically
+         ! meaningless, products of sea ice thicknesses/snow depths and sea ice
+         ! concentration
          DO jk = 1, nlay_i
-            WHERE( h_i_1d(1:npti)>0._wp ) e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i
+            WHERE( (h_i_1d(1:npti) * a_i_1d(1:npti)) > epsi20 )
+               e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i
+            ELSEWHERE
+               e_i_1d(1:npti,jk) = 0._wp
+            ENDWHERE
          END DO
          DO jk = 1, nlay_s
-            WHERE( h_s_1d(1:npti)>0._wp ) e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s
+            WHERE( (h_s_1d(1:npti) * a_i_1d(1:npti)) > epsi20 )
+               e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s
+            ELSEWHERE
+               e_s_1d(1:npti,jk) = 0._wp
+            ENDWHERE
          END DO
          !
          !                 !---------------------!
diff --git a/src/NST/vremap.F90 b/src/NST/vremap.F90
index 6e3853bff00e1117653109e5c839ae2d31aab04e..a16b1f83d21dd9e6666c5ebe9f6e0faf26b9c0fb 100644
--- a/src/NST/vremap.F90
+++ b/src/NST/vremap.F90
@@ -292,8 +292,8 @@ CONTAINS
       !
       INTEGER, PARAMETER :: ndof = 1
       INTEGER  :: jk, jn
-      REAL(wp) ::  zwin(kjpk_in+1) ,  ztin(ndof, kn_var, kjpk_in)
-      REAL(wp) :: zwout(kjpk_out+1), ztout(ndof, kn_var, kjpk_out)
+      REAL(dp) ::  zwin(kjpk_in+1) ,  ztin(ndof, kn_var, kjpk_in)    ! rmap1d uses dp
+      REAL(dp) :: zwout(kjpk_out+1), ztout(ndof, kn_var, kjpk_out)   ! rmap1d uses dp
       TYPE(rmap_work) :: work
       TYPE(rmap_opts) :: opts
       TYPE(rcon_ends) :: bc_l(kn_var)
diff --git a/src/OCE/DOM/domzgr.F90 b/src/OCE/DOM/domzgr.F90
index 9a5b6fc364b48758b08b32d0214c084a3d808f6a..0d63bbb168b69d9059e9691181e5c13e49b7da69 100644
--- a/src/OCE/DOM/domzgr.F90
+++ b/src/OCE/DOM/domzgr.F90
@@ -139,7 +139,7 @@ CONTAINS
       IF( .NOT. ( l_Jperio .OR. l_NFold ) ) THEN                   ! N closed:
          zmsk(:,mj0(jpjglo-nn_hls):mj1(jpjglo-nn_hls)  ) = 0._wp   ! last    line of inner global domain at 0
       ENDIF
-      CALL lbc_lnk( 'usrdef_zgr', zmsk, 'T', 1. )             ! set halos
+      CALL lbc_lnk( 'usrdef_zgr', zmsk, 'T', 1._wp )             ! set halos
       k_top(:,:) = k_top(:,:) * NINT( zmsk(:,:) )
       !
 !!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears
diff --git a/src/OCE/IOM/iom.F90 b/src/OCE/IOM/iom.F90
index 7c336684c5d5048f1adc4191e3a70c20e48dbea2..1c78c9021fea692a2d1e337a2b70666483190543 100644
--- a/src/OCE/IOM/iom.F90
+++ b/src/OCE/IOM/iom.F90
@@ -1044,19 +1044,13 @@ CONTAINS
       INTEGER         , INTENT(in   )                         ::   kdom      ! Type of domain to be read
       CHARACTER(len=*), INTENT(in   )                         ::   cdvar     ! Name of the variable
       REAL(sp)        , INTENT(  out), DIMENSION(:)           ::   pvar      ! read field
-      REAL(dp)        , ALLOCATABLE  , DIMENSION(:)           ::   ztmp_pvar ! tmp var to read field
       INTEGER         , INTENT(in   )              , OPTIONAL ::   ktime     ! record number
       INTEGER         , INTENT(in   ), DIMENSION(1), OPTIONAL ::   kstart    ! start axis position of the reading
       INTEGER         , INTENT(in   ), DIMENSION(1), OPTIONAL ::   kcount    ! number of points in each axis
       !
       IF( kiomid > 0 ) THEN
-         IF( iom_file(kiomid)%nfid > 0 ) THEN
-            ALLOCATE(ztmp_pvar(size(pvar,1)))
-            CALL iom_get_123d( kiomid, kdom       , cdvar        , pv_r1d=ztmp_pvar,   &
-              &                                                     ktime=ktime, kstart=kstart, kcount=kcount )
-            pvar = ztmp_pvar
-            DEALLOCATE(ztmp_pvar)
-         END IF
+         IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar, pvsp1d = pvar,   &
+               &                                            ktime = ktime, kstart = kstart, kcount = kcount )
       ENDIF
    END SUBROUTINE iom_g1d_sp
 
@@ -1071,8 +1065,8 @@ CONTAINS
       INTEGER         , INTENT(in   ), DIMENSION(1), OPTIONAL ::   kcount    ! number of points in each axis
       !
       IF( kiomid > 0 ) THEN
-         IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom       , cdvar        , pv_r1d=pvar,   &
-              &                                                     ktime=ktime, kstart=kstart, kcount=kcount)
+         IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar, pvdp1d = pvar,   &
+            &                                               ktime = ktime, kstart = kstart, kcount = kcount)
       ENDIF
    END SUBROUTINE iom_g1d_dp
 
@@ -1081,23 +1075,17 @@ CONTAINS
       INTEGER         , INTENT(in   )                         ::   kdom      ! Type of domain to be read
       CHARACTER(len=*), INTENT(in   )                         ::   cdvar     ! Name of the variable
       REAL(sp)        , INTENT(  out), DIMENSION(:,:)         ::   pvar      ! read field
-      REAL(dp)        , ALLOCATABLE  , DIMENSION(:,:)         ::   ztmp_pvar ! tmp var to read field
       INTEGER         , INTENT(in   )              , OPTIONAL ::   ktime     ! record number
       CHARACTER(len=1), INTENT(in   )              , OPTIONAL ::   cd_type   ! nature of grid-points (T, U, V, F, W)
-      REAL(dp)        , INTENT(in   )              , OPTIONAL ::   psgn      ! -1.(1.): (not) change sign across the north fold
+      REAL(sp)        , INTENT(in   )              , OPTIONAL ::   psgn      ! -1.(1.): (not) change sign across the north fold
       INTEGER         , INTENT(in   )              , OPTIONAL ::   kfill     ! value of kfillmode in lbc_lbk
       INTEGER         , INTENT(in   ), DIMENSION(2), OPTIONAL ::   kstart    ! start axis position of the reading
       INTEGER         , INTENT(in   ), DIMENSION(2), OPTIONAL ::   kcount    ! number of points in each axis
       !
       IF( kiomid > 0 ) THEN
-         IF( iom_file(kiomid)%nfid > 0 ) THEN
-            ALLOCATE(ztmp_pvar(size(pvar,1), size(pvar,2)))
-            CALL iom_get_123d( kiomid, kdom, cdvar      , pv_r2d = ztmp_pvar  , ktime = ktime,   &
-             &                                                      cd_type = cd_type, psgn   = psgn  , kfill = kfill,   &
-             &                                                      kstart  = kstart , kcount = kcount  )
-            pvar = ztmp_pvar
-            DEALLOCATE(ztmp_pvar)
-         ENDIF
+         IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar, pvsp2d = pvar,                 &
+            &                                               cd_type = cd_type, psgn_sp = psgn, kfill = kfill,   &
+            &                                               ktime = ktime, kstart = kstart, kcount = kcount )
       ENDIF
    END SUBROUTINE iom_g2d_sp
 
@@ -1114,9 +1102,9 @@ CONTAINS
       INTEGER         , INTENT(in   ), DIMENSION(2), OPTIONAL ::   kcount    ! number of points in each axis
       !
       IF( kiomid > 0 ) THEN
-         IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar      , pv_r2d = pvar  , ktime = ktime,   &
-            &                                                       cd_type = cd_type, psgn   = psgn  , kfill = kfill,   &
-            &                                                       kstart  = kstart , kcount = kcount                )
+         IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar, pvdp2d = pvar,                 &
+            &                                               cd_type = cd_type, psgn_dp = psgn, kfill = kfill,   &
+            &                                               ktime = ktime, kstart = kstart, kcount = kcount )
       ENDIF
    END SUBROUTINE iom_g2d_dp
 
@@ -1125,23 +1113,17 @@ CONTAINS
       INTEGER         , INTENT(in   )                         ::   kdom      ! Type of domain to be read
       CHARACTER(len=*), INTENT(in   )                         ::   cdvar     ! Name of the variable
       REAL(sp)        , INTENT(  out), DIMENSION(:,:,:)       ::   pvar      ! read field
-      REAL(dp)        , ALLOCATABLE  , DIMENSION(:,:,:)       ::   ztmp_pvar ! tmp var to read field
       INTEGER         , INTENT(in   )              , OPTIONAL ::   ktime     ! record number
       CHARACTER(len=1), INTENT(in   )              , OPTIONAL ::   cd_type   ! nature of grid-points (T, U, V, F, W)
-      REAL(dp)        , INTENT(in   )              , OPTIONAL ::   psgn      ! -1.(1.) : (not) change sign across the north fold
+      REAL(sp)        , INTENT(in   )              , OPTIONAL ::   psgn      ! -1.(1.) : (not) change sign across the north fold
       INTEGER         , INTENT(in   )              , OPTIONAL ::   kfill     ! value of kfillmode in lbc_lbk
       INTEGER         , INTENT(in   ), DIMENSION(3), OPTIONAL ::   kstart    ! start axis position of the reading
       INTEGER         , INTENT(in   ), DIMENSION(3), OPTIONAL ::   kcount    ! number of points in each axis
       !
       IF( kiomid > 0 ) THEN
-         IF( iom_file(kiomid)%nfid > 0 ) THEN
-            ALLOCATE(ztmp_pvar(size(pvar,1), size(pvar,2), size(pvar,3)))
-            CALL iom_get_123d( kiomid, kdom, cdvar      , pv_r3d = ztmp_pvar  , ktime = ktime,   &
-            &                                                       cd_type = cd_type, psgn   = psgn  , kfill = kfill,   &
-            &                                                       kstart  = kstart , kcount = kcount                )
-            pvar = ztmp_pvar
-            DEALLOCATE(ztmp_pvar)
-         END IF
+         IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar, pvsp3d = pvar,                 &
+            &                                               cd_type = cd_type, psgn_sp = psgn, kfill = kfill,   &
+            &                                               ktime = ktime, kstart = kstart, kcount = kcount )
       ENDIF
    END SUBROUTINE iom_g3d_sp
 
@@ -1158,18 +1140,16 @@ CONTAINS
       INTEGER         , INTENT(in   ), DIMENSION(3), OPTIONAL ::   kcount    ! number of points in each axis
       !
       IF( kiomid > 0 ) THEN
-         IF( iom_file(kiomid)%nfid > 0 ) THEN
-            CALL iom_get_123d( kiomid, kdom, cdvar      , pv_r3d = pvar  , ktime = ktime,   &
-            &                                                       cd_type = cd_type, psgn   = psgn  , kfill = kfill,   &
-            &                                                       kstart  = kstart , kcount = kcount                )
-         END IF
+         IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar, pvdp3d = pvar,                 &
+            &                                               cd_type = cd_type, psgn_dp = psgn, kfill = kfill,   &
+            &                                               ktime = ktime, kstart = kstart, kcount = kcount )
       ENDIF
    END SUBROUTINE iom_g3d_dp
 
    !!----------------------------------------------------------------------
 
-   SUBROUTINE iom_get_123d( kiomid , kdom, cdvar, pv_r1d, pv_r2d, pv_r3d, ktime ,   &
-         &                  cd_type, psgn, kfill, kstart, kcount )
+   SUBROUTINE iom_get_123d( kiomid , kdom, cdvar, pvsp1d, pvsp2d, pvsp3d, pvdp1d, pvdp2d, pvdp3d,   &
+         &                  ktime , cd_type, psgn_sp, psgn_dp, kfill, kstart, kcount )
       !!-----------------------------------------------------------------------
       !!                  ***  ROUTINE  iom_get_123d  ***
       !!
@@ -1180,12 +1160,16 @@ CONTAINS
       INTEGER                    , INTENT(in   )           ::   kiomid    ! Identifier of the file
       INTEGER                    , INTENT(in   )           ::   kdom      ! Type of domain to be read
       CHARACTER(len=*)           , INTENT(in   )           ::   cdvar     ! Name of the variable
-      REAL(dp), DIMENSION(:)     , INTENT(  out), OPTIONAL ::   pv_r1d    ! read field (1D case)
-      REAL(dp), DIMENSION(:,:)   , INTENT(  out), OPTIONAL ::   pv_r2d    ! read field (2D case)
-      REAL(dp), DIMENSION(:,:,:) , INTENT(  out), OPTIONAL ::   pv_r3d    ! read field (3D case)
+      REAL(sp), DIMENSION(:)     , INTENT(  out), OPTIONAL ::   pvsp1d    ! read field (1D case), single precision
+      REAL(sp), DIMENSION(:,:)   , INTENT(  out), OPTIONAL ::   pvsp2d    ! read field (2D case), single precision
+      REAL(sp), DIMENSION(:,:,:) , INTENT(  out), OPTIONAL ::   pvsp3d    ! read field (3D case), single precision
+      REAL(dp), DIMENSION(:)     , INTENT(  out), OPTIONAL ::   pvdp1d    ! read field (1D case), double precision
+      REAL(dp), DIMENSION(:,:)   , INTENT(  out), OPTIONAL ::   pvdp2d    ! read field (2D case), double precision
+      REAL(dp), DIMENSION(:,:,:) , INTENT(  out), OPTIONAL ::   pvdp3d    ! read field (3D case), double precision
       INTEGER                    , INTENT(in   ), OPTIONAL ::   ktime     ! record number
       CHARACTER(len=1)           , INTENT(in   ), OPTIONAL ::   cd_type   ! nature of grid-points (T, U, V, F, W)
-      REAL(dp)                   , INTENT(in   ), OPTIONAL ::   psgn      ! -1.(1.) : (not) change sign across the north fold
+      REAL(sp)                   , INTENT(in   ), OPTIONAL ::   psgn_sp   ! -1.(1.) : (not) change sign across the north fold
+      REAL(dp)                   , INTENT(in   ), OPTIONAL ::   psgn_dp   ! -1.(1.) : (not) change sign across the north fold
       INTEGER                    , INTENT(in   ), OPTIONAL ::   kfill     ! value of kfillmode in lbc_lbk
       INTEGER , DIMENSION(:)     , INTENT(in   ), OPTIONAL ::   kstart    ! start position of the reading in each axis
       INTEGER , DIMENSION(:)     , INTENT(in   ), OPTIONAL ::   kcount    ! number of points to be read in each axis
@@ -1206,22 +1190,36 @@ CONTAINS
       INTEGER, DIMENSION(jpmax_dims) ::   icnt        ! number of value to read along each axis
       INTEGER, DIMENSION(jpmax_dims) ::   idimsz      ! size of the dimensions of the variable
       INTEGER, DIMENSION(jpmax_dims) ::   ishape      ! size of the dimensions of the variable
-      REAL(dp)                       ::   zscf, zofs  ! sacle_factor and add_offset
-      REAL(wp)                       ::   zsgn        ! local value of psgn
+      REAL(sp)                       ::   zscf_sp, zofs_sp  ! sacle_factor and add_offset, single precision
+      REAL(dp)                       ::   zscf_dp, zofs_dp  ! sacle_factor and add_offset, double precision
+      REAL(sp)                       ::   zsgn_sp     ! local value of psgn, single precision
+      REAL(dp)                       ::   zsgn_dp     ! local value of psgn, double precision
       INTEGER                        ::   itmp        ! temporary integer
       CHARACTER(LEN=256)             ::   clinfo      ! info character
       CHARACTER(LEN=256)             ::   clname      ! file name
       CHARACTER(LEN=1)               ::   clrankpv, cldmspc      !
       CHARACTER(LEN=1)               ::   cl_type     ! local value of cd_type
       LOGICAL                        ::   ll_only3rd  ! T => if kstart, kcount present then *only* use values for 3rd spatial dimension.
+      LOGICAL                        ::   llis1d, llis2d, llis3d
+      LOGICAL                        ::   llsp        ! use single precision
       INTEGER                        ::   inlev       ! number of levels for 3D data
-      REAL(dp)                       ::   gma, gmi
       !---------------------------------------------------------------------
       CHARACTER(LEN=lc)                               ::   context
       !
       CALL set_xios_context(kiomid, context)
+      !
+      llsp = PRESENT(pvsp1d) .OR. PRESENT(pvsp2d) .OR. PRESENT(pvsp3d)
+      IF( llsp ) THEN 
+         llis1d = PRESENT(pvsp1d)   ;   IF( llis1d )   ishape(1:1) = SHAPE(pvsp1d)
+         llis2d = PRESENT(pvsp2d)   ;   IF( llis2d )   ishape(1:2) = SHAPE(pvsp2d)
+         llis3d = PRESENT(pvsp3d)   ;   IF( llis3d )   ishape(1:3) = SHAPE(pvsp3d)
+      ELSE
+         llis1d = PRESENT(pvdp1d)   ;   IF( llis1d )   ishape(1:1) = SHAPE(pvdp1d)
+         llis2d = PRESENT(pvdp2d)   ;   IF( llis2d )   ishape(1:2) = SHAPE(pvdp2d)
+         llis3d = PRESENT(pvdp3d)   ;   IF( llis3d )   ishape(1:3) = SHAPE(pvdp3d)
+      ENDIF
       inlev = -1
-      IF( PRESENT(pv_r3d) )   inlev = SIZE(pv_r3d, 3)
+      IF( llis3d )   inlev = ishape(3)
       !
       idom = kdom
       istop = nstop
@@ -1263,27 +1261,27 @@ CONTAINS
             itime = 1
             IF( PRESENT(ktime) ) itime = ktime
             !
-            irankpv = 1 * COUNT( (/PRESENT(pv_r1d)/) ) + 2 * COUNT( (/PRESENT(pv_r2d)/) ) + 3 * COUNT( (/PRESENT(pv_r3d)/) )
+            irankpv = 1 * COUNT( (/ llis1d /) ) + 2 * COUNT( (/ llis2d /) ) + 3 * COUNT( (/ llis3d /) )
             WRITE(clrankpv, fmt='(i1)') irankpv
             WRITE(cldmspc , fmt='(i1)') idmspc
             !
             IF(     idmspc <  irankpv ) THEN                     ! it seems we want to read more than we can...
-               IF(     irankpv == 3 .AND. idmspc == 2 ) THEN     !   3D input array from 2D spatial data in the file:
+               IF(     llis3d .AND. idmspc == 2 ) THEN           !   3D input array from 2D spatial data in the file:
                   llok = inlev == 1                              !     -> 3rd dimension must be equal to 1
-               ELSEIF( irankpv == 3 .AND. idmspc == 1 ) THEN     !   3D input array from 1D spatial data in the file:
-                  llok = inlev == 1 .AND. SIZE(pv_r3d, 2) == 1   !     -> 2nd and 3rd dimensions must be equal to 1
-               ELSEIF( irankpv == 2 .AND. idmspc == 2 ) THEN     !   2D input array from 1D spatial data in the file:
-                  llok = SIZE(pv_r2d, 2) == 1                    !     -> 2nd dimension must be equal to 1
+               ELSEIF( llis3d .AND. idmspc == 1 ) THEN           !   3D input array from 1D spatial data in the file:
+                  llok = inlev == 1 .AND. ishape(2) == 1         !     -> 2nd and 3rd dimensions must be equal to 1
+               ELSEIF( llis3d .AND. idmspc == 2 ) THEN           !   2D input array from 1D spatial data in the file:
+                  llok = ishape(2) == 1                          !     -> 2nd dimension must be equal to 1
                ELSE
                   llok = .FALSE.
                ENDIF
                IF( .NOT. llok )   CALL ctl_stop( TRIM(clinfo), 'The file has only '//cldmspc//' spatial dimension',   &
                   &                                            '=> cannot read a true '//clrankpv//'D array from this file...' )
             ELSEIF( idmspc == irankpv ) THEN
-               IF( PRESENT(pv_r1d) .AND. idom /= jpdom_unknown )   &
+               IF( llis1d .AND. idom /= jpdom_unknown )   &
                   &   CALL ctl_stop( TRIM(clinfo), 'case not coded...You must use jpdom_unknown' )
             ELSEIF( idmspc >  irankpv ) THEN                     ! it seems we want to read less than we should...
-                  IF( PRESENT(pv_r2d) .AND. itime == 1 .AND. idimsz(3) == 1 .AND. idmspc == 3 ) THEN
+                  IF( llis2d .AND. itime == 1 .AND. idimsz(3) == 1 .AND. idmspc == 3 ) THEN
                      CALL ctl_warn( TRIM(clinfo), '2D array input but 3 spatial dimensions in the file...'              ,   &
                            &         'As the size of the z dimension is 1 and as we try to read the first record, ',   &
                            &         'we accept this case, even if there is a possible mix-up between z and time dimension' )
@@ -1309,11 +1307,11 @@ CONTAINS
                ELSE
                   icnt  (1:idmspc) = idimsz(1:idmspc)
                ENDIF
-            ELSE   !   not a 1D array as pv_r1d requires jpdom_unknown
+            ELSE   !   not a 1D array as pv(sd)p1d requires jpdom_unknown
                ! we do not read the overlap and the extra-halos -> from Nis0 to Nie0 and from Njs0 to Nje0
                IF( idom == jpdom_global )   istart(1:2) = (/ mig0(Nis0), mjg0(Njs0) /)
                icnt(1:2) = (/ Ni_0, Nj_0 /)
-               IF( PRESENT(pv_r3d) ) THEN
+               IF( llis3d ) THEN
                   IF( idom == jpdom_auto_xy ) THEN
                      istart(3) = kstart(3)
                      icnt  (3) = kcount(3)
@@ -1337,21 +1335,28 @@ CONTAINS
             ! check that icnt matches the input array
             !-
             IF( idom == jpdom_unknown ) THEN
-               IF( irankpv == 1 )        ishape(1:1) = SHAPE(pv_r1d)
-               IF( irankpv == 2 )        ishape(1:2) = SHAPE(pv_r2d)
-               IF( irankpv == 3 )        ishape(1:3) = SHAPE(pv_r3d)
                ctmp1 = 'd'
-            ELSE
-               IF( irankpv == 2 ) THEN
-                  ishape(1:2) = SHAPE(pv_r2d(Nis0:Nie0,Njs0:Nje0  ))   ;   ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0)'
+            ELSE                    ! we must redefine ishape as we don't read the full array
+               IF( llis2d ) THEN
+                  IF( llsp ) THEN   ;   ishape(1:2) = SHAPE(pvsp2d(Nis0:Nie0,Njs0:Nje0  ))
+                  ELSE              ;   ishape(1:2) = SHAPE(pvdp2d(Nis0:Nie0,Njs0:Nje0  ))
+                  ENDIF
+                  ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0)'
                ENDIF
-               IF( irankpv == 3 ) THEN
-                  ishape(1:3) = SHAPE(pv_r3d(Nis0:Nie0,Njs0:Nje0,:))   ;   ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0,:)'
+               IF( llis3d ) THEN
+                  IF( llsp ) THEN   ;   ishape(1:3) = SHAPE(pvsp3d(Nis0:Nie0,Njs0:Nje0,:))
+                  ELSE              ;   ishape(1:3) = SHAPE(pvdp3d(Nis0:Nie0,Njs0:Nje0,:))
+                  ENDIF
+                  ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0,:)'
                ENDIF
             ENDIF
             DO jl = 1, irankpv
                WRITE( ctmp2, FMT="(', ', i1,'): ', i5,' /= icnt(', i1,'):', i5)" ) jl, ishape(jl), jl, icnt(jl)
-               IF( ishape(jl) /= icnt(jl) )   CALL ctl_stop( TRIM(clinfo), 'size(pv_r'//clrankpv//TRIM(ctmp1)//TRIM(ctmp2) )
+               IF( llsp ) THEN
+                  IF( ishape(jl) /= icnt(jl) )   CALL ctl_stop( TRIM(clinfo), 'size(pvsp'//clrankpv//TRIM(ctmp1)//TRIM(ctmp2) )
+               ELSE
+                  IF( ishape(jl) /= icnt(jl) )   CALL ctl_stop( TRIM(clinfo), 'size(pvdp'//clrankpv//TRIM(ctmp1)//TRIM(ctmp2) )
+               ENDIF
             END DO
 
          ENDIF
@@ -1365,21 +1370,32 @@ CONTAINS
             ELSE                               ;   ix1 = 1      ;   ix2 = icnt(1)   ;   iy1 = 1      ;   iy2 = icnt(2)
             ENDIF
 
-            CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, pv_r3d )
+            CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2,   &
+               &               pvsp1d, pvsp2d, pvsp3d, pvdp1d, pvdp2d, pvdp3d )
 
             IF( istop == nstop ) THEN   ! no additional errors until this point...
                IF(lwp) WRITE(numout,"(10x,' read ',a,' (rec: ',i6,') in ',a,' ok')") TRIM(cdvar), itime, TRIM(iom_file(kiomid)%name)
 
                cl_type = 'T'
                IF( PRESENT(cd_type) )   cl_type = cd_type
-               zsgn = 1._wp
-               IF( PRESENT(psgn   ) )   zsgn    = psgn
-               !--- overlap areas and extra hallows (mpp)
-               IF(     PRESENT(pv_r2d) .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
-                  CALL lbc_lnk( 'iom', pv_r2d, cl_type, zsgn, kfillmode = kfill )
-               ELSEIF( PRESENT(pv_r3d) .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
-                  CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill )
+               IF( llsp ) THEN
+                  zsgn_sp = 1._sp
+                  IF( PRESENT(psgn_sp) )   zsgn_sp = psgn_sp
+                  IF(     llis2d .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+                     CALL lbc_lnk( 'iom', pvsp2d, cl_type, zsgn_sp, kfillmode = kfill )
+                  ELSEIF( llis3d .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+                     CALL lbc_lnk( 'iom', pvsp3d, cl_type, zsgn_sp, kfillmode = kfill )
+                  ENDIF
+               ELSE
+                  zsgn_dp = 1._dp
+                  IF( PRESENT(psgn_dp) )   zsgn_dp = psgn_dp
+                  IF(     llis2d .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+                     CALL lbc_lnk( 'iom', pvdp2d, cl_type, zsgn_dp, kfillmode = kfill )
+                  ELSEIF( llis3d .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+                     CALL lbc_lnk( 'iom', pvdp3d, cl_type, zsgn_dp, kfillmode = kfill )
+                  ENDIF
                ENDIF
+               !--- overlap areas and extra hallows (mpp)
                !
             ELSE
                ! return if istop == nstop is false
@@ -1395,26 +1411,39 @@ CONTAINS
 !would be good to be able to check which context is active and swap only if current is not restart
          idvar = iom_varid( kiomid, cdvar )
          CALL iom_swap(context)
-         zsgn = 1._wp
-         IF( PRESENT(psgn   ) )   zsgn    = psgn
+         IF( llsp ) THEN
+            zsgn_sp = 1._sp   ;   IF( PRESENT(psgn_sp) )   zsgn_sp = psgn_sp
+         ELSE
+            zsgn_dp = 1._dp   ;   IF( PRESENT(psgn_dp) )   zsgn_dp = psgn_dp
+         ENDIF
          cl_type = 'T'
          IF( PRESENT(cd_type) )   cl_type = cd_type
 
-         IF( PRESENT(pv_r3d) ) THEN
+         IF(     llis3d ) THEN
             IF(lwp) WRITE(numout,*) 'XIOS RST READ (3D): ',TRIM(cdvar)
-            CALL xios_recv_field( TRIM(cdvar), pv_r3d(:, :, :))
+            IF( llsp ) THEN   ;   CALL xios_recv_field( TRIM(cdvar), pvsp3d )
+            ELSE              ;   CALL xios_recv_field( TRIM(cdvar), pvdp3d )
+            ENDIF
             IF(idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
-               CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill)
+               IF( llsp ) THEN   ;   CALL lbc_lnk( 'iom', pvsp3d, cl_type, zsgn_sp, kfillmode = kfill)
+               ELSE              ;   CALL lbc_lnk( 'iom', pvdp3d, cl_type, zsgn_dp, kfillmode = kfill)
+               ENDIF
             ENDIF
-         ELSEIF( PRESENT(pv_r2d) ) THEN
+         ELSEIF( llis2d ) THEN
             IF(lwp) WRITE(numout,*) 'XIOS RST READ (2D): ', TRIM(cdvar)
-            CALL xios_recv_field( TRIM(cdvar), pv_r2d(:, :))
+            IF( llsp ) THEN   ;   CALL xios_recv_field( TRIM(cdvar), pvsp2d )
+            ELSE              ;   CALL xios_recv_field( TRIM(cdvar), pvdp2d )
+            ENDIF
             IF(idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
-               CALL lbc_lnk('iom', pv_r2d, cl_type, zsgn, kfillmode = kfill)
+               IF( llsp ) THEN   ;   CALL lbc_lnk('iom', pvsp2d, cl_type, zsgn_sp, kfillmode = kfill)
+               ELSE              ;   CALL lbc_lnk('iom', pvdp2d, cl_type, zsgn_dp, kfillmode = kfill)
+               ENDIF
             ENDIF
-         ELSEIF( PRESENT(pv_r1d) ) THEN
+         ELSEIF( llis1d ) THEN
             IF(lwp) WRITE(numout,*) 'XIOS RST READ (1D): ', TRIM(cdvar)
-            CALL xios_recv_field( TRIM(cdvar), pv_r1d)
+            IF( llsp ) THEN   ;   CALL xios_recv_field( TRIM(cdvar), pvsp1d )
+            ELSE              ;   CALL xios_recv_field( TRIM(cdvar), pvdp1d )
+            ENDIF
          ENDIF
          CALL iom_swap(cxios_context)
 #else
@@ -1424,17 +1453,32 @@ CONTAINS
       ENDIF
 
       !--- Apply scale_factor and offset
-      zscf = iom_file(kiomid)%scf(idvar)      ! scale factor
-      zofs = iom_file(kiomid)%ofs(idvar)      ! offset
-      IF(     PRESENT(pv_r1d) ) THEN
-         IF( zscf /= 1._wp )   pv_r1d(:) = pv_r1d(:) * zscf
-         IF( zofs /= 0._wp )   pv_r1d(:) = pv_r1d(:) + zofs
-      ELSEIF( PRESENT(pv_r2d) ) THEN
-         IF( zscf /= 1._wp)   pv_r2d(:,:) = pv_r2d(:,:) * zscf
-         IF( zofs /= 0._wp)   pv_r2d(:,:) = pv_r2d(:,:) + zofs
-      ELSEIF( PRESENT(pv_r3d) ) THEN
-         IF( zscf /= 1._wp)   pv_r3d(:,:,:) = pv_r3d(:,:,:) * zscf
-         IF( zofs /= 0._wp)   pv_r3d(:,:,:) = pv_r3d(:,:,:) + zofs
+      IF( llsp ) THEN
+         zscf_sp = iom_file(kiomid)%scf(idvar)      ! scale factor
+         zofs_sp = iom_file(kiomid)%ofs(idvar)      ! offset
+         IF(     llis1d ) THEN
+            IF( zscf_sp /= 1._sp )   pvsp1d(:    ) = pvsp1d(:    ) * zscf_sp
+            IF( zofs_sp /= 0._sp )   pvsp1d(:    ) = pvsp1d(:    ) + zofs_sp
+         ELSEIF( llis2d ) THEN
+            IF( zscf_sp /= 1._sp )   pvsp2d(:,:  ) = pvsp2d(:,:  ) * zscf_sp
+            IF( zofs_sp /= 0._sp )   pvsp2d(:,:  ) = pvsp2d(:,:  ) + zofs_sp
+         ELSEIF( llis3d ) THEN
+            IF( zscf_sp /= 1._sp )   pvsp3d(:,:,:) = pvsp3d(:,:,:) * zscf_sp
+            IF( zofs_sp /= 0._sp )   pvsp3d(:,:,:) = pvsp3d(:,:,:) + zofs_sp
+         ENDIF
+      ELSE
+         zscf_dp = iom_file(kiomid)%scf(idvar)      ! scale factor
+         zofs_dp = iom_file(kiomid)%ofs(idvar)      ! offset
+         IF(     llis1d ) THEN
+            IF( zscf_dp /= 1._dp )   pvdp1d(:    ) = pvdp1d(:    ) * zscf_dp
+            IF( zofs_dp /= 0._dp )   pvdp1d(:    ) = pvdp1d(:    ) + zofs_dp
+         ELSEIF( llis2d ) THEN
+            IF( zscf_dp /= 1._dp )   pvdp2d(:,:  ) = pvdp2d(:,:  ) * zscf_dp
+            IF( zofs_dp /= 0._dp )   pvdp2d(:,:  ) = pvdp2d(:,:  ) + zofs_dp
+         ELSEIF( llis3d ) THEN
+            IF( zscf_dp /= 1._dp )   pvdp3d(:,:,:) = pvdp3d(:,:,:) * zscf_dp
+            IF( zofs_dp /= 0._dp )   pvdp3d(:,:,:) = pvdp3d(:,:,:) + zofs_dp
+         ENDIF
       ENDIF
       !
    END SUBROUTINE iom_get_123d
@@ -1620,32 +1664,8 @@ CONTAINS
       INTEGER           :: ivid   ! variable id
       CHARACTER(LEN=lc) :: context
       !
-      CALL set_xios_context(kiomid, context)
-
-      llx = .NOT. (context == "NONE")
-
-      IF( llx ) THEN
-#ifdef key_xios
-         IF( kt == kwrite ) THEN
-            IF(lwp) write(numout,*) 'RESTART: write (XIOS 0D) ',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_put(TRIM(cdvar), pvar)
-            CALL iom_swap(cxios_context)
-         ELSE
-            IF(lwp) write(numout,*) 'RESTART: define (XIOS 0D) ',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_set_rstw_active( TRIM(cdvar), rs0 = pvar )
-            CALL iom_swap(cxios_context)
-         ENDIF
-#endif
-      ELSE
-         IF( kiomid > 0 ) THEN
-            IF( iom_file(kiomid)%nfid > 0 ) THEN
-               ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
-               CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r0d = real(pvar, dp) )
-            ENDIF
-         ENDIF
-      ENDIF
+      CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvsp0d = pvar )
+      !
    END SUBROUTINE iom_rp0d_sp
 
    SUBROUTINE iom_rp0d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype )
@@ -1660,32 +1680,8 @@ CONTAINS
       INTEGER           :: ivid   ! variable id
       CHARACTER(LEN=lc) :: context
       !
-      CALL set_xios_context(kiomid, context)
-
-      llx = .NOT. (context == "NONE")
-
-      IF( llx ) THEN
-#ifdef key_xios
-         IF( kt == kwrite ) THEN
-            IF(lwp) write(numout,*) 'RESTART: write (XIOS 0D) ',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_put(TRIM(cdvar), pvar)
-            CALL iom_swap(cxios_context)
-         ELSE
-            IF(lwp) write(numout,*) 'RESTART: define (XIOS 0D) ',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_set_rstw_active( TRIM(cdvar), rd0 = pvar )
-            CALL iom_swap(cxios_context)
-         ENDIF
-#endif
-      ELSE
-         IF( kiomid > 0 ) THEN
-            IF( iom_file(kiomid)%nfid > 0 ) THEN
-               ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
-               CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r0d = pvar )
-            ENDIF
-         ENDIF
-      ENDIF
+      CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvdp0d = pvar )
+      !
    END SUBROUTINE iom_rp0d_dp
 
 
@@ -1701,32 +1697,8 @@ CONTAINS
       INTEGER           :: ivid   ! variable id
       CHARACTER(LEN=lc) :: context
       !
-      CALL set_xios_context(kiomid, context)
-
-      llx = .NOT. (context == "NONE")
-
-      IF( llx ) THEN
-#ifdef key_xios
-         IF( kt == kwrite ) THEN
-            IF(lwp) write(numout,*) 'RESTART: write (XIOS 1D) ',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_put(TRIM(cdvar), pvar)
-            CALL iom_swap(cxios_context)
-         ELSE
-            IF(lwp) write(numout,*) 'RESTART: define (XIOS 1D)',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_set_rstw_active( TRIM(cdvar), rs1 = pvar )
-            CALL iom_swap(cxios_context)
-         ENDIF
-#endif
-      ELSE
-         IF( kiomid > 0 ) THEN
-            IF( iom_file(kiomid)%nfid > 0 ) THEN
-               ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
-               CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r1d = real(pvar, dp) )
-            ENDIF
-         ENDIF
-      ENDIF
+      CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvsp1d = pvar )
+      !
    END SUBROUTINE iom_rp1d_sp
 
    SUBROUTINE iom_rp1d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype )
@@ -1741,35 +1713,11 @@ CONTAINS
       INTEGER           :: ivid   ! variable id
       CHARACTER(LEN=lc) :: context
       !
-      CALL set_xios_context(kiomid, context)
-
-      llx = .NOT. (context == "NONE")
-
-      IF( llx ) THEN
-#ifdef key_xios
-         IF( kt == kwrite ) THEN
-            IF(lwp) write(numout,*) 'RESTART: write (XIOS 1D) ',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_put(TRIM(cdvar), pvar)
-            CALL iom_swap(cxios_context)
-         ELSE
-            IF(lwp) write(numout,*) 'RESTART: define (XIOS 1D)',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_set_rstw_active( TRIM(cdvar), rd1 = pvar )
-            CALL iom_swap(cxios_context)
-         ENDIF
-#endif
-      ELSE
-         IF( kiomid > 0 ) THEN
-            IF( iom_file(kiomid)%nfid > 0 ) THEN
-               ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
-               CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r1d = pvar )
-            ENDIF
-         ENDIF
-      ENDIF
+      CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvdp1d = pvar )
+      !
    END SUBROUTINE iom_rp1d_dp
 
-
+   
    SUBROUTINE iom_rp2d_sp( kt, kwrite, kiomid, cdvar, pvar, ktype )
       INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
       INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
@@ -1782,32 +1730,8 @@ CONTAINS
       INTEGER            :: ivid   ! variable id
       CHARACTER(LEN=lc)  :: context
       !
-      CALL set_xios_context(kiomid, context)
-
-      llx = .NOT. (context == "NONE")
-
-      IF( llx ) THEN
-#ifdef key_xios
-         IF( kt == kwrite ) THEN
-            IF(lwp) write(numout,*) 'RESTART: write (XIOS 2D) ',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_put(TRIM(cdvar), pvar)
-            CALL iom_swap(cxios_context)
-         ELSE
-            IF(lwp) write(numout,*) 'RESTART: define (XIOS 2D)',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_set_rstw_active( TRIM(cdvar), rs2 = pvar )
-            CALL iom_swap(cxios_context)
-         ENDIF
-#endif
-      ELSE
-         IF( kiomid > 0 ) THEN
-            IF( iom_file(kiomid)%nfid > 0 ) THEN
-               ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
-               CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r2d = real(pvar, dp) )
-            ENDIF
-         ENDIF
-      ENDIF
+      CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvsp2d = pvar )
+      !
    END SUBROUTINE iom_rp2d_sp
 
    SUBROUTINE iom_rp2d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype )
@@ -1822,32 +1746,8 @@ CONTAINS
       INTEGER           :: ivid   ! variable id
       CHARACTER(LEN=lc) :: context
       !
-      CALL set_xios_context(kiomid, context)
-
-      llx = .NOT. (context == "NONE")
-
-      IF( llx ) THEN
-#ifdef key_xios
-         IF( kt == kwrite ) THEN
-            IF(lwp) write(numout,*) 'RESTART: write (XIOS 2D) ',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_put(TRIM(cdvar), pvar)
-            CALL iom_swap(cxios_context)
-         ELSE
-            IF(lwp) write(numout,*) 'RESTART: define (XIOS 2D)',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_set_rstw_active( TRIM(cdvar), rd2 = pvar )
-            CALL iom_swap(cxios_context)
-         ENDIF
-#endif
-      ELSE
-         IF( kiomid > 0 ) THEN
-            IF( iom_file(kiomid)%nfid > 0 ) THEN
-               ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
-               CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r2d = pvar )
-            ENDIF
-         ENDIF
-      ENDIF
+      CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvdp2d = pvar )
+      !
    END SUBROUTINE iom_rp2d_dp
 
 
@@ -1863,32 +1763,8 @@ CONTAINS
       INTEGER           :: ivid   ! variable id
       CHARACTER(LEN=lc) :: context
       !
-      CALL set_xios_context(kiomid, context)
-
-      llx = .NOT. (context == "NONE")
-
-      IF( llx ) THEN
-#ifdef key_xios
-         IF( kt == kwrite ) THEN
-            IF(lwp) write(numout,*) 'RESTART: write (XIOS 3D) ',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_put(TRIM(cdvar), pvar)
-            CALL iom_swap(cxios_context)
-         ELSE
-            IF(lwp) write(numout,*) 'RESTART: define (XIOS 3D)',TRIM(cdvar)
-            CALL iom_swap(context)
-            CALL iom_set_rstw_active( TRIM(cdvar), rs3 = pvar )
-            CALL iom_swap(cxios_context)
-         ENDIF
-#endif
-      ELSE
-         IF( kiomid > 0 ) THEN
-            IF( iom_file(kiomid)%nfid > 0 ) THEN
-               ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
-               CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r3d = real(pvar, dp) )
-            ENDIF
-         ENDIF
-      ENDIF
+      CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvsp3d = pvar )
+      !
    END SUBROUTINE iom_rp3d_sp
 
    SUBROUTINE iom_rp3d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype )
@@ -1903,6 +1779,29 @@ CONTAINS
       INTEGER           :: ivid   ! variable id
       CHARACTER(LEN=lc) :: context
       !
+      CALL iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvdp3d = pvar )
+      !
+   END SUBROUTINE iom_rp3d_dp
+
+   SUBROUTINE iom_rp0123d( kt, kwrite, kiomid, cdvar, ktype, pvsp0d, pvsp1d, pvsp2d, pvsp3d, pvdp0d, pvdp1d, pvdp2d, pvdp3d )
+      INTEGER                    , INTENT(in)           ::   kt       ! ocean time-step
+      INTEGER                    , INTENT(in)           ::   kwrite   ! writing time-step
+      INTEGER                    , INTENT(in)           ::   kiomid   ! Identifier of the file
+      CHARACTER(len=*)           , INTENT(in)           ::   cdvar    ! time axis name
+      INTEGER                    , INTENT(in), OPTIONAL ::   ktype    ! variable external type
+      REAL(sp)                   , INTENT(in), OPTIONAL ::   pvsp0d    ! read field (0D case), single precision
+      REAL(sp), DIMENSION(:)     , INTENT(in), OPTIONAL ::   pvsp1d    ! read field (1D case), single precision
+      REAL(sp), DIMENSION(:,:)   , INTENT(in), OPTIONAL ::   pvsp2d    ! read field (2D case), single precision
+      REAL(sp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL ::   pvsp3d    ! read field (3D case), single precision
+      REAL(dp)                   , INTENT(in), OPTIONAL ::   pvdp0d    ! read field (0D case), double precision
+      REAL(dp), DIMENSION(:)     , INTENT(in), OPTIONAL ::   pvdp1d    ! read field (1D case), double precision
+      REAL(dp), DIMENSION(:,:)   , INTENT(in), OPTIONAL ::   pvdp2d    ! read field (2D case), double precision
+      REAL(dp), DIMENSION(:,:,:) , INTENT(in), OPTIONAL ::   pvdp3d    ! read field (3D case), double precision
+      !
+      LOGICAL           :: llx                 ! local xios write flag
+      INTEGER           :: ivid   ! variable id
+      CHARACTER(LEN=lc) :: context
+      !
       CALL set_xios_context(kiomid, context)
 
       llx = .NOT. (context == "NONE")
@@ -1910,14 +1809,22 @@ CONTAINS
       IF( llx ) THEN
 #ifdef key_xios
          IF( kt == kwrite ) THEN
-            IF(lwp) write(numout,*) 'RESTART: write (XIOS 3D) ',TRIM(cdvar)
+            IF(lwp) write(numout,*) 'RESTART: write (XIOS) ',TRIM(cdvar)
             CALL iom_swap(context)
-            CALL iom_put(TRIM(cdvar), pvar)
+            IF( PRESENT(pvsp0d) )   CALL iom_put(TRIM(cdvar), pvsp0d)
+            IF( PRESENT(pvsp1d) )   CALL iom_put(TRIM(cdvar), pvsp1d)
+            IF( PRESENT(pvsp2d) )   CALL iom_put(TRIM(cdvar), pvsp2d)
+            IF( PRESENT(pvsp3d) )   CALL iom_put(TRIM(cdvar), pvsp3d)
+            IF( PRESENT(pvdp0d) )   CALL iom_put(TRIM(cdvar), pvdp0d)
+            IF( PRESENT(pvdp1d) )   CALL iom_put(TRIM(cdvar), pvdp1d)
+            IF( PRESENT(pvdp2d) )   CALL iom_put(TRIM(cdvar), pvdp2d)
+            IF( PRESENT(pvdp3d) )   CALL iom_put(TRIM(cdvar), pvdp3d)
             CALL iom_swap(cxios_context)
          ELSE
-            IF(lwp) write(numout,*) 'RESTART: define (XIOS 3D)',TRIM(cdvar)
+            IF(lwp) write(numout,*) 'RESTART: define (XIOS)',TRIM(cdvar)
             CALL iom_swap(context)
-            CALL iom_set_rstw_active( TRIM(cdvar), rd3 = pvar )
+            CALL iom_set_rstw_active( TRIM(cdvar), rs0 = pvsp0d, rs1 = pvsp1d, rs2 = pvsp2d, rs3 = pvsp3d   &
+               &                                 , rd0 = pvdp0d, rd1 = pvdp1d, rd2 = pvdp2d, rd3 = pvdp3d )
             CALL iom_swap(cxios_context)
          ENDIF
 #endif
@@ -1925,12 +1832,12 @@ CONTAINS
          IF( kiomid > 0 ) THEN
             IF( iom_file(kiomid)%nfid > 0 ) THEN
                ivid = iom_varid( kiomid, cdvar, ldstop = .FALSE. )
-               CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pv_r3d = pvar )
+               CALL iom_nf90_rstput( kt, kwrite, kiomid, cdvar, ivid, ktype, pvsp0d, pvsp1d, pvsp2d, pvsp3d,   &
+                  &                                                          pvdp0d, pvdp1d, pvdp2d, pvdp3d )
             ENDIF
          ENDIF
       ENDIF
-   END SUBROUTINE iom_rp3d_dp
-
+   END SUBROUTINE iom_rp0123d
 
 
   SUBROUTINE iom_delay_rst( cdaction, cdcpnt, kncid )
diff --git a/src/OCE/IOM/iom_nf90.F90 b/src/OCE/IOM/iom_nf90.F90
index c49d9538e007c312f57e69b2fcb05786be34a42f..f6a156d737c0b9dc61afb4e29afac2b42a157208 100644
--- a/src/OCE/IOM/iom_nf90.F90
+++ b/src/OCE/IOM/iom_nf90.F90
@@ -34,10 +34,7 @@ MODULE iom_nf90
 
    INTERFACE iom_nf90_get
       MODULE PROCEDURE iom_nf90_g0d_sp
-      MODULE PROCEDURE iom_nf90_g0d_dp, iom_nf90_g123d_dp
-   END INTERFACE
-   INTERFACE iom_nf90_rstput
-      MODULE PROCEDURE iom_nf90_rp0123d_dp
+      MODULE PROCEDURE iom_nf90_g0d_dp, iom_nf90_g123d
    END INTERFACE
 
    !!----------------------------------------------------------------------
@@ -310,8 +307,8 @@ CONTAINS
       CALL iom_nf90_check(NF90_GET_VAR(iom_file(kiomid)%nfid, iom_file(kiomid)%nvid(kvid), pvar, start = kstart), clinfo )
    END SUBROUTINE iom_nf90_g0d_dp
 
-   SUBROUTINE iom_nf90_g123d_dp( kiomid, kvid, knbdim, kstart, kcount, kx1, kx2, ky1, ky2,   &
-         &                    pv_r1d, pv_r2d, pv_r3d )
+   SUBROUTINE iom_nf90_g123d( kiomid, kvid, knbdim, kstart, kcount, kx1, kx2, ky1, ky2,   &
+         &                    pvsp1d, pvsp2d, pvsp3d, pvdp1d, pvdp2d, pvdp3d )
       !!-----------------------------------------------------------------------
       !!                  ***  ROUTINE  iom_nf90_g123d  ***
       !!
@@ -325,9 +322,12 @@ CONTAINS
       INTEGER , DIMENSION(:)     , INTENT(in   )           ::   kstart    ! start position of the reading in each axis
       INTEGER , DIMENSION(:)     , INTENT(in   )           ::   kcount    ! number of points to be read in each axis
       INTEGER ,                    INTENT(in   )           ::   kx1, kx2, ky1, ky2   ! subdomain indexes
-      REAL(dp), DIMENSION(:)     , INTENT(  out), OPTIONAL ::   pv_r1d    ! read field (1D case)
-      REAL(dp), DIMENSION(:,:)   , INTENT(  out), OPTIONAL ::   pv_r2d    ! read field (2D case)
-      REAL(dp), DIMENSION(:,:,:) , INTENT(  out), OPTIONAL ::   pv_r3d    ! read field (3D case)
+      REAL(sp), DIMENSION(:)     , INTENT(  out), OPTIONAL ::   pvsp1d    ! read field (1D case), single precision
+      REAL(sp), DIMENSION(:,:)   , INTENT(  out), OPTIONAL ::   pvsp2d    ! read field (2D case), single precision
+      REAL(sp), DIMENSION(:,:,:) , INTENT(  out), OPTIONAL ::   pvsp3d    ! read field (3D case), single precision
+      REAL(dp), DIMENSION(:)     , INTENT(  out), OPTIONAL ::   pvdp1d    ! read field (1D case), double precision
+      REAL(dp), DIMENSION(:,:)   , INTENT(  out), OPTIONAL ::   pvdp2d    ! read field (2D case), double precision
+      REAL(dp), DIMENSION(:,:,:) , INTENT(  out), OPTIONAL ::   pvdp3d    ! read field (3D case), double precision
       !
       CHARACTER(LEN=100) ::   clinfo               ! info character
       INTEGER            ::   if90id               ! nf90 identifier of the opened file
@@ -337,18 +337,27 @@ CONTAINS
       if90id = iom_file(kiomid)%nfid         ! get back NetCDF file id
       ivid   = iom_file(kiomid)%nvid(kvid)   ! get back NetCDF var id
       !
-      IF(     PRESENT(pv_r1d) ) THEN
-         CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pv_r1d(:                ), start = kstart(1:knbdim),   &
+      IF(     PRESENT(pvsp1d) ) THEN
+         CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pvsp1d(:                ), start = kstart(1:knbdim),   &
+            &                                                                       count = kcount(1:knbdim)), clinfo )
+      ELSEIF( PRESENT(pvdp1d) ) THEN
+         CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pvdp1d(:                ), start = kstart(1:knbdim),   &
+            &                                                                       count = kcount(1:knbdim)), clinfo )
+      ELSEIF( PRESENT(pvsp2d) ) THEN
+         CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pvsp2d(kx1:kx2,ky1:ky2  ), start = kstart(1:knbdim),   &
+            &                                                                       count = kcount(1:knbdim)), clinfo )
+      ELSEIF( PRESENT(pvdp2d) ) THEN
+         CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pvdp2d(kx1:kx2,ky1:ky2  ), start = kstart(1:knbdim),   &
             &                                                                       count = kcount(1:knbdim)), clinfo )
-      ELSEIF( PRESENT(pv_r2d) ) THEN
-         CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pv_r2d(kx1:kx2,ky1:ky2  ), start = kstart(1:knbdim),   &
+      ELSEIF( PRESENT(pvsp3d) ) THEN
+         CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pvsp3d(kx1:kx2,ky1:ky2,:), start = kstart(1:knbdim),   &
             &                                                                       count = kcount(1:knbdim)), clinfo )
-      ELSEIF( PRESENT(pv_r3d) ) THEN
-         CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pv_r3d(kx1:kx2,ky1:ky2,:), start = kstart(1:knbdim),   &
+      ELSEIF( PRESENT(pvdp3d) ) THEN
+         CALL iom_nf90_check( NF90_GET_VAR(if90id, ivid, pvdp3d(kx1:kx2,ky1:ky2,:), start = kstart(1:knbdim),   &
             &                                                                       count = kcount(1:knbdim)), clinfo )
       ENDIF
       !
-   END SUBROUTINE iom_nf90_g123d_dp
+   END SUBROUTINE iom_nf90_g123d
 
 
 
@@ -522,8 +531,9 @@ CONTAINS
       !
    END SUBROUTINE iom_nf90_putatt
 
-   SUBROUTINE iom_nf90_rp0123d_dp( kt, kwrite, kiomid, cdvar , kvid  , ktype,   &
-         &                                  pv_r0d, pv_r1d, pv_r2d, pv_r3d )
+   SUBROUTINE iom_nf90_rstput( kt, kwrite, kiomid, cdvar , kvid  , ktype ,   &
+         &                                 pvsp0d, pvsp1d, pvsp2d, pvsp3d,   &
+         &                                 pvdp0d, pvdp1d, pvdp2d, pvdp3d )
       !!--------------------------------------------------------------------
       !!                   ***  SUBROUTINE  iom_nf90_rstput  ***
       !!
@@ -535,10 +545,14 @@ CONTAINS
       CHARACTER(len=*)            , INTENT(in)           ::   cdvar    ! variable name
       INTEGER                     , INTENT(in)           ::   kvid     ! variable id
       INTEGER                     , INTENT(in), OPTIONAL ::   ktype    ! variable type (default R8)
-      REAL(dp)                    , INTENT(in), OPTIONAL ::   pv_r0d   ! written Od field
-      REAL(dp), DIMENSION(      :), INTENT(in), OPTIONAL ::   pv_r1d   ! written 1d field
-      REAL(dp), DIMENSION(:, :   ), INTENT(in), OPTIONAL ::   pv_r2d   ! written 2d field
-      REAL(dp), DIMENSION(:, :, :), INTENT(in), OPTIONAL ::   pv_r3d   ! written 3d field
+      REAL(sp)                    , INTENT(in), OPTIONAL ::   pvsp0d   ! written Od field
+      REAL(sp), DIMENSION(      :), INTENT(in), OPTIONAL ::   pvsp1d   ! written 1d field
+      REAL(sp), DIMENSION(:, :   ), INTENT(in), OPTIONAL ::   pvsp2d   ! written 2d field
+      REAL(sp), DIMENSION(:, :, :), INTENT(in), OPTIONAL ::   pvsp3d   ! written 3d field
+      REAL(dp)                    , INTENT(in), OPTIONAL ::   pvdp0d   ! written Od field
+      REAL(dp), DIMENSION(      :), INTENT(in), OPTIONAL ::   pvdp1d   ! written 1d field
+      REAL(dp), DIMENSION(:, :   ), INTENT(in), OPTIONAL ::   pvdp2d   ! written 2d field
+      REAL(dp), DIMENSION(:, :, :), INTENT(in), OPTIONAL ::   pvdp3d   ! written 3d field
       !
       INTEGER               :: idims                ! number of dimension
       INTEGER               :: idvar                ! variable id
@@ -553,12 +567,18 @@ CONTAINS
       !                                             ! nn_nchunks_[i,j,k,t] namelist parameters
       INTEGER               :: ichunkalg, ishuffle, ideflate, ideflate_level
       !                                             ! NetCDF4 internally fixed parameters
+      INTEGER               :: idlv                 ! local variable
       LOGICAL               :: lchunk               ! logical switch to activate chunking and compression
       !                                             ! when appropriate (currently chunking is applied to 4d fields only)
-      INTEGER               :: idlv                 ! local variable
+      LOGICAL               :: llis0d, llis1d, llis2d, llis3d
       CHARACTER(LEN=256)    :: ccname               ! local variable
       !---------------------------------------------------------------------
       !
+      llis0d = PRESENT(pvsp0d) .OR. PRESENT(pvdp0d)
+      llis1d = PRESENT(pvsp1d) .OR. PRESENT(pvdp1d)
+      llis2d = PRESENT(pvsp2d) .OR. PRESENT(pvdp2d)
+      llis3d = PRESENT(pvsp3d) .OR. PRESENT(pvdp3d)
+      !
       clinfo = '          iom_nf90_rp0123d, file: '//TRIM(iom_file(kiomid)%name)//', var: '//TRIM(cdvar)
       if90id = iom_file(kiomid)%nfid
       !
@@ -607,12 +627,12 @@ CONTAINS
             CALL iom_nf90_check(NF90_REDEF( if90id ), clinfo)   ;   iom_file(kiomid)%irec = -1
          ENDIF
          ! variable definition
-         IF(     PRESENT(pv_r0d) ) THEN   ;   idims = 0
-         ELSEIF( PRESENT(pv_r1d) ) THEN
-                                              idims = 2   ;   idimid(1:idims) = (/3,4/)
-         ELSEIF( PRESENT(pv_r2d) ) THEN   ;   idims = 3   ;   idimid(1:idims) = (/1,2,4/)
-         ELSEIF( PRESENT(pv_r3d) ) THEN
-                                              idims = 4   ;   idimid(1:idims) = (/1,2,3,4/)
+         IF(     llis0d ) THEN   ;   idims = 0
+         ELSEIF( llis1d ) THEN
+                                     idims = 2   ;   idimid(1:idims) = (/3,4/)
+         ELSEIF( llis2d ) THEN   ;   idims = 3   ;   idimid(1:idims) = (/1,2,4/)
+         ELSEIF( llis3d ) THEN
+                                     idims = 4   ;   idimid(1:idims) = (/1,2,3,4/)
          ENDIF
          IF( PRESENT(ktype) ) THEN   ! variable external type
             SELECT CASE (ktype)
@@ -626,7 +646,7 @@ CONTAINS
          ELSE
             itype = NF90_DOUBLE
          ENDIF
-         IF( PRESENT(pv_r0d) ) THEN
+         IF( llis0d ) THEN
             CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype,                    &
                &                              iom_file(kiomid)%nvid(idvar) ), clinfo )
          ELSE
@@ -641,9 +661,7 @@ CONTAINS
          iom_file(kiomid)%scf(idvar)    = 1.
          iom_file(kiomid)%ofs(idvar)    = 0.
          iom_file(kiomid)%ndims(idvar)  = idims
-         IF( .NOT. PRESENT(pv_r0d) ) THEN   ;   iom_file(kiomid)%luld(idvar) = .TRUE.
-         ELSE                               ;   iom_file(kiomid)%luld(idvar) = .FALSE.
-         ENDIF
+         iom_file(kiomid)%luld(idvar)   = .NOT. llis0d
          DO jd = 1, idims
             CALL iom_nf90_check(NF90_INQUIRE_DIMENSION( if90id, idimid(jd), len = iom_file(kiomid)%dimsz(jd,idvar) ), clinfo)
             IF ( lchunk ) ichunksz(jd) = iom_file(kiomid)%dimsz(jd,idvar)
@@ -672,7 +690,7 @@ CONTAINS
             CALL iom_nf90_check(NF90_ENDDEF( if90id ), clinfo)   ;   iom_file(kiomid)%irec = 0
          ENDIF
          ! on what kind of domain must the data be written?
-         IF( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) THEN
+         IF( llis2d .OR. llis3d ) THEN
             idimsz(1:2) = iom_file(kiomid)%dimsz(1:2,idvar)
             IF(     idimsz(1) == Ni_0 .AND. idimsz(2) == Nj_0 ) THEN
                ix1 = Nis0   ;   ix2 = Nie0   ;   iy1 = Njs0   ;   iy2 = Nje0
@@ -708,14 +726,22 @@ CONTAINS
 
          ! write the data
          ! =============
-         IF(     PRESENT(pv_r0d) ) THEN
-            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r0d                    ), clinfo )
-         ELSEIF( PRESENT(pv_r1d) ) THEN
-            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r1d(:)                 ), clinfo )
-         ELSEIF( PRESENT(pv_r2d) ) THEN
-            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r2d(ix1:ix2,iy1:iy2)   ), clinfo )
-         ELSEIF( PRESENT(pv_r3d) ) THEN
-            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r3d(ix1:ix2,iy1:iy2,:) ), clinfo )
+         IF(     PRESENT(pvsp0d) ) THEN
+            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvsp0d                    ), clinfo )
+         ELSEIF( PRESENT(pvdp0d) ) THEN
+            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvdp0d                    ), clinfo )
+         ELSEIF( PRESENT(pvsp1d) ) THEN
+            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvsp1d(:)                 ), clinfo )
+         ELSEIF( PRESENT(pvdp1d) ) THEN
+            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvdp1d(:)                 ), clinfo )
+         ELSEIF( PRESENT(pvsp2d) ) THEN
+            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvsp2d(ix1:ix2,iy1:iy2)   ), clinfo )
+         ELSEIF( PRESENT(pvdp2d) ) THEN
+            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvdp2d(ix1:ix2,iy1:iy2)   ), clinfo )
+         ELSEIF( PRESENT(pvsp3d) ) THEN
+            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvsp3d(ix1:ix2,iy1:iy2,:) ), clinfo )
+         ELSEIF( PRESENT(pvdp3d) ) THEN
+            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pvdp3d(ix1:ix2,iy1:iy2,:) ), clinfo )
          ENDIF
          ! add 1 to the size of the temporal dimension (not really useful...)
          IF( iom_file(kiomid)%luld(idvar) )   iom_file(kiomid)%dimsz(iom_file(kiomid)%ndims(idvar), idvar)    &
@@ -723,7 +749,7 @@ CONTAINS
          IF(lwp) WRITE(numout,*) TRIM(clinfo)//' written ok'
       ENDIF
       !
-   END SUBROUTINE iom_nf90_rp0123d_dp
+   END SUBROUTINE iom_nf90_rstput
 
 
    SUBROUTINE iom_nf90_check( kstatus, cdinfo )
diff --git a/src/OCE/LBC/lib_mpp.F90 b/src/OCE/LBC/lib_mpp.F90
index 0bcde50b8b6b0f5ec11e9d5d80bdee7fed706e1d..81a100a1bd5ba6316299fa8418653387bd1f7744 100644
--- a/src/OCE/LBC/lib_mpp.F90
+++ b/src/OCE/LBC/lib_mpp.F90
@@ -1237,7 +1237,7 @@ CONTAINS
                 + aimag(ydda(ji)) + aimag(yddb(ji))
 
          ! The result is zt1 + zt2, after normalization.
-         yddb(ji) = cmplx ( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1),wp )
+         yddb(ji) = cmplx ( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), dp )
       END DO
       !
    END SUBROUTINE DDPDD_MPI
diff --git a/src/OCE/LDF/ldfslp.F90 b/src/OCE/LDF/ldfslp.F90
index 4a85f6991d3fdfa7cb16b5cad538b18191e43d7d..1b5b09ea5a0ae83a5205a5cac4a56a7b9f87f492 100644
--- a/src/OCE/LDF/ldfslp.F90
+++ b/src/OCE/LDF/ldfslp.F90
@@ -726,7 +726,7 @@ CONTAINS
 !                  END DO
 !               END DO
 !            END DO
-!            CALL lbc_lnk( 'ldfslp', uslp , 'U', -1. ; CALL lbc_lnk( 'ldfslp', vslp , 'V', -1.,  wslpi, 'W', -1.,  wslpj, 'W', -1. )
+!            CALL lbc_lnk( 'ldfslp', uslp , 'U', -1. ; CALL lbc_lnk( 'ldfslp', vslp , 'V', -1._wp,  wslpi, 'W', -1._wp,  wslpj, 'W', -1._wp )
 !!gm         ENDIF
       ENDIF
       !
diff --git a/src/OCE/SBC/sbcblk_algo_ecmwf.F90 b/src/OCE/SBC/sbcblk_algo_ecmwf.F90
index d86c997565a355c17a2eee34b0f6b1ab72c06b92..e4075308d8b329c03c61e516479c6323938221b5 100644
--- a/src/OCE/SBC/sbcblk_algo_ecmwf.F90
+++ b/src/OCE/SBC/sbcblk_algo_ecmwf.F90
@@ -417,27 +417,23 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) :: zta, zx2, zx, ztmp, zpsi_unst, zpsi_stab, zstab, zc
+      REAL(wp) :: zta, zx2, zx, ztmp, zc
       !!----------------------------------------------------------------------------------
       zc = 5._wp/0.35_wp
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             !
             zta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
-
-            ! *** Unstable (Paulson 1970)    [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
-            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 - 16z)^0.5
-            zx  = SQRT(zx2)                          ! (1 - 16z)^0.25
-            ztmp = 1._wp + zx
-            zpsi_unst = LOG( 0.125_wp*ztmp*ztmp*(1._wp + zx2) ) - 2._wp*ATAN( zx ) + 0.5_wp*rpi
-
-            ! *** Stable                   [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
-            zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) &
-               &       - zta - 2._wp/3._wp*zc
             !
-            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
-            !
-            psi_m_ecmwf(ji,jj) =         zstab  * zpsi_stab &  ! (zta > 0) Stable
-               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable
+            IF (zta < 0.0_wp) THEN
+               ! *** Unstable (Paulson 1970)    [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
+               zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 - 16z)^0.5
+               zx  = SQRT(zx2)                        ! (1 - 16z)^0.25
+               ztmp = 1._wp + zx
+               psi_m_ecmwf(ji,jj) = LOG(0.125_wp*ztmp*ztmp*(1._wp + zx2)) - 2._wp*ATAN(zx) + 0.5_wp*rpi
+            ELSE
+               ! *** Stable                   [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
+               psi_m_ecmwf(ji,jj) = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) - zta - 2._wp/3._wp*zc
+            ENDIF
             !
       END_2D
    END FUNCTION psi_m_ecmwf
@@ -458,7 +454,7 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
       !
       INTEGER  ::   ji, jj     ! dummy loop indices
-      REAL(wp) ::  zta, zx2, zpsi_unst, zpsi_stab, zstab, zc
+      REAL(wp) ::  zta, zx2, zc
       !!----------------------------------------------------------------------------------
       zc = 5._wp/0.35_wp
       !
@@ -466,20 +462,15 @@ CONTAINS
             !
             zta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!):
             !
-            ! *** Unstable (Paulson 1970)   [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
-            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 -16z)^0.5
-            zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) )
-            !
-            ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
-            zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) &
+            IF (zta < 0.0_wp) THEN
+               ! *** Unstable (Paulson 1970)   [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
+               zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 -16z)^0.5
+               psi_h_ecmwf(ji,jj) = 2._wp*LOG( 0.5_wp*(1._wp + zx2) )
+            ELSE
+               ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
+               psi_h_ecmwf(ji,jj) = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) &
                &       - ABS(1._wp + 2._wp/3._wp*zta)**1.5_wp - 2._wp/3._wp*zc + 1._wp
-            !
-            ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
-            !
-            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
-            !
-            psi_h_ecmwf(ji,jj) =         zstab  * zpsi_stab &  ! (zta > 0) Stable
-               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable
+            ENDIF
             !
       END_2D
    END FUNCTION psi_h_ecmwf
diff --git a/src/OCE/SBC/sbccpl.F90 b/src/OCE/SBC/sbccpl.F90
index eddf626a522f3479f787e155a6dc76ed02900634..79602758f20f23e2ab645ed1e224474d44d0c1aa 100644
--- a/src/OCE/SBC/sbccpl.F90
+++ b/src/OCE/SBC/sbccpl.F90
@@ -1672,7 +1672,7 @@ CONTAINS
                p_taui(ji,jj) = zztmp1 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
                p_tauj(ji,jj) = zztmp2 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
             END_2D
-            CALL lbc_lnk( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. )
+            CALL lbc_lnk( 'sbccpl', p_taui, 'U',  -1._wp, p_tauj, 'V',  -1._wp )
          END SELECT
 
       ENDIF
diff --git a/src/OCE/SBC/sbcice_cice.F90 b/src/OCE/SBC/sbcice_cice.F90
index 5f7f75770b3ec439806d3c4bfa2a4727652eaae5..07e7c902d014e64660f3056099c23597997983d7 100644
--- a/src/OCE/SBC/sbcice_cice.F90
+++ b/src/OCE/SBC/sbcice_cice.F90
@@ -175,7 +175,7 @@ CONTAINS
       ! there is no restart file.
       ! Values from a CICE restart file would overwrite this
       IF( .NOT. ln_rstart ) THEN    
-         CALL nemo2cice( ts(:,:,1,jp_tem,Kmm) , sst , 'T' , 1.) 
+         CALL nemo2cice( ts(:,:,1,jp_tem,Kmm) , sst , 'T' , 1._wp) 
       ENDIF  
 #endif
 
@@ -207,10 +207,10 @@ CONTAINS
       fr_iu(:,:)=0.0
       fr_iv(:,:)=0.0
 
-      CALL cice2nemo(aice,fr_i, 'T', 1. )
+      CALL cice2nemo(aice,fr_i, 'T', 1._wp )
       IF( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
          DO jl=1,ncat
-            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
+            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1._wp )
          ENDDO
       ENDIF
 
@@ -224,8 +224,8 @@ CONTAINS
       CALL lbc_lnk( 'sbcice_cice', fr_iu , 'U', 1.0_wp,  fr_iv , 'V', 1.0_wp )
 
       ! set the snow+ice mass
-      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
-      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
+      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1._wp )
+      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1._wp )
       snwice_mass  (:,:) = ( rhos * ztmp1(:,:) + rhoi * ztmp2(:,:)  )
       snwice_mass_b(:,:) = snwice_mass(:,:)
 
@@ -315,7 +315,7 @@ CONTAINS
             ztmp(ji,jj) = 0.5 * (  fr_iu(ji,jj) * utau(ji,jj)      &
                                  + fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
          END_2D
-         CALL nemo2cice(ztmp,strax,'F', -1. )
+         CALL nemo2cice(ztmp,strax,'F', -1._wp )
 
 ! y comp of wind stress (CI_2)
 ! V point to F point
@@ -323,7 +323,7 @@ CONTAINS
             ztmp(ji,jj) = 0.5 * (  fr_iv(ji,jj) * vtau(ji,jj)      &
                                  + fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
          END_2D
-         CALL nemo2cice(ztmp,stray,'F', -1. )
+         CALL nemo2cice(ztmp,stray,'F', -1._wp )
 
 ! Surface downward latent heat flux (CI_5)
          IF(ksbc == jp_flx) THEN
@@ -349,7 +349,7 @@ CONTAINS
             END_2D
          ENDIF
          DO jl=1,ncat
-            CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
+            CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1._wp )
 
 ! GBM conductive flux through ice (CI_6)
 !  Convert to GBM
@@ -358,7 +358,7 @@ CONTAINS
             ELSE
                ztmp(:,:) = botmelt(:,:,jl)
             ENDIF
-            CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
+            CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1._wp )
 
 ! GBM surface heat flux (CI_7)
 !  Convert to GBM
@@ -367,7 +367,7 @@ CONTAINS
             ELSE
                ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
             ENDIF
-            CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
+            CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1._wp )
          ENDDO
 
       ELSE IF(ksbc == jp_blk) THEN
@@ -375,39 +375,39 @@ CONTAINS
 ! Pass bulk forcing fields to CICE (which will calculate heat fluxes etc itself)
 ! x comp and y comp of atmosphere surface wind (CICE expects on T points)
          ztmp(:,:) = wndi_ice(:,:)
-         CALL nemo2cice(ztmp,uatm,'T', -1. )
+         CALL nemo2cice(ztmp,uatm,'T', -1._wp )
          ztmp(:,:) = wndj_ice(:,:)
-         CALL nemo2cice(ztmp,vatm,'T', -1. )
+         CALL nemo2cice(ztmp,vatm,'T', -1._wp )
          ztmp(:,:) = SQRT ( wndi_ice(:,:)**2 + wndj_ice(:,:)**2 )
-         CALL nemo2cice(ztmp,wind,'T', 1. )    ! Wind speed (m/s)
+         CALL nemo2cice(ztmp,wind,'T', 1._wp )    ! Wind speed (m/s)
          ztmp(:,:) = qsr_ice(:,:,1)
-         CALL nemo2cice(ztmp,fsw,'T', 1. )     ! Incoming short-wave (W/m^2)
+         CALL nemo2cice(ztmp,fsw,'T', 1._wp )     ! Incoming short-wave (W/m^2)
          ztmp(:,:) = qlw_ice(:,:,1)
-         CALL nemo2cice(ztmp,flw,'T', 1. )     ! Incoming long-wave (W/m^2)
+         CALL nemo2cice(ztmp,flw,'T', 1._wp )     ! Incoming long-wave (W/m^2)
          ztmp(:,:) = tatm_ice(:,:)
-         CALL nemo2cice(ztmp,Tair,'T', 1. )    ! Air temperature (K)
-         CALL nemo2cice(ztmp,potT,'T', 1. )    ! Potential temp (K)
+         CALL nemo2cice(ztmp,Tair,'T', 1._wp )    ! Air temperature (K)
+         CALL nemo2cice(ztmp,potT,'T', 1._wp )    ! Potential temp (K)
 ! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows  
          ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )    
                                                ! Constant (101000.) atm pressure assumed
-         CALL nemo2cice(ztmp,rhoa,'T', 1. )    ! Air density (kg/m^3)
+         CALL nemo2cice(ztmp,rhoa,'T', 1._wp )    ! Air density (kg/m^3)
          ztmp(:,:) = qatm_ice(:,:)
-         CALL nemo2cice(ztmp,Qa,'T', 1. )      ! Specific humidity (kg/kg)
+         CALL nemo2cice(ztmp,Qa,'T', 1._wp )      ! Specific humidity (kg/kg)
          ztmp(:,:)=10.0
-         CALL nemo2cice(ztmp,zlvl,'T', 1. )    ! Atmos level height (m)
+         CALL nemo2cice(ztmp,zlvl,'T', 1._wp )    ! Atmos level height (m)
 
 ! May want to check all values are physically realistic (as in CICE routine 
 ! prepare_forcing)?
 
 ! Divide shortwave into spectral bands (as in prepare_forcing)
          ztmp(:,:)=qsr_ice(:,:,1)*frcvdr       ! visible direct
-         CALL nemo2cice(ztmp,swvdr,'T', 1. )             
+         CALL nemo2cice(ztmp,swvdr,'T', 1._wp )             
          ztmp(:,:)=qsr_ice(:,:,1)*frcvdf       ! visible diffuse
-         CALL nemo2cice(ztmp,swvdf,'T', 1. )              
+         CALL nemo2cice(ztmp,swvdf,'T', 1._wp )              
          ztmp(:,:)=qsr_ice(:,:,1)*frcidr       ! near IR direct
-         CALL nemo2cice(ztmp,swidr,'T', 1. )
+         CALL nemo2cice(ztmp,swidr,'T', 1._wp )
          ztmp(:,:)=qsr_ice(:,:,1)*frcidf       ! near IR diffuse
-         CALL nemo2cice(ztmp,swidf,'T', 1. )
+         CALL nemo2cice(ztmp,swidf,'T', 1._wp )
 
       ENDIF
 
@@ -415,37 +415,37 @@ CONTAINS
 ! Ensure fsnow is positive (as in CICE routine prepare_forcing)
       IF( iom_use('snowpre') )   CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit  
       ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)  
-      CALL nemo2cice(ztmp,fsnow,'T', 1. ) 
+      CALL nemo2cice(ztmp,fsnow,'T', 1._wp ) 
 
 ! Rainfall
       IF( iom_use('precip') )   CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
       ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
-      CALL nemo2cice(ztmp,frain,'T', 1. ) 
+      CALL nemo2cice(ztmp,frain,'T', 1._wp ) 
 
 ! Freezing/melting potential
 ! Calculated over NEMO leapfrog timestep (hence 2*dt)
       nfrzmlt(:,:) = rho0 * rcp * e3t_m(:,:) * ( Tocnfrz-sst_m(:,:) ) / ( 2.0*dt )
 
       ztmp(:,:) = nfrzmlt(:,:)
-      CALL nemo2cice(ztmp,frzmlt,'T', 1. )
+      CALL nemo2cice(ztmp,frzmlt,'T', 1._wp )
 
 ! SST  and SSS
 
-      CALL nemo2cice(sst_m,sst,'T', 1. )
-      CALL nemo2cice(sss_m,sss,'T', 1. )
+      CALL nemo2cice(sst_m,sst,'T', 1._wp )
+      CALL nemo2cice(sss_m,sss,'T', 1._wp )
 
 ! x comp and y comp of surface ocean current
 ! U point to F point
       DO_2D( 1, 1, 1, 0 )
          ztmp(ji,jj)=0.5*(ssu_m(ji,jj)+ssu_m(ji,jj+1))*fmask(ji,jj,1)
       END_2D
-      CALL nemo2cice(ztmp,uocn,'F', -1. )
+      CALL nemo2cice(ztmp,uocn,'F', -1._wp )
 
 ! V point to F point
       DO_2D( 1, 0, 1, 1 )
          ztmp(ji,jj)=0.5*(ssv_m(ji,jj)+ssv_m(ji+1,jj))*fmask(ji,jj,1)
       END_2D
-      CALL nemo2cice(ztmp,vocn,'F', -1. )
+      CALL nemo2cice(ztmp,vocn,'F', -1._wp )
 
       IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
           !
@@ -470,14 +470,14 @@ CONTAINS
          ztmp(ji,jj)=0.5 * (  (zpice(ji+1,jj  )-zpice(ji,jj  )) * r1_e1u(ji,jj  )    &
             &               + (zpice(ji+1,jj+1)-zpice(ji,jj+1)) * r1_e1u(ji,jj+1)  ) * fmask(ji,jj,1)
       END_2D
-      CALL nemo2cice( ztmp,ss_tltx,'F', -1. )
+      CALL nemo2cice( ztmp,ss_tltx,'F', -1._wp )
 
 ! T point to F point
       DO_2D( 1, 0, 1, 0 )
          ztmp(ji,jj)=0.5 * (  (zpice(ji  ,jj+1)-zpice(ji  ,jj)) * r1_e2v(ji  ,jj)    &
             &               + (zpice(ji+1,jj+1)-zpice(ji+1,jj)) * r1_e2v(ji+1,jj)  ) *  fmask(ji,jj,1)
       END_2D
-      CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
+      CALL nemo2cice(ztmp,ss_tlty,'F', -1._wp )
       !
    END SUBROUTINE cice_sbc_in
 
@@ -499,7 +499,7 @@ CONTAINS
       ENDIF
       
 ! x comp of ocean-ice stress 
-      CALL cice2nemo(strocnx,ztmp1,'F', -1. )
+      CALL cice2nemo(strocnx,ztmp1,'F', -1._wp )
       ss_iou(:,:)=0.0
 ! F point to U point
       DO_2D( 0, 0, 0, 0 )
@@ -508,7 +508,7 @@ CONTAINS
       CALL lbc_lnk( 'sbcice_cice', ss_iou , 'U', -1.0_wp )
 
 ! y comp of ocean-ice stress 
-      CALL cice2nemo(strocny,ztmp1,'F', -1. )
+      CALL cice2nemo(strocny,ztmp1,'F', -1._wp )
       ss_iov(:,:)=0.0
 ! F point to V point
 
@@ -527,8 +527,8 @@ CONTAINS
  
 ! Also need ice/ocean stress on T points so that taum can be updated 
 ! This interpolation is already done in CICE so best to use those values 
-      CALL cice2nemo(strocnxT,ztmp1,'T',-1.) 
-      CALL cice2nemo(strocnyT,ztmp2,'T',-1.) 
+      CALL cice2nemo(strocnxT,ztmp1,'T',-1._wp) 
+      CALL cice2nemo(strocnyT,ztmp2,'T',-1._wp) 
  
 ! Update taum with modulus of ice-ocean stress 
 ! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here 
@@ -551,11 +551,11 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
       ENDIF
 
 #if defined key_cice4
-      CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
-      CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
+      CALL cice2nemo(fresh_gbm,ztmp1,'T', 1._wp )
+      CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1._wp )
 #else
-      CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
-      CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
+      CALL cice2nemo(fresh_ai,ztmp1,'T', 1._wp )
+      CALL cice2nemo(fsalt_ai,ztmp2,'T', 1._wp )
 #endif
 
 ! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
@@ -589,9 +589,9 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
 ! Now add in ice / snow related terms
 ! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
 #if defined key_cice4
-      CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
+      CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1._wp )
 #else
-      CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
+      CALL cice2nemo(fswthru_ai,ztmp1,'T', 1._wp )
 #endif
       qsr(:,:)=qsr(:,:)+ztmp1(:,:)
       CALL lbc_lnk( 'sbcice_cice', qsr , 'T', 1.0_wp )
@@ -601,9 +601,9 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
       END_2D
 
 #if defined key_cice4
-      CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
+      CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1._wp )
 #else
-      CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
+      CALL cice2nemo(fhocn_ai,ztmp1,'T', 1._wp )
 #endif
       qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
 
@@ -611,10 +611,10 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
 
 ! Prepare for the following CICE time-step
 
-      CALL cice2nemo(aice,fr_i,'T', 1. )
+      CALL cice2nemo(aice,fr_i,'T', 1._wp )
       IF( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
          DO jl=1,ncat
-            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
+            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1._wp )
          ENDDO
       ENDIF
 
@@ -628,8 +628,8 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
       CALL lbc_lnk( 'sbcice_cice', fr_iu , 'U', 1.0_wp, fr_iv , 'V', 1.0_wp )
 
       ! set the snow+ice mass
-      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
-      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
+      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1._wp )
+      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1._wp )
       snwice_mass  (:,:) = ( rhos * ztmp1(:,:) + rhoi * ztmp2(:,:)  )
       snwice_mass_b(:,:) = snwice_mass(:,:)
       snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
@@ -661,16 +661,16 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
       !
       ! x and y comp of ice velocity
       !
-      CALL cice2nemo(uvel,u_ice,'F', -1. )
-      CALL cice2nemo(vvel,v_ice,'F', -1. )
+      CALL cice2nemo(uvel,u_ice,'F', -1._wp )
+      CALL cice2nemo(vvel,v_ice,'F', -1._wp )
       !
       ! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out  
       !
       ! Snow and ice thicknesses (CO_2 and CO_3)
       !
       DO jl = 1, ncat
-         CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1. )
-         CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1. )
+         CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1._wp )
+         CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1._wp )
       END DO
       !
    END SUBROUTINE cice_sbc_hadgam
diff --git a/src/OCE/SBC/sbcmod.F90 b/src/OCE/SBC/sbcmod.F90
index 29e30de16382c58de2dc09329f82415da1706be3..473c68e977d9d91debf4c2b51ed6583e918b6c77 100644
--- a/src/OCE/SBC/sbcmod.F90
+++ b/src/OCE/SBC/sbcmod.F90
@@ -443,8 +443,8 @@ CONTAINS
             vtau(ji,jj) = vtau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji,jj-1) ) * 0.5_wp
          END_2D
          !
-         CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
-         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
+         CALL lbc_lnk( 'sbcwave', utau, 'U', -1._wp )
+         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1._wp )
          !
          taum(:,:) = taum(:,:)*tauoc_wave(:,:)
          !
@@ -454,8 +454,8 @@ CONTAINS
       ELSEIF( ln_wave .AND. ln_taw ) THEN                  ! Wave stress reduction
          utau(:,:) = utau(:,:) - tawx(:,:) + twox(:,:)
          vtau(:,:) = vtau(:,:) - tawy(:,:) + twoy(:,:)
-         CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
-         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
+         CALL lbc_lnk( 'sbcwave', utau, 'U', -1._wp )
+         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1._wp )
          !
          DO_2D( 0, 0, 0, 0)
              taum(ji,jj) = sqrt((.5*(utau(ji-1,jj)+utau(ji,jj)))**2 + (.5*(vtau(ji,jj-1)+vtau(ji,jj)))**2)
@@ -465,7 +465,7 @@ CONTAINS
             &                                'If not requested select ln_taw=.false.' )
          !
       ENDIF
-      CALL lbc_lnk( 'sbcmod', taum(:,:), 'T', 1. )
+      CALL lbc_lnk( 'sbcmod', taum(:,:), 'T', 1._wp )
       !
       IF( ln_icebergs ) THEN  ! save pure stresses (with no ice-ocean stress) for use by icebergs
          utau_icb(:,:) = utau(:,:) ; vtau_icb(:,:) = vtau(:,:) 
diff --git a/src/OCE/TRA/traadv_cen.F90 b/src/OCE/TRA/traadv_cen.F90
index 3351e1935a5503648b2e9e0ca590bb5a7813c037..8f42bb87a14d2bb9654b8182b85900260358679e 100644
--- a/src/OCE/TRA/traadv_cen.F90
+++ b/src/OCE/TRA/traadv_cen.F90
@@ -136,7 +136,7 @@ CONTAINS
                zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u
                zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v
             END_3D
-            IF (nn_hls==1) CALL lbc_lnk( 'traadv_cen', zwx, 'U', -1. , zwy, 'V', -1. )
+            IF (nn_hls==1) CALL lbc_lnk( 'traadv_cen', zwx, 'U', -1._wp , zwy, 'V', -1._wp )
             !
          CASE DEFAULT
             CALL ctl_stop( 'traadv_cen: wrong value for nn_cen' )
diff --git a/src/OCE/TRA/traadv_fct.F90 b/src/OCE/TRA/traadv_fct.F90
index 500a46d80f8dadc42df2911f9f0e3814acf81b61..bfc80b893c91c2349d8ae39911195b062df5c6cb 100644
--- a/src/OCE/TRA/traadv_fct.F90
+++ b/src/OCE/TRA/traadv_fct.F90
@@ -405,23 +405,25 @@ CONTAINS
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ikm1         ! local integer
-      REAL(dp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
-      REAL(dp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
-      REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zbetup, zbetdo, zbup, zbdo
+      REAL(wp) ::   zpos, zneg, zbt, zbig                 ! local scalars
+      REAL(wp) ::   zup, zdo                              !   -      -
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zbetup, zbetdo, zbup, zbdo
       !!----------------------------------------------------------------------
       !
-      zbig  = 1.e+40_dp
-      zrtrn = 1.e-15_dp
-      zbetup(:,:,:) = 0._dp   ;   zbetdo(:,:,:) = 0._dp
+      zbig = HUGE(1._wp)
+      zbetup(:,:,jpk) = zbig   ;   zbetdo(:,:,jpk) = zbig
 
       ! Search local extrema
       ! --------------------
       ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
-         zbup(ji,jj,jk) = MAX( pbef(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ),   &
-            &                  paft(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) )  )
-         zbdo(ji,jj,jk) = MIN( pbef(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ),   &
-            &                  paft(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) )  )
+         IF( tmask(ji,jj,jk) == 1._wp ) THEN
+            zbup(ji,jj,jk) = MAX( pbef(ji,jj,jk), paft(ji,jj,jk) )
+            zbdo(ji,jj,jk) = MIN( pbef(ji,jj,jk), paft(ji,jj,jk) )
+         ELSE
+            zbup(ji,jj,jk) = -zbig
+            zbdo(ji,jj,jk) =  zbig
+         ENDIF
       END_3D
 
       DO jk = 1, jpkm1
@@ -441,42 +443,39 @@ CONTAINS
                &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  )
 
             ! positive part of the flux
-            zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
-               & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   &
-               & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
+            zpos = MAX( 0._wp, paa(ji-1,jj  ,jk  ) ) - MIN( 0._wp, paa(ji  ,jj  ,jk  ) )   &
+               & + MAX( 0._wp, pbb(ji  ,jj-1,jk  ) ) - MIN( 0._wp, pbb(ji  ,jj  ,jk  ) )   &
+               & + MAX( 0._wp, pcc(ji  ,jj  ,jk+1) ) - MIN( 0._wp, pcc(ji  ,jj  ,jk  ) )
 
             ! negative part of the flux
-            zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
-               & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   &
-               & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
+            zneg = MAX( 0._wp, paa(ji  ,jj  ,jk  ) ) - MIN( 0._wp, paa(ji-1,jj  ,jk  ) )   &
+               & + MAX( 0._wp, pbb(ji  ,jj  ,jk  ) ) - MIN( 0._wp, pbb(ji  ,jj-1,jk  ) )   &
+               & + MAX( 0._wp, pcc(ji  ,jj  ,jk  ) ) - MIN( 0._wp, pcc(ji  ,jj  ,jk+1) )
 
             ! up & down beta terms
             zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt
-            zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
-            zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
+            IF( zup /= -zbig .AND. zpos /= 0._wp ) THEN   ;   zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / zpos * zbt
+            ELSE                                          ;   zbetup(ji,jj,jk) = zbig
+            ENDIF
+            IF( zdo /=  zbig .AND. zneg /= 0._wp ) THEN   ;   zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / zneg * zbt
+            ELSE                                          ;   zbetdo(ji,jj,jk) = zbig
+            ENDIF
          END_2D
       END DO
       IF (nn_hls==1) CALL lbc_lnk( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp, ld4only= .TRUE. )   ! lateral boundary cond. (unchanged sign)
 
-      ! 3. monotonic flux in the i & j direction (paa & pbb)
-      ! ----------------------------------------
+      ! 3. monotonic flux in the i, j and k direction (paa, pbb and pcc)
+      ! ----------------------------------------------------------------
       DO_3D( 1, 0, 1, 0, 1, jpkm1 )
-         zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
-         zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
-         zcu =       ( 0.5  + SIGN( 0.5_wp , paa(ji,jj,jk) ) )
-         paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
-
-         zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
-         zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
-         zcv =       ( 0.5  + SIGN( 0.5_wp , pbb(ji,jj,jk) ) )
-         pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
-
-      ! monotonic flux in the k direction, i.e. pcc
-      ! -------------------------------------------
-         za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
-         zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
-         zc =       ( 0.5  + SIGN( 0.5_wp , pcc(ji,jj,jk+1) ) )
-         pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
+         IF( paa(ji,jj,jk) > 0._wp ) THEN   ;   paa(ji,jj,jk) = paa(ji,jj,jk) * MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
+         ELSE                               ;   paa(ji,jj,jk) = paa(ji,jj,jk) * MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
+         ENDIF
+         IF( pbb(ji,jj,jk) > 0._wp ) THEN   ;   pbb(ji,jj,jk) = pbb(ji,jj,jk) * MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
+         ELSE                               ;   pbb(ji,jj,jk) = pbb(ji,jj,jk) * MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
+         ENDIF
+         IF( pcc(ji,jj,jk+1) > 0._wp ) THEN   ;   pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * MIN( 1._wp, zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
+         ELSE                                 ;   pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * MIN( 1._wp, zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
+         ENDIF
       END_3D
       !
    END SUBROUTINE nonosc
diff --git a/src/OCE/TRA/traadv_ubs.F90 b/src/OCE/TRA/traadv_ubs.F90
index 6d65af219fcdad959749963efcf98ad20d3c42e1..6a254d775c23c7277794f72b0553d4bbf93c1f4f 100644
--- a/src/OCE/TRA/traadv_ubs.F90
+++ b/src/OCE/TRA/traadv_ubs.F90
@@ -294,7 +294,7 @@ CONTAINS
       REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zbetup, zbetdo         ! 3D workspace
       !!----------------------------------------------------------------------
       !
-      zbig  = 1.e+38_wp
+      zbig  = 1.e+20_wp   ! works ok with simple/double precison
       zrtrn = 1.e-15_wp
       zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
       !
diff --git a/src/OCE/TRA/traadv_ubs_lf.F90 b/src/OCE/TRA/traadv_ubs_lf.F90
index ee8acfe2922ef27a36be5c6a0372ae0ea083b20a..034512641488d99d8ebf4e78638b1ca792941a84 100644
--- a/src/OCE/TRA/traadv_ubs_lf.F90
+++ b/src/OCE/TRA/traadv_ubs_lf.F90
@@ -299,7 +299,7 @@ CONTAINS
       REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zbetup, zbetdo         ! 3D workspace
       !!----------------------------------------------------------------------
       !
-      zbig  = 1.e+38_wp
+      zbig  = 1.e+20_wp   ! works ok with simple/double precison
       zrtrn = 1.e-15_wp
       zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
       !
diff --git a/src/OCE/TRA/traqsr.F90 b/src/OCE/TRA/traqsr.F90
index 5f9c8e94e44f98c29f0cbe76fa1d50aa42332f65..0bde7adce5679affb5d610bbf82171cc27ea23d9 100644
--- a/src/OCE/TRA/traqsr.F90
+++ b/src/OCE/TRA/traqsr.F90
@@ -110,7 +110,6 @@ CONTAINS
       INTEGER  ::   irgb                     ! local integers
       REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars
       REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         -
-      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
       REAL(wp) ::   zz0 , zz1 , ze3t, zlui   !    -         -
       REAL(wp) ::   zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze
       REAL(wp) ::   zlogc, zlogze, zlogCtot, zlogCze
@@ -287,8 +286,9 @@ CONTAINS
       !
       ! sea-ice: store the 1st ocean level attenuation coefficient
       DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
-         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
-         ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
+         zz0 = r1_rho0_rcp * qsr(ji,jj)   ! test zz0 and not qsr for rounding errors in single precision
+         IF( zz0 /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / zz0
+         ELSE                      ;   fraqsr_1lev(ji,jj) = 1._wp
          ENDIF
       END_2D
       !
diff --git a/src/OCE/ZDF/zdfddm.F90 b/src/OCE/ZDF/zdfddm.F90
index a08cf3439b254cff6f6864f7eee3f7052aed3ae6..4cfc6ead2360d0b6ed209be7e13c8c347bfb7596 100644
--- a/src/OCE/ZDF/zdfddm.F90
+++ b/src/OCE/ZDF/zdfddm.F90
@@ -142,7 +142,7 @@ CONTAINS
             ! salt fingering
             zrr = zrau(ji,jj) / rn_hsbfr
             zrr = zrr * zrr
-            zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
+            zavfs = rn_avts / ( 1. + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
             zavft = 0.7 * zavfs * zinr
             ! diffusive layering
             zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
diff --git a/src/OCE/ZDF/zdftke.F90 b/src/OCE/ZDF/zdftke.F90
index 8b6d915587952429fdf3237b4bb8d2718b7fa85a..5680ec1e20176e01de75677426a40208fdf571f2 100644
--- a/src/OCE/ZDF/zdftke.F90
+++ b/src/OCE/ZDF/zdftke.F90
@@ -219,6 +219,7 @@ CONTAINS
       REAL(wp), DIMENSION(A2D(nn_hls))     ::   zice_fra, zhlc, zus3, zWlc2
       REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zpelc, zdiag, zd_up, zd_lw
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztmp ! for diags
+      REAL(wp) :: zdiv
       !!--------------------------------------------------------------------
       !
       zbbrau  = rn_ebb / rho0       ! Local constant initialisation
@@ -370,7 +371,16 @@ CONTAINS
             IF (rn2b(ji,jj,jk) <= 0.0_wp) then
                 zri = 0.0_wp
             ELSE
-                zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
+                ! This logic is to avoid divide-by-zero errors which can occur for single-precision
+                ! The actual value you choose for the denominator instead of zero doesn't really
+                ! matter, as long as it is very small and so triggers the same logic below with the
+                ! inverse Prandtl number
+                zdiv = p_sh2(ji,jj,jk) + rn_bshear
+                IF (zdiv == 0.0_wp) THEN
+                   zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / rn_bshear
+                ELSE
+                   zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / zdiv
+                ENDIF
             ENDIF
             !                             ! inverse of Prandtl number
             apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
diff --git a/src/OCE/lib_fortran.F90 b/src/OCE/lib_fortran.F90
index b621373e6eb7b7356db171bbc204f30793c73063..9929defb47e3b42cb645dae8d2488fff8d1e0f97 100644
--- a/src/OCE/lib_fortran.F90
+++ b/src/OCE/lib_fortran.F90
@@ -113,8 +113,8 @@ CONTAINS
 
    FUNCTION local_sum_2d( ptab )
       !!----------------------------------------------------------------------
-      REAL(wp),  INTENT(in   ) ::   ptab(:,:) ! array on which operation is applied
-      COMPLEX(dp)              ::  local_sum_2d
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   ptab ! array on which operation is applied
+      COMPLEX(dp)                          ::   local_sum_2d
       !
       !!-----------------------------------------------------------------------
       !
@@ -127,12 +127,12 @@ CONTAINS
       ipi = SIZE(ptab,1)   ! 1st dimension
       ipj = SIZE(ptab,2)   ! 2nd dimension
       !
-      ctmp = CMPLX( 0.e0, 0.e0, wp )   ! warning ctmp is cumulated
+      ctmp = CMPLX( 0._dp, 0._dp, dp )   ! warning ctmp is cumulated
 
       DO jj = 1, ipj
          DO ji = 1, ipi
             ztmp =  ptab(ji,jj) * tmask_i(ji,jj)
-            CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp )
+            CALL DDPDD( CMPLX( ztmp, 0._dp, dp ), ctmp )
          END DO
       END DO
       !
@@ -142,8 +142,8 @@ CONTAINS
 
    FUNCTION local_sum_3d( ptab )
       !!----------------------------------------------------------------------
-      REAL(wp),  INTENT(in   ) ::   ptab(:,:,:) ! array on which operation is applied
-      COMPLEX(dp)              ::  local_sum_3d
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   ptab ! array on which operation is applied
+      COMPLEX(dp)                            ::   local_sum_3d
       !
       !!-----------------------------------------------------------------------
       !
@@ -157,13 +157,13 @@ CONTAINS
       ipj = SIZE(ptab,2)   ! 2nd dimension
       ipk = SIZE(ptab,3)   ! 3rd dimension
       !
-      ctmp = CMPLX( 0.e0, 0.e0, wp )   ! warning ctmp is cumulated
+      ctmp = CMPLX( 0._dp, 0._dp, dp )   ! warning ctmp is cumulated
 
       DO jk = 1, ipk
         DO jj = 1, ipj
           DO ji = 1, ipi
              ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
-             CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp )
+             CALL DDPDD( CMPLX( ztmp, 0._dp, dp ), ctmp )
           END DO
         END DO
       END DO
@@ -286,9 +286,9 @@ CONTAINS
 
    FUNCTION glob_sum_vec_3d( cdname, ptab ) RESULT( ptmp )
       !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in) ::   cdname      ! name of the calling subroutine
-      REAL(wp),          INTENT(in) ::   ptab(:,:,:) ! array on which operation is applied
-      REAL(wp), DIMENSION(SIZE(ptab,3)) ::   ptmp
+      CHARACTER(len=*),           INTENT(in) ::   cdname   ! name of the calling subroutine
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   ptab     ! array on which operation is applied
+      REAL(wp), DIMENSION(SIZE(ptab,3))      ::   ptmp
       !
       COMPLEX(dp), DIMENSION(:), ALLOCATABLE ::   ctmp
       REAL(wp)    ::   ztmp
@@ -312,11 +312,11 @@ CONTAINS
       ALLOCATE( ctmp(ipk) )
       !
       DO jk = 1, ipk
-         ctmp(jk) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
+         ctmp(jk) = CMPLX( 0._dp, 0._dp, dp )   ! warning ctmp is cumulated
          DO jj = ijs, ije
             DO ji = iis, iie
                ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
-               CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jk) )
+               CALL DDPDD( CMPLX( ztmp, 0._dp, dp ), ctmp(jk) )
             END DO
          END DO
       END DO
@@ -330,9 +330,9 @@ CONTAINS
 
    FUNCTION glob_sum_vec_4d( cdname, ptab ) RESULT( ptmp )
       !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in) ::   cdname        ! name of the calling subroutine
-      REAL(wp),          INTENT(in) ::   ptab(:,:,:,:) ! array on which operation is applied
-      REAL(wp), DIMENSION(SIZE(ptab,4)) ::   ptmp
+      CHARACTER(len=*),             INTENT(in) ::   cdname   ! name of the calling subroutine
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(in) ::   ptab     ! array on which operation is applied
+      REAL(wp), DIMENSION(SIZE(ptab,4))        ::   ptmp
       !
       COMPLEX(dp), DIMENSION(:), ALLOCATABLE ::   ctmp
       REAL(wp)    ::   ztmp
@@ -362,7 +362,7 @@ CONTAINS
             DO jj = ijs, ije
                DO ji = iis, iie
                   ztmp =  ptab(ji,jj,jk,jl) * tmask_i(ji,jj)
-                  CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jl) )
+                  CALL DDPDD( CMPLX( ztmp, 0._dp, dp ), ctmp(jl) )
                END DO
             END DO
          END DO
@@ -377,9 +377,9 @@ CONTAINS
 
    FUNCTION glob_min_vec_3d( cdname, ptab ) RESULT( ptmp )
       !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in) ::   cdname        ! name of the calling subroutine
-      REAL(wp),          INTENT(in) ::   ptab(:,:,:)   ! array on which operation is applied
-      REAL(wp), DIMENSION(SIZE(ptab,3)) ::   ptmp
+      CHARACTER(len=*),           INTENT(in) ::   cdname   ! name of the calling subroutine
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   ptab     ! array on which operation is applied
+      REAL(wp), DIMENSION(SIZE(ptab,3))      ::   ptmp
       !
       INTEGER     ::   jk    ! dummy loop indice & dimension
       INTEGER     ::   ipk   ! dimension
@@ -396,9 +396,9 @@ CONTAINS
 
    FUNCTION glob_min_vec_4d( cdname, ptab ) RESULT( ptmp )
       !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in) ::   cdname          ! name of the calling subroutine
-      REAL(wp),          INTENT(in) ::   ptab(:,:,:,:)   ! array on which operation is applied
-      REAL(wp), DIMENSION(SIZE(ptab,4)) ::   ptmp
+      CHARACTER(len=*),             INTENT(in) ::   cdname   ! name of the calling subroutine
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(in) ::   ptab     ! array on which operation is applied
+      REAL(wp), DIMENSION(SIZE(ptab,4))        ::   ptmp
       !
       INTEGER     ::   jk , jl    ! dummy loop indice & dimension
       INTEGER     ::   ipk, ipl   ! dimension
@@ -419,9 +419,9 @@ CONTAINS
    
    FUNCTION glob_max_vec_3d( cdname, ptab ) RESULT( ptmp )
       !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in) ::   cdname        ! name of the calling subroutine
-      REAL(wp),          INTENT(in) ::   ptab(:,:,:)   ! array on which operation is applied
-      REAL(wp), DIMENSION(SIZE(ptab,3)) ::   ptmp
+      CHARACTER(len=*),           INTENT(in) ::   cdname   ! name of the calling subroutine
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   ptab     ! array on which operation is applied
+      REAL(wp), DIMENSION(SIZE(ptab,3))      ::   ptmp
       !
       INTEGER     ::   jk    ! dummy loop indice & dimension
       INTEGER     ::   ipk   ! dimension
@@ -438,9 +438,9 @@ CONTAINS
 
    FUNCTION glob_max_vec_4d( cdname, ptab ) RESULT( ptmp )
       !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in) ::   cdname          ! name of the calling subroutine
-      REAL(wp),          INTENT(in) ::   ptab(:,:,:,:)   ! array on which operation is applied
-      REAL(wp), DIMENSION(SIZE(ptab,4)) ::   ptmp
+      CHARACTER(len=*),             INTENT(in) ::   cdname   ! name of the calling subroutine
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(in) ::   ptab     ! array on which operation is applied
+      REAL(wp), DIMENSION(SIZE(ptab,4))        ::   ptmp
       !
       INTEGER     ::   jk , jl    ! dummy loop indice & dimension
       INTEGER     ::   ipk, ipl   ! dimension
@@ -489,7 +489,7 @@ CONTAINS
          &   + AIMAG(ydda)         + AIMAG(yddb)
       !
       ! The result is t1 + t2, after normalization.
-      yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), wp )
+      yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), dp )
       !
    END SUBROUTINE DDPDD
 
diff --git a/src/OCE/lib_fortran_generic.h90 b/src/OCE/lib_fortran_generic.h90
index 5201ad6d014dda8259deaeca0042ab9f79a9a581..e09d498522c066fc408a7d5e6d1c22a0aa3cc270 100644
--- a/src/OCE/lib_fortran_generic.h90
+++ b/src/OCE/lib_fortran_generic.h90
@@ -55,12 +55,12 @@
          ijs = 1   ;   ije = jpj
       ENDIF
       !
-      ctmp = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
+      ctmp = CMPLX( 0._dp, 0._dp, dp )   ! warning ctmp is cumulated
       DO jk = 1, ipk
         DO jj = ijs, ije
           DO ji = iis, iie
              ztmp =  ARRAY_IN(ji,jj,jk) * MASK_ARRAY(ji,jj)
-             CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp )
+             CALL DDPDD( CMPLX( ztmp, 0._wp, dp ), ctmp )
           END DO
         END DO
       END DO
@@ -110,7 +110,6 @@
       !
       !!-----------------------------------------------------------------------
       !
-      COMPLEX(dp)::   ctmp
       REAL(wp)   ::   ztmp
       INTEGER    ::   jk       ! dummy loop indices
       INTEGER    ::   ipk      ! dimensions
diff --git a/src/OCE/module_example.F90 b/src/OCE/module_example.F90
index 80c10ce05a5e6e3f1cad08c4502abea7878e6f78..76a25caeb4c6be4e1ecf2535d3ab383e5a059f5c 100644
--- a/src/OCE/module_example.F90
+++ b/src/OCE/module_example.F90
@@ -126,7 +126,7 @@ CONTAINS
       !
       ! WARNING! the lbc_lnk call could not be compatible with the tiling approach
       ! please refer to the manual for how to adapt your code
-      CALL lbc_lnk( 'module_example', avm, 'T', 1. )     ! Lateral boundary conditions (unchanged sign)
+      CALL lbc_lnk( 'module_example', avm, 'T', 1._wp )     ! Lateral boundary conditions (unchanged sign)
       !
    END SUBROUTINE exa_mpl
 
diff --git a/src/OCE/stpmlf.F90 b/src/OCE/stpmlf.F90
index 81952fc6aa5532f132d891dfba2187d1c7627d8d..21bea529fe1880adb237ed2e837645c189559e56 100644
--- a/src/OCE/stpmlf.F90
+++ b/src/OCE/stpmlf.F90
@@ -546,8 +546,8 @@ CONTAINS
             CALL Agrif_dyn( kt )
 # endif
       !                                        ! local domain boundaries  (T-point, unchanged sign)
-      CALL lbc_lnk( 'finalize_lbc', puu(:,:,:,       Kaa), 'U', -1., pvv(:,:,:       ,Kaa), 'V', -1.   &
-                       &          , pts(:,:,:,jp_tem,Kaa), 'T',  1., pts(:,:,:,jp_sal,Kaa), 'T',  1. )
+      CALL lbc_lnk( 'finalize_lbc', puu(:,:,:,       Kaa), 'U', -1._wp, pvv(:,:,:       ,Kaa), 'V', -1._wp   &
+                       &          , pts(:,:,:,jp_tem,Kaa), 'T',  1._wp, pts(:,:,:,jp_sal,Kaa), 'T',  1._wp )
       !
       ! lbc_lnk needed for zdf_sh2 when using nn_hls = 2, moved here to allow tiling in zdf_phy
       IF( nn_hls == 2 .AND. l_zdfsh2 ) CALL lbc_lnk( 'stp', avm_k, 'W', 1.0_wp )
diff --git a/src/OCE/timing.F90 b/src/OCE/timing.F90
index c63003d8b66e5ffb901c471502339be6d98c2eeb..4c1becfe14fdbb9d97515b9e884015d8842ccfa4 100644
--- a/src/OCE/timing.F90
+++ b/src/OCE/timing.F90
@@ -39,7 +39,7 @@ MODULE timing
       CHARACTER(LEN=20)  :: cname
       CHARACTER(LEN=20)  :: surname
       INTEGER :: rank
-      REAL(wp)  :: t_cpu, t_clock, tsum_cpu, tsum_clock, tmax_cpu, tmax_clock, tmin_cpu, tmin_clock, tsub_cpu, tsub_clock
+      REAL(dp)  :: t_cpu, t_clock, tsum_cpu, tsum_clock, tmax_cpu, tmax_clock, tmin_cpu, tmin_clock, tsub_cpu, tsub_clock
       INTEGER :: ncount, ncount_max, ncount_rate
       INTEGER :: niter
       LOGICAL :: l_tdone
@@ -50,8 +50,8 @@ MODULE timing
 
    TYPE alltimer
       CHARACTER(LEN=20), DIMENSION(:), POINTER :: cname => NULL()
-      REAL(wp), DIMENSION(:), POINTER :: tsum_cpu   => NULL()
-      REAL(wp), DIMENSION(:), POINTER :: tsum_clock => NULL()
+      REAL(dp), DIMENSION(:), POINTER :: tsum_cpu   => NULL()
+      REAL(dp), DIMENSION(:), POINTER :: tsum_clock => NULL()
       INTEGER, DIMENSION(:), POINTER :: niter => NULL()
       TYPE(alltimer), POINTER :: next => NULL()
       TYPE(alltimer), POINTER :: prev => NULL()
@@ -62,14 +62,14 @@ MODULE timing
    TYPE(timer), POINTER :: s_timer_old      => NULL()
 
    TYPE(timer), POINTER :: s_wrk        => NULL()
-   REAL(wp) :: t_overclock, t_overcpu
+   REAL(dp) :: t_overclock, t_overcpu
    LOGICAL :: l_initdone = .FALSE.
    INTEGER :: nsize
 
    ! Variables for coarse grain timing
-   REAL(wp) :: tot_etime, tot_ctime
-   REAL(kind=wp), DIMENSION(2)     :: t_elaps, t_cpu
-   REAL(wp), ALLOCATABLE, DIMENSION(:) :: all_etime, all_ctime
+   REAL(dp) :: tot_etime, tot_ctime
+   REAL(kind=dp), DIMENSION(2)     :: t_elaps, t_cpu
+   REAL(dp), ALLOCATABLE, DIMENSION(:) :: all_etime, all_ctime
    INTEGER :: nfinal_count, ncount, ncount_rate, ncount_max
    INTEGER, DIMENSION(8)           :: nvalues
    CHARACTER(LEN=8), DIMENSION(2)  :: cdate
@@ -136,7 +136,7 @@ CONTAINS
       CHARACTER(len=*), INTENT(in), OPTIONAL :: csection
       !
       INTEGER  :: ifinal_count, iperiods
-      REAL(wp) :: zcpu_end, zmpitime,zcpu_raw,zclock_raw
+      REAL(dp) :: zcpu_end, zmpitime,zcpu_raw,zclock_raw
       !
       s_wrk => NULL()
 
@@ -218,7 +218,7 @@ CONTAINS
       !! ** Purpose :   open timing output file
       !!----------------------------------------------------------------------
       INTEGER :: iperiods, istart_count, ifinal_count
-      REAL(wp) :: zdum
+      REAL(dp) :: zdum
       LOGICAL :: ll_f
       CHARACTER(len=*), INTENT(in), OPTIONAL :: clname
       CHARACTER(len=20)                      :: cln
@@ -295,9 +295,9 @@ CONTAINS
       INTEGER :: ji
       LOGICAL :: ll_ord, ll_averep
       CHARACTER(len=120) :: clfmt
-      REAL(wp), DIMENSION(:), ALLOCATABLE ::   timing_glob
-      REAL(wp) ::   zsypd   ! simulated years per day (Balaji 2017)
-      REAL(wp) ::   zperc, ztot
+      REAL(dp), DIMENSION(:), ALLOCATABLE ::   timing_glob
+      REAL(dp) ::   zsypd   ! simulated years per day (Balaji 2017)
+      REAL(dp) ::   zperc, ztot
 
       ll_averep = .TRUE.
 
@@ -387,16 +387,16 @@ CONTAINS
          WRITE(numtime,*) '    warning: includes restarts writing time if output before nitend... '
          WRITE(numtime,*) ' '
          DO ji = 1, jpnij
-            zperc = 0._wp ; zsypd = 0._wp
+            zperc = 0._dp ; zsypd = 0._dp
             ztot = SUM( timing_glob(4*ji-3:4*ji-1) )
             WRITE(numtime,'(A28,F11.6,            A34,I8)') 'Computing       time : ',timing_glob(4*ji-3), ' on MPI rank : ', ji
-            IF ( ztot /= 0._wp ) zperc = timing_glob(4*ji-2) / ztot * 100.
+            IF ( ztot /= 0._dp ) zperc = timing_glob(4*ji-2) / ztot * 100.
             WRITE(numtime,'(A28,F11.6,A2, F4.1,A3,A25,I8)') 'Waiting lbc_lnk time : ',timing_glob(4*ji-2)   &
                &                                                         , ' (',      zperc,' %)',   ' on MPI rank : ', ji
-            IF ( ztot /= 0._wp ) zperc = timing_glob(4*ji-1) / ztot * 100.
+            IF ( ztot /= 0._dp ) zperc = timing_glob(4*ji-1) / ztot * 100.
             WRITE(numtime,'(A28,F11.6,A2, F4.1,A3,A25,I8)') 'Waiting  global time : ',timing_glob(4*ji-1)   &
                &                                                         , ' (',      zperc,' %)',   ' on MPI rank : ', ji
-            IF ( timing_glob(4*ji) /= 0._wp ) zsypd = rn_Dt * REAL(nitend-nit000-1, wp) / (timing_glob(4*ji) * 365.)
+            IF ( timing_glob(4*ji) /= 0._dp ) zsypd = rn_Dt * REAL(nitend-nit000-1, dp) / (timing_glob(4*ji) * 365.)
             WRITE(numtime,'(A28,F11.6,A7,F10.3,A2,A15,I8)') 'Total           time : ',timing_glob(4*ji  )   &
                &                                                         , ' (SYPD: ', zsypd, ')',   ' on MPI rank : ', ji
          END DO
@@ -447,7 +447,7 @@ CONTAINS
       s_timer => s_timer_root
       clfmt = '(1x,a,4x,f12.3,6x,f12.3,x,f12.3,2x,f12.3,6x,f7.3,2x,i9)'
       DO WHILE ( ASSOCIATED(s_timer) )
-         IF( s_timer%tsum_clock > 0._wp )                                &
+         IF( s_timer%tsum_clock > 0._dp )                                &
             WRITE(numtime,TRIM(clfmt))   s_timer%cname,                  &
             &   s_timer%tsum_clock,s_timer%tsum_clock*100./t_elaps(2),   &
             &   s_timer%tsum_cpu  ,s_timer%tsum_cpu*100./t_cpu(2)    ,   &
@@ -489,8 +489,8 @@ CONTAINS
          RETURN
       END IF
       sl_timer_glob_root%cname(:)       = ''
-      sl_timer_glob_root%tsum_cpu(:)   = 0._wp
-      sl_timer_glob_root%tsum_clock(:) = 0._wp
+      sl_timer_glob_root%tsum_cpu(:)   = 0._dp
+      sl_timer_glob_root%tsum_clock(:) = 0._dp
       sl_timer_glob_root%niter(:)      = 0
       sl_timer_glob_root%next => NULL()
       sl_timer_glob_root%prev => NULL()
@@ -506,16 +506,16 @@ CONTAINS
       IF( narea .EQ. 1 ) THEN
          ALLOCATE(sl_timer_ave_root)
          sl_timer_ave_root%cname       = ''
-         sl_timer_ave_root%t_cpu      = 0._wp
-         sl_timer_ave_root%t_clock    = 0._wp
-         sl_timer_ave_root%tsum_cpu   = 0._wp
-         sl_timer_ave_root%tsum_clock = 0._wp
-         sl_timer_ave_root%tmax_cpu   = 0._wp
-         sl_timer_ave_root%tmax_clock = 0._wp
-         sl_timer_ave_root%tmin_cpu   = 0._wp
-         sl_timer_ave_root%tmin_clock = 0._wp
-         sl_timer_ave_root%tsub_cpu   = 0._wp
-         sl_timer_ave_root%tsub_clock = 0._wp
+         sl_timer_ave_root%t_cpu      = 0._dp
+         sl_timer_ave_root%t_clock    = 0._dp
+         sl_timer_ave_root%tsum_cpu   = 0._dp
+         sl_timer_ave_root%tsum_clock = 0._dp
+         sl_timer_ave_root%tmax_cpu   = 0._dp
+         sl_timer_ave_root%tmax_clock = 0._dp
+         sl_timer_ave_root%tmin_cpu   = 0._dp
+         sl_timer_ave_root%tmin_clock = 0._dp
+         sl_timer_ave_root%tsub_cpu   = 0._dp
+         sl_timer_ave_root%tsub_clock = 0._dp
          sl_timer_ave_root%ncount      = 0
          sl_timer_ave_root%ncount_rate = 0
          sl_timer_ave_root%ncount_max  = 0
@@ -642,10 +642,10 @@ CONTAINS
       !
       INTEGER                            :: idum, icode
       INTEGER, ALLOCATABLE, DIMENSION(:) :: iall_rank
-      REAL(wp) :: ztot_ratio
-      REAL(wp) :: zmax_etime, zmax_ctime, zmax_ratio, zmin_etime, zmin_ctime, zmin_ratio
-      REAL(wp) :: zavg_etime, zavg_ctime, zavg_ratio
-      REAL(wp), ALLOCATABLE, DIMENSION(:) :: zall_ratio
+      REAL(dp) :: ztot_ratio
+      REAL(dp) :: zmax_etime, zmax_ctime, zmax_ratio, zmin_etime, zmin_ctime, zmin_ratio
+      REAL(dp) :: zavg_etime, zavg_ctime, zavg_ratio
+      REAL(dp), ALLOCATABLE, DIMENSION(:) :: zall_ratio
       CHARACTER(LEN=128), dimension(8) :: cllignes
       CHARACTER(LEN=128)               :: clhline, clstart_date, clfinal_date
       CHARACTER(LEN=2048)              :: clfmt
@@ -656,19 +656,19 @@ CONTAINS
          iall_rank(:) = (/ (idum,idum=0,jpnij-1) /)
 
          ! Compute elapse user time
-         zavg_etime = tot_etime/REAL(jpnij,wp)
+         zavg_etime = tot_etime/REAL(jpnij,dp)
          zmax_etime = MAXVAL(all_etime(:))
          zmin_etime = MINVAL(all_etime(:))
 
          ! Compute CPU user time
-         zavg_ctime = tot_ctime/REAL(jpnij,wp)
+         zavg_ctime = tot_ctime/REAL(jpnij,dp)
          zmax_ctime = MAXVAL(all_ctime(:))
          zmin_ctime = MINVAL(all_ctime(:))
 
          ! Compute cpu/elapsed ratio
          zall_ratio(:) = all_ctime(:) / all_etime(:)
          ztot_ratio    = SUM(all_ctime(:))/SUM(all_etime(:))
-         zavg_ratio    = SUM(zall_ratio(:))/REAL(jpnij,wp)
+         zavg_ratio    = SUM(zall_ratio(:))/REAL(jpnij,dp)
          zmax_ratio    = MAXVAL(zall_ratio(:))
          zmin_ratio    = MINVAL(zall_ratio(:))
 
@@ -714,16 +714,16 @@ CONTAINS
       IF( .NOT. ASSOCIATED(s_timer_root) ) THEN
          ALLOCATE(s_timer_root)
          s_timer_root%cname       = cdinfo
-         s_timer_root%t_cpu      = 0._wp
-         s_timer_root%t_clock    = 0._wp
-         s_timer_root%tsum_cpu   = 0._wp
-         s_timer_root%tsum_clock = 0._wp
-         s_timer_root%tmax_cpu   = 0._wp
-         s_timer_root%tmax_clock = 0._wp
-         s_timer_root%tmin_cpu   = 0._wp
-         s_timer_root%tmin_clock = 0._wp
-         s_timer_root%tsub_cpu   = 0._wp
-         s_timer_root%tsub_clock = 0._wp
+         s_timer_root%t_cpu      = 0._dp
+         s_timer_root%t_clock    = 0._dp
+         s_timer_root%tsum_cpu   = 0._dp
+         s_timer_root%tsum_clock = 0._dp
+         s_timer_root%tmax_cpu   = 0._dp
+         s_timer_root%tmax_clock = 0._dp
+         s_timer_root%tmin_cpu   = 0._dp
+         s_timer_root%tmin_clock = 0._dp
+         s_timer_root%tsub_cpu   = 0._dp
+         s_timer_root%tsub_clock = 0._dp
          s_timer_root%ncount      = 0
          s_timer_root%ncount_rate = 0
          s_timer_root%ncount_max  = 0
@@ -738,16 +738,16 @@ CONTAINS
          !
          ALLOCATE(s_timer_old)
          s_timer_old%cname       = cdinfo
-         s_timer_old%t_cpu      = 0._wp
-         s_timer_old%t_clock    = 0._wp
-         s_timer_old%tsum_cpu   = 0._wp
-         s_timer_old%tsum_clock = 0._wp
-         s_timer_old%tmax_cpu   = 0._wp
-         s_timer_old%tmax_clock = 0._wp
-         s_timer_old%tmin_cpu   = 0._wp
-         s_timer_old%tmin_clock = 0._wp
-         s_timer_old%tsub_cpu   = 0._wp
-         s_timer_old%tsub_clock = 0._wp
+         s_timer_old%t_cpu      = 0._dp
+         s_timer_old%t_clock    = 0._dp
+         s_timer_old%tsum_cpu   = 0._dp
+         s_timer_old%tsum_clock = 0._dp
+         s_timer_old%tmax_cpu   = 0._dp
+         s_timer_old%tmax_clock = 0._dp
+         s_timer_old%tmin_cpu   = 0._dp
+         s_timer_old%tmin_clock = 0._dp
+         s_timer_old%tsub_cpu   = 0._dp
+         s_timer_old%tsub_clock = 0._dp
          s_timer_old%ncount      = 0
          s_timer_old%ncount_rate = 0
          s_timer_old%ncount_max  = 0
@@ -780,16 +780,16 @@ CONTAINS
     !     write(*,*) 'after allocation of next'
 
          s_timer%next%cname       = cdinfo
-         s_timer%next%t_cpu      = 0._wp
-         s_timer%next%t_clock    = 0._wp
-         s_timer%next%tsum_cpu   = 0._wp
-         s_timer%next%tsum_clock = 0._wp
-         s_timer%next%tmax_cpu   = 0._wp
-         s_timer%next%tmax_clock = 0._wp
-         s_timer%next%tmin_cpu   = 0._wp
-         s_timer%next%tmin_clock = 0._wp
-         s_timer%next%tsub_cpu   = 0._wp
-         s_timer%next%tsub_clock = 0._wp
+         s_timer%next%t_cpu      = 0._dp
+         s_timer%next%t_clock    = 0._dp
+         s_timer%next%tsum_cpu   = 0._dp
+         s_timer%next%tsum_clock = 0._dp
+         s_timer%next%tmax_cpu   = 0._dp
+         s_timer%next%tmax_clock = 0._dp
+         s_timer%next%tmin_cpu   = 0._dp
+         s_timer%next%tmin_clock = 0._dp
+         s_timer%next%tsub_cpu   = 0._dp
+         s_timer%next%tsub_clock = 0._dp
          s_timer%next%ncount      = 0
          s_timer%next%ncount_rate = 0
          s_timer%next%ncount_max  = 0
diff --git a/src/SWE/stpmlf.F90 b/src/SWE/stpmlf.F90
index 5494248c2da3d7d583c2377f3228ce1d73cbae4c..5fca79cb9e66596b8ca32c20f06d804f15d7e03f 100644
--- a/src/SWE/stpmlf.F90
+++ b/src/SWE/stpmlf.F90
@@ -196,9 +196,9 @@ CONTAINS
          ENDIF
       ENDIF
 
-      CALL lbc_lnk( 'stp_MLF', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
-         &                     uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
-      IF (nn_hls==2) CALL lbc_lnk( 'stp_MLF', r3u(:,:,Naa), 'U', 1., r3v(:,:,Naa), 'V', 1.)
+      CALL lbc_lnk( 'stp_MLF', uu(:,:,:,Nnn), 'U', -1._wp, vv(:,:,:,Nnn), 'V', -1._wp,   &   !* local domain boundaries
+         &                     uu(:,:,:,Naa), 'U', -1._wp, vv(:,:,:,Naa), 'V', -1._wp    )     
+      IF (nn_hls==2) CALL lbc_lnk( 'stp_MLF', r3u(:,:,Naa), 'U', 1._wp, r3v(:,:,Naa), 'V', 1._wp)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Set boundary conditions, time filter and swap time levels
diff --git a/src/SWE/stprk3.F90 b/src/SWE/stprk3.F90
index 0e1a1910fc971a67837aafea9430c53efa272d2e..c7e1ee5306056f56eb84d1bc63b9f05ec1e1d9a0 100644
--- a/src/SWE/stprk3.F90
+++ b/src/SWE/stprk3.F90
@@ -170,8 +170,8 @@ CONTAINS
          END_3D
       ENDIF
       !
-      CALL lbc_lnk( 'stp_RK3', uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
-      IF (nn_hls==2) CALL lbc_lnk( 'stp_RK3', r3u(:,:,Naa), 'U', 1., r3v(:,:,Naa), 'V', 1.)
+      CALL lbc_lnk( 'stp_RK3', uu(:,:,:,Naa), 'U', -1._wp, vv(:,:,:,Naa), 'V', -1._wp )
+      IF (nn_hls==2) CALL lbc_lnk( 'stp_RK3', r3u(:,:,Naa), 'U', 1._wp, r3v(:,:,Naa), 'V', 1._wp)
       !
       !                                 !==  Swap time levels  ==!
       Nrhs= Nnn
@@ -236,8 +236,8 @@ CONTAINS
          END_3D
       ENDIF
       !
-      CALL lbc_lnk( 'stp_RK3', uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
-      IF (nn_hls==2) CALL lbc_lnk( 'stp_RK3', r3u(:,:,Naa), 'U', 1., r3v(:,:,Naa), 'V', 1.)
+      CALL lbc_lnk( 'stp_RK3', uu(:,:,:,Naa), 'U', -1._wp, vv(:,:,:,Naa), 'V', -1._wp )
+      IF (nn_hls==2) CALL lbc_lnk( 'stp_RK3', r3u(:,:,Naa), 'U', 1._wp, r3v(:,:,Naa), 'V', 1._wp)
       !
       !                                 !==  Swap time levels  ==!
       Nrhs= Nnn
@@ -300,8 +300,8 @@ CONTAINS
          END_3D
       ENDIF
       !
-      CALL lbc_lnk( 'stp_RK3', uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
-      IF (nn_hls==2) CALL lbc_lnk( 'stp_RK3', r3u(:,:,Naa), 'U', 1., r3v(:,:,Naa), 'V', 1.)
+      CALL lbc_lnk( 'stp_RK3', uu(:,:,:,Naa), 'U', -1._wp, vv(:,:,:,Naa), 'V', -1._wp )
+      IF (nn_hls==2) CALL lbc_lnk( 'stp_RK3', r3u(:,:,Naa), 'U', 1._wp, r3v(:,:,Naa), 'V', 1._wp)
       !
       !                                 !==  Swap time levels  ==!
       !
diff --git a/src/TOP/TRP/trdmxl_trc.F90 b/src/TOP/TRP/trdmxl_trc.F90
index 1525d1d72618af2c559c1f126a17fac215c45ebd..59374e4876ca11b90676de62e049852a595e1a7d 100644
--- a/src/TOP/TRP/trdmxl_trc.F90
+++ b/src/TOP/TRP/trdmxl_trc.F90
@@ -283,7 +283,7 @@ CONTAINS
          DO jn = 1, jptra
             IF( ln_trdtrc(jn) ) THEN
                DO jl = 1, jpltrd_trc
-                  CALL lbc_lnk( 'trdmxl_trc', tmltrd_trc(:,:,jl,jn), 'T', 1. )        ! lateral boundary conditions
+                  CALL lbc_lnk( 'trdmxl_trc', tmltrd_trc(:,:,jl,jn), 'T', 1._wp )        ! lateral boundary conditions
                END DO
             ENDIF
          END DO
@@ -418,8 +418,8 @@ CONTAINS
                
          !-- Lateral boundary conditions
                IF ( cn_cfg .NE. 'gyre' ) THEN
-                  CALL lbc_lnk( 'trdmxl_trc', ztmltot(:,:,jn) , 'T', 1. , ztmlres(:,:,jn) , 'T', 1., &
-                     &                        ztmlatf(:,:,jn) , 'T', 1. , ztmlrad(:,:,jn) , 'T', 1. )
+                  CALL lbc_lnk( 'trdmxl_trc', ztmltot(:,:,jn) , 'T', 1._wp , ztmlres(:,:,jn) , 'T', 1._wp, &
+                     &                        ztmlatf(:,:,jn) , 'T', 1._wp , ztmlrad(:,:,jn) , 'T', 1._wp )
                ENDIF
 
 
@@ -469,9 +469,9 @@ CONTAINS
 
          !-- Lateral boundary conditions 
                IF ( cn_cfg .NE. 'gyre' ) THEN            ! other than GYRE configuration    
-                  CALL lbc_lnk( 'trdmxl_trc', ztmltot2(:,:,jn), 'T', 1., ztmlres2(:,:,jn), 'T', 1. )
+                  CALL lbc_lnk( 'trdmxl_trc', ztmltot2(:,:,jn), 'T', 1._wp, ztmlres2(:,:,jn), 'T', 1._wp )
                   DO jl = 1, jpltrd_trc
-                     CALL lbc_lnk( 'trdmxl_trc', ztmltrd2(:,:,jl,jn), 'T', 1. )       ! will be output in the NetCDF trends file
+                     CALL lbc_lnk( 'trdmxl_trc', ztmltrd2(:,:,jl,jn), 'T', 1._wp )       ! will be output in the NetCDF trends file
                   END DO
                ENDIF
 
diff --git a/tests/ADIAB_WAVE/MY_SRC/usrdef_zgr.F90 b/tests/ADIAB_WAVE/MY_SRC/usrdef_zgr.F90
index 5c7fac8382fe52b9c044774838d5f5641b28c06a..8e26909d398a7503158c8b713465bf88630a277e 100644
--- a/tests/ADIAB_WAVE/MY_SRC/usrdef_zgr.F90
+++ b/tests/ADIAB_WAVE/MY_SRC/usrdef_zgr.F90
@@ -103,7 +103,7 @@ CONTAINS
       DO_2D(0,0,0,0)
          zhu(ji,jj) = 0.5_wp * ( zht(ji,jj) + zht(ji+1,jj) )
       END_2D
-      CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1. )     ! boundary condition: this mask the surrouding grid-points
+      CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1._wp )     ! boundary condition: this mask the surrouding grid-points
       !                                ! ==>>>  set by hand non-zero value on first/last columns & rows 
       DO ji = mi0(1), mi1(1)              ! first row of global domain only
          zhu(ji,2) = zht(ji,2)
@@ -122,7 +122,7 @@ CONTAINS
       ! no ocean cavities : top ocean level is ONE, except over land
       ! the ocean basin surrounded by land (1+nn_hls grid-points) set through lbc_lnk call
       z2d(:,:) = 1._wp                    ! surface ocean is the 1st level
-      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )        ! closed basin, see userdef_nam.F90
+      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp )        ! closed basin, see userdef_nam.F90
       k_top(:,:) = NINT( z2d(:,:) )
       !
       !                              
diff --git a/tests/BENCH/MY_SRC/usrdef_istate.F90 b/tests/BENCH/MY_SRC/usrdef_istate.F90
index 69da90f1b8425eec665b26e6d8b8c2fc29628abe..a800d13ecefd9b0add39412a78e0e33ff254c5e4 100644
--- a/tests/BENCH/MY_SRC/usrdef_istate.F90
+++ b/tests/BENCH/MY_SRC/usrdef_istate.F90
@@ -82,9 +82,9 @@ CONTAINS
       pu( :,:,jpk  ) = 0._wp
       pv( :,:,jpk  ) = 0._wp
       !
-      CALL lbc_lnk('usrdef_istate',  pts, 'T',  1. )            ! apply boundary conditions
-      CALL lbc_lnk('usrdef_istate',   pu, 'U', -1. )            ! apply boundary conditions
-      CALL lbc_lnk('usrdef_istate',   pv, 'V', -1. )            ! apply boundary conditions
+      CALL lbc_lnk('usrdef_istate',  pts, 'T',  1._wp )            ! apply boundary conditions
+      CALL lbc_lnk('usrdef_istate',   pu, 'U', -1._wp )            ! apply boundary conditions
+      CALL lbc_lnk('usrdef_istate',   pv, 'V', -1._wp )            ! apply boundary conditions
       
    END SUBROUTINE usr_def_istate
 
@@ -111,7 +111,7 @@ CONTAINS
          pssh(ji,jj) = 0.1 * ( 0.5 - REAL( mig0(ji) + (mjg0(jj)-1) * Ni0glo, wp ) / REAL( Ni0glo * Nj0glo, wp ) )
       END_2D
       !
-      CALL lbc_lnk('usrdef_istate', pssh, 'T',  1. )   ! apply boundary conditions
+      CALL lbc_lnk('usrdef_istate', pssh, 'T',  1._wp )   ! apply boundary conditions
       !
    END SUBROUTINE usr_def_istate_ssh
    
diff --git a/tests/BENCH/MY_SRC/usrdef_sbc.F90 b/tests/BENCH/MY_SRC/usrdef_sbc.F90
index e9ac757b669a92d58a27241281a7abf2db635ba8..93877cbd06bef9271a1bd6263d727a1c8aa0e56a 100644
--- a/tests/BENCH/MY_SRC/usrdef_sbc.F90
+++ b/tests/BENCH/MY_SRC/usrdef_sbc.F90
@@ -109,7 +109,7 @@ CONTAINS
          vtau_ice(ji,jj) = 0.1_wp + zztmp
       END_2D
 
-      CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
+      CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'U', -1._wp, vtau_ice, 'V', -1._wp )
 #endif
       !
    END SUBROUTINE usrdef_sbc_ice_tau
diff --git a/tests/CANAL/MY_SRC/usrdef_istate.F90 b/tests/CANAL/MY_SRC/usrdef_istate.F90
index fca0099da7637f80e53b9425d0112546370a13cb..498f8c49962883e5d31ead56f74cbc0542642bf4 100644
--- a/tests/CANAL/MY_SRC/usrdef_istate.F90
+++ b/tests/CANAL/MY_SRC/usrdef_istate.F90
@@ -237,8 +237,8 @@ CONTAINS
          !            
       END SELECT
       !
-      CALL lbc_lnk( 'usrdef_istate', pts , 'T',  1. )
-      CALL lbc_lnk( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. )
+      CALL lbc_lnk( 'usrdef_istate', pts , 'T',  1._wp )
+      CALL lbc_lnk( 'usrdef_istate', pu, 'U', -1._wp, pv, 'V', -1._wp )
 
    END SUBROUTINE usr_def_istate
 
@@ -322,7 +322,7 @@ CONTAINS
             DO_2D( 0, 0, 0, 0 )
                pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * rn_uzonal * EXP( - 0.5 * gphit(ji,jj)**2 / rn_lambda**2 ) * e2t(ji,jj)
             END_2D
-            CALL lbc_lnk( 'usrdef_istate_ssh', pssh, 'T',  1. )
+            CALL lbc_lnk( 'usrdef_istate_ssh', pssh, 'T',  1._wp )
          END DO
          !            
       CASE(4)                      !==  geostrophic zonal pulse !!st need to implement a way to separate ssh properly  ==!
@@ -374,7 +374,7 @@ CONTAINS
          CALL RANDOM_NUMBER(zrandom)
          pssh(:,:) = pssh(:,:) + ( 0.1  * zrandom(:,:) - 0.05 )
       ENDIF
-      CALL lbc_lnk( 'usrdef_istate_ssh', pssh, 'T',  1. )
+      CALL lbc_lnk( 'usrdef_istate_ssh', pssh, 'T',  1._wp )
       !
    END SUBROUTINE usr_def_istate_ssh
    
diff --git a/tests/CANAL/MY_SRC/usrdef_zgr.F90 b/tests/CANAL/MY_SRC/usrdef_zgr.F90
index 1f9da256991333df8dbc50228df7dfddcdd00f62..b4490cf5724c6ffa172b770948680f38ba0f6cf2 100644
--- a/tests/CANAL/MY_SRC/usrdef_zgr.F90
+++ b/tests/CANAL/MY_SRC/usrdef_zgr.F90
@@ -201,7 +201,7 @@ CONTAINS
          z2d(:,:) = REAL(jpkm1 - NINT( 0.5 * REAL(jpkm1,wp) * z2d(:,:) ), wp)
       END SELECT
       !
-      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (closed boundaries)
+      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp )           ! set surrounding land to zero (closed boundaries)
       !
       k_bot(:,:) = NINT( z2d(:,:) )          ! =jpkm1 over the ocean point, =0 elsewhere
       !
diff --git a/tests/DOME/MY_SRC/usrdef_zgr.F90 b/tests/DOME/MY_SRC/usrdef_zgr.F90
index 797006a772fd6df5b35684cf7bc4f997fe1dc69c..cf112efa4bef75735fbbb93e1fe30ebafc352c96 100644
--- a/tests/DOME/MY_SRC/usrdef_zgr.F90
+++ b/tests/DOME/MY_SRC/usrdef_zgr.F90
@@ -104,7 +104,7 @@ CONTAINS
       WHERE (gphit(:,:)>0._wp) z2d(:,:) = 0._wp
       ! Dig inlet:
       WHERE ((gphit(:,:)>0._wp).AND.(glamt(:,:)>-50._wp).AND.(glamt(:,:)<50._wp)) z2d(:,:) = 1._wp
-      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )        ! closed basin, see userdef_nam.F90
+      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp )        ! closed basin, see userdef_nam.F90
       k_top(:,:) = NINT( z2d(:,:) )
       !
       !                              
diff --git a/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90 b/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
index 8d5fe4eb615de8df2477cec3d0c65b87c3c64b02..4fe72a986d89687e271fff3179732e030cc4bec0 100644
--- a/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
+++ b/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
@@ -125,7 +125,7 @@ CONTAINS
          windu(ji,jj) = Umax/sqrt(d*1000)*(d-2*mig(ji)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*min(kt*30./21600,1.)
          windv(ji,jj) = Umax/sqrt(d*1000)*(d-2*mjg(jj)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*Rwind*min(kt*30./21600,1.)
       END_2D
-      CALL lbc_lnk( 'usrdef_sbc', windu, 'U', -1., windv, 'V', -1. )
+      CALL lbc_lnk( 'usrdef_sbc', windu, 'U', -1._wp, windv, 'V', -1._wp )
 
       wndm_ice(:,:) = 0._wp      !!gm brutal....
 
@@ -138,7 +138,7 @@ CONTAINS
          zwndj_t = (  windv(ji,jj) - r_vfac * 0.5 * ( v_ice(ji,jj-1) + v_ice(ji,jj) )  )
          wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
       END_2D
-      CALL lbc_lnk( 'usrdef_sbc', wndm_ice, 'T',  1. )
+      CALL lbc_lnk( 'usrdef_sbc', wndm_ice, 'T',  1._wp )
 
       !!gm brutal....
       utau_ice  (:,:) = 0._wp
@@ -155,7 +155,7 @@ CONTAINS
          vtau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
             &          * ( 0.5 * (windv(ji,jj+1) + windv(ji,jj) ) - r_vfac * v_ice(ji,jj) )
       END_2D
-      CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
+      CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'U', -1._wp, vtau_ice, 'V', -1._wp )
       !
    END SUBROUTINE usrdef_sbc_ice_tau
 
diff --git a/tests/LOCK_EXCHANGE/MY_SRC/usrdef_zgr.F90 b/tests/LOCK_EXCHANGE/MY_SRC/usrdef_zgr.F90
index a6b980b7f512a4de82052b9700c96964f0b138dc..5daa0f23f4c1153cccb8db68e7b8ba530ff0271c 100644
--- a/tests/LOCK_EXCHANGE/MY_SRC/usrdef_zgr.F90
+++ b/tests/LOCK_EXCHANGE/MY_SRC/usrdef_zgr.F90
@@ -85,7 +85,7 @@ CONTAINS
       ! no ocean cavities : top ocean level is ONE, except over land
       ! the ocean basin surrounded by land (1+nn_hls grid-points) set through lbc_lnk call
       z2d(:,:) = 1._wp                    ! surface ocean is the 1st level
-      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )        ! closed basin, see userdef_nam.F90
+      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp )        ! closed basin, see userdef_nam.F90
       k_top(:,:) = NINT( z2d(:,:) )
       !
       !                              
diff --git a/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90 b/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90
index b39e22ab155ee6eaf7ba10d463f3edfc0205a74d..e0c4b1d0aaf9e02d57ccf1841868999043f28fa2 100644
--- a/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90
+++ b/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90
@@ -92,7 +92,7 @@ CONTAINS
       DO_2D( 0, 0, 0, 0 )
          zhu(ji,jj) = 0.5_wp * ( zht(ji,jj) + zht(ji+1,jj) )
       END_2D
-      CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1. )     ! boundary condition: this mask the surrouding grid-points
+      CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1._wp )     ! boundary condition: this mask the surrouding grid-points
       !                                ! ==>>>  set by hand non-zero value on first/last columns & rows 
       DO ji = mi0(1), mi1(1)              ! first row of global domain only
          zhu(ji,2) = zht(ji,2)
@@ -111,7 +111,7 @@ CONTAINS
       ! no ocean cavities : top ocean level is ONE, except over land
       ! the ocean basin surrounded by land (1+nn_hls grid-points) set through lbc_lnk call
       z2d(:,:) = 1._wp                    ! surface ocean is the 1st level
-      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )        ! closed basin, see userdef_nam.F90
+      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp )        ! closed basin, see userdef_nam.F90
       k_top(:,:) = NINT( z2d(:,:) )
       !
       !                              
diff --git a/tests/SWG/MY_SRC/usrdef_sbc.F90 b/tests/SWG/MY_SRC/usrdef_sbc.F90
index a3069954721d86aa464b7bb297e5d233f808c8b3..9ee4a9aaa6192a536774b15160cf2253170a9bda 100644
--- a/tests/SWG/MY_SRC/usrdef_sbc.F90
+++ b/tests/SWG/MY_SRC/usrdef_sbc.F90
@@ -101,7 +101,7 @@ CONTAINS
          taum(ji,jj) = zmod
          wndm(ji,jj) = SQRT( zmod * zcoef )
       END_2D
-      CALL lbc_lnk( 'usrdef_sbc', taum(:,:), 'T', 1. , wndm(:,:), 'T', 1. )
+      CALL lbc_lnk( 'usrdef_sbc', taum(:,:), 'T', 1._wp , wndm(:,:), 'T', 1._wp )
       !
    END SUBROUTINE usrdef_sbc_oce
 
diff --git a/tests/SWG/MY_SRC/usrdef_zgr.F90 b/tests/SWG/MY_SRC/usrdef_zgr.F90
index acd295e3dd6197b329b8749b6cc37fe6b1d78a17..a96c15cdf6ee1262b6b1e31ccfb9c1fb0cbf3771 100644
--- a/tests/SWG/MY_SRC/usrdef_zgr.F90
+++ b/tests/SWG/MY_SRC/usrdef_zgr.F90
@@ -191,7 +191,7 @@ CONTAINS
       !
       z2d(:,:) = REAL( jpkm1 , wp )          ! flat bottom
       !
-      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (closed boundaries)
+      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp )           ! set surrounding land to zero (closed boundaries)
       !
       zylim0 =   10000._wp    ! +10km 
       zylim1 = 2010000._wp    ! 2010km
diff --git a/tests/TSUNAMI/MY_SRC/usrdef_istate.F90 b/tests/TSUNAMI/MY_SRC/usrdef_istate.F90
index 0ca6c3492fd66504612e2c8f224d4c330d8eee84..69202f700203415cb384af1d7f562d3486fc8c27 100644
--- a/tests/TSUNAMI/MY_SRC/usrdef_istate.F90
+++ b/tests/TSUNAMI/MY_SRC/usrdef_istate.F90
@@ -103,7 +103,7 @@ CONTAINS
       END_2D
 
       !
-      CALL lbc_lnk('usrdef_istate', pssh, 'T',  1. )            ! apply boundary conditions
+      CALL lbc_lnk('usrdef_istate', pssh, 'T',  1._wp )            ! apply boundary conditions
       !
    END SUBROUTINE usr_def_istate_ssh
    
diff --git a/tests/VORTEX/MY_SRC/usrdef_istate.F90 b/tests/VORTEX/MY_SRC/usrdef_istate.F90
index 6de0b19191cffc5940d9cc6763bd90c93f90ec48..7694c43c9c5ec788f0a943a2615b431b2e0a1b4c 100644
--- a/tests/VORTEX/MY_SRC/usrdef_istate.F90
+++ b/tests/VORTEX/MY_SRC/usrdef_istate.F90
@@ -122,7 +122,7 @@ CONTAINS
          END DO
       END_2D
       !
-      CALL lbc_lnk( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. )
+      CALL lbc_lnk( 'usrdef_istate', pu, 'U', -1._wp, pv, 'V', -1._wp )
       !   
    END SUBROUTINE usr_def_istate
 
diff --git a/tests/VORTEX/MY_SRC/usrdef_zgr.F90 b/tests/VORTEX/MY_SRC/usrdef_zgr.F90
index 6b02fbcb9e6c72a41b1f35cd97a1f79e9a18c94f..575bc330527ddd4cb4579696557eece89d7e2bc2 100644
--- a/tests/VORTEX/MY_SRC/usrdef_zgr.F90
+++ b/tests/VORTEX/MY_SRC/usrdef_zgr.F90
@@ -189,7 +189,7 @@ CONTAINS
       !
       z2d(:,:) = REAL( jpkm1 , wp )          ! flat bottom
       !
-      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (closed boundaries)
+      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp )           ! set surrounding land to zero (closed boundaries)
       !
       k_bot(:,:) = NINT( z2d(:,:) )          ! =jpkm1 over the ocean point, =0 elsewhere
       !
diff --git a/tests/WAD/MY_SRC/usrdef_zgr.F90 b/tests/WAD/MY_SRC/usrdef_zgr.F90
index eab2b9a325da1789b58db614c2e3bc6bc181996e..2543e30e5ce0e8e838397b6af4246897d1ae7cfe 100644
--- a/tests/WAD/MY_SRC/usrdef_zgr.F90
+++ b/tests/WAD/MY_SRC/usrdef_zgr.F90
@@ -232,7 +232,7 @@ CONTAINS
       DO_2D( 0, 0, 0, 0 )
          zhu(ji,jj) = 0.5_wp * ( zht(ji,jj) + zht(ji+1,jj) )
       END_2D
-      CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1. )     ! boundary condition: this mask the surrounding grid-points
+      CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1._wp )     ! boundary condition: this mask the surrounding grid-points
       !                                ! ==>>>  set by hand non-zero value on first/last columns & rows 
       DO ji = mi0(1), mi1(1)              ! first row of global domain only
          zhu(ji,:) = zht(1,:)
@@ -245,7 +245,7 @@ CONTAINS
       DO_2D( 0, 0, 0, 0 )
          zhv(ji,jj) = 0.5_wp * ( zht(ji,jj) + zht(ji,jj+1) )
       END_2D
-      CALL lbc_lnk( 'usrdef_zgr', zhv, 'V', 1. )     ! boundary condition: this mask the surrounding grid-points
+      CALL lbc_lnk( 'usrdef_zgr', zhv, 'V', 1._wp )     ! boundary condition: this mask the surrounding grid-points
       DO jj = mj0(1), mj1(1)   ! first  row of global domain only
          zhv(:,jj) = zht(:,jj)
       END DO
@@ -266,7 +266,7 @@ CONTAINS
       z2d(:,mj0(1):mj1(1)) = 0._wp
       z2d(:,mj0(jpjglo):mj1(jpjglo)) = 0._wp
 
-      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )        ! closed basin, see userdef_nam.F90
+      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1._wp )        ! closed basin, see userdef_nam.F90
       k_top(:,:) = NINT( z2d(:,:) )
       !
       !                              
@@ -301,15 +301,15 @@ CONTAINS
                pe3vw(ji,jj,jk) = zwet * z1_jpkm1
             END DO
          END_2D     
-         CALL lbc_lnk( 'usrdef_zgr', pdept, 'T', 1. )
-         CALL lbc_lnk( 'usrdef_zgr', pdepw, 'T', 1. )
-         CALL lbc_lnk( 'usrdef_zgr', pe3t , 'T', 1. )
-         CALL lbc_lnk( 'usrdef_zgr', pe3w , 'T', 1. )
-         CALL lbc_lnk( 'usrdef_zgr', pe3u , 'U', 1. )
-         CALL lbc_lnk( 'usrdef_zgr', pe3uw, 'U', 1. )
-         CALL lbc_lnk( 'usrdef_zgr', pe3f , 'F', 1. )
-         CALL lbc_lnk( 'usrdef_zgr', pe3v , 'V', 1. )
-         CALL lbc_lnk( 'usrdef_zgr', pe3vw, 'V', 1. )
+         CALL lbc_lnk( 'usrdef_zgr', pdept, 'T', 1._wp )
+         CALL lbc_lnk( 'usrdef_zgr', pdepw, 'T', 1._wp )
+         CALL lbc_lnk( 'usrdef_zgr', pe3t , 'T', 1._wp )
+         CALL lbc_lnk( 'usrdef_zgr', pe3w , 'T', 1._wp )
+         CALL lbc_lnk( 'usrdef_zgr', pe3u , 'U', 1._wp )
+         CALL lbc_lnk( 'usrdef_zgr', pe3uw, 'U', 1._wp )
+         CALL lbc_lnk( 'usrdef_zgr', pe3f , 'F', 1._wp )
+         CALL lbc_lnk( 'usrdef_zgr', pe3v , 'V', 1._wp )
+         CALL lbc_lnk( 'usrdef_zgr', pe3vw, 'V', 1._wp )
          WHERE( pe3t (:,:,:) == 0._wp )   pe3t (:,:,:) = 1._wp
          WHERE( pe3u (:,:,:) == 0._wp )   pe3u (:,:,:) = 1._wp
          WHERE( pe3v (:,:,:) == 0._wp )   pe3v (:,:,:) = 1._wp
diff --git a/tools/MISCELLANEOUS/add_wp_in_call_lbc_lnk.sh b/tools/MISCELLANEOUS/add_wp_in_call_lbc_lnk.sh
new file mode 100755
index 0000000000000000000000000000000000000000..5739b9283aa5b0ef6db462c6358f4c3f0ddc0ff5
--- /dev/null
+++ b/tools/MISCELLANEOUS/add_wp_in_call_lbc_lnk.sh
@@ -0,0 +1,15 @@
+#!/bin/bash
+#
+# this sed makes sure that we alwas specify the _wp in the calls to lbc_lnk
+#
+for ff in $( find src -type f ) */*/MY_SRC/*  # for all nemo files
+do
+    # look for something like ", 'X', 1. )" and replace it by ", 'X', 1.wp )"
+    #   - X can be TUVFW
+    #   - the "1" can be followed by a "." or not. It must be followed by "," or ")" or "&"
+    #   - any number of blanks is acceptable.
+    sed -e "s/\(, *'[TUVFW]' *, *-*1\)\.*\( *[,)&]\)/\1._wp\2/g" $ff > tmp
+    mv tmp $ff
+done
+
+