diff --git a/src/OCE/DIA/diahth.F90 b/src/OCE/DIA/diahth.F90
index 2cee430d390d650c650df5e2745bb32a5cfa2ce3..0b8650acf6c31eacf0ef56eb870ebf898a7aee5f 100644
--- a/src/OCE/DIA/diahth.F90
+++ b/src/OCE/DIA/diahth.F90
@@ -106,12 +106,13 @@ CONTAINS
IF( ln_timing ) CALL timing_start('dia_hth')
IF( kt == nit000 ) THEN
- l_hth = .FALSE.
- IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) .OR. &
- & iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' ) .OR. &
- & iom_use( '20d' ) .OR. iom_use( '26d' ) .OR. iom_use( '28d' ) .OR. &
- & iom_use( 'hc300' ) .OR. iom_use( 'hc700' ) .OR. iom_use( 'hc2000' ) .OR. &
- & iom_use( 'pycndep' ) .OR. iom_use( 'tinv' ) .OR. iom_use( 'depti' ) ) l_hth = .TRUE.
+ !
+ l_hth = iom_use( 'mlddzt' ) .OR. iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) .OR. &
+ & iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' ) .OR. &
+ & iom_use( '20d' ) .OR. iom_use( '26d' ) .OR. iom_use( '28d' ) .OR. &
+ & iom_use( 'hc300' ) .OR. iom_use( 'hc700' ) .OR. iom_use( 'hc2000' ) .OR. &
+ & iom_use( 'pycndep' ) .OR. iom_use( 'tinv' ) .OR. iom_use( 'depti' )
+ !
! ! allocate dia_hth array
IF( l_hth ) THEN
IF( dia_hth_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' )
@@ -124,11 +125,12 @@ CONTAINS
IF( l_hth ) THEN
!
- IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN
- ! initialization
- ztinv (:,:) = 0._wp
- zdepinv(:,:) = 0._wp
- zmaxdzT(:,:) = 0._wp
+ ! initialization
+ IF( iom_use( 'tinv' ) ) ztinv (:,:) = 0._wp
+ IF( iom_use( 'depti' ) ) zdepinv(:,:) = 0._wp
+ IF( iom_use( 'mlddzt' ) ) zmaxdzT(:,:) = 0._wp
+ IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) &
+ & .OR. iom_use( 'mldr10_3' ) .OR. iom_use( 'pycndep' ) ) THEN
DO_2D( 1, 1, 1, 1 )
zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)
hth (ji,jj) = zztmp
@@ -137,6 +139,8 @@ CONTAINS
zrho10_3(ji,jj) = zztmp
zpycn (ji,jj) = zztmp
END_2D
+ ENDIF
+ IF( iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN
IF( nla10 > 1 ) THEN
DO_2D( 1, 1, 1, 1 )
zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)
@@ -144,25 +148,9 @@ CONTAINS
zrho0_1(ji,jj) = zztmp
END_2D
ENDIF
-
- ! Preliminary computation
- ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
- DO_2D( 1, 1, 1, 1 )
- IF( tmask(ji,jj,nla10) == 1. ) THEN
- zu = 1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80 * ts(ji,jj,nla10,jp_sal,Kmm) &
- & - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) &
- & - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm)
- zv = 5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00 * ts(ji,jj,nla10,jp_sal,Kmm) &
- & - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)
- zut = 11.25 - 0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01 * ts(ji,jj,nla10,jp_sal,Kmm)
- zvt = 38.00 - 0.750 * ts(ji,jj,nla10,jp_tem,Kmm)
- zw = (zu + 0.698*zv) * (zu + 0.698*zv)
- zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
- ELSE
- zdelr(ji,jj) = 0._wp
- ENDIF
- END_2D
-
+ ENDIF
+
+ IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN
! ------------------------------------------------------------- !
! thermocline depth: strongest vertical gradient of temperature !
! turbocline depth (mixing layer depth): avt = zavt5 !
@@ -198,6 +186,25 @@ CONTAINS
!
IF( iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' ) .OR. &
& iom_use( 'pycndep' ) .OR. iom_use( 'tinv' ) .OR. iom_use( 'depti' ) ) THEN
+ !
+ ! Preliminary computation
+ ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
+ DO_2D( 1, 1, 1, 1 )
+ IF( tmask(ji,jj,nla10) == 1. ) THEN
+ zu = 1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80 * ts(ji,jj,nla10,jp_sal,Kmm) &
+ & - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) &
+ & - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm)
+ zv = 5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00 * ts(ji,jj,nla10,jp_sal,Kmm) &
+ & - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)
+ zut = 11.25 - 0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01 * ts(ji,jj,nla10,jp_sal,Kmm)
+ zvt = 38.00 - 0.750 * ts(ji,jj,nla10,jp_tem,Kmm)
+ zw = (zu + 0.698*zv) * (zu + 0.698*zv)
+ zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
+ ELSE
+ zdelr(ji,jj) = 0._wp
+ ENDIF
+ END_2D
+ !
! ------------------------------------------------------------- !
! MLD: abs( tn - tn(10m) ) = ztem2 !
! Top of thermocline: tn = tn(10m) - ztem2 !
diff --git a/src/OCE/DOM/dtatsd.F90 b/src/OCE/DOM/dtatsd.F90
index 7ffccc4cfbc85e95a31aa1518138450e0fc5329a..5863789cc4d8e4c1f09657a00610197ef9c7c599 100644
--- a/src/OCE/DOM/dtatsd.F90
+++ b/src/OCE/DOM/dtatsd.F90
@@ -199,8 +199,7 @@ CONTAINS
ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk)
END_3D
!
-! JC I think it's more convenient to consider the general sco case as the rule
-! IF( ln_sco ) THEN !== s- or mixed s-zps-coordinate ==!
+ IF( ln_sco ) THEN !== s- or mixed s-zps-coordinate ==!
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 .AND. lwp )THEN
@@ -236,31 +235,35 @@ CONTAINS
ptsd(ji,jj,jpk,jp_sal) = 0._wp
END_2D
!
-! ELSE !== z- or zps- coordinate ==!
-! !
-! DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
-! ptsd(ji,jj,jk,jp_tem) = ptsd(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) ! Mask
-! ptsd(ji,jj,jk,jp_sal) = ptsd(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
-! END_3D
-! !
-! IF( ln_zps ) THEN ! zps-coordinate (partial steps) interpolation at the last ocean level
-! DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-! ik = mbkt(ji,jj)
-! IF( ik > 1 ) THEN
-! zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
-! ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
-! ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
-! ENDIF
-! ik = mikt(ji,jj)
-! IF( ik > 1 ) THEN
-! zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
-! ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
-! ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
-! END IF
-! END_2D
-! ENDIF
-! !
-! ENDIF
+ ELSE !== z- or zps- coordinate ==!
+ !
+ ! We must keep this definition in a case different from the general case of s-coordinate as we don't
+ ! want to use "underground" values (levels below ocean bottom) to be able to start the model from
+ ! masked temp and sal (read for example in a restart or in output.init)
+ !
+ DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+ ptsd(ji,jj,jk,jp_tem) = ptsd(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) ! Mask
+ ptsd(ji,jj,jk,jp_sal) = ptsd(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
+ END_3D
+ !
+ IF( ln_zps ) THEN ! zps-coordinate (partial steps) interpolation at the last ocean level
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ ik = mbkt(ji,jj)
+ IF( ik > 1 ) THEN
+ zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+ ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
+ ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
+ ENDIF
+ ik = mikt(ji,jj)
+ IF( ik > 1 ) THEN
+ zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
+ ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
+ ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
+ END IF
+ END_2D
+ ENDIF
+ !
+ ENDIF
!
IF( .NOT.ln_tsd_dmp ) THEN !== deallocate T & S structure ==!
! (data used only for initialisation)
diff --git a/src/OCE/SBC/sbcssm.F90 b/src/OCE/SBC/sbcssm.F90
index 1c108df091e2a0b9d59535414d041ad4ea5d0669..95e64875363eb85952e4b5b8c5a1dbd72963b9e8 100644
--- a/src/OCE/SBC/sbcssm.F90
+++ b/src/OCE/SBC/sbcssm.F90
@@ -222,7 +222,7 @@ CONTAINS
!
IF( zf_sbc /= REAL( nn_fsbc, wp ) ) THEN ! nn_fsbc has changed between 2 runs
IF(lwp) WRITE(numout,*) ' restart with a change in the frequency of mean from ', zf_sbc, ' to ', nn_fsbc
- zcoef = REAL( nn_fsbc - 1, wp ) / zf_sbc
+ zcoef = REAL( nn_fsbc - 1, wp ) / ( zf_sbc - 1._wp ) ! zf_sbc /= 1 as it was written in the restart
ssu_m(:,:) = zcoef * ssu_m(:,:)
ssv_m(:,:) = zcoef * ssv_m(:,:)
sst_m(:,:) = zcoef * sst_m(:,:)