Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_oce_interp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_oce_interp.F90	(revision 11099)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_oce_interp.F90	(revision 11463)
@@ -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
