diff --git a/src/OCE/DIA/diahth.F90 b/src/OCE/DIA/diahth.F90
index 26c97b7dff6b7fffa247c568d1fc887c811d6741..cdb657b234083caeb59fa2a360ea155fcee00902 100644
--- a/src/OCE/DIA/diahth.F90
+++ b/src/OCE/DIA/diahth.F90
@@ -326,8 +326,8 @@ CONTAINS
          iid = iktem(ji,jj)
          IF( iid /= 1 ) THEN 
              zztmp =     gdept(ji,jj,iid  ,Kmm)   &                     ! linear interpolation
-               &  + (    gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm)                       )   &
-               &  * ( 20.*tmask(ji,jj,iid+1) - ts(ji,jj,iid,jp_tem,Kmm)                       )   &
+               &  + (    gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm)                            )   &
+               &  * ( ptem * tmask(ji,jj,iid+1)  - ts(ji,jj,iid,jp_tem,Kmm)                       )   &
                &  / ( ts(ji,jj,iid+1,jp_tem,Kmm) - ts(ji,jj,iid,jp_tem,Kmm) + (1.-tmask(ji,jj,1)) )
             pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
          ELSE