Index: NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcpl.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo_TEST_MERGE/src/OCE/ISF/isfcpl.F90	(revision 11970)
+++ NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcpl.F90	(revision 12068)
@@ -46,7 +46,9 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE isfcpl_init()
+   SUBROUTINE isfcpl_init(Kbb, Kmm, Kaa)
       !!---------------------------------------------------------------------
       !!
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa      ! ocean time level indices
       !!---------------------------------------------------------------------
       INTEGER :: id
@@ -74,29 +76,29 @@
       ELSE
          ! extrapolation ssh
-         CALL isfcpl_ssh()
+         CALL isfcpl_ssh(Kbb, Kmm, Kaa)
          !
          ! extrapolation tracer properties
-         CALL isfcpl_tra()
+         CALL isfcpl_tra(Kmm)
          !
          ! correction of the horizontal divergence and associated temp. and salt content flux
          ! Need to : - include in the cpl cons the risfcpl_vol/tsc contribution
          !           - decide how to manage thickness level change in conservation
-         CALL isfcpl_vol()
+         CALL isfcpl_vol(Kmm)
          !
          ! apply the 'conservation' method
-         IF ( ln_isfcpl_cons ) CALL isfcpl_cons()
+         IF ( ln_isfcpl_cons ) CALL isfcpl_cons(Kmm)
          !
       END IF
       !
       ! mask velocity properly (mask used in restart not compatible with new mask)
-      un(:,:,:) = un(:,:,:) * umask(:,:,:)
-      vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)
+      uu(:,:,:,Kmm) = uu(:,:,:,Kmm) * umask(:,:,:)
+      vv(:,:,:,Kmm) = vv(:,:,:,Kmm) * vmask(:,:,:)
       !
       ! all before fields set to now values
-      tsb  (:,:,:,:) = tsn  (:,:,:,:)
-      ub   (:,:,:)   = un   (:,:,:)
-      vb   (:,:,:)   = vn   (:,:,:)
-      sshb (:,:)     = sshn (:,:)
-      e3t_b(:,:,:)   = e3t_n(:,:,:)
+      ts  (:,:,:,:,Kbb) = ts  (:,:,:,:,Kmm)
+      uu   (:,:,:,Kbb)   = uu   (:,:,:,Kmm)
+      vv   (:,:,:,Kbb)   = vv   (:,:,:,Kmm)
+      ssh (:,:,Kbb)     = ssh (:,:,Kmm)
+      e3t(:,:,:,Kbb)   = e3t(:,:,:,Kmm)
  
       ! prepare writing restart
@@ -111,5 +113,5 @@
    END SUBROUTINE isfcpl_init
    ! 
-   SUBROUTINE isfcpl_rst_write(kt)
+   SUBROUTINE isfcpl_rst_write(kt, Kmm)
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE iscpl_rst_write  ***
@@ -119,4 +121,5 @@
       !!-------------------------- IN  --------------------------------------
       INTEGER, INTENT(in) :: kt
+      INTEGER, INTENT(in) :: Kmm    ! ocean time level index
       !!----------------------------------------------------------------------
       !
@@ -124,13 +127,13 @@
       CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask , ldxios = lwxios )
       CALL iom_rstput( kt, nitrst, numrow, 'ssmask' , ssmask, ldxios = lwxios )
-      CALL iom_rstput( kt, nitrst, numrow, 'e3t_n'  , e3t_n , ldxios = lwxios )
-      CALL iom_rstput( kt, nitrst, numrow, 'e3u_n'  , e3u_n , ldxios = lwxios )
-      CALL iom_rstput( kt, nitrst, numrow, 'e3v_n'  , e3v_n , ldxios = lwxios )
-      CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw_n , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'e3t_n'  , e3t(:,:,:,Kmm) , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'e3u_n'  , e3u(:,:,:,Kmm) , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'e3v_n'  , e3v(:,:,:,Kmm) , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw(:,:,:,Kmm) , ldxios = lwxios )
       IF( lwxios ) CALL iom_swap( cxios_context )
       !
    END SUBROUTINE isfcpl_rst_write
 
-   SUBROUTINE isfcpl_ssh()
+   SUBROUTINE isfcpl_ssh(Kbb, Kmm, Kaa)
       !!---------------------------------------------------------------------- 
       !!                   ***  ROUTINE iscpl_ssh  ***
@@ -142,4 +145,6 @@
       !!----------------------------------------------------------------------
       !!
+      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa    ! ocean time level indices
+      !!----------------------------------------------------------------------
       INTEGER :: ji, jj, jd, jk      !! loop index
       INTEGER :: jip1, jim1, jjp1, jjm1
@@ -154,5 +159,5 @@
       ! rude average of the closest neigbourgs (e1e2t not taking into account)
       !
-      zssh(:,:)     = sshn(:,:)
+      zssh(:,:)     = ssh(:,:,Kmm)
       zssmask0(:,:) = zssmask_b(:,:)
       !
@@ -168,5 +173,5 @@
                !
                IF (zdssmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp) THEN
-                  sshn(ji,jj)=( zssh(jip1,jj)*zssmask0(jip1,jj)     &
+                  ssh(ji,jj,Kmm)=( zssh(jip1,jj)*zssmask0(jip1,jj)     &
                   &           + zssh(jim1,jj)*zssmask0(jim1,jj)     &
                   &           + zssh(ji,jjp1)*zssmask0(ji,jjp1)     &
@@ -177,5 +182,5 @@
          END DO
          !
-         zssh(:,:) = sshn(:,:)
+         zssh(:,:) = ssh(:,:,Kmm)
          zssmask0(:,:) = zssmask_b(:,:)
          !
@@ -184,25 +189,25 @@
       END DO
       !
-      ! update sshn
-      sshn(:,:) = zssh(:,:) * ssmask(:,:)
-      !
-      sshb(:,:) = sshn(:,:)
-      !
-      IF ( ln_isfdebug ) CALL debug('isfcpl_ssh: sshn',sshn(:,:))
+      ! update ssh(:,:,Kmm)
+      ssh(:,:,Kmm) = zssh(:,:) * ssmask(:,:)
+      !
+      ssh(:,:,Kbb) = ssh(:,:,Kmm)
+      !
+      IF ( ln_isfdebug ) CALL debug('isfcpl_ssh: sshn',ssh(:,:,Kmm))
       !
       ! recompute the vertical scale factor, depth and water thickness
-      IF(lwp) write(numout,*) 'isfcpl_ssh : recompute scale factor from sshn (new wet cell)'
+      IF(lwp) write(numout,*) 'isfcpl_ssh : recompute scale factor from ssh (new wet cell,Kmm)'
       IF(lwp) write(numout,*) '~~~~~~~~~~~'
       DO jk = 1, jpk
-         e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
+         e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
              &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
              &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
       END DO
-      e3t_b(:,:,:) = e3t_n(:,:,:)
-      CALL dom_vvl_zgr()
+      e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
+      CALL dom_vvl_zgr(Kbb, Kmm, Kaa)
       !
    END SUBROUTINE isfcpl_ssh
 
-   SUBROUTINE isfcpl_tra()
+   SUBROUTINE isfcpl_tra(Kmm)
       !!---------------------------------------------------------------------- 
       !!                   ***  ROUTINE iscpl_tra  ***
@@ -213,4 +218,6 @@
       !!
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kmm    ! ocean time level index
+      !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b
       !REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
@@ -237,5 +244,5 @@
       !bugged : to be corrected (PM)
       ! back up original t/s/mask
-      !tsb (:,:,:,:) = tsn(:,:,:,:)
+      !tsb (:,:,:,:) = ts(:,:,:,:,Kmm)
       !
      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
@@ -248,7 +255,7 @@
 !
 !                     !compute weight
-!                     zdzp1 = MAX(0._wp,pdepw_b(ji,jj,jk+1) - gdepw_n(ji,jj,jk+1))
-!                     zdzm1 = MAX(0._wp,gdepw_n(ji,jj,jk  ) - pdepw_b(ji,jj,jk  ))
-!                     zdz   = e3t_n(ji,jj,jk) - zdzp1 - zdzm1 ! if isf : e3t = gdepw_n(ji,jj,jk+1)- gdepw_n(ji,jj,jk)
+!                     zdzp1 = MAX(0._wp,pdepw_b(ji,jj,jk+1) - gdepw(ji,jj,jk+1,Kmm))
+!                     zdzm1 = MAX(0._wp,gdepw(ji,jj,jk  ,Kmm) - pdepw_b(ji,jj,jk  ))
+!                     zdz   = e3t(ji,jj,jk,Kmm) - zdzp1 - zdzm1 ! if isf : e3t = gdepw(ji,jj,jk+1,Kmm)- gdepw(ji,jj,jk,Kmm)
 !
 !                     IF (zdz .LT. 0._wp) THEN
@@ -256,11 +263,11 @@
 !                     END IF
 !
-!                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem) &
-!                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem) &
-!                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem) )/e3t_n(ji,jj,jk)
+!                     ts(ji,jj,jk,jp_tem,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_tem,Kbb) &
+!                        &                   + zdz  *ts(ji,jj,jk  ,jp_tem,Kbb) &
+!                        &                   + zdzm1*ts(ji,jj,jk-1,jp_tem,Kbb) )/e3t(ji,jj,jk,Kmm)
 !
-!                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal) &
-!                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal) &
-!                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal) )/e3t_n(ji,jj,jk)
+!                     ts(ji,jj,jk,jp_sal,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_sal,Kbb) &
+!                        &                   + zdz  *ts(ji,jj,jk  ,jp_sal,Kbb) &
+!                        &                   + zdzm1*ts(ji,jj,jk-1,jp_sal,Kbb) )/e3t(ji,jj,jk,Kmm)
 !
 !                  END IF
@@ -270,5 +277,5 @@
 !      END IF
 
-      zts0(:,:,:,:)  = tsn(:,:,:,:)
+      zts0(:,:,:,:)  = ts(:,:,:,:,Kmm)
       ztmask0(:,:,:) = ztmask_b(:,:,:)
       ztmask1(:,:,:) = ztmask_b(:,:,:)
@@ -294,9 +301,9 @@
                      !
                      ! horizontal basic extrapolation
-                     tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1) * ztmask0(jip1,jj  ,jk) &
+                     ts(ji,jj,jk,1,Kmm)=( zts0(jip1,jj  ,jk,1) * ztmask0(jip1,jj  ,jk) &
                      &               + zts0(jim1,jj  ,jk,1) * ztmask0(jim1,jj  ,jk) &
                      &               + zts0(ji  ,jjp1,jk,1) * ztmask0(ji  ,jjp1,jk) &
                      &               + zts0(ji  ,jjm1,jk,1) * ztmask0(ji  ,jjm1,jk) ) / zsummsk
-                     tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2) * ztmask0(jip1,jj  ,jk) &
+                     ts(ji,jj,jk,2,Kmm)=( zts0(jip1,jj  ,jk,2) * ztmask0(jip1,jj  ,jk) &
                      &               + zts0(jim1,jj  ,jk,2) * ztmask0(jim1,jj  ,jk) &
                      &               + zts0(ji  ,jjp1,jk,2) * ztmask0(ji  ,jjp1,jk) &
@@ -316,7 +323,7 @@
                      zsummsk = ztmask0(ji,jj,jkm1) + ztmask0(ji,jj,jkp1)
                      IF (zdmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp ) THEN
-                        tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
+                        ts(ji,jj,jk,1,Kmm)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
                         &               + zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / zsummsk
-                        tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
+                        ts(ji,jj,jk,2,Kmm)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
                         &               + zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / zsummsk
                         !
@@ -330,5 +337,5 @@
          !
          ! update temperature and salinity and mask
-         zts0(:,:,:,:)  = tsn(:,:,:,:)
+         zts0(:,:,:,:)  = ts(:,:,:,:,Kmm)
          ztmask0(:,:,:) = ztmask1(:,:,:)
          !
@@ -337,7 +344,7 @@
       END DO  ! nn_drown
       !
-      ! mask new tsn field
-      tsn(:,:,:,jp_tem) = zts0(:,:,:,jp_tem) * tmask(:,:,:)
-      tsn(:,:,:,jp_sal) = zts0(:,:,:,jp_sal) * tmask(:,:,:)
+      ! mask new ts(:,:,:,:,Kmm) field
+      ts(:,:,:,jp_tem,Kmm) = zts0(:,:,:,jp_tem) * tmask(:,:,:)
+      ts(:,:,:,jp_sal,Kmm) = zts0(:,:,:,jp_sal) * tmask(:,:,:)
       !
       ! sanity check
@@ -347,5 +354,5 @@
          DO jj = 1,jpj
             DO ji = 1,jpi
-               IF (tmask(ji,jj,jk) == 1._wp .AND. tsn(ji,jj,jk,2) == 0._wp)              &
+               IF (tmask(ji,jj,jk) == 1._wp .AND. ts(ji,jj,jk,2,Kmm) == 0._wp)              &
                   &   CALL ctl_stop('STOP', 'failing to fill all new weet cell,     &
                   &                          try increase nn_drown or activate XXXX &
@@ -357,5 +364,5 @@
    END SUBROUTINE isfcpl_tra
 
-   SUBROUTINE isfcpl_vol()
+   SUBROUTINE isfcpl_vol(Kmm)
       !!---------------------------------------------------------------------- 
       !!                   ***  ROUTINE iscpl_vol  ***
@@ -370,4 +377,6 @@
       !!----------------------------------------------------------------------
       !!
+      INTEGER, INTENT(in) :: Kmm    ! ocean time level index
+      !!----------------------------------------------------------------------
       INTEGER :: ji, jj, jk 
       INTEGER :: ikb, ikt
@@ -388,6 +397,6 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               zqvolb(ji,jj,jk) =  (   e2u(ji,jj) * ze3u_b(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  ) * ze3u_b(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
-                  &                  + e1v(ji,jj) * ze3v_b(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1) * ze3v_b(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  ) &
+               zqvolb(ji,jj,jk) =  (   e2u(ji,jj) * ze3u_b(ji,jj,jk) * uu(ji,jj,jk,Kmm) - e2u(ji-1,jj  ) * ze3u_b(ji-1,jj  ,jk) * uu(ji-1,jj  ,jk,Kmm)    &
+                  &                  + e1v(ji,jj) * ze3v_b(ji,jj,jk) * vv(ji,jj,jk,Kmm) - e1v(ji  ,jj-1) * ze3v_b(ji  ,jj-1,jk) * vv(ji  ,jj-1,jk,Kmm)  ) &
                   &                * ztmask_b(ji,jj,jk)
             END DO
@@ -397,11 +406,11 @@
          ! properly mask velocity 
          ! (velocity are still mask with old mask at this stage)
-         un(:,:,jk) = un(:,:,jk) * umask(:,:,jk)
-         vn(:,:,jk) = vn(:,:,jk) * vmask(:,:,jk)
+         uu(:,:,jk,Kmm) = uu(:,:,jk,Kmm) * umask(:,:,jk)
+         vv(:,:,jk,Kmm) = vv(:,:,jk,Kmm) * vmask(:,:,jk)
          ! compute volume flux divergence after coupling
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               zqvoln(ji,jj,jk) = (   e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  ) * e3u_n(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
-                  &                 + e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1) * e3v_n(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  ) &
+               zqvoln(ji,jj,jk) = (   e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) - e2u(ji-1,jj  ) * e3u(ji-1,jj  ,jk,Kmm) * uu(ji-1,jj  ,jk,Kmm)    &
+                  &                 + e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) - e1v(ji  ,jj-1) * e3v(ji  ,jj-1,jk,Kmm) * vv(ji  ,jj-1,jk,Kmm)  ) &
                   &               * tmask(ji,jj,jk)
             END DO
@@ -428,15 +437,15 @@
       CALL lbc_lnk( 'iscpl', risfcpl_vol, 'T', 1. )
       !
-      ! 3.0: set total correction (div, tra, ssh)
+      ! 3.0: set total correction (div, tr(:,:,:,:,Krhs), ssh)
       !
       ! 3.1: mask volume flux divergence correction
       risfcpl_vol(:,:,:) = risfcpl_vol(:,:,:) * tmask(:,:,:)
       !
-      ! 3.2: get 3d tra increment to apply at the first time step
-      ! temperature and salt content flux computed using local tsn 
+      ! 3.2: get 3d tr(:,:,:,:,Krhs) increment to apply at the first time step
+      ! temperature and salt content flux computed using local ts(:,:,:,:,Kmm) 
       ! (very simple advection scheme)
       ! (>0 out)
-      risfcpl_tsc(:,:,:,jp_tem) = -risfcpl_vol(:,:,:) * tsn(:,:,:,jp_tem)
-      risfcpl_tsc(:,:,:,jp_sal) = -risfcpl_vol(:,:,:) * tsn(:,:,:,jp_sal)
+      risfcpl_tsc(:,:,:,jp_tem) = -risfcpl_vol(:,:,:) * ts(:,:,:,jp_tem,Kmm)
+      risfcpl_tsc(:,:,:,jp_sal) = -risfcpl_vol(:,:,:) * ts(:,:,:,jp_sal,Kmm)
       !
       ! 3.3: ssh correction (for dynspg_ts)
@@ -448,5 +457,5 @@
    END SUBROUTINE isfcpl_vol
 
-   SUBROUTINE isfcpl_cons()
+   SUBROUTINE isfcpl_cons(Kmm)
       !!---------------------------------------------------------------------- 
       !!                   ***  ROUTINE iscpl_cons  ***
@@ -463,4 +472,7 @@
       TYPE(isfcons), DIMENSION(:),ALLOCATABLE :: zisfpts ! list of point receiving a correction
       !
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kmm    ! ocean time level index
+      !!----------------------------------------------------------------------
       INTEGER  ::   ji   , jj  , jk  , jproc          ! loop index
       INTEGER  ::   jip1 , jim1, jjp1, jjm1           ! dummy indices
@@ -513,12 +525,12 @@
 
                ! volume diff
-               zdvol = e3t_n(ji,jj,jk) * tmask(ji,jj,jk) - ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
+               zdvol = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) - ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
 
                ! heat diff
-               zdtem = tsn (ji,jj,jk,jp_tem) *  e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
+               zdtem = ts (ji,jj,jk,jp_tem,Kmm) *  e3t(ji,jj,jk,Kmm) *  tmask  (ji,jj,jk)   &
                      - zt_b(ji,jj,jk)        * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
 
                ! salt diff
-               zdsal = tsn(ji,jj,jk,jp_sal) *  e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
+               zdsal = ts(ji,jj,jk,jp_sal,Kmm) *  e3t(ji,jj,jk,Kmm) *  tmask  (ji,jj,jk)   &
                      - zs_b(ji,jj,jk)       * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
             
