Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90	(revision 10919)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90	(revision 10921)
@@ -94,5 +94,5 @@
 CONTAINS
 
-   SUBROUTINE dta_dyn( kt )
+   SUBROUTINE dta_dyn( kt, Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dta_dyn  ***
@@ -105,5 +105,6 @@
       !!             - interpolates data if needed
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt             ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa  ! ocean time level indices
       !
       INTEGER             ::   ji, jj, jk
@@ -123,6 +124,6 @@
       IF( l_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt )    ! Computation of slopes
       !
-      tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
-      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity
+      ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
+      ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity
       wndm(:,:)         = sf_dyn(jf_wnd)%fnow(:,:,1)  * tmask(:,:,1)    ! wind speed - needed for gas exchange
       fmmflx(:,:)       = sf_dyn(jf_fmf)%fnow(:,:,1)  * tmask(:,:,1)    ! downward salt flux (v3.5+)
@@ -132,10 +133,10 @@
       IF( ln_dynrnf ) THEN 
          rnf (:,:)      = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
-         IF( ln_dynrnf_depth .AND. .NOT. ln_linssh )    CALL  dta_dyn_hrnf
-      ENDIF
-      !
-      un(:,:,:)        = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport
-      vn(:,:,:)        = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport
-      wn(:,:,:)        = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport
+         IF( ln_dynrnf_depth .AND. .NOT. ln_linssh )    CALL  dta_dyn_hrnf(Kmm)
+      ENDIF
+      !
+      uu(:,:,:,Kmm)        = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport
+      vv(:,:,:,Kmm)        = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport
+      ww(:,:,:)        = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport
       !
       IF( .NOT.ln_linssh ) THEN
@@ -144,5 +145,5 @@
          emp_b  (:,:)   = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
          zemp   (:,:)   = ( 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr ) * tmask(:,:,1)
-         CALL dta_dyn_ssh( kt, zhdivtr, sshb, zemp, ssha, e3t_a(:,:,:) )  !=  ssh, vertical scale factor & vertical transport
+         CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Kbb), zemp, ssh(:,:,Kaa), e3t(:,:,:,Kaa) )  !=  ssh, vertical scale factor & vertical transport
          DEALLOCATE( zemp , zhdivtr )
          !                                           Write in the tracer restart file
@@ -152,12 +153,12 @@
             IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file at it= ', kt,' date= ', ndastp
             IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha )
-            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn )
+            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssh(:,:,Kaa) )
+            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', ssh(:,:,Kmm) )
          ENDIF
       ENDIF
       !
-      CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
-      CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points
-      CALL bn2    ( tsn, rab_n, rn2 ) ! before Brunt-Vaisala frequency need for zdfmxl
+      CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
+      CALL eos_rab( ts(:,:,:,:,Kmm), rab_n )       ! now    local thermal/haline expension ratio at T-points
+      CALL bn2    ( ts(:,:,:,:,Kmm), rab_n, rn2 ) ! before Brunt-Vaisala frequency need for zdfmxl
 
       rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
@@ -174,12 +175,12 @@
       !
       !
-      CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
+      CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
       !
       IF(ln_ctl) THEN                  ! print control
-         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   )
-         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   )
-         CALL prt_ctl(tab3d_1=un               , clinfo1=' un      - : ', mask1=umask,  kdim=jpk   )
-         CALL prt_ctl(tab3d_1=vn               , clinfo1=' vn      - : ', mask1=vmask,  kdim=jpk   )
-         CALL prt_ctl(tab3d_1=wn               , clinfo1=' wn      - : ', mask1=tmask,  kdim=jpk   )
+         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   )
+         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   )
+         CALL prt_ctl(tab3d_1=uu(:,:,:,Kmm)               , clinfo1=' uu(:,:,:,Kmm)      - : ', mask1=umask,  kdim=jpk   )
+         CALL prt_ctl(tab3d_1=vv(:,:,:,Kmm)               , clinfo1=' vv(:,:,:,Kmm)      - : ', mask1=vmask,  kdim=jpk   )
+         CALL prt_ctl(tab3d_1=wn               , clinfo1=' ww      - : ', mask1=tmask,  kdim=jpk   )
          CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask,  kdim=jpk   )
          CALL prt_ctl(tab3d_1=uslp             , clinfo1=' slp  - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk)
@@ -192,5 +193,5 @@
 
 
-   SUBROUTINE dta_dyn_init
+   SUBROUTINE dta_dyn_init( Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dta_dyn_init  ***
@@ -199,4 +200,6 @@
       !! ** Method  : - read the data namdta_dyn namelist
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa  ! ocean time level indices
+      !
       INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
       INTEGER  :: ifpr                               ! dummy loop indice
@@ -312,37 +315,37 @@
         IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file
            iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
-           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation'
-           CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:)   )
-           CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:)   )
+           IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation'
+           CALL iom_get( numrtr, jpdom_autoglo, 'sshn', ssh(:,:,Kmm)   )
+           CALL iom_get( numrtr, jpdom_autoglo, 'sshb', ssh(:,:,Kbb)   )
         ELSE
-           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation'
+           IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation'
            CALL iom_open( 'restart', inum )
-           CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:)   )
-           CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:)   )
+           CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh(:,:,Kmm)   )
+           CALL iom_get( inum, jpdom_autoglo, 'sshb', ssh(:,:,Kbb)   )
            CALL iom_close( inum )                                        ! close file
         ENDIF
         !
         DO jk = 1, jpkm1
-           e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
+           e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
         ENDDO
-        e3t_a(:,:,jpk) = e3t_0(:,:,jpk)
+        e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)
 
         ! Horizontal scale factor interpolations
         ! --------------------------------------
-        CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-        CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
+        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
+        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
 
         ! Vertical scale factor interpolations
         ! ------------------------------------
-        CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' )
+        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' )
   
-        e3t_b(:,:,:)  = e3t_n(:,:,:)
-        e3u_b(:,:,:)  = e3u_n(:,:,:)
-        e3v_b(:,:,:)  = e3v_n(:,:,:)
+        e3t(:,:,:,Kbb)  = e3t(:,:,:,Kmm)
+        e3u(:,:,:,Kbb)  = e3u(:,:,:,Kmm)
+        e3v(:,:,:,Kbb)  = e3v(:,:,:,Kmm)
 
         ! t- and w- points depth
         ! ----------------------
-        gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-        gdepw_n(:,:,1) = 0.0_wp
+        gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
+        gdepw(:,:,1,Kmm) = 0.0_wp
 
         DO jk = 2, jpk
@@ -354,13 +357,13 @@
                                                                    ! mikt
                  zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-                 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-                 gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
-                     &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk))
+                 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
+                 gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
+                     &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm))
               END DO
            END DO
         END DO
 
-        gdept_b(:,:,:) = gdept_n(:,:,:)
-        gdepw_b(:,:,:) = gdepw_n(:,:,:)
+        gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
+        gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
         !
       ENDIF
@@ -393,5 +396,5 @@
                h_rnf(ji,jj) = 0._wp
                DO jk = 1, nk_rnf(ji,jj)
-                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)
+                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -399,5 +402,5 @@
       ELSE                                       ! runoffs applied at the surface
          nk_rnf(:,:) = 1
-         h_rnf (:,:) = e3t_n(:,:,1)
+         h_rnf (:,:) = e3t(:,:,1,Kmm)
       ENDIF
       nkrnf_max = MAXVAL( nk_rnf(:,:) )
@@ -411,9 +414,9 @@
       IF(lwp) WRITE(numout,*) ' '
       !
-      CALL dta_dyn( nit000 )
+      CALL dta_dyn( nit000, Kbb, Kmm, Kaa )
       !
    END SUBROUTINE dta_dyn_init
 
-   SUBROUTINE dta_dyn_sed( kt )
+   SUBROUTINE dta_dyn_sed( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dta_dyn  ***
@@ -427,4 +430,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
       !
       !!----------------------------------------------------------------------
@@ -439,12 +443,12 @@
       CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==!
       !
-      tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
-      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity
-      !
-      CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
+      ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
+      ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity
+      !
+      CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
 
       IF(ln_ctl) THEN                  ! print control
-         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   )
-         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   )
+         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   )
+         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   )
       ENDIF
       !
@@ -454,5 +458,5 @@
 
 
-   SUBROUTINE dta_dyn_sed_init
+   SUBROUTINE dta_dyn_sed_init( Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dta_dyn_init  ***
@@ -461,4 +465,6 @@
       !! ** Method  : - read the data namdta_dyn namelist
       !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: Kmm                   ! ocean time level index
+      !
       INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
       INTEGER  :: ifpr                               ! dummy loop indice
@@ -525,9 +531,9 @@
       END DO
       !
-      CALL dta_dyn_sed( nit000 )
+      CALL dta_dyn_sed( nit000, Kmm )
       !
    END SUBROUTINE dta_dyn_sed_init
 
-   SUBROUTINE dta_dyn_swp( kt )
+   SUBROUTINE dta_dyn_swp( kt, Kbb, Kmm, Kaa )
      !!---------------------------------------------------------------------
       !!                    ***  ROUTINE dta_dyn_swp  ***
@@ -536,5 +542,6 @@
       !!              at U/V/W pointand the depht
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) :: kt       ! time step
+      INTEGER, INTENT(in) :: kt             ! time step
+      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa  ! ocean time level indices
       !
       INTEGER             :: ji, jj, jk
@@ -548,8 +555,8 @@
       ENDIF
 
-      sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:))  ! before <-- now filtered
-      sshn(:,:) = ssha(:,:)
-
-      e3t_n(:,:,:) = e3t_a(:,:,:)
+      ssh(:,:,Kbb) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa))  ! before <-- now filtered
+      ssh(:,:,Kmm) = ssh(:,:,Kaa)
+
+      e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa)
 
       ! Reconstruction of all vertical scale factors at now and before time steps
@@ -558,19 +565,19 @@
       ! Horizontal scale factor interpolations
       ! --------------------------------------
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
 
       ! Vertical scale factor interpolations
       ! ------------------------------------
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
-
-      e3t_b(:,:,:)  = e3t_n(:,:,:)
-      e3u_b(:,:,:)  = e3u_n(:,:,:)
-      e3v_b(:,:,:)  = e3v_n(:,:,:)
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' )
+
+      e3t(:,:,:,Kbb)  = e3t(:,:,:,Kmm)
+      e3u(:,:,:,Kbb)  = e3u(:,:,:,Kmm)
+      e3v(:,:,:,Kbb)  = e3v(:,:,:,Kmm)
 
       ! t- and w- points depth
       ! ----------------------
-      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-      gdepw_n(:,:,1) = 0.0_wp
+      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
+      gdepw(:,:,1,Kmm) = 0.0_wp
       !
       DO jk = 2, jpk
@@ -578,13 +585,13 @@
             DO ji = 1,jpi
                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
-                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk))
+               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
+               gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
+                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm))
             END DO
          END DO
       END DO
       !
-      gdept_b(:,:,:) = gdept_n(:,:,:)
-      gdepw_b(:,:,:) = gdepw_n(:,:,:)
+      gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
+      gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
       !
    END SUBROUTINE dta_dyn_swp
@@ -595,5 +602,5 @@
       !!                ***  ROUTINE dta_dyn_wzv  ***
       !!                   
-      !! ** Purpose :   compute the after ssh (ssha) and the now vertical velocity
+      !! ** Purpose :   compute the after ssh (ssh(:,:,Kaa)) and the now vertical velocity
       !!
       !! ** Method  : Using the incompressibility hypothesis, 
@@ -608,5 +615,5 @@
       !!          The boundary conditions are w=0 at the bottom (no flux)
       !!
-      !! ** action  :   ssha / e3t_a / wn
+      !! ** action  :   ssh(:,:,Kaa) / e3t(:,:,:,Kaa) / ww
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
@@ -641,5 +648,5 @@
 
 
-   SUBROUTINE dta_dyn_hrnf
+   SUBROUTINE dta_dyn_hrnf( Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_rnf  ***
@@ -654,5 +661,7 @@
       !!----------------------------------------------------------------------
       !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER, INTENT(in) :: Kmm  ! ocean time level index
+      !
+      INTEGER  ::   ji, jj, jk    ! dummy loop indices
       !!----------------------------------------------------------------------
       !
@@ -661,5 +670,5 @@
             h_rnf(ji,jj) = 0._wp
             DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres
-                h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)   ! to the bottom of the relevant grid box
+                h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)   ! to the bottom of the relevant grid box
             END DO
         END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90	(revision 10919)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90	(revision 10921)
@@ -117,10 +117,10 @@
                                 CALL iom_setkt  ( istp - nit000 + 1, cxios_context )   ! say to iom that we are at time step kstp
 #if defined key_sed_off
-                                CALL dta_dyn_sed( istp )         ! Interpolation of the dynamical fields
+                                CALL dta_dyn_sed( istp,      Nnn      )       ! Interpolation of the dynamical fields
 #else
-                                CALL dta_dyn    ( istp )         ! Interpolation of the dynamical fields
-         IF( .NOT.ln_linssh )   CALL dta_dyn_swp( istp )         ! swap of sea  surface height and vertical scale factors
+                                CALL dta_dyn    ( istp, Nbb, Nnn, Naa )       ! Interpolation of the dynamical fields
+         IF( .NOT.ln_linssh )   CALL dta_dyn_swp( istp, Nbb, Nnn, Naa )       ! swap of sea  surface height and vertical scale factors
 #endif
-                                CALL trc_stp    ( istp, Nbb, Nnn, Nrhs, Naa )         ! time-stepping
+                                CALL trc_stp    ( istp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping
                                 CALL stp_ctl    ( istp, indic )  ! Time loop: control and print
          istp = istp + 1
@@ -331,7 +331,7 @@
                            CALL trc_rst_cal( nit000, 'READ' )   ! calendar
 #if defined key_sed_off
-                           CALL dta_dyn_sed_init ! Initialization for the dynamics
+                           CALL dta_dyn_sed_init(  Nnn      )        ! Initialization for the dynamics
 #else
-                           CALL dta_dyn_init   ! Initialization for the dynamics
+                           CALL dta_dyn_init( Nbb, Nnn, Naa )        ! Initialization for the dynamics
 #endif
 
