Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_interp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_interp.F90	(revision 10988)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_interp.F90	(revision 10989)
@@ -107,13 +107,13 @@
          !
          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            ua_b(ibdy1:ibdy2,:) = 0._wp
+            uu_b(ibdy1:ibdy2,:,Krhs) = 0._wp
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) & 
-                      & + e3u_a(ibdy1:ibdy2,jj,jk) * ua(ibdy1:ibdy2,jj,jk) * umask(ibdy1:ibdy2,jj,jk)
+                  uu_b(ibdy1:ibdy2,jj,Krhs) = uu_b(ibdy1:ibdy2,jj,Krhs) & 
+                      & + e3u(ibdy1:ibdy2,jj,jk,Krhs) * uu(ibdy1:ibdy2,jj,jk,Krhs) * umask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
             DO jj = 1, jpj
-               ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
+               uu_b(ibdy1:ibdy2,jj,Krhs) = uu_b(ibdy1:ibdy2,jj,Krhs) * r1_hu_a(ibdy1:ibdy2,jj)
             END DO
          ENDIF
@@ -122,5 +122,5 @@
             DO jk=1,jpkm1              ! Smooth
                DO jj=j1,j2
-                  ua(ibdy2,jj,jk) = 0.25_wp*(ua(ibdy2-1,jj,jk)+2._wp*ua(ibdy2,jj,jk)+ua(ibdy2+1,jj,jk))
+                  uu(ibdy2,jj,jk,Krhs) = 0.25_wp*(uu(ibdy2-1,jj,jk,Krhs)+2._wp*uu(ibdy2,jj,jk,Krhs)+uu(ibdy2+1,jj,jk,Krhs))
                END DO
             END DO
@@ -131,5 +131,5 @@
             DO jj = 1, jpj
                zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) & 
-                  & + e3u_a(ibdy1:ibdy2,jj,jk)  * ua(ibdy1:ibdy2,jj,jk)*umask(ibdy1:ibdy2,jj,jk)
+                  & + e3u(ibdy1:ibdy2,jj,jk,Krhs)  * uu(ibdy1:ibdy2,jj,jk,Krhs)*umask(ibdy1:ibdy2,jj,jk)
             END DO
          END DO
@@ -140,6 +140,6 @@
          DO jk = 1, jpkm1
             DO jj = 1, jpj
-               ua(ibdy1:ibdy2,jj,jk) = ( ua(ibdy1:ibdy2,jj,jk) &
-                 & + ua_b(ibdy1:ibdy2,jj)-zub(ibdy1:ibdy2,jj)) * umask(ibdy1:ibdy2,jj,jk)
+               uu(ibdy1:ibdy2,jj,jk,Krhs) = ( uu(ibdy1:ibdy2,jj,jk,Krhs) &
+                 & + uu_b(ibdy1:ibdy2,jj,Krhs)-zub(ibdy1:ibdy2,jj)) * umask(ibdy1:ibdy2,jj,jk)
             END DO
          END DO
@@ -150,5 +150,5 @@
                DO jj = 1, jpj
                   zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) & 
-                     & + e3v_a(ibdy1:ibdy2,jj,jk) * va(ibdy1:ibdy2,jj,jk) * vmask(ibdy1:ibdy2,jj,jk)
+                     & + e3v(ibdy1:ibdy2,jj,jk,Krhs) * vv(ibdy1:ibdy2,jj,jk,Krhs) * vmask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -158,6 +158,6 @@
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  va(ibdy1:ibdy2,jj,jk) = ( va(ibdy1:ibdy2,jj,jk) & 
-                    & + va_b(ibdy1:ibdy2,jj)-zvb(ibdy1:ibdy2,jj))*vmask(ibdy1:ibdy2,jj,jk)
+                  vv(ibdy1:ibdy2,jj,jk,Krhs) = ( vv(ibdy1:ibdy2,jj,jk,Krhs) & 
+                    & + vv_b(ibdy1:ibdy2,jj,Krhs)-zvb(ibdy1:ibdy2,jj))*vmask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -166,6 +166,6 @@
          DO jk = 1, jpkm1              ! Mask domain edges
             DO jj = 1, jpj
-               ua(1,jj,jk) = 0._wp
-               va(1,jj,jk) = 0._wp
+               uu(1,jj,jk,Krhs) = 0._wp
+               vv(1,jj,jk,Krhs) = 0._wp
             END DO
          END DO 
@@ -178,13 +178,13 @@
          !
          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            ua_b(ibdy1:ibdy2,:) = 0._wp
+            uu_b(ibdy1:ibdy2,:,Krhs) = 0._wp
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) & 
-                      & + e3u_a(ibdy1:ibdy2,jj,jk) * ua(ibdy1:ibdy2,jj,jk) * umask(ibdy1:ibdy2,jj,jk)
+                  uu_b(ibdy1:ibdy2,jj,Krhs) = uu_b(ibdy1:ibdy2,jj,Krhs) & 
+                      & + e3u(ibdy1:ibdy2,jj,jk,Krhs) * uu(ibdy1:ibdy2,jj,jk,Krhs) * umask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
             DO jj = 1, jpj
-               ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
+               uu_b(ibdy1:ibdy2,jj,Krhs) = uu_b(ibdy1:ibdy2,jj,Krhs) * r1_hu_a(ibdy1:ibdy2,jj)
             END DO
          ENDIF
@@ -193,5 +193,5 @@
             DO jk=1,jpkm1              ! Smooth
                DO jj=j1,j2
-                  ua(ibdy1,jj,jk) = 0.25_wp*(ua(ibdy1-1,jj,jk)+2._wp*ua(ibdy1,jj,jk)+ua(ibdy1+1,jj,jk))
+                  uu(ibdy1,jj,jk,Krhs) = 0.25_wp*(uu(ibdy1-1,jj,jk,Krhs)+2._wp*uu(ibdy1,jj,jk,Krhs)+uu(ibdy1+1,jj,jk,Krhs))
                END DO
             END DO
@@ -202,5 +202,5 @@
             DO jj = 1, jpj
                zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) & 
-                  & + e3u_a(ibdy1:ibdy2,jj,jk)  * ua(ibdy1:ibdy2,jj,jk) * umask(ibdy1:ibdy2,jj,jk)
+                  & + e3u(ibdy1:ibdy2,jj,jk,Krhs)  * uu(ibdy1:ibdy2,jj,jk,Krhs) * umask(ibdy1:ibdy2,jj,jk)
             END DO
          END DO
@@ -211,6 +211,6 @@
          DO jk = 1, jpkm1
             DO jj = 1, jpj
-               ua(ibdy1:ibdy2,jj,jk) = ( ua(ibdy1:ibdy2,jj,jk) & 
-                 & + ua_b(ibdy1:ibdy2,jj)-zub(ibdy1:ibdy2,jj))*umask(ibdy1:ibdy2,jj,jk)
+               uu(ibdy1:ibdy2,jj,jk,Krhs) = ( uu(ibdy1:ibdy2,jj,jk,Krhs) & 
+                 & + uu_b(ibdy1:ibdy2,jj,Krhs)-zub(ibdy1:ibdy2,jj))*umask(ibdy1:ibdy2,jj,jk)
             END DO
          END DO
@@ -223,5 +223,5 @@
                DO jj = 1, jpj
                   zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) &
-                     & + e3v_a(ibdy1:ibdy2,jj,jk) * va(ibdy1:ibdy2,jj,jk) * vmask(ibdy1:ibdy2,jj,jk)
+                     & + e3v(ibdy1:ibdy2,jj,jk,Krhs) * vv(ibdy1:ibdy2,jj,jk,Krhs) * vmask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -231,6 +231,6 @@
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  va(ibdy1:ibdy2,jj,jk) = ( va(ibdy1:ibdy2,jj,jk) & 
-                      & + va_b(ibdy1:ibdy2,jj)-zvb(ibdy1:ibdy2,jj)) * vmask(ibdy1:ibdy2,jj,jk)
+                  vv(ibdy1:ibdy2,jj,jk,Krhs) = ( vv(ibdy1:ibdy2,jj,jk,Krhs) & 
+                      & + vv_b(ibdy1:ibdy2,jj,Krhs)-zvb(ibdy1:ibdy2,jj)) * vmask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -239,6 +239,6 @@
          DO jk = 1, jpkm1              ! Mask domain edges
             DO jj = 1, jpj
-               ua(nlci-1,jj,jk) = 0._wp
-               va(nlci  ,jj,jk) = 0._wp
+               uu(nlci-1,jj,jk,Krhs) = 0._wp
+               vv(nlci  ,jj,jk,Krhs) = 0._wp
             END DO
          END DO 
@@ -251,13 +251,13 @@
          !
          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            va_b(:,jbdy1:jbdy2) = 0._wp
+            vv_b(:,jbdy1:jbdy2,Krhs) = 0._wp
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) & 
-                      & + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)
+                  vv_b(ji,jbdy1:jbdy2,Krhs) = vv_b(ji,jbdy1:jbdy2,Krhs) & 
+                      & + e3v(ji,jbdy1:jbdy2,jk,Krhs) * vv(ji,jbdy1:jbdy2,jk,Krhs) * vmask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
             DO ji=1,jpi
-               va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
+               vv_b(ji,jbdy1:jbdy2,Krhs) = vv_b(ji,jbdy1:jbdy2,Krhs) * r1_hv_a(ji,jbdy1:jbdy2)
             END DO
          ENDIF
@@ -266,5 +266,5 @@
             DO jk = 1, jpkm1           ! Smooth
                DO ji = i1, i2
-                  va(ji,jbdy2,jk) = 0.25_wp*(va(ji,jbdy2-1,jk)+2._wp*va(ji,jbdy2,jk)+va(ji,jbdy2+1,jk))
+                  vv(ji,jbdy2,jk,Krhs) = 0.25_wp*(vv(ji,jbdy2-1,jk,Krhs)+2._wp*vv(ji,jbdy2,jk,Krhs)+vv(ji,jbdy2+1,jk,Krhs))
                END DO
             END DO
@@ -275,5 +275,5 @@
             DO ji=1,jpi
                zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) & 
-                  & + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)
+                  & + e3v(ji,jbdy1:jbdy2,jk,Krhs) * vv(ji,jbdy1:jbdy2,jk,Krhs) * vmask(ji,jbdy1:jbdy2,jk)
             END DO
          END DO
@@ -284,6 +284,6 @@
          DO jk = 1, jpkm1
             DO ji = 1, jpi
-               va(ji,jbdy1:jbdy2,jk) = ( va(ji,jbdy1:jbdy2,jk) & 
-                 & + va_b(ji,jbdy1:jbdy2) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
+               vv(ji,jbdy1:jbdy2,jk,Krhs) = ( vv(ji,jbdy1:jbdy2,jk,Krhs) & 
+                 & + vv_b(ji,jbdy1:jbdy2,Krhs) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
             END DO
          END DO
@@ -294,5 +294,5 @@
                DO ji = 1, jpi
                   zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) & 
-                     & + e3u_a(ji,jbdy1:jbdy2,jk) * ua(ji,jbdy1:jbdy2,jk) * umask(ji,jbdy1:jbdy2,jk)
+                     & + e3u(ji,jbdy1:jbdy2,jk,Krhs) * uu(ji,jbdy1:jbdy2,jk,Krhs) * umask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -303,6 +303,6 @@
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  ua(ji,jbdy1:jbdy2,jk) = ( ua(ji,jbdy1:jbdy2,jk) & 
-                    & + ua_b(ji,jbdy1:jbdy2) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
+                  uu(ji,jbdy1:jbdy2,jk,Krhs) = ( uu(ji,jbdy1:jbdy2,jk,Krhs) & 
+                    & + uu_b(ji,jbdy1:jbdy2,Krhs) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -311,6 +311,6 @@
          DO jk = 1, jpkm1              ! Mask domain edges
             DO ji = 1, jpi
-               ua(ji,1,jk) = 0._wp
-               va(ji,1,jk) = 0._wp
+               uu(ji,1,jk,Krhs) = 0._wp
+               vv(ji,1,jk,Krhs) = 0._wp
             END DO
          END DO 
@@ -323,13 +323,13 @@
          !
          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            va_b(:,jbdy1:jbdy2) = 0._wp
+            vv_b(:,jbdy1:jbdy2,Krhs) = 0._wp
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) & 
-                      & + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)
+                  vv_b(ji,jbdy1:jbdy2,Krhs) = vv_b(ji,jbdy1:jbdy2,Krhs) & 
+                      & + e3v(ji,jbdy1:jbdy2,jk,Krhs) * vv(ji,jbdy1:jbdy2,jk,Krhs) * vmask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
             DO ji=1,jpi
-               va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
+               vv_b(ji,jbdy1:jbdy2,Krhs) = vv_b(ji,jbdy1:jbdy2,Krhs) * r1_hv_a(ji,jbdy1:jbdy2)
             END DO
          ENDIF
@@ -338,5 +338,5 @@
             DO jk = 1, jpkm1           ! Smooth
                DO ji = i1, i2
-                  va(ji,jbdy1,jk) = 0.25_wp*(va(ji,jbdy1-1,jk)+2._wp*va(ji,jbdy1,jk)+va(ji,jbdy1+1,jk))
+                  vv(ji,jbdy1,jk,Krhs) = 0.25_wp*(vv(ji,jbdy1-1,jk,Krhs)+2._wp*vv(ji,jbdy1,jk,Krhs)+vv(ji,jbdy1+1,jk,Krhs))
                END DO
             END DO
@@ -347,5 +347,5 @@
             DO ji=1,jpi
                zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) & 
-                  & + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)
+                  & + e3v(ji,jbdy1:jbdy2,jk,Krhs) * vv(ji,jbdy1:jbdy2,jk,Krhs) * vmask(ji,jbdy1:jbdy2,jk)
             END DO
          END DO
@@ -356,6 +356,6 @@
          DO jk = 1, jpkm1
             DO ji = 1, jpi
-               va(ji,jbdy1:jbdy2,jk) = ( va(ji,jbdy1:jbdy2,jk) & 
-                 & + va_b(ji,jbdy1:jbdy2) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
+               vv(ji,jbdy1:jbdy2,jk,Krhs) = ( vv(ji,jbdy1:jbdy2,jk,Krhs) & 
+                 & + vv_b(ji,jbdy1:jbdy2,Krhs) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
             END DO
          END DO
@@ -368,5 +368,5 @@
                DO ji = 1, jpi
                   zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) & 
-                     & + e3u_a(ji,jbdy1:jbdy2,jk) * ua(ji,jbdy1:jbdy2,jk) * umask(ji,jbdy1:jbdy2,jk)
+                     & + e3u(ji,jbdy1:jbdy2,jk,Krhs) * uu(ji,jbdy1:jbdy2,jk,Krhs) * umask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -377,6 +377,6 @@
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  ua(ji,jbdy1:jbdy2,jk) = ( ua(ji,jbdy1:jbdy2,jk) & 
-                    & + ua_b(ji,jbdy1:jbdy2) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
+                  uu(ji,jbdy1:jbdy2,jk,Krhs) = ( uu(ji,jbdy1:jbdy2,jk,Krhs) & 
+                    & + uu_b(ji,jbdy1:jbdy2,Krhs) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -385,6 +385,6 @@
          DO jk = 1, jpkm1              ! Mask domain edges
             DO ji = 1, jpi
-               ua(ji,nlcj  ,jk) = 0._wp
-               va(ji,nlcj-1,jk) = 0._wp
+               uu(ji,nlcj  ,jk,Krhs) = 0._wp
+               vv(ji,nlcj-1,jk,Krhs) = 0._wp
             END DO
          END DO 
@@ -520,5 +520,5 @@
          DO jj = 1, jpj
             DO ji = 2, indx
-               ssha(ji,jj) = hbdy_w(ji-1,jj)
+               ssh(ji,jj,Krhs) = hbdy_w(ji-1,jj)
             ENDDO
          ENDDO
@@ -530,5 +530,5 @@
          DO jj = 1, jpj
             DO ji = indx, nlci-1
-               ssha(ji,jj) = hbdy_e(ji-indx+1,jj)
+               ssh(ji,jj,Krhs) = hbdy_e(ji-indx+1,jj)
             ENDDO
          ENDDO
@@ -540,5 +540,5 @@
          DO jj = 2, indy
             DO ji = 1, jpi
-               ssha(ji,jj) = hbdy_s(ji,jj-1)
+               ssh(ji,jj,Krhs) = hbdy_s(ji,jj-1)
             ENDDO
          ENDDO
@@ -550,5 +550,5 @@
          DO jj = indy, nlcj-1
             DO ji = 1, jpi
-               ssha(ji,jj) = hbdy_n(ji,jj-indy+1)
+               ssh(ji,jj,Krhs) = hbdy_n(ji,jj-indy+1)
             ENDDO
          ENDDO
@@ -659,5 +659,5 @@
                DO jj=j1,j2
                  DO ji=i1,i2
-                       ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
+                       ptab(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm)
                  END DO
               END DO
@@ -669,5 +669,5 @@
            DO jj=j1,j2
               DO ji=i1,i2
-                 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
+                 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 
               END DO
            END DO
@@ -699,5 +699,5 @@
                   IF (tmask(iref,jref,jk) == 0) EXIT 
                   N_out = N_out + 1
-                  h_out(jk) = e3t_n(iref,jref,jk)
+                  h_out(jk) = e3t(iref,jref,jk,Kmm)
                ENDDO
                IF (N_in > 0) THEN
@@ -713,5 +713,5 @@
          !
          DO jn=1, jpts
-            tsa(i1:i2,j1:j2,1:jpk,jn)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
+            ts(i1:i2,j1:j2,1:jpk,jn,Krhs)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
          END DO
 
@@ -737,14 +737,14 @@
                ibdy = nlci-nbghostcells
                DO jn = 1, jpts
-                  tsa(ibdy+1,jmin:jmax,1:jpkm1,jn) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
+                  ts(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
                   DO jk = 1, jpkm1
                      DO jj = jmin,jmax
                         IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN
-                           tsa(ibdy,jj,jk,jn) = tsa(ibdy+1,jj,jk,jn) * tmask(ibdy,jj,jk)
+                           ts(ibdy,jj,jk,jn,Krhs) = ts(ibdy+1,jj,jk,jn,Krhs) * tmask(ibdy,jj,jk)
                         ELSE
-                           tsa(ibdy,jj,jk,jn)=(z4*tsa(ibdy+1,jj,jk,jn)+z3*tsa(ibdy-1,jj,jk,jn))*tmask(ibdy,jj,jk)
-                           IF( un(ibdy-1,jj,jk) > 0._wp ) THEN
-                              tsa(ibdy,jj,jk,jn)=( z6*tsa(ibdy-1,jj,jk,jn)+z5*tsa(ibdy+1,jj,jk,jn) & 
-                                                 + z7*tsa(ibdy-2,jj,jk,jn) ) * tmask(ibdy,jj,jk)
+                           ts(ibdy,jj,jk,jn,Krhs)=(z4*ts(ibdy+1,jj,jk,jn,Krhs)+z3*ts(ibdy-1,jj,jk,jn,Krhs))*tmask(ibdy,jj,jk)
+                           IF( uu(ibdy-1,jj,jk,Kmm) > 0._wp ) THEN
+                              ts(ibdy,jj,jk,jn,Krhs)=( z6*ts(ibdy-1,jj,jk,jn,Krhs)+z5*ts(ibdy+1,jj,jk,jn,Krhs) & 
+                                                 + z7*ts(ibdy-2,jj,jk,jn,Krhs) ) * tmask(ibdy,jj,jk)
                            ENDIF
                         ENDIF
@@ -752,5 +752,5 @@
                   END DO
                   ! Restore ghost points:
-                  tsa(ibdy+1,jmin:jmax,1:jpkm1,jn) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1)
+                  ts(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1)
                END DO
             ENDIF
@@ -766,14 +766,14 @@
                jbdy = nlcj-nbghostcells         
                DO jn = 1, jpts
-                  tsa(imin:imax,jbdy+1,1:jpkm1,jn) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
+                  ts(imin:imax,jbdy+1,1:jpkm1,jn,Krhs) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
                   DO jk = 1, jpkm1
                      DO ji = imin,imax
                         IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN
-                           tsa(ji,jbdy,jk,jn) = tsa(ji,jbdy+1,jk,jn) * tmask(ji,jbdy,jk)
+                           ts(ji,jbdy,jk,jn,Krhs) = ts(ji,jbdy+1,jk,jn,Krhs) * tmask(ji,jbdy,jk)
                         ELSE
-                           tsa(ji,jbdy,jk,jn)=(z4*tsa(ji,jbdy+1,jk,jn)+z3*tsa(ji,jbdy-1,jk,jn))*tmask(ji,jbdy,jk)        
-                           IF (vn(ji,jbdy-1,jk) > 0._wp ) THEN
-                              tsa(ji,jbdy,jk,jn)=( z6*tsa(ji,jbdy-1,jk,jn)+z5*tsa(ji,jbdy+1,jk,jn)  &
-                                                 + z7*tsa(ji,jbdy-2,jk,jn) ) * tmask(ji,jbdy,jk)
+                           ts(ji,jbdy,jk,jn,Krhs)=(z4*ts(ji,jbdy+1,jk,jn,Krhs)+z3*ts(ji,jbdy-1,jk,jn,Krhs))*tmask(ji,jbdy,jk)        
+                           IF (vv(ji,jbdy-1,jk,Kmm) > 0._wp ) THEN
+                              ts(ji,jbdy,jk,jn,Krhs)=( z6*ts(ji,jbdy-1,jk,jn,Krhs)+z5*ts(ji,jbdy+1,jk,jn,Krhs)  &
+                                                 + z7*ts(ji,jbdy-2,jk,jn,Krhs) ) * tmask(ji,jbdy,jk)
                            ENDIF
                         ENDIF
@@ -781,5 +781,5 @@
                   END DO
                   ! Restore ghost points:
-                  tsa(imin:imax,jbdy+1,1:jpkm1,jn) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1)
+                  ts(imin:imax,jbdy+1,1:jpkm1,jn,Krhs) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1)
                END DO
             ENDIF
@@ -795,14 +795,14 @@
                ibdy = 1+nbghostcells       
                DO jn = 1, jpts
-                  tsa(ibdy-1,jmin:jmax,1:jpkm1,jn) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
+                  ts(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
                   DO jk = 1, jpkm1
                      DO jj = jmin,jmax
                         IF( umask(ibdy,jj,jk) == 0._wp ) THEN
-                           tsa(ibdy,jj,jk,jn) = tsa(ibdy-1,jj,jk,jn) * tmask(ibdy,jj,jk)
+                           ts(ibdy,jj,jk,jn,Krhs) = ts(ibdy-1,jj,jk,jn,Krhs) * tmask(ibdy,jj,jk)
                         ELSE
-                           tsa(ibdy,jj,jk,jn)=(z4*tsa(ibdy-1,jj,jk,jn)+z3*tsa(ibdy+1,jj,jk,jn))*tmask(ibdy,jj,jk)        
-                           IF( un(ibdy,jj,jk) < 0._wp ) THEN
-                              tsa(ibdy,jj,jk,jn)=( z6*tsa(ibdy+1,jj,jk,jn)+z5*tsa(ibdy-1,jj,jk,jn) &
-                                                 + z7*tsa(ibdy+2,jj,jk,jn) ) * tmask(ibdy,jj,jk)
+                           ts(ibdy,jj,jk,jn,Krhs)=(z4*ts(ibdy-1,jj,jk,jn,Krhs)+z3*ts(ibdy+1,jj,jk,jn,Krhs))*tmask(ibdy,jj,jk)        
+                           IF( uu(ibdy,jj,jk,Kmm) < 0._wp ) THEN
+                              ts(ibdy,jj,jk,jn,Krhs)=( z6*ts(ibdy+1,jj,jk,jn,Krhs)+z5*ts(ibdy-1,jj,jk,jn,Krhs) &
+                                                 + z7*ts(ibdy+2,jj,jk,jn,Krhs) ) * tmask(ibdy,jj,jk)
                            ENDIF
                         ENDIF
@@ -810,5 +810,5 @@
                   END DO
                   ! Restore ghost points:
-                  tsa(ibdy-1,jmin:jmax,1:jpkm1,jn) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1)
+                  ts(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1)
                END DO
             ENDIF
@@ -824,14 +824,14 @@
                jbdy=1+nbghostcells        
                DO jn = 1, jpts
-                  tsa(imin:imax,jbdy-1,1:jpkm1,jn) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
+                  ts(imin:imax,jbdy-1,1:jpkm1,jn,Krhs) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
                   DO jk = 1, jpkm1      
                      DO ji = imin,imax
                         IF( vmask(ji,jbdy,jk) == 0._wp ) THEN
-                           tsa(ji,jbdy,jk,jn)=tsa(ji,jbdy-1,jk,jn) * tmask(ji,jbdy,jk)
+                           ts(ji,jbdy,jk,jn,Krhs)=ts(ji,jbdy-1,jk,jn,Krhs) * tmask(ji,jbdy,jk)
                         ELSE
-                           tsa(ji,jbdy,jk,jn)=(z4*tsa(ji,jbdy-1,jk,jn)+z3*tsa(ji,jbdy+1,jk,jn))*tmask(ji,jbdy,jk)
-                           IF( vn(ji,jbdy,jk) < 0._wp ) THEN
-                              tsa(ji,jbdy,jk,jn)=( z6*tsa(ji,jbdy+1,jk,jn)+z5*tsa(ji,jbdy-1,jk,jn) & 
-                                                 + z7*tsa(ji,jbdy+2,jk,jn) ) * tmask(ji,jbdy,jk)
+                           ts(ji,jbdy,jk,jn,Krhs)=(z4*ts(ji,jbdy-1,jk,jn,Krhs)+z3*ts(ji,jbdy+1,jk,jn,Krhs))*tmask(ji,jbdy,jk)
+                           IF( vv(ji,jbdy,jk,Kmm) < 0._wp ) THEN
+                              ts(ji,jbdy,jk,jn,Krhs)=( z6*ts(ji,jbdy+1,jk,jn,Krhs)+z5*ts(ji,jbdy-1,jk,jn,Krhs) & 
+                                                 + z7*ts(ji,jbdy+2,jk,jn,Krhs) ) * tmask(ji,jbdy,jk)
                            ENDIF
                         ENDIF
@@ -839,5 +839,5 @@
                   END DO
                   ! Restore ghost points:
-                  tsa(imin:imax,jbdy-1,1:jpkm1,jn) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1)
+                  ts(imin:imax,jbdy-1,1:jpkm1,jn,Krhs) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1)
                END DO
             ENDIF
@@ -861,5 +861,5 @@
       !
       IF( before) THEN
-         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
+         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm)
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -900,7 +900,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) 
+                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm)*umask(ji,jj,jk)) 
 # if defined key_vertical
-                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk))
+                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm))
 # endif
                END DO
@@ -928,5 +928,5 @@
          
               IF (N_in == 0) THEN
-                 ua(ji,jj,:) = 0._wp
+                 uu(ji,jj,:,Krhs) = 0._wp
                  CYCLE
               ENDIF
@@ -936,9 +936,9 @@
                  if (umask(iref,jj,jk) == 0) EXIT
                  N_out = N_out + 1
-                 h_out(N_out) = e3u_a(iref,jj,jk)
+                 h_out(N_out) = e3u(iref,jj,jk,Krhs)
               ENDDO
          
               IF (N_out == 0) THEN
-                 ua(ji,jj,:) = 0._wp
+                 uu(ji,jj,:,Krhs) = 0._wp
                  CYCLE
               ENDIF
@@ -952,5 +952,5 @@
                  endif
               ENDIF
-              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
+              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs),h_out(1:N_out),N_in,N_out)
             ENDDO
          ENDDO
@@ -959,5 +959,5 @@
          DO jk = 1, jpkm1
             DO jj=j1,j2
-               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) )
+               uu(i1:i2,jj,jk,Krhs) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs) )
             END DO
          END DO
@@ -992,7 +992,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk))
+                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm)*vmask(ji,jj,jk))
 # if defined key_vertical
-                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
+                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
 # endif
                END DO
@@ -1019,5 +1019,5 @@
                END DO
                IF (N_in == 0) THEN
-                  va(ji,jj,:) = 0._wp
+                  vv(ji,jj,:,Krhs) = 0._wp
                   CYCLE
                ENDIF
@@ -1027,16 +1027,16 @@
                   if (vmask(ji,jref,jk) == 0) EXIT
                   N_out = N_out + 1
-                  h_out(N_out) = e3v_a(ji,jref,jk)
+                  h_out(N_out) = e3v(ji,jref,jk,Krhs)
                END DO
                IF (N_out == 0) THEN
-                 va(ji,jj,:) = 0._wp
+                 vv(ji,jj,:,Krhs) = 0._wp
                  CYCLE
                ENDIF
-               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
+               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs),h_out(1:N_out),N_in,N_out)
             END DO
          END DO
 # else
          DO jk = 1, jpkm1
-            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) )
+            vv(i1:i2,j1:j2,jk,Krhs) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs) )
          END DO
 # endif
@@ -1060,5 +1060,5 @@
       !
       IF( before ) THEN 
-         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
+         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * uu_b(i1:i2,j1:j2,Kmm)
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -1113,5 +1113,5 @@
       ! 
       IF( before ) THEN 
-         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
+         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vv_b(i1:i2,j1:j2,Kmm)
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -1394,5 +1394,5 @@
            DO jj=j1,j2
               DO ji=i1,i2
-                 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w_n(ji,jj,jk) 
+                 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w(ji,jj,jk,Kmm) 
               END DO
            END DO
@@ -1415,5 +1415,5 @@
                   IF (wmask(ji,jj,jk) == 0) EXIT 
                   N_out = N_out + 1
-                  h_out(jk) = e3t_n(ji,jj,jk)
+                  h_out(jk) = e3t(ji,jj,jk,Kmm)
                ENDDO
                IF (N_in > 0) THEN
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_sponge.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_sponge.F90	(revision 10988)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_sponge.F90	(revision 10989)
@@ -218,5 +218,5 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     tabres(ji,jj,jk,jn) = tsb(ji,jj,jk,jn)
+                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb)
                   END DO
                END DO
@@ -228,5 +228,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
+                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 
                END DO
             END DO
@@ -251,5 +251,5 @@
                   IF (tmask(ji,jj,jk) == 0) EXIT 
                   N_out = N_out + 1
-                  h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above
+                  h_out(jk) = e3t(ji,jj,jk,Kmm) !Child grid scale factors. Could multiply by e1e2t here instead of division above
                ENDDO
                IF (N_in > 0) THEN
@@ -268,7 +268,7 @@
                DO jk=1,jpkm1
 # if defined key_vertical
-                  tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)
+                  tsbdiff(ji,jj,jk,1:jpts) = ts(ji,jj,jk,1:jpts,Kbb) - tabres_child(ji,jj,jk,1:jpts)
 # else
-                  tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts)
+                  tsbdiff(ji,jj,jk,1:jpts) = ts(ji,jj,jk,1:jpts,Kbb) - tabres(ji,jj,jk,1:jpts)
 # endif
                ENDDO
@@ -281,5 +281,5 @@
                DO jj = j1,j2
                   DO ji = i1,i2-1
-                     zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
+                     zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)
                      ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
                   END DO
@@ -288,5 +288,5 @@
                DO ji = i1,i2
                   DO jj = j1,j2-1
-                     zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
+                     zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
                      ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
                   END DO
@@ -310,9 +310,9 @@
                   DO ji = i1+1,i2-1
                      IF (.NOT. tabspongedone_tsn(ji,jj)) THEN 
-                        zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                         ! horizontal diffusive trends
                         ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
                         ! add it to the general tracer trends
-                        tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
+                        ts(ji,jj,jk,jn,Krhs) = ts(ji,jj,jk,jn,Krhs) + ztsa
                      ENDIF
                   END DO
@@ -353,7 +353,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,m1) = ub(ji,jj,jk)
+                  tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb)
 # if defined key_vertical
-                  tabres(ji,jj,jk,m2) = e3u_n(ji,jj,jk)*umask(ji,jj,jk)
+                  tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kmm)*umask(ji,jj,jk)
 # endif
                END DO
@@ -384,5 +384,5 @@
                  if (umask(ji,jj,jk) == 0) EXIT
                  N_out = N_out + 1
-                 h_out(N_out) = e3u_n(ji,jj,jk)
+                 h_out(N_out) = e3u(ji,jj,jk,Kmm)
               ENDDO
          
@@ -403,7 +403,7 @@
          ENDDO
 
-         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
+         ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
 #else
-         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
+         ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
 #endif
          !
@@ -416,7 +416,7 @@
             DO jj = j1,j2
                DO ji = i1+1,i2   ! vector opt.
-                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
-                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u_n(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) &
-                                     &   -e2u(ji-1,jj)*e3u_n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr
+                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * fsahm_spt(ji,jj)
+                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm) * ubdiff(ji  ,jj,jk) &
+                                     &   -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm) * ubdiff(ji-1,jj,jk) ) * zbtr
                END DO
             END DO
@@ -424,5 +424,5 @@
             DO jj = j1,j2-1
                DO ji = i1,i2   ! vector opt.
-                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
+                  zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * fsahm_spf(ji,jj)
                   rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk)   &
                                     &   +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) ) * fmask(ji,jj,jk) * zbtr 
@@ -439,9 +439,9 @@
                      ze1v = hdivdiff(ji,jj,jk)
                      ! horizontal diffusive trends
-                     zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )   &
+                     zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm) )   &
                            + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj)
 
                      ! add it to the general momentum trends
-                     ua(ji,jj,jk) = ua(ji,jj,jk) + zua
+                     uu(ji,jj,jk,Krhs) = uu(ji,jj,jk,Krhs) + zua
 
                   END DO
@@ -465,9 +465,9 @@
 
                      ! horizontal diffusive trends
-                     zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
+                     zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm) )   &
                            + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
 
                      ! add it to the general momentum trends
-                     va(ji,jj,jk) = va(ji,jj,jk) + zva
+                     vv(ji,jj,jk,Krhs) = vv(ji,jj,jk,Krhs) + zva
                   END DO
                ENDIF
@@ -506,7 +506,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,m1) = vb(ji,jj,jk)
+                  tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb)
 # if defined key_vertical
-                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_n(ji,jj,jk)
+                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kmm)
 # endif
                END DO
@@ -536,5 +536,5 @@
                  if (vmask(ji,jj,jk) == 0) EXIT
                  N_out = N_out + 1
-                 h_out(N_out) = e3v_n(ji,jj,jk)
+                 h_out(N_out) = e3v(ji,jj,jk,Kmm)
               ENDDO
          
@@ -549,7 +549,7 @@
          ENDDO
 
-         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)  
+         vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)  
 # else
-         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
+         vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
 # endif
          !
@@ -562,12 +562,12 @@
             DO jj = j1+1,j2
                DO ji = i1,i2   ! vector opt.
-                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
-                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  &
-                                     &  -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr
+                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * fsahm_spt(ji,jj)
+                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vbdiff(ji,jj  ,jk)  &
+                                     &  -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vbdiff(ji,jj-1,jk)  ) * zbtr
                END DO
             END DO
             DO jj = j1,j2
                DO ji = i1,i2-1   ! vector opt.
-                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
+                  zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * fsahm_spf(ji,jj)
                   rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
                                     &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk)  ) * fmask(ji,jj,jk) * zbtr
@@ -586,6 +586,6 @@
                IF( .NOT. tabspongedone_u(ji,jj) ) THEN
                   DO jk = 1, jpkm1
-                     ua(ji,jj,jk) = ua(ji,jj,jk)                                                               &
-                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )  &
+                     uu(ji,jj,jk,Krhs) = uu(ji,jj,jk,Krhs)                                                               &
+                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm) )  &
                         & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
                   END DO
@@ -600,6 +600,6 @@
                IF( .NOT. tabspongedone_v(ji,jj) ) THEN
                   DO jk = 1, jpkm1
-                     va(ji,jj,jk) = va(ji,jj,jk)                                                                  &
-                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
+                     vv(ji,jj,jk,Krhs) = vv(ji,jj,jk,Krhs)                                                                  &
+                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm) )   &
                         &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)
                   END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_update.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_update.F90	(revision 10988)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_update.F90	(revision 10989)
@@ -230,8 +230,8 @@
       ! -----------------------
       !
-      e3u_a(:,:,:) = e3u_n(:,:,:)
-      e3v_a(:,:,:) = e3v_n(:,:,:)
-!      ua(:,:,:) = e3u_b(:,:,:)
-!      va(:,:,:) = e3v_b(:,:,:)
+      e3u(:,:,:,Krhs) = e3u(:,:,:,Kmm)
+      e3v(:,:,:,Krhs) = e3v(:,:,:,Kmm)
+!      uu(:,:,:,Krhs) = e3u(:,:,:,Kbb)
+!      vv(:,:,:,Krhs) = e3v(:,:,:,Kbb)
       hu_a(:,:) = hu_n(:,:)
       hv_a(:,:) = hv_n(:,:)
@@ -242,10 +242,10 @@
          ! Vertical scale factor interpolations
          ! ------------------------------------
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) ,  'U' )
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) ,  'V' )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) ,  'F' )
-
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm) ,  'U' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm) ,  'V' )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:) ,  'F' )
+
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
 
          ! Update total depths:
@@ -254,6 +254,6 @@
       hv_n(:,:) = 0._wp                        ! Ocean depth at V-points
       DO jk = 1, jpkm1
-         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
+         hu_n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
+         hv_n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
       END DO
       !                                        ! Inverse of the local depth
@@ -268,9 +268,9 @@
          ! Vertical scale factor interpolations
          ! ------------------------------------
-         CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:),  'U'  )
-         CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:),  'V'  )
-
-         CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-         CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+         CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb),  'U'  )
+         CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb),  'V'  )
+
+         CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+         CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
 
          ! Update total depths:
@@ -279,6 +279,6 @@
          hv_b(:,:) = 0._wp                     ! Ocean depth at V-points
          DO jk = 1, jpkm1
-            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
-            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
+            hu_b(:,:) = hu_b(:,:) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
+            hv_b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
          END DO
          !                                     ! Inverse of the local depth
@@ -315,5 +315,5 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) &
+                     tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Kmm) ) &
                                            * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp
                   END DO
@@ -324,5 +324,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) &
+                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) &
                                            + (tmask(ji,jj,jk)-1)*999._wp
                END DO
@@ -345,5 +345,5 @@
                   IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF
                   N_out = N_out + 1
-                  h_out(N_out) = e3t_n(ji,jj,jk) 
+                  h_out(N_out) = e3t(ji,jj,jk,Kmm) 
                ENDDO
                IF (N_in > 0) THEN !Remove this?
@@ -369,7 +369,7 @@
                      DO ji=i1,i2
                         IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
-                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
+                           ts(ji,jj,jk,jn,Kbb) = ts(ji,jj,jk,jn,Kbb) & 
                                  & + atfp * ( tabres_child(ji,jj,jk,jn) &
-                                 &          - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
+                                 &          - ts(ji,jj,jk,jn,Kmm) ) * tmask(ji,jj,jk)
                         ENDIF
                      ENDDO
@@ -383,5 +383,5 @@
                   DO ji=i1,i2
                      IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
-                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
+                        ts(ji,jj,jk,jn,Kmm) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
                      END IF
                   END DO
@@ -413,6 +413,6 @@
                   DO ji=i1,i2
 !> jc tmp
-                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
-!                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk)
+                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm)  * e3t(ji,jj,jk,Kmm) / e3t_0(ji,jj,jk)
+!                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm)  * e3t(ji,jj,jk,Kmm)
 !< jc tmp
                   END DO
@@ -434,9 +434,9 @@
                      DO ji = i1, i2
                         IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
-                           ztb  = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used
+                           ztb  = ts(ji,jj,jk,jn,Kbb) * e3t(ji,jj,jk,Kbb) ! fse3t_b prior update should be used
                            ztnu = tabres(ji,jj,jk,jn)
-                           ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
-                           tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
-                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
+                           ztno = ts(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Krhs)
+                           ts(ji,jj,jk,jn,Kbb) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb)
                         ENDIF
                      END DO
@@ -450,5 +450,5 @@
                   DO ji=i1,i2
                      IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 
-                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
+                        ts(ji,jj,jk,jn,Kmm) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm)
                      END IF
                   END DO
@@ -458,5 +458,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            tsb(i1:i2,j1:j2,k1:k2,1:jpts)  = tsn(i1:i2,j1:j2,k1:k2,1:jpts)
+            ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb)  = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm)
          ENDIF
          !
@@ -495,7 +495,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  &
+                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm)  &
                                        + (umask(ji,jj,jk)-1)*999._wp
-                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  &
+                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)  &
                                        + (umask(ji,jj,jk)-1)*999._wp
                END DO
@@ -520,5 +520,5 @@
                   IF (umask(ji,jj,jk) == 0) EXIT
                   N_out = N_out + 1
-                  h_out(N_out) = e3u_n(ji,jj,jk)
+                  h_out(N_out) = e3u(ji,jj,jk,Kmm)
                ENDDO
                IF (N_in * N_out > 0) THEN
@@ -550,9 +550,9 @@
                DO ji=i1,i2
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     ub(ji,jj,jk) = ub(ji,jj,jk) & 
-                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
+                     uu(ji,jj,jk,Kbb) = uu(ji,jj,jk,Kbb) & 
+                           & + atfp * ( tabres_child(ji,jj,jk) - uu(ji,jj,jk,Kmm) ) * umask(ji,jj,jk)
                   ENDIF
                   !
-                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
+                  uu(ji,jj,jk,Kmm) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
                END DO
             END DO
@@ -579,5 +579,5 @@
          zrhoy = Agrif_Rhoy()
          DO jk = k1, k2
-            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
+            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm) * uu(i1:i2,j1:j2,jk,Kmm)
          END DO
       ELSE
@@ -588,12 +588,12 @@
                   !
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)  ! fse3t_b prior update should be used
-                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
+                     zub  = uu(ji,jj,jk,Kbb) * e3u(ji,jj,jk,Kbb)  ! fse3t_b prior update should be used
+                     zuno = uu(ji,jj,jk,Kmm) * e3u(ji,jj,jk,Krhs)
                      zunu = tabres(ji,jj,jk,1)
-                     ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &      
-                                    & * umask(ji,jj,jk) / e3u_b(ji,jj,jk)
+                     uu(ji,jj,jk,Kbb) = ( zub + atfp * ( zunu - zuno) ) &      
+                                    & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb)
                   ENDIF
                   !
-                  un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)
+                  uu(ji,jj,jk,Kmm) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -601,5 +601,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2)
+            uu(i1:i2,j1:j2,k1:k2,Kbb)  = uu(i1:i2,j1:j2,k1:k2,Kmm)
          ENDIF
          !
@@ -632,8 +632,8 @@
          IF (western_side) THEN
             DO jj=j1,j2
-               zcor = un_b(i1-1,jj) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - un_b(i1-1,jj)
-               un_b(i1-1,jj) = un_b(i1-1,jj) + zcor
+               zcor = uu_b(i1-1,jj,Kmm) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - uu_b(i1-1,jj,Kmm)
+               uu_b(i1-1,jj,Kmm) = uu_b(i1-1,jj,Kmm) + zcor
                DO jk=1,jpkm1
-                  un(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk)
+                  uu(i1-1,jj,jk,Kmm) = uu(i1-1,jj,jk,Kmm) + zcor * umask(i1-1,jj,jk)
                END DO 
             END DO
@@ -642,8 +642,8 @@
          IF (eastern_side) THEN
             DO jj=j1,j2
-               zcor = un_b(i2+1,jj) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - un_b(i2+1,jj)
-               un_b(i2+1,jj) = un_b(i2+1,jj) + zcor
+               zcor = uu_b(i2+1,jj,Kmm) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - uu_b(i2+1,jj,Kmm)
+               uu_b(i2+1,jj,Kmm) = uu_b(i2+1,jj,Kmm) + zcor
                DO jk=1,jpkm1
-                  un(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk)
+                  uu(i2+1,jj,jk,Kmm) = uu(i2+1,jj,jk,Kmm) + zcor * umask(i2+1,jj,jk)
                END DO 
             END DO
@@ -682,7 +682,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) &
+                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm) &
                                        + (vmask(ji,jj,jk)-1)*999._wp
-                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) &
+                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) &
                                        + (vmask(ji,jj,jk)-1)*999._wp
                END DO
@@ -705,5 +705,5 @@
                   IF (vmask(ji,jj,jk) == 0) EXIT
                   N_out = N_out + 1
-                  h_out(N_out) = e3v_n(ji,jj,jk)
+                  h_out(N_out) = e3v(ji,jj,jk,Kmm)
                ENDDO
                IF (N_in * N_out > 0) THEN
@@ -736,9 +736,9 @@
                   !
                   IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part
-                     vb(ji,jj,jk) = vb(ji,jj,jk) & 
-                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
+                     vv(ji,jj,jk,Kbb) = vv(ji,jj,jk,Kbb) & 
+                           & + atfp * ( tabres_child(ji,jj,jk) - vv(ji,jj,jk,Kmm) ) * vmask(ji,jj,jk)
                   ENDIF
                   !
-                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
+                  vv(ji,jj,jk,Kmm) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
                END DO
             END DO
@@ -767,5 +767,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
+                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -778,12 +778,12 @@
                   !
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used
-                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
+                     zvb  = vv(ji,jj,jk,Kbb) * e3v(ji,jj,jk,Kbb) ! fse3t_b prior update should be used
+                     zvno = vv(ji,jj,jk,Kmm) * e3v(ji,jj,jk,Krhs)
                      zvnu = tabres(ji,jj,jk,1)
-                     vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &      
-                                    & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk)
+                     vv(ji,jj,jk,Kbb) = ( zvb + atfp * ( zvnu - zvno) ) &      
+                                    & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb)
                   ENDIF
                   !
-                  vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)
+                  vv(ji,jj,jk,Kmm) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -791,5 +791,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2)
+            vv(i1:i2,j1:j2,k1:k2,Kbb)  = vv(i1:i2,j1:j2,k1:k2,Kmm)
          ENDIF
          !
@@ -822,8 +822,8 @@
          IF (southern_side) THEN
             DO ji=i1,i2
-               zcor = vn_b(ji,j1-1) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vn_b(ji,j1-1)
-               vn_b(ji,j1-1) = vn_b(ji,j1-1) + zcor
+               zcor = vv_b(ji,j1-1,Kmm) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vv_b(ji,j1-1,Kmm)
+               vv_b(ji,j1-1,Kmm) = vv_b(ji,j1-1,Kmm) + zcor
                DO jk=1,jpkm1
-                  vn(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk)
+                  vv(ji,j1-1,jk,Kmm) = vv(ji,j1-1,jk,Kmm) + zcor * vmask(ji,j1-1,jk)
                END DO 
             END DO
@@ -832,8 +832,8 @@
          IF (northern_side) THEN
             DO ji=i1,i2
-               zcor = vn_b(ji,j2+1) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vn_b(ji,j2+1)
-               vn_b(ji,j2+1) = vn_b(ji,j2+1) + zcor
+               zcor = vv_b(ji,j2+1,Kmm) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vv_b(ji,j2+1,Kmm)
+               vv_b(ji,j2+1,Kmm) = vv_b(ji,j2+1,Kmm) + zcor
                DO jk=1,jpkm1
-                  vn(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk)
+                  vv(ji,j2+1,jk,Kmm) = vv(ji,j2+1,jk,Kmm) + zcor * vmask(ji,j2+1,jk)
                END DO 
             END DO
@@ -862,5 +862,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
+               tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm) * hu_n(ji,jj) * e2u(ji,jj)
             END DO
          END DO
@@ -873,10 +873,10 @@
                spgu(ji,jj) = 0._wp
                DO jk=1,jpkm1
-                  spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
+                  spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm)
                END DO
                !
                zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj)
                DO jk=1,jpkm1              
-                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
+                  uu(ji,jj,jk,Kmm) = uu(ji,jj,jk,Kmm) + zcorr * umask(ji,jj,jk)           
                END DO
                !
@@ -884,19 +884,19 @@
                IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
-                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
+                     zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
+                     uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + atfp * zcorr * umask(ji,jj,1)
                   END IF
                ENDIF    
-               un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
+               uu_b(ji,jj,Kmm) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
                !       
                ! Correct "before" velocities to hold correct bt component:
                spgu(ji,jj) = 0.e0
                DO jk=1,jpkm1
-                  spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
+                  spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb)
                END DO
                !
-               zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)
+               zcorr = uu_b(ji,jj,Kbb) - spgu(ji,jj) * r1_hu_b(ji,jj)
                DO jk=1,jpkm1              
-                  ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
+                  uu(ji,jj,jk,Kbb) = uu(ji,jj,jk,Kbb) + zcorr * umask(ji,jj,jk)           
                END DO
                !
@@ -905,5 +905,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2)
+            uu_b(i1:i2,j1:j2,Kbb)  = uu_b(i1:i2,j1:j2,Kmm)
          ENDIF
       ENDIF
@@ -928,5 +928,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj) 
+               tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm) * hv_n(ji,jj) * e1v(ji,jj) 
             END DO
          END DO
@@ -939,10 +939,10 @@
                spgv(ji,jj) = 0.e0
                DO jk=1,jpkm1
-                  spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
+                  spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm)
                END DO
                !
                zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj)
                DO jk=1,jpkm1              
-                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
+                  vv(ji,jj,jk,Kmm) = vv(ji,jj,jk,Kmm) + zcorr * vmask(ji,jj,jk)           
                END DO
                !
@@ -950,19 +950,19 @@
                IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
-                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
+                     zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
+                     vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + atfp * zcorr * vmask(ji,jj,1)
                   END IF
                ENDIF              
-               vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
+               vv_b(ji,jj,Kmm) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
                !       
                ! Correct "before" velocities to hold correct bt component:
                spgv(ji,jj) = 0.e0
                DO jk=1,jpkm1
-                  spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
+                  spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb)
                END DO
                !
-               zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)
+               zcorr = vv_b(ji,jj,Kbb) - spgv(ji,jj) * r1_hv_b(ji,jj)
                DO jk=1,jpkm1              
-                  vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
+                  vv(ji,jj,jk,Kbb) = vv(ji,jj,jk,Kbb) + zcorr * vmask(ji,jj,jk)           
                END DO
                !
@@ -971,5 +971,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2)
+            vv_b(i1:i2,j1:j2,Kbb)  = vv_b(i1:i2,j1:j2,Kmm)
          ENDIF
          !
@@ -993,5 +993,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) = sshn(ji,jj)
+               tabres(ji,jj) = ssh(ji,jj,Kmm)
             END DO
          END DO
@@ -1000,6 +1000,6 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  sshb(ji,jj) =   sshb(ji,jj) &
-                        & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
+                  ssh(ji,jj,Kbb) =   ssh(ji,jj,Kbb) &
+                        & + atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm) ) * tmask(ji,jj,1)
                END DO
             END DO
@@ -1008,10 +1008,10 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
+               ssh(ji,jj,Kmm) = tabres(ji,jj) * tmask(ji,jj,1)
             END DO
          END DO
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2)
+            ssh(i1:i2,j1:j2,Kbb)  = ssh(i1:i2,j1:j2,Kmm)
          ENDIF
          !
@@ -1094,6 +1094,6 @@
             DO jj=j1,j2
                zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 
-               sshn(i1  ,jj) = sshn(i1  ,jj) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor
+               ssh(i1  ,jj,Kmm) = ssh(i1  ,jj,Kmm) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i1  ,jj,Kbb) = ssh(i1  ,jj,Kbb) + atfp * zcor
             END DO
          ENDIF
@@ -1101,6 +1101,6 @@
             DO jj=j1,j2
                zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
-               sshn(i2+1,jj) = sshn(i2+1,jj) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor
+               ssh(i2+1,jj,Kmm) = ssh(i2+1,jj,Kmm) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i2+1,jj,Kbb) = ssh(i2+1,jj,Kbb) + atfp * zcor
             END DO
          ENDIF
@@ -1182,6 +1182,6 @@
             DO ji=i1,i2
                zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
-               sshn(ji,j1  ) = sshn(ji,j1  ) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor
+               ssh(ji,j1  ,Kmm) = ssh(ji,j1  ,Kmm) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j1  ,Kbb) = ssh(ji,j1,Kbb) + atfp * zcor
             END DO
          ENDIF
@@ -1189,6 +1189,6 @@
             DO ji=i1,i2
                zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
-               sshn(ji,j2+1) = sshn(ji,j2+1) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor
+               ssh(ji,j2+1,Kmm) = ssh(ji,j2+1,Kmm) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j2+1,Kbb) = ssh(ji,j2+1,Kbb) + atfp * zcor
             END DO
          ENDIF
@@ -1319,5 +1319,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + sshn(ji,jj) &
+                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm) &
                                      & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
                END DO
@@ -1330,8 +1330,8 @@
          ! Save "old" scale factor (prior update) for subsequent asselin correction
          ! of prognostic variables
-         e3t_a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1)
-
-         ! One should also save e3t_b, but lacking of workspace...
-!         hdivn(i1:i2,j1:j2,1:jpkm1)   = e3t_b(i1:i2,j1:j2,1:jpkm1)
+         e3t(i1:i2,j1:j2,1:jpkm1,Krhs) = e3t(i1:i2,j1:j2,1:jpkm1,Kmm)
+
+         ! One should also save e3t(:,:,:,Kbb), but lacking of workspace...
+!         hdiv(i1:i2,j1:j2,1:jpkm1)   = e3t(i1:i2,j1:j2,1:jpkm1,Kbb)
 
          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
@@ -1339,13 +1339,13 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) &
-                           & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) )
+                     e3t(ji,jj,jk,Kbb) =  e3t(ji,jj,jk,Kbb) &
+                           & + atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm) )
                   END DO
                END DO
             END DO
             !
-            e3w_b  (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_b(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
-            gdepw_b(i1:i2,j1:j2,1) = 0.0_wp
-            gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1)
+            e3w  (i1:i2,j1:j2,1,Kbb) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kbb) - e3t_0(i1:i2,j1:j2,1)
+            gdepw(i1:i2,j1:j2,1,Kbb) = 0.0_wp
+            gdept(i1:i2,j1:j2,1,Kbb) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kbb)
             !
             DO jk = 2, jpk
@@ -1353,11 +1353,11 @@
                   DO ji = i1,i2            
                      zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-                     e3w_b(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
-                     &                                        ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )  &
+                     e3w(ji,jj,jk,Kbb)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
+                     &                                        ( e3t(ji,jj,jk-1,Kbb) - e3t_0(ji,jj,jk-1) )  &
                      &                                  +            0.5_wp * tmask(ji,jj,jk)   *        &
-                     &                                        ( e3t_b(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
-                     gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
-                     gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
-                         &               + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
+                     &                                        ( e3t(ji,jj,jk  ,Kbb) - e3t_0(ji,jj,jk  ) )
+                     gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb)
+                     gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  &
+                         &               + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb)) 
                   END DO
                END DO
@@ -1370,17 +1370,17 @@
          !
          ! Update vertical scale factor at T-points:
-         e3t_n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1)
+         e3t(i1:i2,j1:j2,1:jpkm1,Kmm) = ptab(i1:i2,j1:j2,1:jpkm1)
          !
          ! Update total depth:
          ht_n(i1:i2,j1:j2) = 0._wp
          DO jk = 1, jpkm1
-            ht_n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t_n(i1:i2,j1:j2,jk) * tmask(i1:i2,j1:j2,jk)
+            ht_n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t(i1:i2,j1:j2,jk,Kmm) * tmask(i1:i2,j1:j2,jk)
          END DO
          !
          ! Update vertical scale factor at W-points and depths:
-         e3w_n (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_n(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
-         gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1)
-         gdepw_n(i1:i2,j1:j2,1) = 0.0_wp
-         gde3w_n(i1:i2,j1:j2,1) = gdept_n(i1:i2,j1:j2,1) - (ht_n(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh
+         e3w (i1:i2,j1:j2,1,Kmm) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm) - e3t_0(i1:i2,j1:j2,1)
+         gdept(i1:i2,j1:j2,1,Kmm) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm)
+         gdepw(i1:i2,j1:j2,1,Kmm) = 0.0_wp
+         gde3w(i1:i2,j1:j2,1) = gdept(i1:i2,j1:j2,1,Kmm) - (ht_n(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh
          !
          DO jk = 2, jpk
@@ -1388,10 +1388,10 @@
                DO ji = i1,i2            
                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               e3w_n(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t_n(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )   &
-               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t_n(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
-               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
-                   &               + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
-               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
+               e3w(ji,jj,jk,Kmm)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t(ji,jj,jk-1,Kmm) - e3t_0(ji,jj,jk-1) )   &
+               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t(ji,jj,jk  ,Kmm) - e3t_0(ji,jj,jk  ) )
+               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
+               gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
+                   &               + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm)) 
+               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
                END DO
             END DO
@@ -1399,8 +1399,8 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk)
-            e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk)
-            gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk)
-            gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk)
+            e3t (i1:i2,j1:j2,1:jpk,Kbb)  = e3t (i1:i2,j1:j2,1:jpk,Kmm)
+            e3w (i1:i2,j1:j2,1:jpk,Kbb)  = e3w (i1:i2,j1:j2,1:jpk,Kmm)
+            gdepw(i1:i2,j1:j2,1:jpk,Kbb) = gdepw(i1:i2,j1:j2,1:jpk,Kmm)
+            gdept(i1:i2,j1:j2,1:jpk,Kbb) = gdept(i1:i2,j1:j2,1:jpk,Kmm)
          ENDIF
          !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_interp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_interp.F90	(revision 10988)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_interp.F90	(revision 10989)
@@ -73,5 +73,5 @@
                DO jj=j1,j2
                  DO ji=i1,i2
-                       ptab(ji,jj,jk,jn) = trn(ji,jj,jk,jn)
+                       ptab(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm)
                  END DO
               END DO
@@ -83,5 +83,5 @@
            DO jj=j1,j2
               DO ji=i1,i2
-                 ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
+                 ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 
               END DO
            END DO
@@ -113,5 +113,5 @@
                   IF (tmask(iref,jref,jk) == 0) EXIT 
                   N_out = N_out + 1
-                  h_out(jk) = e3t_n(iref,jref,jk)
+                  h_out(jk) = e3t(iref,jref,jk,Kmm)
                ENDDO
                IF (N_in > 0) THEN
@@ -127,5 +127,5 @@
          !
          DO jn=1, jptra
-            tra(i1:i2,j1:j2,1:jpk,jn)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
+            tr(i1:i2,j1:j2,1:jpk,jn,Krhs)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
          END DO
 
@@ -151,20 +151,20 @@
                ibdy = nlci-nbghostcells
                DO jn = 1, jptra
-                  tra(ibdy+1,jmin:jmax,1:jpkm1,jn) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
+                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
                   DO jk = 1, jpkm1
                      DO jj = jmin,jmax
                         IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN
-                           tra(ibdy,jj,jk,jn) = tra(ibdy+1,jj,jk,jn) * tmask(ibdy,jj,jk)
-                        ELSE
-                           tra(ibdy,jj,jk,jn)=(z4*tra(ibdy+1,jj,jk,jn)+z3*tra(ibdy-1,jj,jk,jn))*tmask(ibdy,jj,jk)
-                           IF( un(ibdy-1,jj,jk) > 0._wp ) THEN
-                              tra(ibdy,jj,jk,jn)=( z6*tra(ibdy-1,jj,jk,jn)+z5*tra(ibdy+1,jj,jk,jn) & 
-                                                 + z7*tra(ibdy-2,jj,jk,jn) ) * tmask(ibdy,jj,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tra(ibdy+1,jmin:jmax,1:jpkm1,jn) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1)
+                           tr(ibdy,jj,jk,jn,Krhs) = tr(ibdy+1,jj,jk,jn,Krhs) * tmask(ibdy,jj,jk)
+                        ELSE
+                           tr(ibdy,jj,jk,jn,Krhs)=(z4*tr(ibdy+1,jj,jk,jn,Krhs)+z3*tr(ibdy-1,jj,jk,jn,Krhs))*tmask(ibdy,jj,jk)
+                           IF( uu(ibdy-1,jj,jk,Kmm) > 0._wp ) THEN
+                              tr(ibdy,jj,jk,jn,Krhs)=( z6*tr(ibdy-1,jj,jk,jn,Krhs)+z5*tr(ibdy+1,jj,jk,jn,Krhs) & 
+                                                 + z7*tr(ibdy-2,jj,jk,jn,Krhs) ) * tmask(ibdy,jj,jk)
+                           ENDIF
+                        ENDIF
+                     END DO
+                  END DO
+                  ! Restore ghost points:
+                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1)
                END DO
             ENDIF
@@ -180,20 +180,20 @@
                jbdy = nlcj-nbghostcells         
                DO jn = 1, jptra
-                  tra(imin:imax,jbdy+1,1:jpkm1,jn) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
+                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
                   DO jk = 1, jpkm1
                      DO ji = imin,imax
                         IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN
-                           tra(ji,jbdy,jk,jn) = tra(ji,jbdy+1,jk,jn) * tmask(ji,jbdy,jk)
-                        ELSE
-                           tra(ji,jbdy,jk,jn)=(z4*tra(ji,jbdy+1,jk,jn)+z3*tra(ji,jbdy-1,jk,jn))*tmask(ji,jbdy,jk)        
-                           IF (vn(ji,jbdy-1,jk) > 0._wp ) THEN
-                              tra(ji,jbdy,jk,jn)=( z6*tra(ji,jbdy-1,jk,jn)+z5*tra(ji,jbdy+1,jk,jn)  &
-                                                 + z7*tra(ji,jbdy-2,jk,jn) ) * tmask(ji,jbdy,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tra(imin:imax,jbdy+1,1:jpkm1,jn) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1)
+                           tr(ji,jbdy,jk,jn,Krhs) = tr(ji,jbdy+1,jk,jn,Krhs) * tmask(ji,jbdy,jk)
+                        ELSE
+                           tr(ji,jbdy,jk,jn,Krhs)=(z4*tr(ji,jbdy+1,jk,jn,Krhs)+z3*tr(ji,jbdy-1,jk,jn,Krhs))*tmask(ji,jbdy,jk)        
+                           IF (vv(ji,jbdy-1,jk,Kmm) > 0._wp ) THEN
+                              tr(ji,jbdy,jk,jn,Krhs)=( z6*tr(ji,jbdy-1,jk,jn,Krhs)+z5*tr(ji,jbdy+1,jk,jn,Krhs)  &
+                                                 + z7*tr(ji,jbdy-2,jk,jn,Krhs) ) * tmask(ji,jbdy,jk)
+                           ENDIF
+                        ENDIF
+                     END DO
+                  END DO
+                  ! Restore ghost points:
+                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1)
                END DO
             ENDIF
@@ -209,20 +209,20 @@
                ibdy = 1+nbghostcells       
                DO jn = 1, jptra
-                  tra(ibdy-1,jmin:jmax,1:jpkm1,jn) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
+                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
                   DO jk = 1, jpkm1
                      DO jj = jmin,jmax
                         IF( umask(ibdy,jj,jk) == 0._wp ) THEN
-                           tra(ibdy,jj,jk,jn) = tra(ibdy-1,jj,jk,jn) * tmask(ibdy,jj,jk)
-                        ELSE
-                           tra(ibdy,jj,jk,jn)=(z4*tra(ibdy-1,jj,jk,jn)+z3*tra(ibdy+1,jj,jk,jn))*tmask(ibdy,jj,jk)        
-                           IF( un(ibdy,jj,jk) < 0._wp ) THEN
-                              tra(ibdy,jj,jk,jn)=( z6*tra(ibdy+1,jj,jk,jn)+z5*tra(ibdy-1,jj,jk,jn) &
-                                                 + z7*tra(ibdy+2,jj,jk,jn) ) * tmask(ibdy,jj,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tra(ibdy-1,jmin:jmax,1:jpkm1,jn) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1)
+                           tr(ibdy,jj,jk,jn,Krhs) = tr(ibdy-1,jj,jk,jn,Krhs) * tmask(ibdy,jj,jk)
+                        ELSE
+                           tr(ibdy,jj,jk,jn,Krhs)=(z4*tr(ibdy-1,jj,jk,jn,Krhs)+z3*tr(ibdy+1,jj,jk,jn,Krhs))*tmask(ibdy,jj,jk)        
+                           IF( uu(ibdy,jj,jk,Kmm) < 0._wp ) THEN
+                              tr(ibdy,jj,jk,jn,Krhs)=( z6*tr(ibdy+1,jj,jk,jn,Krhs)+z5*tr(ibdy-1,jj,jk,jn,Krhs) &
+                                                 + z7*tr(ibdy+2,jj,jk,jn,Krhs) ) * tmask(ibdy,jj,jk)
+                           ENDIF
+                        ENDIF
+                     END DO
+                  END DO
+                  ! Restore ghost points:
+                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1)
                END DO
             ENDIF
@@ -238,20 +238,20 @@
                jbdy=1+nbghostcells        
                DO jn = 1, jptra
-                  tra(imin:imax,jbdy-1,1:jpkm1,jn) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
+                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
                   DO jk = 1, jpkm1      
                      DO ji = imin,imax
                         IF( vmask(ji,jbdy,jk) == 0._wp ) THEN
-                           tra(ji,jbdy,jk,jn)=tra(ji,jbdy-1,jk,jn) * tmask(ji,jbdy,jk)
-                        ELSE
-                           tra(ji,jbdy,jk,jn)=(z4*tra(ji,jbdy-1,jk,jn)+z3*tra(ji,jbdy+1,jk,jn))*tmask(ji,jbdy,jk)
-                           IF( vn(ji,jbdy,jk) < 0._wp ) THEN
-                              tra(ji,jbdy,jk,jn)=( z6*tra(ji,jbdy+1,jk,jn)+z5*tra(ji,jbdy-1,jk,jn) & 
-                                                 + z7*tra(ji,jbdy+2,jk,jn) ) * tmask(ji,jbdy,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tra(imin:imax,jbdy-1,1:jpkm1,jn) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1)
+                           tr(ji,jbdy,jk,jn,Krhs)=tr(ji,jbdy-1,jk,jn,Krhs) * tmask(ji,jbdy,jk)
+                        ELSE
+                           tr(ji,jbdy,jk,jn,Krhs)=(z4*tr(ji,jbdy-1,jk,jn,Krhs)+z3*tr(ji,jbdy+1,jk,jn,Krhs))*tmask(ji,jbdy,jk)
+                           IF( vv(ji,jbdy,jk,Kmm) < 0._wp ) THEN
+                              tr(ji,jbdy,jk,jn,Krhs)=( z6*tr(ji,jbdy+1,jk,jn,Krhs)+z5*tr(ji,jbdy-1,jk,jn,Krhs) & 
+                                                 + z7*tr(ji,jbdy+2,jk,jn,Krhs) ) * tmask(ji,jbdy,jk)
+                           ENDIF
+                        ENDIF
+                     END DO
+                  END DO
+                  ! Restore ghost points:
+                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1)
                END DO
             ENDIF
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_sponge.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_sponge.F90	(revision 10988)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_sponge.F90	(revision 10989)
@@ -83,5 +83,5 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     tabres(ji,jj,jk,jn) = trb(ji,jj,jk,jn)
+                     tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kbb)
                   END DO
                END DO
@@ -93,5 +93,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
+                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 
                END DO
             END DO
@@ -114,5 +114,5 @@
                   IF (tmask(ji,jj,jk) == 0) EXIT 
                   N_out = N_out + 1
-                  h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above
+                  h_out(jk) = e3t(ji,jj,jk,Kmm) !Child grid scale factors. Could multiply by e1e2t here instead of division above
                ENDDO
                IF (N_in > 0) THEN
@@ -131,7 +131,7 @@
                DO jk=1,jpkm1
 # if defined key_vertical
-                  trbdiff(ji,jj,jk,1:jptra) = trb(ji,jj,jk,1:jptra) - tabres_child(ji,jj,jk,1:jptra)
+                  trbdiff(ji,jj,jk,1:jptra) = tr(ji,jj,jk,1:jptra,Kbb) - tabres_child(ji,jj,jk,1:jptra)
 # else
-                  trbdiff(ji,jj,jk,1:jptra) = trb(ji,jj,jk,1:jptra) - tabres(ji,jj,jk,1:jptra)
+                  trbdiff(ji,jj,jk,1:jptra) = tr(ji,jj,jk,1:jptra,Kbb) - tabres(ji,jj,jk,1:jptra)
 # endif
                ENDDO
@@ -143,6 +143,6 @@
                DO jj = j1,j2-1
                   DO ji = i1,i2-1
-                     zabe1 = fsaht_spu(ji,jj) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk)
-                     zabe2 = fsaht_spv(ji,jj) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk)
+                     zabe1 = fsaht_spu(ji,jj) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
+                     zabe2 = fsaht_spv(ji,jj) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
                      ztu(ji,jj) = zabe1 * ( trbdiff(ji+1,jj  ,jk,jn) - trbdiff(ji,jj,jk,jn) )
                      ztv(ji,jj) = zabe2 * ( trbdiff(ji  ,jj+1,jk,jn) - trbdiff(ji,jj,jk,jn) )
@@ -153,7 +153,7 @@
                   DO ji = i1+1,i2-1
                      IF( .NOT. tabspongedone_trn(ji,jj) ) THEN 
-                        tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + (  ztu(ji,jj) - ztu(ji-1,jj  )     &
+                        tr(ji,jj,jk,jn,Krhs) = tr(ji,jj,jk,jn,Krhs) + (  ztu(ji,jj) - ztu(ji-1,jj  )     &
                            &                                   + ztv(ji,jj) - ztv(ji  ,jj-1)  )  &
-                           &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                           &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                      ENDIF
                   END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_update.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_update.F90	(revision 10988)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_update.F90	(revision 10989)
@@ -84,5 +84,5 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     tabres(ji,jj,jk,jn) = (trn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) &
+                     tabres(ji,jj,jk,jn) = (tr(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Kmm) ) &
                                            * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp
                   END DO
@@ -93,5 +93,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) &
+                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) &
                                            + (tmask(ji,jj,jk)-1)*999._wp
                END DO
@@ -114,5 +114,5 @@
                   IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF
                   N_out = N_out + 1
-                  h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above
+                  h_out(N_out) = e3t(ji,jj,jk,Kmm) !Parent grid scale factors. Could multiply by e1e2t here instead of division above
                ENDDO
                IF (N_in > 0) THEN !Remove this?
@@ -138,7 +138,7 @@
                      DO ji=i1,i2
                         IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
-                           trb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
+                           tr(ji,jj,jk,jn,Kbb) = ts(ji,jj,jk,jn,Kbb) & 
                                  & + atfp * ( tabres_child(ji,jj,jk,jn) &
-                                 &          - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
+                                 &          - tr(ji,jj,jk,jn,Kmm) ) * tmask(ji,jj,jk)
                         ENDIF
                      ENDDO
@@ -152,5 +152,5 @@
                   DO ji=i1,i2
                      IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
-                        trn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
+                        tr(ji,jj,jk,jn,Kmm) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
                      END IF
                   END DO
@@ -183,6 +183,6 @@
                   DO ji=i1,i2
 !> jc tmp
-                     tabres(ji,jj,jk,jn) = trn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
-!                     tabres(ji,jj,jk,jn) = trn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk)
+                     tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm)  * e3t(ji,jj,jk,Kmm) / e3t_0(ji,jj,jk)
+!                     tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm)  * e3t(ji,jj,jk,Kmm)
 !< jc tmp
                   END DO
@@ -204,9 +204,9 @@
                      DO ji=i1,i2
                         IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
-                           ztb  = trb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used
+                           ztb  = tr(ji,jj,jk,jn,Kbb) * e3t(ji,jj,jk,Kbb) ! fse3t_b prior update should be used
                            ztnu = tabres(ji,jj,jk,jn)
-                           ztno = trn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
-                           trb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
-                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
+                           ztno = tr(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Krhs)
+                           tr(ji,jj,jk,jn,Kbb) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb)
                         ENDIF
                      ENDDO
@@ -220,5 +220,5 @@
                   DO ji=i1,i2
                      IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 
-                        trn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
+                        tr(ji,jj,jk,jn,Kmm) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm)
                      END IF
                   END DO
@@ -228,5 +228,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            trb(i1:i2,j1:j2,k1:k2,n1:n2)  = trn(i1:i2,j1:j2,k1:k2,n1:n2)
+            tr(i1:i2,j1:j2,k1:k2,n1:n2,Kbb)  = tr(i1:i2,j1:j2,k1:k2,n1:n2,Kmm)
          ENDIF
          !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_user.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_user.F90	(revision 10988)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_user.F90	(revision 10989)
@@ -175,6 +175,6 @@
    tabspongedone_tsn = .FALSE.
    CALL Agrif_Bc_variable(tsn_sponge_id,calledweight=1.,procname=interptsn_sponge)
-   ! reset tsa to zero
-   tsa(:,:,:,:) = 0.
+   ! reset ts(:,:,:,:,Krhs) to zero
+   ts(:,:,:,:,Krhs) = 0.
 
    Agrif_UseSpecialValue = ln_spc_dyn
@@ -191,5 +191,5 @@
    CALL Agrif_Bc_variable(sshn_id,calledweight=1., procname=interpsshn )
    hbdy_w(:,:) = 0.e0 ; hbdy_e(:,:) = 0.e0 ; hbdy_n(:,:) = 0.e0 ; hbdy_s(:,:) = 0.e0
-   ssha(:,:) = 0.e0
+   ssh(:,:,Krhs) = 0.e0
 
    IF ( ln_dynspg_ts ) THEN
@@ -207,6 +207,6 @@
    Agrif_UseSpecialValue = .FALSE. 
    ! reset velocities to zero
-   ua(:,:,:) = 0.
-   va(:,:,:) = 0.
+   uu(:,:,:,Krhs) = 0.
+   vv(:,:,:,Krhs) = 0.
 
    ! 3. Some controls
@@ -591,6 +591,6 @@
    tabspongedone_trn = .FALSE.
    CALL Agrif_Bc_variable(trn_sponge_id,calledweight=1.,procname=interptrn_sponge)
-   ! reset tsa to zero
-   tra(:,:,:,:) = 0.
+   ! reset ts(:,:,:,:,Krhs) to zero
+   tr(:,:,:,:,Krhs) = 0.
 
 
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diawri.F90	(revision 10988)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diawri.F90	(revision 10989)
@@ -57,6 +57,6 @@
    USE lib_mpp         ! MPP library
    USE timing          ! preformance summary
-   USE diurnal_bulk    ! diurnal warm layer
-   USE cool_skin       ! Cool skin
+   USE diu_bulk        ! diurnal warm layer
+   USE diu_coolskin    ! Cool skin
 
    IMPLICIT NONE
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diu_bulk.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diu_bulk.F90	(revision 10989)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diu_bulk.F90	(revision 10989)
@@ -0,0 +1,267 @@
+MODULE diu_bulk
+   !!======================================================================
+   !!                    ***  MODULE  diu_bulk  ***
+   !!     Takaya model of diurnal warming (Takaya, 2010)
+   !!=====================================================================
+   !! History :        !  11-10  (J. While)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   diurnal_sst_bulk_init  : initialise diurnal model
+   !!   diurnal_sst_bulk_step  : time-step the diurnal model
+   !!----------------------------------------------------------------------
+   USE par_kind
+   USE phycst
+   USE dom_oce
+   USE lib_mpp
+   USE solfrac_mod
+   USE in_out_manager
+   
+   IMPLICIT NONE
+   PRIVATE
+   
+   ! Namelist parameters
+   LOGICAL, PUBLIC :: ln_diurnal
+   LOGICAL, PUBLIC :: ln_diurnal_only
+
+   ! Parameters
+   REAL(wp), PRIVATE, PARAMETER :: pp_alpha = 2.0e-4_wp
+   REAL(wp), PRIVATE, PARAMETER :: pp_veltol = 0._wp
+   REAL(wp), PRIVATE, PARAMETER :: pp_min_fvel = 1.e-10_wp 
+   
+   ! Key variables
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_dsst     ! Delta SST
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_solfrac  ! Fraction of 
+   !                                                           ! absorbed radiation
+
+   PUBLIC diurnal_sst_bulk_init, diurnal_sst_takaya_step
+   
+   !!----------------------------------------------------------------------
+CONTAINS 
+   
+   SUBROUTINE diurnal_sst_bulk_init
+      !!----------------------------------------------------------------------
+      !! *** ROUTINE diurnal_sst_init ***
+      !!
+      !! ** Purpose : Initialise the Takaya diurnal model   
+      !!----------------------------------------------------------------------      
+      INTEGER ::   ios   ! local integer
+      !!
+      NAMELIST /namdiu/ ln_diurnal, ln_diurnal_only
+      !!----------------------------------------------------------------------      
+
+      ! Read the namelist
+      REWIND( numnam_ref )
+      READ  ( numnam_ref, namdiu, IOSTAT = ios, ERR = 901 )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdiu in reference namelist', lwp )
+      REWIND( numnam_cfg )
+      READ  ( numnam_cfg, namdiu, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdiu in configuration namelist', lwp )      
+      !
+      IF( ln_diurnal_only .AND. ( .NOT. ln_diurnal ) ) THEN
+         CALL ctl_stop( "ln_diurnal_only set, but ln_diurnal = FALSE !" )
+      ENDIF
+      
+      IF( ln_diurnal ) THEN      
+         !         
+         ALLOCATE( x_dsst(jpi,jpj), x_solfrac(jpi,jpj) )
+         !
+         x_solfrac = 0._wp         ! Initialise the solar fraction
+         x_dsst    = 0._wp
+         !
+         IF( ln_diurnal_only ) THEN
+            CALL ctl_warn( "ln_diurnal_only set; only the diurnal component of SST will be calculated" )
+         ENDIF
+      ENDIF
+      
+   END SUBROUTINE diurnal_sst_bulk_init
+
+
+   SUBROUTINE diurnal_sst_takaya_step(kt, psolflux, pqflux, ptauflux, prho, p_rdt,   &
+            &                  pla, pthick, pcoolthick, pmu, &
+            &                  p_fvel_bkginc, p_hflux_bkginc)
+      !!----------------------------------------------------------------------
+      !! *** ROUTINE diurnal_sst_takaya_step ***
+      !!
+      !! ** Purpose :   Time-step the Takaya diurnal model
+      !!
+      !! ** Method :    1) Calculate the Obukhov length
+      !!                2) Calculate the Similarity function
+      !!                2) Calculate the increment to dsst
+      !!                3) Apply the increment
+      !! ** Reference : Refinements to a prognostic scheme of skin sea surface
+      !!                temperature, Takaya et al, JGR, 2010 
+      !!----------------------------------------------------------------------
+      INTEGER                               , INTENT(in) ::   kt             ! time step
+      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   psolflux       ! solar flux (Watts)
+      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   pqflux         ! heat (non-solar) flux (Watts)
+      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   ptauflux       ! wind stress  (kg/ m s^2)
+      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   prho           ! water density  (kg/m^3)
+      REAL(wp)                              , INTENT(in) ::   p_rdt          ! time-step
+      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pLa            ! Langmuir number
+      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pthick         ! warm layer thickness (m)
+      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pcoolthick     ! cool skin thickness (m)
+      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pmu            ! mu parameter
+      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   p_hflux_bkginc ! increment to the heat flux
+      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   p_fvel_bkginc  ! increment to the friction velocity
+      !
+      INTEGER :: ji,jj
+      LOGICAL  :: ll_calcfrac
+      REAL(wp), DIMENSION(jpi,jpj) :: z_fvel              ! friction velocity     
+      REAL(wp), DIMENSION(jpi,jpj) :: zthick, zcoolthick, zmu, zla
+      REAL(wp), DIMENSION(jpi,jpj) :: z_abflux            ! absorbed flux           
+      REAL(wp), DIMENSION(jpi,jpj) :: z_fla               ! Langmuir function value 
+      !!----------------------------------------------------------------------
+
+      ! Set optional arguments to their defaults
+      IF( .NOT. PRESENT( pthick )   ) THEN   ;   zthick(:,:) = 3._wp
+      ELSE                                   ;   zthick(:,:) = pthick(:,:)
+      ENDIF
+      IF( .NOT. PRESENT(pcoolthick) ) THEN   ;   zcoolthick(:,:) = 0._wp
+      ELSE                                   ;   zcoolthick(:,:) = pcoolthick(:,:)
+      ENDIF
+      IF( .NOT. PRESENT( pmu )      ) THEN   ;   zmu(:,:) = 0.3_wp
+      ELSE                                   ;   zmu(:,:) = pmu(:,:)
+      ENDIF
+      IF( .NOT. PRESENT(pla) ) THEN          ;   zla(:,:) = 0.3_wp
+      ELSE                                   ;   zla(:,:) = pla(:,:)
+      ENDIF
+      
+      ! If not done already, calculate the solar fraction
+      IF ( kt==nit000 ) THEN
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               IF(  ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) &
+                  &   x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
+            END DO
+         END DO   
+      ENDIF   
+
+      ! convert solar flux and heat flux to absorbed flux   
+      WHERE ( tmask(:,:,1) == 1._wp) 
+         z_abflux(:,:) = (  x_solfrac(:,:) *  psolflux (:,:)) + pqflux(:,:)     
+      ELSEWHERE
+         z_abflux(:,:) = 0._wp
+      ENDWHERE
+      IF( PRESENT(p_hflux_bkginc) ) z_abflux(:,:) = z_abflux(:,:) + p_hflux_bkginc   ! Optional increment
+      WHERE ( ABS( z_abflux(:,:) ) < rsmall )
+         z_abflux(:,:) = rsmall
+      ENDWHERE 
+      
+      ! Calculate the friction velocity
+      WHERE ( (ptauflux /= 0) .AND. ( tmask(:,:,1) == 1.) )   
+         z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(:,:) )
+      ELSEWHERE
+         z_fvel(:,:) = 0._wp
+      ENDWHERE
+      IF( PRESENT(p_fvel_bkginc) ) z_fvel(:,:) = z_fvel(:,:) + p_fvel_bkginc   ! Optional increment
+      
+       
+       
+      ! Calculate the Langmuir function value
+      WHERE ( tmask(:,:,1) == 1.)
+         z_fla(:,:) = MAX( 1._wp, zla(:,:)**( -2._wp / 3._wp ) )  
+      ELSEWHERE
+         z_fla(:,:) = 0._wp
+      ENDWHERE     
+      
+      ! Increment the temperature using the implicit solution
+      x_dsst(:,:) = t_imp( x_dsst(:,:), p_rdt, z_abflux(:,:), z_fvel(:,:),   &
+         &                       z_fla(:,:), zmu(:,:), zthick(:,:), prho(:,:) )
+      !
+   END SUBROUTINE diurnal_sst_takaya_step
+
+   
+   FUNCTION t_imp(p_dsst, p_rdt, p_abflux, p_fvel, &
+                          p_fla, pmu, pthick, prho )
+                          
+      IMPLICIT NONE
+      
+      ! Function definition
+      REAL(wp), DIMENSION(jpi,jpj) :: t_imp
+      ! Dummy variables
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst     ! Delta SST
+      REAL(wp), INTENT(IN)                     :: p_rdt      ! Time-step
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux   ! Heat forcing
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel     ! Friction velocity
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fla      ! Langmuir number
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pmu        ! Structure parameter
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pthick     ! Layer thickness
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: prho       ! Water density
+   
+      ! Local variables
+      REAL(wp) :: z_olength          ! Obukhov length
+      REAL(wp) :: z_sigma, z_sigma2
+      REAL(wp) :: z_term1, z_term2     
+      REAL(wp) :: z_stabfunc         ! stability function value
+      REAL(wp) :: z_fvel      
+      
+      CHARACTER(200) :: warn_string
+       
+      INTEGER :: ji,jj
+                     
+      DO jj = 1, jpj
+         DO ji = 1, jpi   
+            
+            ! Only calculate outside tmask
+            IF ( tmask(ji,jj,1) /= 1._wp ) THEN
+               t_imp(ji,jj) = 0._wp
+               CYCLE   
+            END IF
+            
+            IF (p_fvel(ji,jj) < pp_min_fvel) THEN
+               z_fvel = pp_min_fvel
+               WRITE(warn_string,*) "diurnal_sst_takaya step: "&
+               &//"friction velocity < minimum\n" &
+               &//"Setting friction velocity =",pp_min_fvel
+               CALL ctl_warn(warn_string)
+               
+            ELSE
+               z_fvel = p_fvel(ji,jj)
+            ENDIF
+                 
+            ! Calculate the Obukhov length
+            IF ( (z_fvel < pp_veltol ) .AND. &
+            &    (p_dsst(ji,jj) > 0._wp) ) THEN
+               z_olength =  z_fvel  / &
+               &     SQRT( p_dsst(ji,jj) * vkarmn * grav * &
+               &             pp_alpha / ( 5._wp * pthick(ji,jj) ) )
+            ELSE
+               z_olength = &
+               &   ( prho(ji,jj) * rcp * z_fvel**3._wp ) / &
+               &   ( vkarmn * grav * pp_alpha *&
+               &     p_abflux(ji,jj) )
+            ENDIF
+             
+            ! Calculate the stability function          
+            z_sigma = pthick(ji,jj) / z_olength
+            z_sigma2 = z_sigma * z_sigma
+  
+            IF ( z_sigma >= 0. ) THEN
+               z_stabfunc = 1._wp + &
+               &  ( ( 5._wp * z_sigma + 4._wp * z_sigma2 ) / &
+               &  ( 1._wp + 3._wp * z_sigma + 0.25_wp * &
+               &    z_sigma2 ) )
+            ELSE
+               z_stabfunc = 1._wp / &
+               &                   SQRT( 1._wp - 16._wp * z_sigma )
+            ENDIF
+
+            ! Calculate the T increment
+            z_term1 = ( p_abflux(ji,jj) * ( pmu(ji,jj) + 1._wp)  / &
+            & ( pmu(ji,jj) * pthick(ji,jj) * prho(ji,jj) * rcp ) )
+            
+             
+            z_term2 = -( ( pmu(ji,jj) + 1._wp) * &
+            &                       ( vkarmn * z_fvel * p_fla(ji,jj) ) / &
+            &      ( pthick(ji,jj) * z_stabfunc ) )     
+          
+            t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / &
+                           ( 1._wp - p_rdt * z_term2 )
+
+         END DO
+      END DO
+      
+      END FUNCTION t_imp
+
+END MODULE diu_bulk
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diu_coolskin.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diu_coolskin.F90	(revision 10989)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diu_coolskin.F90	(revision 10989)
@@ -0,0 +1,145 @@
+MODULE diu_coolskin
+   !!======================================================================
+   !!                    ***  MODULE  diu_coolskin  ***
+   !!     Cool skin thickness and delta T correction using Artele et al. (2002)
+   !!     [see also Tu and Tsuang (2005)]
+   !!
+   !!=====================================================================
+   !! History :        !  2012-01  (P. Sykes)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   diurnal_sst_coolskin_init : initialisation of the cool skin
+   !!   diurnal_sst_coolskin_step : time-stepping  of the cool skin corrections
+   !!----------------------------------------------------------------------
+   USE par_kind
+   USE phycst
+   USE dom_oce
+   USE in_out_manager
+   USE sbc_oce
+   USE lib_mpp
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   ! Namelist parameters
+
+   ! Parameters
+   REAL(wp), PRIVATE, PARAMETER :: pp_k = 0.596_wp          ! Thermal conductivity of seawater
+   REAL(wp), PRIVATE, PARAMETER :: pp_v = 1.05e-6_wp        ! Kinematic viscosity of seawater
+   REAL(wp), PRIVATE, PARAMETER :: pp_C = 86400             ! seconds [see Tu and Tsuang (2005)]
+   REAL(wp), PRIVATE, PARAMETER :: pp_cw = 3993._wp         ! specific heat capacity of seawater
+   REAL(wp), PRIVATE, PARAMETER :: pp_h = 10._wp            ! reference depth [using 10m from Artale et al. (2002)]
+   REAL(wp), PRIVATE, PARAMETER :: pp_rhoa = 1.20421_wp     ! density of air (at 20C)
+   REAL(wp), PRIVATE, PARAMETER :: pp_cda = 1.45e-3_wp      ! assumed air-sea drag coefficient for calculating wind speed
+   
+   ! Key variables
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_csdsst    ! Cool skin delta SST
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_csthick   ! Cool skin thickness
+   PUBLIC diurnal_sst_coolskin_step, diurnal_sst_coolskin_init
+
+      !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------   
+   CONTAINS 
+   
+   SUBROUTINE diurnal_sst_coolskin_init
+      !!----------------------------------------------------------------------
+      !! *** ROUTINE diurnal_sst_coolskin_init ***
+      !!
+      !! ** Purpose :   initialise the cool skin model
+      !!
+      !! ** Method : 
+      !!
+      !! ** Reference :
+      !! 
+      !!----------------------------------------------------------------------
+      ALLOCATE( x_csdsst(jpi,jpj), x_csthick(jpi,jpj) )
+      x_csdsst = 0.
+      x_csthick = 0.
+      !
+   END SUBROUTINE diurnal_sst_coolskin_init
+
+
+   SUBROUTINE diurnal_sst_coolskin_step(psqflux, pstauflux, psrho, rdt)
+      !!----------------------------------------------------------------------
+      !! *** ROUTINE diurnal_sst_takaya_step ***
+      !!
+      !! ** Purpose :   Time-step the Artale cool skin model
+      !!
+      !! ** Method : 
+      !!
+      !! ** Reference : 
+      !!----------------------------------------------------------------------
+      ! Dummy variables
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psqflux     ! Heat (non-solar)(Watts)
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pstauflux   ! Wind stress (kg/ m s^2)
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psrho       ! Water density (kg/m^3)
+      REAL(wp), INTENT(IN) :: rdt                             ! Time-step
+     
+      ! Local variables
+      REAL(wp), DIMENSION(jpi,jpj) :: z_fv                    ! Friction velocity     
+      REAL(wp), DIMENSION(jpi,jpj) :: z_gamma                 ! Dimensionless function of wind speed
+      REAL(wp), DIMENSION(jpi,jpj) :: z_lamda                 ! Sauders (dimensionless) proportionality constant
+      REAL(wp), DIMENSION(jpi,jpj) :: z_wspd                  ! Wind speed (m/s)
+      REAL(wp) :: z_ztx                                       ! Temporary u wind stress
+      REAL(wp) :: z_zty                                       ! Temporary v wind stress
+      REAL(wp) :: z_zmod                                      ! Temporary total wind stress
+     
+      INTEGER :: ji,jj
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. ln_blk )   CALL ctl_stop("diu_coolskin.f90: diurnal flux processing only implemented for bulk forcing")
+      !
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            !
+            ! Calcualte wind speed from wind stress and friction velocity
+            IF( tmask(ji,jj,1) == 1. .AND. pstauflux(ji,jj) /= 0 .AND. psrho(ji,jj) /=0 ) THEN
+               z_fv(ji,jj) = SQRT( pstauflux(ji,jj) / psrho(ji,jj) )
+               z_wspd(ji,jj) = SQRT( pstauflux(ji,jj) / ( pp_cda * pp_rhoa ) )
+            ELSE
+               z_fv(ji,jj) = 0.
+               z_wspd(ji,jj) = 0.     
+            ENDIF
+            !
+            ! Calculate gamma function which is dependent upon wind speed
+            IF( tmask(ji,jj,1) == 1. ) THEN
+               IF( ( z_wspd(ji,jj) <= 7.5 ) ) z_gamma(ji,jj) = ( 0.2 * z_wspd(ji,jj) ) + 0.5
+               IF( ( z_wspd(ji,jj) > 7.5 ) .AND. ( z_wspd(ji,jj) < 10. ) ) z_gamma(ji,jj) = ( 1.6 * z_wspd(ji,jj) ) - 10.
+               IF( ( z_wspd(ji,jj) >= 10. ) ) z_gamma(ji,jj) = 6.
+            ENDIF
+            !
+            ! Calculate lamda function
+            IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 ) THEN
+               z_lamda(ji,jj) = ( z_fv(ji,jj) * pp_k * pp_C ) / ( z_gamma(ji,jj) * psrho(ji,jj) * pp_cw * pp_h * pp_v )
+            ELSE
+               z_lamda(ji,jj) = 0.
+            ENDIF
+            !
+            ! Calculate the cool skin thickness - only when heat flux is out of the ocean
+            IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 .AND. psqflux(ji,jj) < 0 ) THEN
+               x_csthick(ji,jj) = ( z_lamda(ji,jj) * pp_v ) / z_fv(ji,jj)
+            ELSE
+               x_csthick(ji,jj) = 0.
+            ENDIF
+            !
+            ! Calculate the cool skin correction - only when the heat flux is out of the ocean
+            IF( tmask(ji,jj,1) == 1. .AND. x_csthick(ji,jj) /= 0. .AND. psqflux(ji,jj) < 0. ) THEN
+               x_csdsst(ji,jj) = ( psqflux(ji,jj) * x_csthick(ji,jj) ) / pp_k
+             ELSE
+               x_csdsst(ji,jj) = 0.
+            ENDIF
+            !
+         END DO
+      END DO
+      !
+   END SUBROUTINE diurnal_sst_coolskin_step
+
+   !!======================================================================
+END MODULE diu_coolskin
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diu_layers.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diu_layers.F90	(revision 10989)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diu_layers.F90	(revision 10989)
@@ -0,0 +1,51 @@
+MODULE diu_layers
+   !!======================================================================
+   !!                       ***  MODULE diu_layers **
+   !! Apply coolskin and warm layer calculations
+   !!======================================================================
+   !! History :  3.7  ! 2015-11  (J. While)  Original code
+
+   USE diu_bulk     ! diurnal SST bulk routines  (diurnal_sst_takaya routine) 
+   USE diu_coolskin ! diurnal cool skin correction (diurnal_sst_coolskin routine)   
+   USE oce
+   USE iom
+   USE sbc_oce
+   USE sbcmod       ! surface boundary condition       (sbc     routine)
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC diurnal_sst_bulk_init      ! called by nemogcm.F90
+   PUBLIC diurnal_sst_coolskin_init  ! called by nemogcm.F90
+   PUBLIC diurnal_layers             ! called by step.F90 or step_diu.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: step_diu.F90 10922 2019-05-02 15:10:39Z acc $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+
+   CONTAINS
+
+   SUBROUTINE diurnal_layers( kstp ) 
+      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index 
+      !!---------------------------------------------------------------------- 
+      !!                     ***  ROUTINE diurnal_layers  *** 
+      !!                       
+      !! ** Purpose : - Apply coolskin and warm layer calculations
+      !!   
+      !!---------------------------------------------------------------------- 
+      
+      ! Cool skin
+
+      CALL diurnal_sst_coolskin_step( qns, taum, rhop(:,:,1), rdt)
+
+      CALL iom_put( "sst_wl"   , x_dsst               )    ! warm layer (write out before update below).
+      CALL iom_put( "sst_cs"   , x_csdsst             )    ! cool skin
+
+      ! Diurnal warm layer model       
+      CALL diurnal_sst_takaya_step( kstp, qsr, qns, taum, rhop(:,:,1), rdt) 
+
+   END SUBROUTINE diurnal_layers  
+   
+END MODULE diu_layers
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diurnal_bulk.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diurnal_bulk.F90	(revision 10988)
+++ 	(revision )
@@ -1,267 +1,0 @@
-MODULE diurnal_bulk
-   !!======================================================================
-   !!                    ***  MODULE  diurnal_bulk  ***
-   !!     Takaya model of diurnal warming (Takaya, 2010)
-   !!=====================================================================
-   !! History :        !  11-10  (J. While)  Original code
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   diurnal_sst_bulk_init  : initialise diurnal model
-   !!   diurnal_sst_bulk_step  : time-step the diurnal model
-   !!----------------------------------------------------------------------
-   USE par_kind
-   USE phycst
-   USE dom_oce
-   USE lib_mpp
-   USE solfrac_mod
-   USE in_out_manager
-   
-   IMPLICIT NONE
-   PRIVATE
-   
-   ! Namelist parameters
-   LOGICAL, PUBLIC :: ln_diurnal
-   LOGICAL, PUBLIC :: ln_diurnal_only
-
-   ! Parameters
-   REAL(wp), PRIVATE, PARAMETER :: pp_alpha = 2.0e-4_wp
-   REAL(wp), PRIVATE, PARAMETER :: pp_veltol = 0._wp
-   REAL(wp), PRIVATE, PARAMETER :: pp_min_fvel = 1.e-10_wp 
-   
-   ! Key variables
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_dsst     ! Delta SST
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_solfrac  ! Fraction of 
-   !                                                           ! absorbed radiation
-
-   PUBLIC diurnal_sst_bulk_init, diurnal_sst_takaya_step
-   
-   !!----------------------------------------------------------------------
-CONTAINS 
-   
-   SUBROUTINE diurnal_sst_bulk_init
-      !!----------------------------------------------------------------------
-      !! *** ROUTINE diurnal_sst_init ***
-      !!
-      !! ** Purpose : Initialise the Takaya diurnal model   
-      !!----------------------------------------------------------------------      
-      INTEGER ::   ios   ! local integer
-      !!
-      NAMELIST /namdiu/ ln_diurnal, ln_diurnal_only
-      !!----------------------------------------------------------------------      
-
-      ! Read the namelist
-      REWIND( numnam_ref )
-      READ  ( numnam_ref, namdiu, IOSTAT = ios, ERR = 901 )
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdiu in reference namelist', lwp )
-      REWIND( numnam_cfg )
-      READ  ( numnam_cfg, namdiu, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdiu in configuration namelist', lwp )      
-      !
-      IF( ln_diurnal_only .AND. ( .NOT. ln_diurnal ) ) THEN
-         CALL ctl_stop( "ln_diurnal_only set, but ln_diurnal = FALSE !" )
-      ENDIF
-      
-      IF( ln_diurnal ) THEN      
-         !         
-         ALLOCATE( x_dsst(jpi,jpj), x_solfrac(jpi,jpj) )
-         !
-         x_solfrac = 0._wp         ! Initialise the solar fraction
-         x_dsst    = 0._wp
-         !
-         IF( ln_diurnal_only ) THEN
-            CALL ctl_warn( "ln_diurnal_only set; only the diurnal component of SST will be calculated" )
-         ENDIF
-      ENDIF
-      
-   END SUBROUTINE diurnal_sst_bulk_init
-
-
-   SUBROUTINE diurnal_sst_takaya_step(kt, psolflux, pqflux, ptauflux, prho, p_rdt,   &
-            &                  pla, pthick, pcoolthick, pmu, &
-            &                  p_fvel_bkginc, p_hflux_bkginc)
-      !!----------------------------------------------------------------------
-      !! *** ROUTINE diurnal_sst_takaya_step ***
-      !!
-      !! ** Purpose :   Time-step the Takaya diurnal model
-      !!
-      !! ** Method :    1) Calculate the Obukhov length
-      !!                2) Calculate the Similarity function
-      !!                2) Calculate the increment to dsst
-      !!                3) Apply the increment
-      !! ** Reference : Refinements to a prognostic scheme of skin sea surface
-      !!                temperature, Takaya et al, JGR, 2010 
-      !!----------------------------------------------------------------------
-      INTEGER                               , INTENT(in) ::   kt             ! time step
-      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   psolflux       ! solar flux (Watts)
-      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   pqflux         ! heat (non-solar) flux (Watts)
-      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   ptauflux       ! wind stress  (kg/ m s^2)
-      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   prho           ! water density  (kg/m^3)
-      REAL(wp)                              , INTENT(in) ::   p_rdt          ! time-step
-      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pLa            ! Langmuir number
-      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pthick         ! warm layer thickness (m)
-      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pcoolthick     ! cool skin thickness (m)
-      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pmu            ! mu parameter
-      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   p_hflux_bkginc ! increment to the heat flux
-      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   p_fvel_bkginc  ! increment to the friction velocity
-      !
-      INTEGER :: ji,jj
-      LOGICAL  :: ll_calcfrac
-      REAL(wp), DIMENSION(jpi,jpj) :: z_fvel              ! friction velocity     
-      REAL(wp), DIMENSION(jpi,jpj) :: zthick, zcoolthick, zmu, zla
-      REAL(wp), DIMENSION(jpi,jpj) :: z_abflux            ! absorbed flux           
-      REAL(wp), DIMENSION(jpi,jpj) :: z_fla               ! Langmuir function value 
-      !!----------------------------------------------------------------------
-
-      ! Set optional arguments to their defaults
-      IF( .NOT. PRESENT( pthick )   ) THEN   ;   zthick(:,:) = 3._wp
-      ELSE                                   ;   zthick(:,:) = pthick(:,:)
-      ENDIF
-      IF( .NOT. PRESENT(pcoolthick) ) THEN   ;   zcoolthick(:,:) = 0._wp
-      ELSE                                   ;   zcoolthick(:,:) = pcoolthick(:,:)
-      ENDIF
-      IF( .NOT. PRESENT( pmu )      ) THEN   ;   zmu(:,:) = 0.3_wp
-      ELSE                                   ;   zmu(:,:) = pmu(:,:)
-      ENDIF
-      IF( .NOT. PRESENT(pla) ) THEN          ;   zla(:,:) = 0.3_wp
-      ELSE                                   ;   zla(:,:) = pla(:,:)
-      ENDIF
-      
-      ! If not done already, calculate the solar fraction
-      IF ( kt==nit000 ) THEN
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               IF(  ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) &
-                  &   x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
-            END DO
-         END DO   
-      ENDIF   
-
-      ! convert solar flux and heat flux to absorbed flux   
-      WHERE ( tmask(:,:,1) == 1._wp) 
-         z_abflux(:,:) = (  x_solfrac(:,:) *  psolflux (:,:)) + pqflux(:,:)     
-      ELSEWHERE
-         z_abflux(:,:) = 0._wp
-      ENDWHERE
-      IF( PRESENT(p_hflux_bkginc) ) z_abflux(:,:) = z_abflux(:,:) + p_hflux_bkginc   ! Optional increment
-      WHERE ( ABS( z_abflux(:,:) ) < rsmall )
-         z_abflux(:,:) = rsmall
-      ENDWHERE 
-      
-      ! Calculate the friction velocity
-      WHERE ( (ptauflux /= 0) .AND. ( tmask(:,:,1) == 1.) )   
-         z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(:,:) )
-      ELSEWHERE
-         z_fvel(:,:) = 0._wp
-      ENDWHERE
-      IF( PRESENT(p_fvel_bkginc) ) z_fvel(:,:) = z_fvel(:,:) + p_fvel_bkginc   ! Optional increment
-      
-       
-       
-      ! Calculate the Langmuir function value
-      WHERE ( tmask(:,:,1) == 1.)
-         z_fla(:,:) = MAX( 1._wp, zla(:,:)**( -2._wp / 3._wp ) )  
-      ELSEWHERE
-         z_fla(:,:) = 0._wp
-      ENDWHERE     
-      
-      ! Increment the temperature using the implicit solution
-      x_dsst(:,:) = t_imp( x_dsst(:,:), p_rdt, z_abflux(:,:), z_fvel(:,:),   &
-         &                       z_fla(:,:), zmu(:,:), zthick(:,:), prho(:,:) )
-      !
-   END SUBROUTINE diurnal_sst_takaya_step
-
-   
-   FUNCTION t_imp(p_dsst, p_rdt, p_abflux, p_fvel, &
-                          p_fla, pmu, pthick, prho )
-                          
-      IMPLICIT NONE
-      
-      ! Function definition
-      REAL(wp), DIMENSION(jpi,jpj) :: t_imp
-      ! Dummy variables
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst     ! Delta SST
-      REAL(wp), INTENT(IN)                     :: p_rdt      ! Time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux   ! Heat forcing
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel     ! Friction velocity
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fla      ! Langmuir number
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pmu        ! Structure parameter
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pthick     ! Layer thickness
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: prho       ! Water density
-   
-      ! Local variables
-      REAL(wp) :: z_olength          ! Obukhov length
-      REAL(wp) :: z_sigma, z_sigma2
-      REAL(wp) :: z_term1, z_term2     
-      REAL(wp) :: z_stabfunc         ! stability function value
-      REAL(wp) :: z_fvel      
-      
-      CHARACTER(200) :: warn_string
-       
-      INTEGER :: ji,jj
-                     
-      DO jj = 1, jpj
-         DO ji = 1, jpi   
-            
-            ! Only calculate outside tmask
-            IF ( tmask(ji,jj,1) /= 1._wp ) THEN
-               t_imp(ji,jj) = 0._wp
-               CYCLE   
-            END IF
-            
-            IF (p_fvel(ji,jj) < pp_min_fvel) THEN
-               z_fvel = pp_min_fvel
-               WRITE(warn_string,*) "diurnal_sst_takaya step: "&
-               &//"friction velocity < minimum\n" &
-               &//"Setting friction velocity =",pp_min_fvel
-               CALL ctl_warn(warn_string)
-               
-            ELSE
-               z_fvel = p_fvel(ji,jj)
-            ENDIF
-                 
-            ! Calculate the Obukhov length
-            IF ( (z_fvel < pp_veltol ) .AND. &
-            &    (p_dsst(ji,jj) > 0._wp) ) THEN
-               z_olength =  z_fvel  / &
-               &     SQRT( p_dsst(ji,jj) * vkarmn * grav * &
-               &             pp_alpha / ( 5._wp * pthick(ji,jj) ) )
-            ELSE
-               z_olength = &
-               &   ( prho(ji,jj) * rcp * z_fvel**3._wp ) / &
-               &   ( vkarmn * grav * pp_alpha *&
-               &     p_abflux(ji,jj) )
-            ENDIF
-             
-            ! Calculate the stability function          
-            z_sigma = pthick(ji,jj) / z_olength
-            z_sigma2 = z_sigma * z_sigma
-  
-            IF ( z_sigma >= 0. ) THEN
-               z_stabfunc = 1._wp + &
-               &  ( ( 5._wp * z_sigma + 4._wp * z_sigma2 ) / &
-               &  ( 1._wp + 3._wp * z_sigma + 0.25_wp * &
-               &    z_sigma2 ) )
-            ELSE
-               z_stabfunc = 1._wp / &
-               &                   SQRT( 1._wp - 16._wp * z_sigma )
-            ENDIF
-
-            ! Calculate the T increment
-            z_term1 = ( p_abflux(ji,jj) * ( pmu(ji,jj) + 1._wp)  / &
-            & ( pmu(ji,jj) * pthick(ji,jj) * prho(ji,jj) * rcp ) )
-            
-             
-            z_term2 = -( ( pmu(ji,jj) + 1._wp) * &
-            &                       ( vkarmn * z_fvel * p_fla(ji,jj) ) / &
-            &      ( pthick(ji,jj) * z_stabfunc ) )     
-          
-            t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / &
-                           ( 1._wp - p_rdt * z_term2 )
-
-         END DO
-      END DO
-      
-      END FUNCTION t_imp
-
-END MODULE diurnal_bulk
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diurnal_coolskin.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diurnal_coolskin.F90	(revision 10988)
+++ 	(revision )
@@ -1,146 +1,0 @@
-MODULE cool_skin
-   !!======================================================================
-   !!                    ***  MODULE  cool_skin  ***
-   !!     Cool skin thickness and delta T correction using Artele et al. (2002)
-   !!     [see also Tu and Tsuang (2005)]
-   !!
-   !!=====================================================================
-   !! History :        !  2012-01  (P. Sykes)  Original code
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   diurnal_sst_coolskin_init : initialisation of the cool skin
-   !!   diurnal_sst_coolskin_step : time-stepping  of the cool skin corrections
-   !!----------------------------------------------------------------------
-   USE par_kind
-   USE phycst
-   USE dom_oce
-   USE in_out_manager
-   USE sbc_oce
-   USE lib_mpp
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   
-   IMPLICIT NONE
-   PRIVATE
-
-   ! Namelist parameters
-
-   ! Parameters
-   REAL(wp), PRIVATE, PARAMETER :: pp_k = 0.596_wp          ! Thermal conductivity of seawater
-   REAL(wp), PRIVATE, PARAMETER :: pp_v = 1.05e-6_wp        ! Kinematic viscosity of seawater
-   REAL(wp), PRIVATE, PARAMETER :: pp_C = 86400             ! seconds [see Tu and Tsuang (2005)]
-   REAL(wp), PRIVATE, PARAMETER :: pp_cw = 3993._wp         ! specific heat capacity of seawater
-   REAL(wp), PRIVATE, PARAMETER :: pp_h = 10._wp            ! reference depth [using 10m from Artale et al. (2002)]
-   REAL(wp), PRIVATE, PARAMETER :: pp_rhoa = 1.20421_wp     ! density of air (at 20C)
-   REAL(wp), PRIVATE, PARAMETER :: pp_cda = 1.45e-3_wp      ! assumed air-sea drag coefficient for calculating wind speed
-   
-   ! Key variables
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_csdsst    ! Cool skin delta SST
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_csthick   ! Cool skin thickness
-
-   PUBLIC diurnal_sst_coolskin_step, diurnal_sst_coolskin_init
-
-      !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------   
-   CONTAINS 
-   
-   SUBROUTINE diurnal_sst_coolskin_init
-      !!----------------------------------------------------------------------
-      !! *** ROUTINE diurnal_sst_coolskin_init ***
-      !!
-      !! ** Purpose :   initialise the cool skin model
-      !!
-      !! ** Method : 
-      !!
-      !! ** Reference :
-      !! 
-      !!----------------------------------------------------------------------
-      ALLOCATE( x_csdsst(jpi,jpj), x_csthick(jpi,jpj) )
-      x_csdsst = 0.
-      x_csthick = 0.
-      !
-   END SUBROUTINE diurnal_sst_coolskin_init
-
-
-   SUBROUTINE diurnal_sst_coolskin_step(psqflux, pstauflux, psrho, rdt)
-      !!----------------------------------------------------------------------
-      !! *** ROUTINE diurnal_sst_takaya_step ***
-      !!
-      !! ** Purpose :   Time-step the Artale cool skin model
-      !!
-      !! ** Method : 
-      !!
-      !! ** Reference : 
-      !!----------------------------------------------------------------------
-      ! Dummy variables
-      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psqflux     ! Heat (non-solar)(Watts)
-      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pstauflux   ! Wind stress (kg/ m s^2)
-      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psrho       ! Water density (kg/m^3)
-      REAL(wp), INTENT(IN) :: rdt                             ! Time-step
-     
-      ! Local variables
-      REAL(wp), DIMENSION(jpi,jpj) :: z_fv                    ! Friction velocity     
-      REAL(wp), DIMENSION(jpi,jpj) :: z_gamma                 ! Dimensionless function of wind speed
-      REAL(wp), DIMENSION(jpi,jpj) :: z_lamda                 ! Sauders (dimensionless) proportionality constant
-      REAL(wp), DIMENSION(jpi,jpj) :: z_wspd                  ! Wind speed (m/s)
-      REAL(wp) :: z_ztx                                       ! Temporary u wind stress
-      REAL(wp) :: z_zty                                       ! Temporary v wind stress
-      REAL(wp) :: z_zmod                                      ! Temporary total wind stress
-     
-      INTEGER :: ji,jj
-      !!----------------------------------------------------------------------
-      !
-      IF( .NOT. ln_blk )   CALL ctl_stop("cool_skin.f90: diurnal flux processing only implemented for bulk forcing")
-      !
-      DO jj = 1,jpj
-         DO ji = 1,jpi
-            !
-            ! Calcualte wind speed from wind stress and friction velocity
-            IF( tmask(ji,jj,1) == 1. .AND. pstauflux(ji,jj) /= 0 .AND. psrho(ji,jj) /=0 ) THEN
-               z_fv(ji,jj) = SQRT( pstauflux(ji,jj) / psrho(ji,jj) )
-               z_wspd(ji,jj) = SQRT( pstauflux(ji,jj) / ( pp_cda * pp_rhoa ) )
-            ELSE
-               z_fv(ji,jj) = 0.
-               z_wspd(ji,jj) = 0.     
-            ENDIF
-            !
-            ! Calculate gamma function which is dependent upon wind speed
-            IF( tmask(ji,jj,1) == 1. ) THEN
-               IF( ( z_wspd(ji,jj) <= 7.5 ) ) z_gamma(ji,jj) = ( 0.2 * z_wspd(ji,jj) ) + 0.5
-               IF( ( z_wspd(ji,jj) > 7.5 ) .AND. ( z_wspd(ji,jj) < 10. ) ) z_gamma(ji,jj) = ( 1.6 * z_wspd(ji,jj) ) - 10.
-               IF( ( z_wspd(ji,jj) >= 10. ) ) z_gamma(ji,jj) = 6.
-            ENDIF
-            !
-            ! Calculate lamda function
-            IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 ) THEN
-               z_lamda(ji,jj) = ( z_fv(ji,jj) * pp_k * pp_C ) / ( z_gamma(ji,jj) * psrho(ji,jj) * pp_cw * pp_h * pp_v )
-            ELSE
-               z_lamda(ji,jj) = 0.
-            ENDIF
-            !
-            ! Calculate the cool skin thickness - only when heat flux is out of the ocean
-            IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 .AND. psqflux(ji,jj) < 0 ) THEN
-               x_csthick(ji,jj) = ( z_lamda(ji,jj) * pp_v ) / z_fv(ji,jj)
-            ELSE
-               x_csthick(ji,jj) = 0.
-            ENDIF
-            !
-            ! Calculate the cool skin correction - only when the heat flux is out of the ocean
-            IF( tmask(ji,jj,1) == 1. .AND. x_csthick(ji,jj) /= 0. .AND. psqflux(ji,jj) < 0. ) THEN
-               x_csdsst(ji,jj) = ( psqflux(ji,jj) * x_csthick(ji,jj) ) / pp_k
-             ELSE
-               x_csdsst(ji,jj) = 0.
-            ENDIF
-            !
-         END DO
-      END DO
-      !
-   END SUBROUTINE diurnal_sst_coolskin_step
-
-   !!======================================================================
-END MODULE cool_skin
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diurnal_layers.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIU/diurnal_layers.F90	(revision 10988)
+++ 	(revision )
@@ -1,51 +1,0 @@
-MODULE diu_layers
-   !!======================================================================
-   !!                       ***  MODULE diu_layers **
-   !! Apply coolskin and warm layer calculations
-   !!======================================================================
-   !! History :  3.7  ! 2015-11  (J. While)  Original code
-
-   USE diurnal_bulk    ! diurnal SST bulk routines  (diurnal_sst_takaya routine) 
-   USE cool_skin       ! diurnal cool skin correction (diurnal_sst_coolskin routine)   
-   USE oce
-   USE iom
-   USE sbc_oce
-   USE sbcmod           ! surface boundary condition       (sbc     routine)
-   
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC diurnal_sst_bulk_init      ! called by nemogcm.F90
-   PUBLIC diurnal_sst_coolskin_init  ! called by nemogcm.F90
-   PUBLIC diurnal_layers             ! called by step.F90 or step_diu.F90
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id: step_diu.F90 10922 2019-05-02 15:10:39Z acc $
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-
-   CONTAINS
-
-   SUBROUTINE diurnal_layers( kstp ) 
-      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index 
-      !!---------------------------------------------------------------------- 
-      !!                     ***  ROUTINE diurnal_layers  *** 
-      !!                       
-      !! ** Purpose : - Apply coolskin and warm layer calculations
-      !!   
-      !!---------------------------------------------------------------------- 
-      
-      ! Cool skin
-
-      CALL diurnal_sst_coolskin_step( qns, taum, rhop(:,:,1), rdt)
-
-      CALL iom_put( "sst_wl"   , x_dsst               )    ! warm layer (write out before update below).
-      CALL iom_put( "sst_cs"   , x_csdsst             )    ! cool skin
-
-      ! Diurnal warm layer model       
-      CALL diurnal_sst_takaya_step( kstp, qsr, qns, taum, rhop(:,:,1), rdt) 
-
-   END SUBROUTINE diurnal_layers  
-   
-END MODULE diu_layers
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/IOM/iom.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/IOM/iom.F90	(revision 10988)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/IOM/iom.F90	(revision 10989)
@@ -46,5 +46,5 @@
 #endif
    USE lib_fortran 
-   USE diurnal_bulk, ONLY : ln_diurnal_only, ln_diurnal
+   USE diu_bulk, ONLY : ln_diurnal_only, ln_diurnal
 
    IMPLICIT NONE
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/IOM/restart.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/IOM/restart.F90	(revision 10988)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/IOM/restart.F90	(revision 10989)
@@ -27,5 +27,5 @@
    USE in_out_manager  ! I/O manager
    USE iom             ! I/O module
-   USE diurnal_bulk
+   USE diu_bulk
    USE lib_mpp         ! distribued memory computing library
 
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcmod.F90	(revision 10988)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcmod.F90	(revision 10989)
@@ -59,5 +59,5 @@
    USE timing         ! Timing
    USE wet_dry
-   USE diurnal_bulk, ONLY:   ln_diurnal_only   ! diurnal SST diagnostic
+   USE diu_bulk, ONLY:   ln_diurnal_only   ! diurnal SST diagnostic
 
    IMPLICIT NONE
