#define TWO_WAY        /* TWO WAY NESTING */
#undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/
 
MODULE agrif_opa_update
#if defined key_agrif 
   USE par_oce
   USE oce
   USE dom_oce
   USE agrif_oce
   USE in_out_manager  ! I/O manager
   USE lib_mpp
   USE wrk_nemo  
   USE zdf_oce        ! vertical physics: ocean variables 
   USE domvvl         ! Need interpolation routines 

   IMPLICIT NONE
   PRIVATE

   PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Update_Scales, Agrif_Update_vvl

# if defined key_zdftke
   PUBLIC Agrif_Update_Tke
# endif
   !!----------------------------------------------------------------------
   !! NEMO/NST 3.6 , NEMO Consortium (2010)
   !! $Id$
   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE Agrif_Update_Tra( )
      !!---------------------------------------------
      !!   *** ROUTINE Agrif_Update_Tra ***
      !!---------------------------------------------
      ! 
      IF (Agrif_Root()) RETURN
      !
#if defined TWO_WAY  
      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed(), 'nbcline', nbcline

      Agrif_UseSpecialValueInUpdate = .TRUE.
      Agrif_SpecialValueFineGrid = 0.
      ! 
      IF (MOD(nbcline,nbclineupdate) == 0) THEN
# if ! defined DECAL_FEEDBACK
         CALL Agrif_Update_Variable(tsn_id, procname=updateTS)
# else
         CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS)
# endif
      ELSE
# if ! defined DECAL_FEEDBACK
         CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS)
# else
         CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS)
# endif
      ENDIF
      !
      Agrif_UseSpecialValueInUpdate = .FALSE.
      !
#endif
      !
   END SUBROUTINE Agrif_Update_Tra

   SUBROUTINE Agrif_Update_Dyn( )
      !!---------------------------------------------
      !!   *** ROUTINE Agrif_Update_Dyn ***
      !!---------------------------------------------
      ! 
      IF (Agrif_Root()) RETURN
      !
#if defined TWO_WAY
      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed(), 'nbcline', nbcline

      Agrif_UseSpecialValueInUpdate = .FALSE.
      Agrif_SpecialValueFineGrid = 0.
      !     
      IF (mod(nbcline,nbclineupdate) == 0) THEN
# if ! defined DECAL_FEEDBACK
         CALL Agrif_Update_Variable(un_update_id,procname = updateU)
         CALL Agrif_Update_Variable(vn_update_id,procname = updateV)
# else
         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU)
         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV)
# endif
      ELSE
# if ! defined DECAL_FEEDBACK
         CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU)
         CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV)         
# else
         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU)
         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV)
# endif
      ENDIF

# if ! defined DECAL_FEEDBACK
      CALL Agrif_Update_Variable(e1u_id,procname = updateU2d)
      CALL Agrif_Update_Variable(e2v_id,procname = updateV2d)  
# else
      CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d)
      CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)  
# endif

      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
         ! Update time integrated transports
         IF (mod(nbcline,nbclineupdate) == 0) THEN
#  if ! defined DECAL_FEEDBACK
            CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)
            CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)
#  else
            CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b)
            CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b)
#  endif
         ELSE
#  if ! defined DECAL_FEEDBACK
            CALL Agrif_Update_Variable(ub2b_update_id,locupdate=(/0,1/),procname = updateub2b)
            CALL Agrif_Update_Variable(vb2b_update_id,locupdate=(/0,1/),procname = updatevb2b)
#  else
            CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateub2b)
            CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updatevb2b)
#  endif
         ENDIF
      END IF
      !
      nbcline = nbcline + 1
      !
      Agrif_UseSpecialValueInUpdate = .TRUE.
      Agrif_SpecialValueFineGrid = 0.
# if ! defined DECAL_FEEDBACK
      CALL Agrif_Update_Variable(sshn_id,procname = updateSSH)
# else
      CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH)
# endif
      Agrif_UseSpecialValueInUpdate = .FALSE.
      ! 
#endif
      !
   END SUBROUTINE Agrif_Update_Dyn

# if defined key_zdftke

   SUBROUTINE Agrif_Update_Tke( kt )
      !!---------------------------------------------
      !!   *** ROUTINE Agrif_Update_Tke ***
      !!---------------------------------------------
      !!
      INTEGER, INTENT(in) :: kt 
      ! 
      IF (Agrif_Root()) RETURN
      !       
      IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN
#  if defined TWO_WAY

      Agrif_UseSpecialValueInUpdate = .TRUE.
      Agrif_SpecialValueFineGrid = 0.

      CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN  )
      CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT )
      CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM )

      Agrif_UseSpecialValueInUpdate = .FALSE.

#  endif
      
   END SUBROUTINE Agrif_Update_Tke
   
# endif /* key_zdftke */

   SUBROUTINE Agrif_Update_vvl( )
      !!---------------------------------------------
      !!   *** ROUTINE Agrif_Update_vvl ***
      !!---------------------------------------------
      !
      IF (Agrif_Root()) RETURN
      !
#if defined TWO_WAY  
      !
      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step()
      !
      Agrif_UseSpecialValueInUpdate = .TRUE.
      Agrif_SpecialValueFineGrid = 0.
      ! 
# if ! defined DECAL_FEEDBACK
      CALL Agrif_Update_Variable(e3t_id, procname=updatee3t)
# else
      CALL Agrif_Update_Variable(e3t_id, locupdate=(/1,0/), procname=updatee3t)
# endif 
      !
      Agrif_UseSpecialValueInUpdate = .FALSE.
      !
      CALL Agrif_ChildGrid_To_ParentGrid()
      CALL dom_vvl_update_UVF
      CALL Agrif_ParentGrid_To_ChildGrid()
      !
#endif
      !
   END SUBROUTINE Agrif_Update_vvl

   SUBROUTINE dom_vvl_update_UVF
      !!---------------------------------------------
      !!       *** ROUTINE dom_vvl_update_UVF ***
      !!---------------------------------------------
      !!
      INTEGER :: jk
      REAL(wp):: zcoef
      !!---------------------------------------------

      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', &
                  & Agrif_Fixed(), 'Step', Agrif_Nb_Step()

      ! Save "old" scale factor (prior update) for subsequent asselin correction
      ! of prognostic variables (needed to update initial state only)
      ! -------------------------------------------------------------
      !
      e3u_a(:,:,:) = e3u_n(:,:,:)
      e3v_a(:,:,:) = e3v_n(:,:,:)
!      ua(:,:,:) = e3u_b(:,:,:)
!      va(:,:,:) = e3v_b(:,:,:)
      hu_a(:,:) = hu_n(:,:)
      hv_a(:,:) = hv_n(:,:)

      ! 1) NOW fields
      !--------------
      
         ! Vertical scale factor interpolations
         ! ------------------------------------
      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) ,  'U' )
      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) ,  'V' )
      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) ,  'F' )

      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )

         ! Update total depths:
         ! --------------------
      hu_n(:,:) = 0._wp                        ! Ocean depth at U-points
      hv_n(:,:) = 0._wp                        ! Ocean depth at V-points
      DO jk = 1, jpkm1
         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
      END DO
      !                                        ! Inverse of the local depth
      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )


      ! 2) BEFORE fields:
      !------------------
      IF (     (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) &
         & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts    &
         & .AND.(.NOT.ln_bt_fw)))) THEN
         !
         ! Vertical scale factor interpolations
         ! ------------------------------------
         CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:),  'U'  )
         CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:),  'V'  )

         CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
         CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )

         ! Update total depths:
         ! --------------------
         hu_b(:,:) = 0._wp                     ! Ocean depth at U-points
         hv_b(:,:) = 0._wp                     ! Ocean depth at V-points
         DO jk = 1, jpkm1
            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
         END DO
         !                                     ! Inverse of the local depth
         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
      ENDIF
      !
   END SUBROUTINE dom_vvl_update_UVF

   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
      !!---------------------------------------------
      !!           *** ROUTINE updateT ***
      !!---------------------------------------------
      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
      LOGICAL, INTENT(in) :: before
      !!
      INTEGER :: ji,jj,jk,jn
      REAL(wp) :: ztb, ztnu, ztno
      !!---------------------------------------------
      !
      IF (before) THEN
         DO jn = n1,n2
            DO jk=k1,k2
               DO jj=j1,j2
                  DO ji=i1,i2
!> jc tmp
                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
!                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk)
!< jc tmp
                  END DO
               END DO
            END DO
         END DO
      ELSE
!> jc tmp
         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)
         ENDDO
!< jc tmp
         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
            ! Add asselin part
            DO jn = n1,n2
               DO jk=k1,k2
                  DO jj=j1,j2
                     DO ji=i1,i2
                        IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
                           ztb  = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used
                           ztnu = tabres(ji,jj,jk,jn)
                           ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
                           tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
                        ENDIF
                     ENDDO
                  ENDDO
               ENDDO
            ENDDO
         ENDIF
         DO jn = n1,n2
            DO jk=k1,k2
               DO jj=j1,j2
                  DO ji=i1,i2
                     IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 
                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
                     END IF
                  END DO
               END DO
            END DO
         END DO
         !
         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
            tsb(i1:i2,j1:j2,k1:k2,n1:n2)  = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
         ENDIF
         !
      ENDIF
      ! 
   END SUBROUTINE updateTS


   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before )
      !!---------------------------------------------
      !!           *** ROUTINE updateu ***
      !!---------------------------------------------
      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
      LOGICAL                               , INTENT(in   ) :: before
      !
      INTEGER  :: ji, jj, jk
      REAL(wp) :: zrhoy, zub, zunu, zuno
      !!---------------------------------------------
      ! 
      IF( before ) THEN
         zrhoy = Agrif_Rhoy()
         DO jk = k1, k2
            tabres(i1:i2,j1:j2,jk) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
         END DO
      ELSE
         DO jk=k1,k2
            DO jj=j1,j2
               DO ji=i1,i2
                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) 
                  !
                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
                     zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)
                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
                     zunu = tabres(ji,jj,jk)
                     ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &      
                                    & * umask(ji,jj,jk) / e3u_b(ji,jj,jk)
                  ENDIF
                  !
                  un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)
               END DO
            END DO
         END DO
         !
         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2)
         ENDIF
         !
      ENDIF
      ! 
   END SUBROUTINE updateu


   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before)
      !!---------------------------------------------
      !!           *** ROUTINE updatev ***
      !!---------------------------------------------
      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
      LOGICAL                               , INTENT(in   ) :: before
      !
      INTEGER  :: ji, jj, jk
      REAL(wp) :: zrhox, zvb, zvnu, zvno
      !!---------------------------------------------      
      !
      IF (before) THEN
         zrhox = Agrif_Rhox()
         DO jk=k1,k2
            DO jj=j1,j2
               DO ji=i1,i2
                  tabres(ji,jj,jk) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
               END DO
            END DO
         END DO
      ELSE
         DO jk=k1,k2
            DO jj=j1,j2
               DO ji=i1,i2
                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj)
                  !
                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
                     zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk)
                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
                     zvnu = tabres(ji,jj,jk)
                     vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &      
                                    & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk)
                  ENDIF
                  !
                  vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)
               END DO
            END DO
         END DO
         !
         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2)
         ENDIF
         !
      ENDIF
      ! 
   END SUBROUTINE updatev


   SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
      !!---------------------------------------------
      !!          *** ROUTINE updateu2d ***
      !!---------------------------------------------
      INTEGER, INTENT(in) :: i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
      LOGICAL, INTENT(in) :: before
      !! 
      INTEGER  :: ji, jj, jk
      REAL(wp) :: zrhoy
      REAL(wp) :: zcorr
      !!---------------------------------------------
      !
      IF (before) THEN
         zrhoy = Agrif_Rhoy()
         DO jj=j1,j2
            DO ji=i1,i2
               tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
            END DO
         END DO
      ELSE
         DO jj=j1,j2
            DO ji=i1,i2
               tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj)  
               !    
               ! Update "now" 3d velocities:
               spgu(ji,jj) = 0._wp
               DO jk=1,jpkm1
                  spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
               END DO
               !
               zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj)
               DO jk=1,jpkm1              
                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
               END DO
               !
               ! 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) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
                  END IF
               ENDIF    
               un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
               !       
               ! Correct "before" velocities to hold correct bt component:
               spgu(ji,jj) = 0.e0
               DO jk=1,jpkm1
                  spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
               END DO
               !
               zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)
               DO jk=1,jpkm1              
                  ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
               END DO
               !
            END DO
         END DO
         !
         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2)
         ENDIF
      ENDIF
      !
   END SUBROUTINE updateu2d


   SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
      !!---------------------------------------------
      !!          *** ROUTINE updatev2d ***
      !!---------------------------------------------
      INTEGER, INTENT(in) :: i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
      LOGICAL, INTENT(in) :: before
      !! 
      INTEGER  :: ji, jj, jk
      REAL(wp) :: zrhox
      REAL(wp) :: zcorr
      !!---------------------------------------------
      !
      IF (before) THEN
         zrhox = Agrif_Rhox()
         DO jj=j1,j2
            DO ji=i1,i2
               tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj) 
            END DO
         END DO
      ELSE
         DO jj=j1,j2
            DO ji=i1,i2
               tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj)  
               !    
               ! Update "now" 3d velocities:
               spgv(ji,jj) = 0.e0
               DO jk=1,jpkm1
                  spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
               END DO
               !
               zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj)
               DO jk=1,jpkm1              
                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
               END DO
               !
               ! 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) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
                  END IF
               ENDIF              
               vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
               !       
               ! Correct "before" velocities to hold correct bt component:
               spgv(ji,jj) = 0.e0
               DO jk=1,jpkm1
                  spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
               END DO
               !
               zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)
               DO jk=1,jpkm1              
                  vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
               END DO
               !
            END DO
         END DO
         !
         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2)
         ENDIF
         !
      ENDIF
      ! 
   END SUBROUTINE updatev2d


   SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
      !!---------------------------------------------
      !!          *** ROUTINE updateSSH ***
      !!---------------------------------------------
      INTEGER, INTENT(in) :: i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
      LOGICAL, INTENT(in) :: before
      !!
      INTEGER :: ji, jj
      !!---------------------------------------------
      ! 
      IF (before) THEN
         DO jj=j1,j2
            DO ji=i1,i2
               tabres(ji,jj) = sshn(ji,jj)
            END DO
         END DO
      ELSE
         IF( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. .NOT.ln_bt_fw ) ) THEN
            IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
               DO jj=j1,j2
                  DO ji=i1,i2
                     sshb(ji,jj) =   sshb(ji,jj) &
                           & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
                  END DO
               END DO
            ENDIF
         ENDIF
         !
         DO jj=j1,j2
            DO ji=i1,i2
               sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
            END DO
         END DO
         !
         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2)
         ENDIF
         !

      ENDIF
      !
   END SUBROUTINE updateSSH


   SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before, nb, ndir )
      !!---------------------------------------------
      !!          *** ROUTINE updateub2b ***
      !!---------------------------------------------
      INTEGER, INTENT(in) :: i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
      LOGICAL, INTENT(in) :: before
      INTEGER, INTENT(in) :: nb, ndir
      !!
      LOGICAL :: western_side, eastern_side 
      INTEGER :: ji, jj
      REAL(wp) :: zrhoy, za1
      !!---------------------------------------------
      !
      IF (before) THEN
         zrhoy = Agrif_Rhoy()
         DO jj=j1,j2
            DO ji=i1,i2
               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
            END DO
         END DO
         tabres = zrhoy * tabres
      ELSE
         ! 
         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
         !
         ! Refluxing here:
#if defined DECAL_FEEDBACK
         western_side  = (nb == 1).AND.(ndir == 1)
         eastern_side  = (nb == 1).AND.(ndir == 2)
         !
         IF (western_side) THEN
            DO jj=j1,j2
               sshn(i1  ,jj) = sshn(i1  ,jj) + rdt * r1_e1e2t(i1  ,jj) &
                     &         * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))
            END DO
         ENDIF
         IF (eastern_side) THEN
            DO jj=j1,j2
               sshn(i2+1,jj) = sshn(i2+1,jj) - rdt * r1_e1e2t(i2+1,jj) &
                     &         * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
            END DO
         ENDIF
         ! 
#endif
         za1 = 1._wp / REAL(Agrif_rhot(), wp)
         !
         DO jj=j1,j2
            DO ji=i1,i2
               ! Update time integrated fluxes also in case of multiply nested grids:
               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) & 
                & + za1 * (tabres(ji,jj) - ub2_b(ji,jj))
               ! Update half step back fluxes:
               ub2_b(ji,jj) = tabres(ji,jj)
            END DO
         END DO
      ENDIF
      !
   END SUBROUTINE updateub2b


   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before, nb, ndir )
      !!---------------------------------------------
      !!          *** ROUTINE updatevb2b ***
      !!---------------------------------------------
      INTEGER, INTENT(in) :: i1, i2, j1, j2
      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
      LOGICAL, INTENT(in) :: before
      INTEGER, INTENT(in) :: nb, ndir
      !!
      LOGICAL :: southern_side, northern_side 
      INTEGER :: ji, jj
      REAL(wp) :: zrhox, za1
      !!---------------------------------------------
      !
      IF (before) THEN
         zrhox = Agrif_Rhox()
         DO jj=j1,j2
            DO ji=i1,i2
               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
            END DO
         END DO
         tabres = zrhox * tabres
      ELSE
         ! 
         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
         !
         ! Refluxing here:
#if defined DECAL_FEEDBACK
         southern_side = (nb == 2).AND.(ndir == 1)
         northern_side = (nb == 2).AND.(ndir == 2)
         !
         IF (southern_side) THEN
            DO ji=i1,i2
               sshn(ji,j1  ) = sshn(ji,j1  ) + rdt * r1_e1e2t(ji,j1  ) &
                     &         * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
            END DO
         ENDIF
         IF (northern_side) THEN               
            DO ji=i1,i2
               sshn(ji,j2+1) = sshn(ji,j2+1) - rdt * r1_e1e2t(ji,j2+1) &
                     &         * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
            END DO
         ENDIF
         ! 
#endif
         za1 = 1._wp / REAL(Agrif_rhot(), wp)
         DO jj=j1,j2
            DO ji=i1,i2
               ! Update time integrated fluxes also in case of multiply nested grids:
               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) & 
                & + za1 * (tabres(ji,jj) - vb2_b(ji,jj))
               ! Update half step back fluxes:
               vb2_b(ji,jj) = tabres(ji,jj)
            END DO
         END DO
      ENDIF
      !
   END SUBROUTINE updatevb2b


   SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
      ! currently not used
      !!---------------------------------------------
      !!           *** ROUTINE updateT ***
      !!---------------------------------------------
      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
      LOGICAL, iNTENT(in) :: before
      !
      INTEGER :: ji,jj,jk
      REAL(wp) :: ztemp
      !!---------------------------------------------

      IF (before) THEN
         DO jk=k1,k2
            DO jj=j1,j2
               DO ji=i1,i2
                  tabres(ji,jj,jk,1) = e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
                  tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
                  tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
               END DO
            END DO
         END DO
         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
            DO jj=j1,j2
               DO ji=i1,i2
                  IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN 
                     print *,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
                     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)))
                     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
                  END IF
               END DO
            END DO
         END DO
      ENDIF
      !
   END SUBROUTINE update_scales

# if defined key_zdftke

   SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
      !!---------------------------------------------
      !!           *** ROUTINE updateen ***
      !!---------------------------------------------
      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
      LOGICAL, INTENT(in) :: before
      !!---------------------------------------------
      !
      IF (before) THEN
         ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
      ELSE
         en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
      ENDIF
      !
   END SUBROUTINE updateEN


   SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
      !!---------------------------------------------
      !!           *** ROUTINE updateavt ***
      !!---------------------------------------------
      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
      LOGICAL, INTENT(in) :: before
      !!---------------------------------------------
      !
      IF (before) THEN
         ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
      ELSE
         avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
      ENDIF
      !
   END SUBROUTINE updateAVT


   SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
      !!---------------------------------------------
      !!           *** ROUTINE updateavm ***
      !!---------------------------------------------
      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
      LOGICAL, INTENT(in) :: before
      !!---------------------------------------------
      !
      IF (before) THEN
         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
      ELSE
         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
      ENDIF
      !
   END SUBROUTINE updateAVM

# endif /* key_zdftke */ 

   SUBROUTINE updatee3t( ptab, i1, i2, j1, j2, k1, k2, before )
      !!---------------------------------------------
      !!           *** ROUTINE updatee3t ***
      !!---------------------------------------------
      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
      LOGICAL, INTENT(in) :: before
      INTEGER :: ji,jj,jk
      REAL(wp) :: zcoef
      !!---------------------------------------------
      !
      IF (before) THEN
!> jc tmp:
!         ptab(i1:i2,j1:j2,k1:k2) = e3t_n(i1:i2,j1:j2,k1:k2)
         ptab(i1:i2,j1:j2,k1:k2) = e3t_n(i1:i2,j1:j2,k1:k2) / e3t_0(i1:i2,j1:j2,k1:k2) * tmask(i1:i2,j1:j2,k1:k2)
!< jc tmp:
      ELSE
         !
         ! 1) Updates at BEFORE time step:
         ! -------------------------------
         !
!> jc tmp:
         DO jk = 1, jpkm1
            DO jj=j1,j2
               DO ji=i1,i2
                  IF (tmask(ji,jj,jk)==1) THEN
                     ptab(ji,jj,jk) = ptab(ji,jj,jk) * e3t_0(ji,jj,jk)
                  ELSE
                     ptab(ji,jj,jk) = e3t_0(ji,jj,jk)
                  ENDIF
               END DO
            END DO
         END DO
!         ptab(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
!< jc tmp:

         ! Save "old" scale factor (prior update) for subsequent asselin correction
         ! of prognostic variables (needed to update initial state only)
         e3t_a(i1:i2,j1:j2,k1:k2) = e3t_n(i1:i2,j1:j2,k1:k2)
!         hdivn(i1:i2,j1:j2,k1:k2)   = e3t_b(i1:i2,j1:j2,k1:k2)

         IF (     (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) &
            & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts    &
            & .AND.(.NOT.ln_bt_fw)))) THEN

            DO jk = 1, jpkm1
               DO jj=j1,j2
                  DO ji=i1,i2
                     e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) &
                           & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) )
                  END DO
               END DO
            END DO
            !
            e3w_b  (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_b(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
            gdepw_b(i1:i2,j1:j2,1) = 0.0_wp
            gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1)
            !
            DO jk = 2, jpk
               DO jj = j1,j2
                  DO ji = i1,i2            
                     zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
                     e3w_b(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
                     &                                        ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )  &
                     &                                  +            0.5_wp * tmask(ji,jj,jk)   *        &
                     &                                        ( e3t_b(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
                     gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
                     gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
                         &               + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
                  END DO
               END DO
            END DO
            !
         ENDIF        
         !
         ! 2) Updates at NOW time step:
         ! ----------------------------
         !
         ! Update vertical scale factor at T-points:
         e3t_n(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2)
         !
         ! Update total depth:
         ht_n(i1:i2,j1:j2) = 0._wp
         DO jk = 1, jpkm1
            ht_n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t_n(i1:i2,j1:j2,jk) * tmask(i1:i2,j1:j2,jk)
         END DO
         !
         ! Update vertical scale factor at W-points and depths:
         e3w_n (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_n(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
         gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1)
         gdepw_n(i1:i2,j1:j2,1) = 0.0_wp
         gde3w_n(i1:i2,j1:j2,1) = gdept_n(i1:i2,j1:j2,1) - (ht_n(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh
         !
         DO jk = 2, jpk
            DO jj = j1,j2
               DO ji = i1,i2            
               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
               e3w_n(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t_n(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )   &
               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t_n(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
                   &               + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
               END DO
            END DO
         END DO
         !
         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
            e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk)
            e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk)
            gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk)
            gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk)
         ENDIF
         !
      ENDIF
      !
   END SUBROUTINE updatee3t

#else
CONTAINS
   SUBROUTINE agrif_opa_update_empty
      !!---------------------------------------------
      !!   *** ROUTINE agrif_opa_update_empty ***
      !!---------------------------------------------
      WRITE(*,*)  'agrif_opa_update : You should not have seen this print! error?'
   END SUBROUTINE agrif_opa_update_empty
#endif
END MODULE agrif_opa_update