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 11427)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_interp.F90	(revision 11428)
@@ -110,6 +110,7 @@
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) & 
-                      & + e3u(ibdy1:ibdy2,jj,jk,Krhs_a) * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)
+                  uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,   Krhs_a) & 
+                      &                        + e3u(ibdy1:ibdy2,jj,jk,Krhs_a) &
+                      &                        *  uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -122,5 +123,6 @@
             DO jk=1,jpkm1              ! Smooth
                DO jj=j1,j2
-                  uu(ibdy2,jj,jk,Krhs_a) = 0.25_wp*(uu(ibdy2-1,jj,jk,Krhs_a)+2._wp*uu(ibdy2,jj,jk,Krhs_a)+uu(ibdy2+1,jj,jk,Krhs_a))
+                  uu(ibdy2,jj,jk,Krhs_a) = 0.25_wp*( uu(ibdy2-1,jj,jk,Krhs_a)+2._wp*uu(ibdy2,jj,jk,Krhs_a) &
+                      &                            + uu(ibdy2+1,jj,jk,Krhs_a) )
                END DO
             END DO
@@ -130,6 +132,6 @@
          DO jk = 1, jpkm1
             DO jj = 1, jpj
-               zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) & 
-                  & + e3u(ibdy1:ibdy2,jj,jk,Krhs_a)  * uu(ibdy1:ibdy2,jj,jk,Krhs_a)*umask(ibdy1:ibdy2,jj,jk)
+               zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj)           +   e3u(ibdy1:ibdy2,jj,jk,Krhs_a)  & 
+                 &                  * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)
             END DO
          END DO
@@ -140,6 +142,8 @@
          DO jk = 1, jpkm1
             DO jj = 1, jpj
-               uu(ibdy1:ibdy2,jj,jk,Krhs_a) = ( uu(ibdy1:ibdy2,jj,jk,Krhs_a) &
-                 & + uu_b(ibdy1:ibdy2,jj,Krhs_a)-zub(ibdy1:ibdy2,jj)) * umask(ibdy1:ibdy2,jj,jk)
+               uu(ibdy1:ibdy2,jj,jk,Krhs_a) = (    uu(ibdy1:ibdy2,jj,jk,Krhs_a)   &
+                 &                             + uu_b(ibdy1:ibdy2,jj,   Krhs_a)   &
+                 &                             -  zub(ibdy1:ibdy2,jj)           ) &
+                 &                            * umask(ibdy1:ibdy2,jj,jk)
             END DO
          END DO
@@ -149,6 +153,6 @@
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) & 
-                     & + e3v(ibdy1:ibdy2,jj,jk,Krhs_a) * vv(ibdy1:ibdy2,jj,jk,Krhs_a) * vmask(ibdy1:ibdy2,jj,jk)
+                  zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj)           +   e3v(ibdy1:ibdy2,jj,jk,Krhs_a)  &
+                    &                  * vv(ibdy1:ibdy2,jj,jk,Krhs_a) * vmask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -158,6 +162,8 @@
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  vv(ibdy1:ibdy2,jj,jk,Krhs_a) = ( vv(ibdy1:ibdy2,jj,jk,Krhs_a) & 
-                    & + vv_b(ibdy1:ibdy2,jj,Krhs_a)-zvb(ibdy1:ibdy2,jj))*vmask(ibdy1:ibdy2,jj,jk)
+                  vv(ibdy1:ibdy2,jj,jk,Krhs_a) = (    vv(ibdy1:ibdy2,jj,jk,Krhs_a)   & 
+                    &                             + vv_b(ibdy1:ibdy2,jj,   Krhs_a)   &
+                    &                             -  zvb(ibdy1:ibdy2,jj)           ) &
+                    &                            * vmask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -181,6 +187,8 @@
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) & 
-                      & + e3u(ibdy1:ibdy2,jj,jk,Krhs_a) * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)
+                  uu_b(ibdy1:ibdy2,jj,Krhs_a) =   uu_b(ibdy1:ibdy2,jj,   Krhs_a) & 
+                      &                        +   e3u(ibdy1:ibdy2,jj,jk,Krhs_a) &
+                      &                        *    uu(ibdy1:ibdy2,jj,jk,Krhs_a) &
+                      &                        * umask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -193,5 +201,7 @@
             DO jk=1,jpkm1              ! Smooth
                DO jj=j1,j2
-                  uu(ibdy1,jj,jk,Krhs_a) = 0.25_wp*(uu(ibdy1-1,jj,jk,Krhs_a)+2._wp*uu(ibdy1,jj,jk,Krhs_a)+uu(ibdy1+1,jj,jk,Krhs_a))
+                  uu(ibdy1,jj,jk,Krhs_a) = 0.25_wp*(        uu(ibdy1-1,jj,jk,Krhs_a)  &
+                    &                               + 2._wp*uu(ibdy1  ,jj,jk,Krhs_a)  &
+                    &                                     + uu(ibdy1+1,jj,jk,Krhs_a) )
                END DO
             END DO
@@ -201,6 +211,7 @@
          DO jk = 1, jpkm1
             DO jj = 1, jpj
-               zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) & 
-                  & + e3u(ibdy1:ibdy2,jj,jk,Krhs_a)  * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)
+               zub(ibdy1:ibdy2,jj) =  zub(ibdy1:ibdy2,jj)                                      & 
+                  &                 + e3u(ibdy1:ibdy2,jj,jk,Krhs_a)                            &
+                  &                 *  uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)
             END DO
          END DO
@@ -211,6 +222,8 @@
          DO jk = 1, jpkm1
             DO jj = 1, jpj
-               uu(ibdy1:ibdy2,jj,jk,Krhs_a) = ( uu(ibdy1:ibdy2,jj,jk,Krhs_a) & 
-                 & + uu_b(ibdy1:ibdy2,jj,Krhs_a)-zub(ibdy1:ibdy2,jj))*umask(ibdy1:ibdy2,jj,jk)
+               uu(ibdy1:ibdy2,jj,jk,Krhs_a) = (      uu(ibdy1:ibdy2,jj,jk,Krhs_a) & 
+                 &                             +   uu_b(ibdy1:ibdy2,jj,   Krhs_a) &
+                 &                             -    zub(ibdy1:ibdy2,jj)           &
+                 &                            ) * umask(ibdy1:ibdy2,jj,jk)
             END DO
          END DO
@@ -222,6 +235,8 @@
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) &
-                     & + e3v(ibdy1:ibdy2,jj,jk,Krhs_a) * vv(ibdy1:ibdy2,jj,jk,Krhs_a) * vmask(ibdy1:ibdy2,jj,jk)
+                  zvb(ibdy1:ibdy2,jj) =    zvb(ibdy1:ibdy2,jj)                     &
+                     &                 +   e3v(ibdy1:ibdy2,jj,jk,Krhs_a)           &
+                     &                 *    vv(ibdy1:ibdy2,jj,jk,Krhs_a)           &
+                     &                 * vmask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -231,6 +246,8 @@
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  vv(ibdy1:ibdy2,jj,jk,Krhs_a) = ( vv(ibdy1:ibdy2,jj,jk,Krhs_a) & 
-                      & + vv_b(ibdy1:ibdy2,jj,Krhs_a)-zvb(ibdy1:ibdy2,jj)) * vmask(ibdy1:ibdy2,jj,jk)
+                  vv(ibdy1:ibdy2,jj,jk,Krhs_a) = (      vv(ibdy1:ibdy2,jj,jk,Krhs_a) & 
+                      &                           +   vv_b(ibdy1:ibdy2,jj,   Krhs_a) &
+                      &                           -    zvb(ibdy1:ibdy2,jj)           &
+                      &                          ) * vmask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -254,6 +271,8 @@
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) & 
-                      & + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) * vv(ji,jbdy1:jbdy2,jk,Krhs_a) * vmask(ji,jbdy1:jbdy2,jk)
+                  vv_b(ji,jbdy1:jbdy2,Krhs_a) =   vv_b(ji,jbdy1:jbdy2,   Krhs_a) & 
+                      &                        +   e3v(ji,jbdy1:jbdy2,jk,Krhs_a) &
+                      &                        *    vv(ji,jbdy1:jbdy2,jk,Krhs_a) &
+                      &                        * vmask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -266,5 +285,7 @@
             DO jk = 1, jpkm1           ! Smooth
                DO ji = i1, i2
-                  vv(ji,jbdy2,jk,Krhs_a) = 0.25_wp*(vv(ji,jbdy2-1,jk,Krhs_a)+2._wp*vv(ji,jbdy2,jk,Krhs_a)+vv(ji,jbdy2+1,jk,Krhs_a))
+                  vv(ji,jbdy2,jk,Krhs_a) = 0.25_wp*(        vv(ji,jbdy2-1,jk,Krhs_a)  &
+                    &                               + 2._wp*vv(ji,jbdy2  ,jk,Krhs_a)  &
+                    &                               +       vv(ji,jbdy2+1,jk,Krhs_a) )
                END DO
             END DO
@@ -274,6 +295,8 @@
          DO jk=1,jpkm1
             DO ji=1,jpi
-               zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) & 
-                  & + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) * vv(ji,jbdy1:jbdy2,jk,Krhs_a) * vmask(ji,jbdy1:jbdy2,jk)
+               zvb(ji,jbdy1:jbdy2) =    zvb(ji,jbdy1:jbdy2)           & 
+                  &                 +   e3v(ji,jbdy1:jbdy2,jk,Krhs_a) &
+                  &                 *    vv(ji,jbdy1:jbdy2,jk,Krhs_a) &
+                  &                 * vmask(ji,jbdy1:jbdy2,jk)
             END DO
          END DO
@@ -284,6 +307,8 @@
          DO jk = 1, jpkm1
             DO ji = 1, jpi
-               vv(ji,jbdy1:jbdy2,jk,Krhs_a) = ( vv(ji,jbdy1:jbdy2,jk,Krhs_a) & 
-                 & + vv_b(ji,jbdy1:jbdy2,Krhs_a) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
+               vv(ji,jbdy1:jbdy2,jk,Krhs_a) = (      vv(ji,jbdy1:jbdy2,jk,Krhs_a) & 
+                 &                             +   vv_b(ji,jbdy1:jbdy2,   Krhs_a) &
+                 &                             -    zvb(ji,jbdy1:jbdy2)           &
+                 &                            ) * vmask(ji,jbdy1:jbdy2,jk)
             END DO
          END DO
@@ -293,6 +318,8 @@
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) & 
-                     & + e3u(ji,jbdy1:jbdy2,jk,Krhs_a) * uu(ji,jbdy1:jbdy2,jk,Krhs_a) * umask(ji,jbdy1:jbdy2,jk)
+                  zub(ji,jbdy1:jbdy2) =    zub(ji,jbdy1:jbdy2)           & 
+                     &                 +   e3u(ji,jbdy1:jbdy2,jk,Krhs_a) &
+                     &                 *    uu(ji,jbdy1:jbdy2,jk,Krhs_a) &
+                     &                 * umask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -303,6 +330,8 @@
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  uu(ji,jbdy1:jbdy2,jk,Krhs_a) = ( uu(ji,jbdy1:jbdy2,jk,Krhs_a) & 
-                    & + uu_b(ji,jbdy1:jbdy2,Krhs_a) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
+                  uu(ji,jbdy1:jbdy2,jk,Krhs_a) = (      uu(ji,jbdy1:jbdy2,jk,Krhs_a) & 
+                    &                             +   uu_b(ji,jbdy1:jbdy2,   Krhs_a) &
+                    &                             -    zub(ji,jbdy1:jbdy2)           &
+                    &                            ) * umask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -326,6 +355,8 @@
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) & 
-                      & + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) * vv(ji,jbdy1:jbdy2,jk,Krhs_a) * vmask(ji,jbdy1:jbdy2,jk)
+                  vv_b(ji,jbdy1:jbdy2,Krhs_a) =   vv_b(ji,jbdy1:jbdy2,   Krhs_a) & 
+                      &                        +   e3v(ji,jbdy1:jbdy2,jk,Krhs_a) &
+                      &                        *    vv(ji,jbdy1:jbdy2,jk,Krhs_a) &
+                      &                        * vmask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -338,5 +369,7 @@
             DO jk = 1, jpkm1           ! Smooth
                DO ji = i1, i2
-                  vv(ji,jbdy1,jk,Krhs_a) = 0.25_wp*(vv(ji,jbdy1-1,jk,Krhs_a)+2._wp*vv(ji,jbdy1,jk,Krhs_a)+vv(ji,jbdy1+1,jk,Krhs_a))
+                  vv(ji,jbdy1,jk,Krhs_a) = 0.25_wp*(        vv(ji,jbdy1-1,jk,Krhs_a)  &
+                    &                               + 2._wp*vv(ji,jbdy1  ,jk,Krhs_a)  &
+                    &                               +       vv(ji,jbdy1+1,jk,Krhs_a) )
                END DO
             END DO
@@ -346,6 +379,8 @@
          DO jk=1,jpkm1
             DO ji=1,jpi
-               zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) & 
-                  & + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) * vv(ji,jbdy1:jbdy2,jk,Krhs_a) * vmask(ji,jbdy1:jbdy2,jk)
+               zvb(ji,jbdy1:jbdy2) =    zvb(ji,jbdy1:jbdy2)           & 
+                  &                 +   e3v(ji,jbdy1:jbdy2,jk,Krhs_a) &
+                  &                 *    vv(ji,jbdy1:jbdy2,jk,Krhs_a) &
+                  &                 * vmask(ji,jbdy1:jbdy2,jk)
             END DO
          END DO
@@ -356,6 +391,8 @@
          DO jk = 1, jpkm1
             DO ji = 1, jpi
-               vv(ji,jbdy1:jbdy2,jk,Krhs_a) = ( vv(ji,jbdy1:jbdy2,jk,Krhs_a) & 
-                 & + vv_b(ji,jbdy1:jbdy2,Krhs_a) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
+               vv(ji,jbdy1:jbdy2,jk,Krhs_a) = (      vv(ji,jbdy1:jbdy2,jk,Krhs_a) & 
+                 &                             +   vv_b(ji,jbdy1:jbdy2,   Krhs_a) &
+                 &                             -    zvb(ji,jbdy1:jbdy2)           &
+                 &                            ) * vmask(ji,jbdy1:jbdy2,jk)
             END DO
          END DO
@@ -367,6 +404,8 @@
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) & 
-                     & + e3u(ji,jbdy1:jbdy2,jk,Krhs_a) * uu(ji,jbdy1:jbdy2,jk,Krhs_a) * umask(ji,jbdy1:jbdy2,jk)
+                  zub(ji,jbdy1:jbdy2) =    zub(ji,jbdy1:jbdy2)           & 
+                     &                 +   e3u(ji,jbdy1:jbdy2,jk,Krhs_a) &
+                     &                 *    uu(ji,jbdy1:jbdy2,jk,Krhs_a) &
+                     &                 * umask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -377,6 +416,8 @@
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  uu(ji,jbdy1:jbdy2,jk,Krhs_a) = ( uu(ji,jbdy1:jbdy2,jk,Krhs_a) & 
-                    & + uu_b(ji,jbdy1:jbdy2,Krhs_a) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
+                  uu(ji,jbdy1:jbdy2,jk,Krhs_a) = (      uu(ji,jbdy1:jbdy2,jk,Krhs_a) & 
+                    &                             +   uu_b(ji,jbdy1:jbdy2,   Krhs_a) &
+                    &                             -    zub(ji,jbdy1:jbdy2)           &
+                    &                            ) * umask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -684,7 +725,7 @@
                iref = ji
                jref = jj
-               if(western_side) iref=MAX(2,ji)
-               if(eastern_side) iref=MIN(nlci-1,ji)
-               if(southern_side) jref=MAX(2,jj)
+               if(western_side)  iref=MAX(2     ,ji)
+               if(eastern_side)  iref=MIN(nlci-1,ji)
+               if(southern_side) jref=MAX(2     ,jj)
                if(northern_side) jref=MIN(nlcj-1,jj)
                N_in = 0
@@ -693,5 +734,5 @@
                   N_in = N_in + 1
                   tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
-                  h_in(N_in) = ptab(ji,jj,jk,n2)
+                  h_in(N_in)  = ptab(ji,jj,jk,n2)
                END DO
                N_out = 0
@@ -722,7 +763,7 @@
             ! 
             ! Remove CORNERS
-            IF((nbondj == -1).OR.(nbondj == 2)) jmin = 2 + nbghostcells
+            IF((nbondj == -1).OR.(nbondj == 2)) jmin = 2    + nbghostcells
             IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj - nbghostcells - 1
-            IF((nbondi == -1).OR.(nbondi == 2)) imin = 2 + nbghostcells
+            IF((nbondi == -1).OR.(nbondi == 2)) imin = 2    + nbghostcells
             IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci - nbghostcells - 1      
             !
@@ -737,5 +778,6 @@
                ibdy = nlci-nbghostcells
                DO jn = 1, jpts
-                  ts(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = 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_a) =  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
@@ -743,8 +785,12 @@
                            ts(ibdy,jj,jk,jn,Krhs_a) = ts(ibdy+1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)
                         ELSE
-                           ts(ibdy,jj,jk,jn,Krhs_a)=(z4*ts(ibdy+1,jj,jk,jn,Krhs_a)+z3*ts(ibdy-1,jj,jk,jn,Krhs_a))*tmask(ibdy,jj,jk)
+                           ts(ibdy,jj,jk,jn,Krhs_a) = (   z4 * ts(ibdy+1,jj,jk,jn,Krhs_a) &
+                    &                                  +  z3 * ts(ibdy-1,jj,jk,jn,Krhs_a) &
+                    &                                 ) *   tmask(ibdy  ,jj,jk)
                            IF( uu(ibdy-1,jj,jk,Kmm_a) > 0._wp ) THEN
-                              ts(ibdy,jj,jk,jn,Krhs_a)=( z6*ts(ibdy-1,jj,jk,jn,Krhs_a)+z5*ts(ibdy+1,jj,jk,jn,Krhs_a) & 
-                                                 + z7*ts(ibdy-2,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)
+                              ts(ibdy,jj,jk,jn,Krhs_a) = (  z6 * ts(ibdy-1,jj,jk,jn,Krhs_a) &
+                    &                                     + z5 * ts(ibdy+1,jj,jk,jn,Krhs_a) & 
+                    &                                     + z7 * ts(ibdy-2,jj,jk,jn,Krhs_a) &
+                    &                                    ) *  tmask(ibdy  ,jj,jk)
                            ENDIF
                         ENDIF
@@ -752,5 +798,6 @@
                   END DO
                   ! Restore ghost points:
-                  ts(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = 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_a) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) &
+                    &                                     *     tmask(ibdy+1,jmin:jmax,1:jpkm1)
                END DO
             ENDIF
@@ -766,5 +813,6 @@
                jbdy = nlcj-nbghostcells         
                DO jn = 1, jpts
-                  ts(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = 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_a) =  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
@@ -772,8 +820,12 @@
                            ts(ji,jbdy,jk,jn,Krhs_a) = ts(ji,jbdy+1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
                         ELSE
-                           ts(ji,jbdy,jk,jn,Krhs_a)=(z4*ts(ji,jbdy+1,jk,jn,Krhs_a)+z3*ts(ji,jbdy-1,jk,jn,Krhs_a))*tmask(ji,jbdy,jk)        
+                           ts(ji,jbdy,jk,jn,Krhs_a)=(  z4 * ts(ji,jbdy+1,jk,jn,Krhs_a)   &
+                    &                                + z3 * ts(ji,jbdy-1,jk,jn,Krhs_a)   &
+                    &                               ) *  tmask(ji,jbdy  ,jk)        
                            IF (vv(ji,jbdy-1,jk,Kmm_a) > 0._wp ) THEN
-                              ts(ji,jbdy,jk,jn,Krhs_a)=( z6*ts(ji,jbdy-1,jk,jn,Krhs_a)+z5*ts(ji,jbdy+1,jk,jn,Krhs_a)  &
-                                                 + z7*ts(ji,jbdy-2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
+                              ts(ji,jbdy,jk,jn,Krhs_a)=(  z6 * ts(ji,jbdy-1,jk,jn,Krhs_a) &
+                    &                                   + z5 * ts(ji,jbdy+1,jk,jn,Krhs_a) &
+                    &                                   + z7 * ts(ji,jbdy-2,jk,jn,Krhs_a) &
+                    &                                  ) *  tmask(ji,jbdy  ,jk)
                            ENDIF
                         ENDIF
@@ -781,5 +833,6 @@
                   END DO
                   ! Restore ghost points:
-                  ts(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = 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_a) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) &
+                    &                                     *     tmask(imin:imax,jbdy+1,1:jpkm1)
                END DO
             ENDIF
@@ -795,5 +848,6 @@
                ibdy = 1+nbghostcells       
                DO jn = 1, jpts
-                  ts(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = 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_a) =  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
@@ -801,8 +855,12 @@
                            ts(ibdy,jj,jk,jn,Krhs_a) = ts(ibdy-1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)
                         ELSE
-                           ts(ibdy,jj,jk,jn,Krhs_a)=(z4*ts(ibdy-1,jj,jk,jn,Krhs_a)+z3*ts(ibdy+1,jj,jk,jn,Krhs_a))*tmask(ibdy,jj,jk)        
+                           ts(ibdy,jj,jk,jn,Krhs_a) = (  z4 * ts(ibdy-1,jj,jk,jn,Krhs_a) &
+                    &                                  + z3 * ts(ibdy+1,jj,jk,jn,Krhs_a) &
+                    &                                 ) *  tmask(ibdy  ,jj,jk)        
                            IF( uu(ibdy,jj,jk,Kmm_a) < 0._wp ) THEN
-                              ts(ibdy,jj,jk,jn,Krhs_a)=( z6*ts(ibdy+1,jj,jk,jn,Krhs_a)+z5*ts(ibdy-1,jj,jk,jn,Krhs_a) &
-                                                 + z7*ts(ibdy+2,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)
+                              ts(ibdy,jj,jk,jn,Krhs_a) = (  z6 * ts(ibdy+1,jj,jk,jn,Krhs_a) &
+                    &                                     + z5 * ts(ibdy-1,jj,jk,jn,Krhs_a) &
+                    &                                     + z7 * ts(ibdy+2,jj,jk,jn,Krhs_a) &
+                    &                                    ) *  tmask(ibdy,jj,jk)
                            ENDIF
                         ENDIF
@@ -810,5 +868,6 @@
                   END DO
                   ! Restore ghost points:
-                  ts(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = 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_a) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) &
+                    &                                     *     tmask(ibdy-1,jmin:jmax,1:jpkm1)
                END DO
             ENDIF
@@ -824,14 +883,19 @@
                jbdy=1+nbghostcells        
                DO jn = 1, jpts
-                  ts(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = 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_a) =  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
-                           ts(ji,jbdy,jk,jn,Krhs_a)=ts(ji,jbdy-1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
+                           ts(ji,jbdy,jk,jn,Krhs_a) = ts(ji,jbdy-1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
                         ELSE
-                           ts(ji,jbdy,jk,jn,Krhs_a)=(z4*ts(ji,jbdy-1,jk,jn,Krhs_a)+z3*ts(ji,jbdy+1,jk,jn,Krhs_a))*tmask(ji,jbdy,jk)
+                           ts(ji,jbdy,jk,jn,Krhs_a) = (  z4 * ts(ji,jbdy-1,jk,jn,Krhs_a) &
+                    &                                  + z3 * ts(ji,jbdy+1,jk,jn,Krhs_a) &
+                    &                                 ) *  tmask(ji,jbdy  ,jk)
                            IF( vv(ji,jbdy,jk,Kmm_a) < 0._wp ) THEN
-                              ts(ji,jbdy,jk,jn,Krhs_a)=( z6*ts(ji,jbdy+1,jk,jn,Krhs_a)+z5*ts(ji,jbdy-1,jk,jn,Krhs_a) & 
-                                                 + z7*ts(ji,jbdy+2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
+                              ts(ji,jbdy,jk,jn,Krhs_a) = (  z6 * ts(ji,jbdy+1,jk,jn,Krhs_a) &
+                    &                                     + z5 * ts(ji,jbdy-1,jk,jn,Krhs_a) & 
+                    &                                     + z7 * ts(ji,jbdy+2,jk,jn,Krhs_a) &
+                    &                                    ) *  tmask(ji,jbdy  ,jk)
                            ENDIF
                         ENDIF
@@ -839,5 +903,6 @@
                   END DO
                   ! Restore ghost points:
-                  ts(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = 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_a) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) &
+                    &                                     *     tmask(imin:imax,jbdy-1,1:jpkm1)
                END DO
             ENDIF
@@ -900,5 +965,6 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)) 
+                  ptab(ji,jj,jk,1) = (  e2u(ji,jj)          *   e3u(ji,jj,jk,Kmm_a)  &
+                    &                 *  uu(ji,jj,jk,Kmm_a) * umask(ji,jj,jk)      ) 
 # if defined key_vertical
                   ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a))
@@ -923,5 +989,5 @@
                   IF (ptab(ji,jj,jk,2) == 0) EXIT
                   N_in = N_in + 1
-                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
+                  tabin(jk)  = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
                   h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
               ENDDO
@@ -945,12 +1011,13 @@
          
               IF (N_in * N_out > 0) THEN
-                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                 h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) )
 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
                  if (h_diff < -1.e4) then
-                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
+                    print *,'CHECK YOUR BATHY ...', h_diff, SUM( h_out(1:N_out) ), SUM( h_in(1:N_in) )
 !                    stop
                  endif
               ENDIF
-              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),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_a), &
+                    &                   h_out(1:N_out), N_in        , N_out                    )
             ENDDO
          ENDDO
@@ -992,5 +1059,6 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk))
+                  ptab(ji,jj,jk,1) = (  e1v(ji,jj)          *   e3v(ji,jj,jk,Kmm_a) &
+                    &                 *  vv(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk)      )
 # if defined key_vertical
                   ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
@@ -1015,5 +1083,5 @@
                   if (ptab(ji,jj,jk,2) == 0) EXIT
                   N_in = N_in + 1
-                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
+                  tabin(jk)  = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
                   h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
                END DO
@@ -1033,5 +1101,6 @@
                  CYCLE
                ENDIF
-               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),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_a), &
+                    &                    h_out(1:N_out), N_in        , N_out                    )
             END DO
          END DO
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 11427)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_sponge.F90	(revision 11428)
@@ -393,7 +393,7 @@
          
               IF (N_in * N_out > 0) THEN
-                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                 h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) )
                  if (h_diff < -1.e4) then
-                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
+                    print *,'CHECK YOUR BATHY ...', h_diff, SUM( h_out(1:N_out) ), SUM( h_in(1:N_in) )
                  endif
               ENDIF
@@ -403,7 +403,7 @@
          ENDDO
 
-         ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
+         ubdiff(i1:i2,j1:j2,:) = ( uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:  ) )*umask(i1:i2,j1:j2,:)
 #else
-         ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
+         ubdiff(i1:i2,j1:j2,:) = ( uu(i1:i2,j1:j2,:,Kbb_a) -       tabres(i1:i2,j1:j2,:,1) )*umask(i1:i2,j1:j2,:)
 #endif
          !
@@ -540,7 +540,7 @@
          
               IF (N_in * N_out > 0) THEN
-                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                 h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) )
                  if (h_diff < -1.e4) then
-                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
+                    print *,'CHECK YOUR BATHY ...', h_diff, SUM( h_out(1:N_out) ), SUM( h_in(1:N_in) )
                  endif
               ENDIF
@@ -549,7 +549,7 @@
          ENDDO
 
-         vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)  
+         vbdiff(i1:i2,j1:j2,:) = ( vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:  ) )*vmask(i1:i2,j1:j2,:)  
 # else
-         vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
+         vbdiff(i1:i2,j1:j2,:) = ( vv(i1:i2,j1:j2,:,Kbb_a) -       tabres(i1:i2,j1:j2,:,1) )*vmask(i1:i2,j1:j2,:)
 # endif
          !
@@ -570,6 +570,6 @@
                DO ji = i1,i2-1   ! vector opt.
                   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
+                  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
                END DO
             END DO
@@ -586,7 +586,8 @@
                IF( .NOT. tabspongedone_u(ji,jj) ) THEN
                   DO jk = 1, jpkm1
-                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a)                                                               &
-                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )  &
-                        & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
+                     uu(ji,jj,jk,Krhs_a) =          uu(ji  ,jj,jk,Krhs_a)                                   &
+                        &                 - (  rotdiff(ji  ,jj,jk) -  rotdiff(ji,jj-1,jk) )                 &
+                        &                   / (    e2u(ji  ,jj)    *      e3u(ji,jj  ,jk,Kmm_a) )           &
+                        &                 + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk) ) * r1_e1u(ji,jj)
                   END DO
                ENDIF
@@ -600,7 +601,8 @@
                IF( .NOT. tabspongedone_v(ji,jj) ) THEN
                   DO jk = 1, jpkm1
-                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a)                                                                  &
-                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   &
-                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)
+                     vv(ji,jj,jk,Krhs_a) =          vv(ji,jj  ,jk,Krhs_a)                                   &
+                        &                 + (  rotdiff(ji,jj  ,jk) -  rotdiff(ji-1,jj,jk) )                 &
+                        &                   / (    e1v(ji,jj  ) *         e3v(ji  ,jj,jk,Kmm_a) )           &
+                        &                 + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)
                   END DO
                ENDIF
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 11427)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_update.F90	(revision 11428)
@@ -315,6 +315,6 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) &
-                                           * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp
+                     tabres(ji,jj,jk,jn) = (    ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) )          &
+                       &                   * tmask(ji,jj,jk)       + (tmask(ji,jj,jk) - 1._wp)*999._wp
                   END DO
                END DO
@@ -325,10 +325,10 @@
                DO ji=i1,i2
                   tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) &
-                                           + (tmask(ji,jj,jk)-1)*999._wp
+                    &                      + (tmask(ji,jj,jk) - 1._wp)*999._wp
                END DO
             END DO
          END DO
       ELSE
-         tabres_child(:,:,:,:) = 0.
+         tabres_child(:,:,:,:) = 0._wp
          AGRIF_SpecialValue = 0._wp
          DO jj=j1,j2
@@ -356,5 +356,6 @@
                   ENDIF
                   DO jn=n1,n2-1
-                     CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out)
+                     CALL reconstructandremap( tabin(1:N_in,jn), h_in(1:N_in), tabres_child(ji,jj,1:N_out,jn),  &
+                       &                       h_out(1:N_out)  , N_in        , N_out                           )
                   ENDDO
                ENDIF
@@ -369,7 +370,8 @@
                      DO ji=i1,i2
                         IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
-                           ts(ji,jj,jk,jn,Kbb_a) = ts(ji,jj,jk,jn,Kbb_a) & 
-                                 & + atfp * ( tabres_child(ji,jj,jk,jn) &
-                                 &          - ts(ji,jj,jk,jn,Kmm_a) ) * tmask(ji,jj,jk)
+                           ts(ji,jj,jk,jn,Kbb_a) =                     ts(ji,jj,jk,jn,Kbb_a) & 
+                                 &                + atfp * ( tabres_child(ji,jj,jk,jn)       &
+                                 &                                   - ts(ji,jj,jk,jn,Kmm_a) &
+                                 &                         ) * tmask(ji,jj,jk)
                         ENDIF
                      ENDDO
@@ -414,5 +416,5 @@
 !> jc tmp
                      tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk)
-!                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a)
+!                    tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a)
 !< jc tmp
                   END DO
@@ -438,5 +440,5 @@
                            ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
                            ts(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) )  & 
-                                     &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
+                                     &            * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
                         ENDIF
                      END DO
@@ -496,7 +498,7 @@
                DO ji=i1,i2
                   tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a)  &
-                                       + (umask(ji,jj,jk)-1)*999._wp
+                    &                  + (umask(ji,jj,jk)-1._wp)*999._wp
                   tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)  &
-                                       + (umask(ji,jj,jk)-1)*999._wp
+                    &                  + (umask(ji,jj,jk)-1._wp)*999._wp
                END DO
             END DO
@@ -513,5 +515,5 @@
                   IF( tabres(ji,jj,jk,2) < -900) EXIT
                   N_in = N_in + 1
-                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
+                  tabin(jk)  = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
                   h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
                ENDDO
@@ -523,5 +525,5 @@
                ENDDO
                IF (N_in * N_out > 0) THEN
-                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                  h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) )
                   IF (h_diff < -1.e-4) THEN
 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 
@@ -540,6 +542,7 @@
                      ENDDO
                   ENDIF
-                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
-                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
+                  CALL reconstructandremap( tabin(1:N_in) , h_in(1:N_in), tabres_child(ji,jj,1:N_out),       &
+                    &                       h_out(1:N_out), N_in        , N_out                            )
+                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/( e2u(ji,jj)*h_out(N_out) )
                ENDIF
             ENDDO
@@ -550,6 +553,6 @@
                DO ji=i1,i2
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) & 
-                           & + atfp * ( tabres_child(ji,jj,jk) - uu(ji,jj,jk,Kmm_a) ) * umask(ji,jj,jk)
+                     uu(ji,jj,jk,Kbb_a) =                     uu(ji,jj,jk,Kbb_a)                                   & 
+                           &             + atfp * ( tabres_child(ji,jj,jk) - uu(ji,jj,jk,Kmm_a) ) * umask(ji,jj,jk)
                   ENDIF
                   !
@@ -579,5 +582,6 @@
          zrhoy = Agrif_Rhoy()
          DO jk = k1, k2
-            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) * uu(i1:i2,j1:j2,jk,Kmm_a)
+            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) &
+              &                                                 *  uu(i1:i2,j1:j2,jk,Kmm_a)
          END DO
       ELSE
@@ -591,6 +595,6 @@
                      zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a)
                      zunu = tabres(ji,jj,jk,1)
-                     uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno) ) &      
-                                    & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
+                     uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno ) )       &      
+                       &                  * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
                   ENDIF
                   !
@@ -682,6 +686,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) &
-                                       + (vmask(ji,jj,jk)-1)*999._wp
+                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)                   &
+                    &                  *  vmask(ji,jj,jk) *  vv(ji,jj,jk,Kmm_a)                   &
+                    &                  + (vmask(ji,jj,jk)-1)*999._wp
                   tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) &
                                        + (vmask(ji,jj,jk)-1)*999._wp
@@ -708,5 +713,5 @@
                ENDDO
                IF (N_in * N_out > 0) THEN
-                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                  h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) )
                   IF (h_diff < -1.e-4) then
 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 
@@ -714,7 +719,7 @@
                      excess = 0._wp
                      DO jk=N_in,1,-1
-                        thick = MIN(-1*h_diff, h_in(jk))
+                        thick = MIN( -1._wp*h_diff, h_in(jk) )
                         excess = excess + tabin(jk)*thick*e2u(ji,jj)
-                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
+                        tabin(jk) = tabin(jk)*(1._wp - thick/h_in(jk))
                         h_diff = h_diff + thick
                         IF ( h_diff == 0) THEN
@@ -725,5 +730,6 @@
                      ENDDO
                   ENDIF
-                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(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), tabres_child(ji,jj,1:N_out),     &
+                    &                       h_out(1:N_out), N_in        , N_out                          )
                   tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
                ENDIF
@@ -736,6 +742,7 @@
                   !
                   IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part
-                     vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) & 
-                           & + atfp * ( tabres_child(ji,jj,jk) - vv(ji,jj,jk,Kmm_a) ) * vmask(ji,jj,jk)
+                     vv(ji,jj,jk,Kbb_a) =                     vv(ji,jj,jk,Kbb_a)                  & 
+                           &             + atfp * ( tabres_child(ji,jj,jk) - vv(ji,jj,jk,Kmm_a) ) &
+                           &                      *        vmask(ji,jj,jk)
                   ENDIF
                   !
@@ -781,6 +788,6 @@
                      zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a)
                      zvnu = tabres(ji,jj,jk,1)
-                     vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) ) &      
-                                    & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
+                     vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) )        &      
+                       &                  * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
                   ENDIF
                   !
@@ -883,6 +890,6 @@
                ! Update barotropic velocities:
                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) - uu_b(ji,jj,Kmm_a) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
+                  IF ( .NOT.( lk_agrif_fstep .AND. (neuler==0) ) ) THEN                          ! Add asselin part
+                     zcorr = ( tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu_a(ji,jj) ) * r1_hu_b(ji,jj)
                      uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + atfp * zcorr * umask(ji,jj,1)
                   END IF
@@ -948,8 +955,8 @@
                !
                ! Update barotropic velocities:
-               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) - vv_b(ji,jj,Kmm_a) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
-                     vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + atfp * zcorr * vmask(ji,jj,1)
+               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) - vv_b(ji,jj,Kmm_a) * hv_a(ji,jj) ) * r1_hv_b(ji,jj)
+                     vv_b(ji,jj,Kbb_a) =       vv_b(ji,jj,Kbb_a) + atfp * zcorr  *   vmask(ji,jj,1)
                   END IF
                ENDIF              
@@ -1000,6 +1007,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ssh(ji,jj,Kbb_a) =   ssh(ji,jj,Kbb_a) &
-                        & + atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1)
+                  ssh(ji,jj,Kbb_a) =              ssh(ji,jj,Kbb_a)                &
+                    &               + atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) &
+                    &                        *  tmask(ji,jj,1)
                END DO
             END DO
@@ -1095,5 +1103,6 @@
                zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 
                ssh(i1  ,jj,Kmm_a) = ssh(i1  ,jj,Kmm_a) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i1  ,jj,Kbb_a) = ssh(i1  ,jj,Kbb_a) + atfp * zcor
+               IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) )  &
+                 &                  ssh(i1  ,jj,Kbb_a) = ssh(i1  ,jj,Kbb_a) + atfp * zcor
             END DO
          ENDIF
@@ -1102,5 +1111,6 @@
                zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
                ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + atfp * zcor
+               IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) )   &
+                 &                  ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + atfp * zcor
             END DO
          ENDIF
@@ -1181,14 +1191,16 @@
          IF (southern_side) THEN
             DO ji=i1,i2
-               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
+               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * ( vb2_b(ji,j1)-tabres(ji,j1) )
                ssh(ji,j1  ,Kmm_a) = ssh(ji,j1  ,Kmm_a) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j1  ,Kbb_a) = ssh(ji,j1,Kbb_a) + atfp * zcor
+               IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) )   &
+                 &                 ssh(ji,j1  ,Kbb_a) = ssh(ji,j1,Kbb_a) + atfp * zcor
             END DO
          ENDIF
          IF (northern_side) THEN               
             DO ji=i1,i2
-               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
+               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * ( vb2_b(ji,j2)-tabres(ji,j2) )
                ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + atfp * zcor
+               IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) )   &
+                 &                  ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + atfp * zcor
             END DO
          ENDIF
@@ -1223,7 +1235,7 @@
             END DO
          END DO
-         tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()
-         tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
-         tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
+         tabres(:,:,:,1) = tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()
+         tabres(:,:,:,2) = tabres(:,:,:,2)*Agrif_Rhox()
+         tabres(:,:,:,3) = tabres(:,:,:,3)*Agrif_Rhoy()
       ELSE
          DO jk=k1,k2
@@ -1234,10 +1246,10 @@
                      print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk)
                      print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk)
-                     ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
+                     ztemp = SQRT( tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)) )
                      print *,'CORR = ',ztemp-1.
                      print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, &
                            tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp
-                     e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
-                     e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
+                     e1t(ji,jj) = tabres(ji,jj,jk,2) * ztemp
+                     e2t(ji,jj) = tabres(ji,jj,jk,3) * ztemp
                   END IF
                END DO
@@ -1319,6 +1331,8 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm_a) &
-                                     & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
+                  ptab(ji,jj,jk) =   e3t_0(ji,jj,jk)                                    &
+                    &                * ( 1._wp + ssh(ji,jj,Kmm_a)                       &
+                    &                          * ssmask(ji,jj)                          &
+                    &                          / ( ht_0(ji,jj)-1._wp + ssmask(ji,jj) ) )
                END DO
             END DO
@@ -1339,6 +1353,6 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     e3t(ji,jj,jk,Kbb_a) =  e3t(ji,jj,jk,Kbb_a) &
-                           & + atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) )
+                     e3t(ji,jj,jk,Kbb_a) =             e3t(ji,jj,jk,Kbb_a)                  &
+                       &                   + atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) )
                   END DO
                END DO
@@ -1353,11 +1367,11 @@
                   DO ji = i1,i2            
                      zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-                     e3w(ji,jj,jk,Kbb_a)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
-                     &                                        ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) )  &
-                     &                                  +            0.5_wp * tmask(ji,jj,jk)   *        &
-                     &                                        ( e3t(ji,jj,jk  ,Kbb_a) - e3t_0(ji,jj,jk  ) )
-                     gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk-1,Kbb_a)
-                     gdept(ji,jj,jk,Kbb_a) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb_a) + 0.5 * e3w(ji,jj,jk,Kbb_a))  &
-                         &               + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) +       e3w(ji,jj,jk,Kbb_a)) 
+                     e3w(ji,jj,jk,Kbb_a)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) )          & 
+                       &                                    * ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) )  &
+                       &                                    + 0.5_wp * tmask(ji,jj,jk)                       &
+                       &                                    * ( e3t(ji,jj,jk  ,Kbb_a) - e3t_0(ji,jj,jk  ) )
+                     gdepw(ji,jj,jk,Kbb_a) =                       gdepw(ji,jj,jk-1,Kbb_a) +       e3t(ji,jj,jk-1,Kbb_a)
+                     gdept(ji,jj,jk,Kbb_a) =            zcoef  * ( gdepw(ji,jj,jk  ,Kbb_a) + 0.5 * e3w(ji,jj,jk  ,Kbb_a) )  &
+                       &                     + (1._wp - zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) +       e3w(ji,jj,jk  ,Kbb_a) ) 
                   END DO
                END DO
@@ -1379,8 +1393,8 @@
          !
          ! Update vertical scale factor at W-points and depths:
-         e3w (i1:i2,j1:j2,1,Kmm_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm_a) - e3t_0(i1:i2,j1:j2,1)
+           e3w(i1:i2,j1:j2,1,Kmm_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm_a) - e3t_0(i1:i2,j1:j2,1)
          gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a)
          gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp
-         gde3w(i1:i2,j1:j2,1) = gdept(i1:i2,j1:j2,1,Kmm_a) - (ht_n(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh
+         gde3w(i1:i2,j1:j2,1      ) = gdept(i1:i2,j1:j2,1,Kmm_a) - ( 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 +1402,12 @@
                DO ji = i1,i2            
                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               e3w(ji,jj,jk,Kmm_a)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t(ji,jj,jk-1,Kmm_a) - e3t_0(ji,jj,jk-1) )   &
-               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t(ji,jj,jk  ,Kmm_a) - e3t_0(ji,jj,jk  ) )
-               gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a)
-               gdept(ji,jj,jk,Kmm_a) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm_a) + 0.5 * e3w(ji,jj,jk,Kmm_a))  &
-                   &               + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm_a) +       e3w(ji,jj,jk,Kmm_a)) 
-               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm_a) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
+                e3w(ji,jj,jk,Kmm_a)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) )          &
+                 &                                     * ( e3t(ji,jj,jk-1,Kmm_a) - e3t_0(ji,jj,jk-1) )  &
+                 &                                     + 0.5_wp * tmask(ji,jj,jk)                       &
+                 &                                     * ( e3t(ji,jj,jk  ,Kmm_a) - e3t_0(ji,jj,jk  ) )
+               gdepw(ji,jj,jk,Kmm_a) =                         gdepw(ji,jj,jk-1,Kmm_a) +       e3t(ji,jj,jk-1,Kmm_a)
+               gdept(ji,jj,jk,Kmm_a) =             zcoef   * ( gdepw(ji,jj,jk  ,Kmm_a) + 0.5 * e3w(ji,jj,jk  ,Kmm_a) )  &
+                 &                     + ( 1._wp - zcoef ) * ( gdept(ji,jj,jk-1,Kmm_a) +       e3w(ji,jj,jk  ,Kmm_a) ) 
+               gde3w(ji,jj,jk      ) =                         gdept(ji,jj,jk  ,Kmm_a) - ( ht_n(ji,jj)-ht_0(ji,jj) )    ! Last term in the rhs is ssh
                END DO
             END DO
@@ -1399,6 +1415,6 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            e3t (i1:i2,j1:j2,1:jpk,Kbb_a)  = e3t (i1:i2,j1:j2,1:jpk,Kmm_a)
-            e3w (i1:i2,j1:j2,1:jpk,Kbb_a)  = e3w (i1:i2,j1:j2,1:jpk,Kmm_a)
+              e3t(i1:i2,j1:j2,1:jpk,Kbb_a) =   e3t(i1:i2,j1:j2,1:jpk,Kmm_a)
+              e3w(i1:i2,j1:j2,1:jpk,Kbb_a) =   e3w(i1:i2,j1:j2,1:jpk,Kmm_a)
             gdepw(i1:i2,j1:j2,1:jpk,Kbb_a) = gdepw(i1:i2,j1:j2,1:jpk,Kmm_a)
             gdept(i1:i2,j1:j2,1:jpk,Kbb_a) = gdept(i1:i2,j1:j2,1:jpk,Kmm_a)
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 11427)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_interp.F90	(revision 11428)
@@ -151,5 +151,6 @@
                ibdy = nlci-nbghostcells
                DO jn = 1, jptra
-                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = 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_a) =  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
@@ -157,14 +158,19 @@
                            tr(ibdy,jj,jk,jn,Krhs_a) = tr(ibdy+1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)
                         ELSE
-                           tr(ibdy,jj,jk,jn,Krhs_a)=(z4*tr(ibdy+1,jj,jk,jn,Krhs_a)+z3*tr(ibdy-1,jj,jk,jn,Krhs_a))*tmask(ibdy,jj,jk)
+                           tr(ibdy,jj,jk,jn,Krhs_a) = (  z4 * tr(ibdy+1,jj,jk,jn,Krhs_a)     &
+                             &                         + z3 * tr(ibdy-1,jj,jk,jn,Krhs_a)     &
+                             &                        ) *  tmask(ibdy  ,jj,jk)
                            IF( uu(ibdy-1,jj,jk,Kmm_a) > 0._wp ) THEN
-                              tr(ibdy,jj,jk,jn,Krhs_a)=( z6*tr(ibdy-1,jj,jk,jn,Krhs_a)+z5*tr(ibdy+1,jj,jk,jn,Krhs_a) & 
-                                                 + z7*tr(ibdy-2,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1)
+                              tr(ibdy,jj,jk,jn,Krhs_a) = (  z6 * tr(ibdy-1,jj,jk,jn,Krhs_a)  &
+                                &                         + z5 * tr(ibdy+1,jj,jk,jn,Krhs_a)  & 
+                                &                         + z7 * tr(ibdy-2,jj,jk,jn,Krhs_a)  &
+                                &                         ) * tmask(ibdy  ,jj,jk)
+                           ENDIF
+                        ENDIF
+                     END DO
+                  END DO
+                  ! Restore ghost points:
+                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) &
+                    &                                     *     tmask(ibdy+1,jmin:jmax,1:jpkm1)
                END DO
             ENDIF
@@ -180,5 +186,6 @@
                jbdy = nlcj-nbghostcells         
                DO jn = 1, jptra
-                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = 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_a) =  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
@@ -186,14 +193,17 @@
                            tr(ji,jbdy,jk,jn,Krhs_a) = tr(ji,jbdy+1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
                         ELSE
-                           tr(ji,jbdy,jk,jn,Krhs_a)=(z4*tr(ji,jbdy+1,jk,jn,Krhs_a)+z3*tr(ji,jbdy-1,jk,jn,Krhs_a))*tmask(ji,jbdy,jk)        
+                           tr(ji,jbdy,jk,jn,Krhs_a) = ( z4 * tr(ji,jbdy+1,jk,jn,Krhs_a) 
+                             &                        + z3 * tr(ji,jbdy-1,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)        
                            IF (vv(ji,jbdy-1,jk,Kmm_a) > 0._wp ) THEN
-                              tr(ji,jbdy,jk,jn,Krhs_a)=( z6*tr(ji,jbdy-1,jk,jn,Krhs_a)+z5*tr(ji,jbdy+1,jk,jn,Krhs_a)  &
-                                                 + z7*tr(ji,jbdy-2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1)
+                              tr(ji,jbdy,jk,jn,Krhs_a) = (  z6 * tr(ji,jbdy-1,jk,jn,Krhs_a)                       &
+                                &                         + z5 * tr(ji,jbdy+1,jk,jn,Krhs_a)                       &
+                                &                         + z7 * tr(ji,jbdy-2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
+                           ENDIF
+                        ENDIF
+                     END DO
+                  END DO
+                  ! Restore ghost points:
+                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn)  &
+                    &                                     *     tmask(imin:imax,jbdy+1,1:jpkm1)
                END DO
             ENDIF
@@ -209,5 +219,6 @@
                ibdy = 1+nbghostcells       
                DO jn = 1, jptra
-                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = 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_a) =  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
@@ -215,14 +226,17 @@
                            tr(ibdy,jj,jk,jn,Krhs_a) = tr(ibdy-1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)
                         ELSE
-                           tr(ibdy,jj,jk,jn,Krhs_a)=(z4*tr(ibdy-1,jj,jk,jn,Krhs_a)+z3*tr(ibdy+1,jj,jk,jn,Krhs_a))*tmask(ibdy,jj,jk)        
+                           tr(ibdy,jj,jk,jn,Krhs_a) = (  z4 * tr(ibdy-1,jj,jk,jn,Krhs_a)                          &
+                             &                         + z3 * tr(ibdy+1,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)        
                            IF( uu(ibdy,jj,jk,Kmm_a) < 0._wp ) THEN
-                              tr(ibdy,jj,jk,jn,Krhs_a)=( z6*tr(ibdy+1,jj,jk,jn,Krhs_a)+z5*tr(ibdy-1,jj,jk,jn,Krhs_a) &
-                                                 + z7*tr(ibdy+2,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1)
+                              tr(ibdy,jj,jk,jn,Krhs_a) = (  z6 * tr(ibdy+1,jj,jk,jn,Krhs_a)                       &
+                                &                         + z5 * tr(ibdy-1,jj,jk,jn,Krhs_a)                       &
+                                &                         + z7 * tr(ibdy+2,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)
+                           ENDIF
+                        ENDIF
+                     END DO
+                  END DO
+                  ! Restore ghost points:
+                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn)                &
+                    &                                     *     tmask(ibdy-1,jmin:jmax,1:jpkm1)
                END DO
             ENDIF
@@ -238,20 +252,24 @@
                jbdy=1+nbghostcells        
                DO jn = 1, jptra
-                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = 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_a) =  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
-                           tr(ji,jbdy,jk,jn,Krhs_a)=tr(ji,jbdy-1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
-                        ELSE
-                           tr(ji,jbdy,jk,jn,Krhs_a)=(z4*tr(ji,jbdy-1,jk,jn,Krhs_a)+z3*tr(ji,jbdy+1,jk,jn,Krhs_a))*tmask(ji,jbdy,jk)
+                           tr(ji,jbdy,jk,jn,Krhs_a) = tr(ji,jbdy-1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
+                        ELSE
+                           tr(ji,jbdy,jk,jn,Krhs_a) = (  z4 * tr(ji,jbdy-1,jk,jn,Krhs_a)                          &
+                             &                         + z3 * tr(ji,jbdy+1,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
                            IF( vv(ji,jbdy,jk,Kmm_a) < 0._wp ) THEN
-                              tr(ji,jbdy,jk,jn,Krhs_a)=( z6*tr(ji,jbdy+1,jk,jn,Krhs_a)+z5*tr(ji,jbdy-1,jk,jn,Krhs_a) & 
-                                                 + z7*tr(ji,jbdy+2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1)
+                              tr(ji,jbdy,jk,jn,Krhs_a) = (  z6 * tr(ji,jbdy+1,jk,jn,Krhs_a)                       &
+                                &                         + z5 * tr(ji,jbdy-1,jk,jn,Krhs_a)                       & 
+                                &                         + z7 * tr(ji,jbdy+2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
+                           ENDIF
+                        ENDIF
+                     END DO
+                  END DO
+                  ! Restore ghost points:
+                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = 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 11427)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_sponge.F90	(revision 11428)
@@ -108,14 +108,14 @@
                   N_in = N_in + 1
                   tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
-                  h_in(N_in) = tabres(ji,jj,jk,n2)
+                  h_in(N_in)  = tabres(ji,jj,jk,n2)
                END DO
                N_out = 0
                DO jk=1,jpk ! jpk of child grid
                   IF (tmask(ji,jj,jk) == 0) EXIT 
-                  N_out = N_out + 1
+                  N_out     = N_out + 1
                   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
-                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                  h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) )
                   tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for?????
                   DO jn=1,jptra
@@ -133,5 +133,5 @@
                   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) = tr(ji,jj,jk,1:jptra,Kbb) - 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
@@ -153,7 +153,7 @@
                   DO ji = i1+1,i2-1
                      IF( .NOT. tabspongedone_trn(ji,jj) ) THEN 
-                        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(ji,jj,jk,Kmm)
+                        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(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 11427)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_update.F90	(revision 11428)
@@ -84,6 +84,6 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     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
+                     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
                END DO
@@ -93,6 +93,6 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) &
-                                           + (tmask(ji,jj,jk)-1)*999._wp
+                  tabres(ji,jj,jk,n2) =   tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) &
+                    &                  + (tmask(ji,jj,jk)-1)*999._wp
                END DO
             END DO
@@ -138,7 +138,8 @@
                      DO ji=i1,i2
                         IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
-                           tr(ji,jj,jk,jn,Kbb) = ts(ji,jj,jk,jn,Kbb) & 
-                                 & + atfp * ( tabres_child(ji,jj,jk,jn) &
-                                 &          - tr(ji,jj,jk,jn,Kmm) ) * tmask(ji,jj,jk)
+                           tr(ji,jj,jk,jn,Kbb) = ts(ji,jj,jk,jn,Kbb)                     & 
+                             &                  + atfp * ( tabres_child(ji,jj,jk,jn)     &
+                             &                                     - tr(ji,jj,jk,jn,Kmm) &
+                             &                           ) *      tmask(ji,jj,jk)
                         ENDIF
                      ENDDO
@@ -183,6 +184,6 @@
                   DO ji=i1,i2
 !> jc tmp
-                     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)
+                     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
@@ -194,5 +195,5 @@
          DO jn = n1,n2
             tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
-                                         & * tmask(i1:i2,j1:j2,k1:k2)
+              &                            * tmask(i1:i2,j1:j2,k1:k2)
          ENDDO
 !< jc tmp
@@ -204,9 +205,9 @@
                      DO ji=i1,i2
                         IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
-                           ztb  = tr(ji,jj,jk,jn,Kbb) * e3t(ji,jj,jk,Kbb) ! 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 = 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)
+                           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
