Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 7761)
@@ -445,4 +445,18 @@
       sxyage (:,:,:)  = 0._wp
 
+      !------------------------------------
+      ! 6) store fields at before time-step
+      !------------------------------------
+      ! it is only necessary for the 1st interpolation by Agrif
+      a_i_b  (:,:,:)   = a_i  (:,:,:)
+      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
+      v_i_b  (:,:,:)   = v_i  (:,:,:)
+      v_s_b  (:,:,:)   = v_s  (:,:,:)
+      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
+      smv_i_b(:,:,:)   = smv_i(:,:,:)
+      oa_i_b (:,:,:)   = oa_i (:,:,:)
+      u_ice_b(:,:)     = u_ice(:,:)
+      v_ice_b(:,:)     = v_ice(:,:)
+
 !!!clem
 !!      ! Output the initial state and forcings
Index: /trunk/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90	(revision 7761)
@@ -54,10 +54,15 @@
       IF( Agrif_Root() )  RETURN
       !
-      IF( PRESENT( kiter ) ) THEN  ! interpolation at the child sub-time step (for ice rheology)
-         zbeta = ( REAL(lim_nbstep) - REAL(kitermax - kiter) / REAL(kitermax) ) /  &
-            &    ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
-      ELSE                         ! interpolation at the child time step
-         zbeta = REAL(lim_nbstep) / ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
-      ENDIF
+      SELECT CASE(cd_type)
+      CASE('U','V')
+         IF( PRESENT( kiter ) ) THEN  ! interpolation at the child sub-time step (only for ice rheology)
+            zbeta = ( REAL(lim_nbstep) - REAL(kitermax - kiter) / REAL(kitermax) ) /  &
+               &    ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
+         ELSE                         ! interpolation at the child time step
+            zbeta = REAL(lim_nbstep) / ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
+         ENDIF
+      CASE('T')
+            zbeta = REAL(lim_nbstep-1) / ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
+      END SELECT
       !
       Agrif_SpecialValue=-9999.
@@ -151,6 +156,6 @@
 
       !!-----------------------------------------------------------------------
-      ! clem: pkoi on n'utilise pas les quantités intégrées ici => before: * e1e2t ; after: * r1_e1e2t / rhox / rhoy
-      ! a priori c'est ok comme ca (cf ce qui est fait dans l'ocean). Je ne sais pas pkoi ceci dit
+      ! tracers are not multiplied by grid cell here => before: * e12t ; after: * r1_e12t / rhox / rhoy
+      ! and it is ok since we conserve tracers (same as in the ocean).
       ALLOCATE( ztab(SIZE(a_i_b,1),SIZE(a_i_b,2),SIZE(ptab,3)) )
      
@@ -177,133 +182,142 @@
       ELSE               ! child grid
 !! ==> The easiest interpolation is the following commented lines
-!!         jm = 1
-!!         DO jl = 1, jpl
-!!            a_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-!!            v_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-!!            v_s  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-!!            smv_i(i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-!!            oa_i (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-!!            DO jk = 1, nlay_s
-!!               e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-!!            ENDDO
-!!            DO jk = 1, nlay_i
-!!               e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-!!            ENDDO
-!!         ENDDO
+         jm = 1
+         DO jl = 1, jpl
+            a_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
+            v_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
+            v_s  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
+            smv_i(i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
+            oa_i (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
+            DO jk = 1, nlay_s
+               e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
+            ENDDO
+            DO jk = 1, nlay_i
+               e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
+            ENDDO
+         ENDDO
 
 !! ==> this is a more complex interpolation since we mix solutions over a couple of grid points
 !!     it is advised to use it for fields modified by high order schemes (e.g. advection UM5...)
-         ! record ztab
-         jm = 1
-         DO jl = 1, jpl
-            ztab(:,:,jm) = a_i_b  (:,:,jl) ; jm = jm + 1
-            ztab(:,:,jm) = v_i_b  (:,:,jl) ; jm = jm + 1
-            ztab(:,:,jm) = v_s_b  (:,:,jl) ; jm = jm + 1
-            ztab(:,:,jm) = smv_i_b(:,:,jl) ; jm = jm + 1
-            ztab(:,:,jm) = oa_i_b (:,:,jl) ; jm = jm + 1
-            DO jk = 1, nlay_s
-               ztab(:,:,jm) = e_s_b(:,:,jk,jl) ; jm = jm + 1
-            ENDDO
-            DO jk = 1, nlay_i
-               ztab(:,:,jm) = e_i_b(:,:,jk,jl) ; jm = jm + 1
-            ENDDO
-         ENDDO
-         !
-         ! borders of the domain
-         western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2)
-         southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2)
-         !
-         ! spatial smoothing
-         zrhox = Agrif_Rhox()
-         z1 =      ( zrhox - 1. ) * 0.5 
-         z3 =      ( zrhox - 1. ) / ( zrhox + 1. )
-         z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
-         z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
-         z2 = 1. - z1
-         z4 = 1. - z3
-         z5 = 1. - z6 - z7
-         !
-         ! Remove corners
-         imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2
-         IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3
-         IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2
-         IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3
-         IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2
-
-         ! smoothed fields
-         IF( eastern_side ) THEN
-            ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:)
-            DO jj = jmin, jmax
-               rswitch = 0.
-               IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1.
-               ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  &
-                  &                +      umask(nlci-2,jj,1)   *  &
-                  &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  &
-                  &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) )
-               ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1)
-            END DO
-         ENDIF
-         ! 
-         IF( northern_side ) THEN
-            ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:)
-            DO ji = imin, imax
-               rswitch = 0.
-               IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1.
-               ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  &
-                  &                +      vmask(ji,nlcj-2,1)   *  &
-                  &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) &
-                  &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) )
-               ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1)
-            END DO
-         END IF
-         !
-         IF( western_side) THEN
-            ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:)
-            DO jj = jmin, jmax
-               rswitch = 0.
-               IF( u_ice(2,jj) > 0._wp ) rswitch = 1.
-               ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  &
-                  &           +      umask(2,jj,1)   *   &
-                  &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) &
-                  &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) )
-               ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1)
-            END DO
-         ENDIF
-         !
-         IF( southern_side ) THEN
-            ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:)
-            DO ji = imin, imax
-               rswitch = 0.
-               IF( v_ice(ji,2) > 0._wp ) rswitch = 1.
-               ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  &
-                  &           +      vmask(ji,2,1)   *  &
-                  &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) &
-                  &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) )
-               ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1)
-            END DO
-         END IF
-         !
-         ! Treatment of corners
-         IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south
-         IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north
-         IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south
-         IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north
-
-         ! retrieve ice tracers
-         jm = 1
-         DO jl = 1, jpl
-            a_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            v_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            v_s  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            smv_i(i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            oa_i (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            DO jk = 1, nlay_s
-               e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            ENDDO
-            DO jk = 1, nlay_i
-               e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            ENDDO
-         ENDDO
-         
+!!        clem: for some reason (I don't know why), the following lines do not work 
+!!              with mpp (or in realistic configurations?). It makes the model crash
+!         ! record ztab
+!         jm = 1
+!         DO jl = 1, jpl
+!            ztab(:,:,jm) = a_i  (:,:,jl) ; jm = jm + 1
+!            ztab(:,:,jm) = v_i  (:,:,jl) ; jm = jm + 1
+!            ztab(:,:,jm) = v_s  (:,:,jl) ; jm = jm + 1
+!            ztab(:,:,jm) = smv_i(:,:,jl) ; jm = jm + 1
+!            ztab(:,:,jm) = oa_i (:,:,jl) ; jm = jm + 1
+!            DO jk = 1, nlay_s
+!               ztab(:,:,jm) = e_s(:,:,jk,jl) ; jm = jm + 1
+!            ENDDO
+!            DO jk = 1, nlay_i
+!               ztab(:,:,jm) = e_i(:,:,jk,jl) ; jm = jm + 1
+!            ENDDO
+!         ENDDO
+!         !
+!         ! borders of the domain
+!         western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2)
+!         southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2)
+!         !
+!         ! spatial smoothing
+!         zrhox = Agrif_Rhox()
+!         z1 =      ( zrhox - 1. ) * 0.5 
+!         z3 =      ( zrhox - 1. ) / ( zrhox + 1. )
+!         z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
+!         z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
+!         z2 = 1. - z1
+!         z4 = 1. - z3
+!         z5 = 1. - z6 - z7
+!         !
+!         ! Remove corners
+!         imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2
+!         IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3
+!         IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2
+!         IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3
+!         IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2
+!
+!         ! smoothed fields
+!         IF( eastern_side ) THEN
+!            ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:)
+!            DO jj = jmin, jmax
+!               rswitch = 0.
+!               IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1.
+!               ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  &
+!                  &                +      umask(nlci-2,jj,1)   *  &
+!                  &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  &
+!                  &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) )
+!               ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1)
+!            END DO
+!         ENDIF
+!         ! 
+!         IF( northern_side ) THEN
+!            ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:)
+!            DO ji = imin, imax
+!               rswitch = 0.
+!               IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1.
+!               ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  &
+!                  &                +      vmask(ji,nlcj-2,1)   *  &
+!                  &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) &
+!                  &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) )
+!               ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1)
+!            END DO
+!         END IF
+!         !
+!         IF( western_side) THEN
+!            ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:)
+!            DO jj = jmin, jmax
+!               rswitch = 0.
+!               IF( u_ice(2,jj) < 0._wp ) rswitch = 1.
+!               ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  &
+!                  &           +      umask(2,jj,1)   *   &
+!                  &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) &
+!                  &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) )
+!               ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1)
+!            END DO
+!         ENDIF
+!         !
+!         IF( southern_side ) THEN
+!            ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:)
+!            DO ji = imin, imax
+!               rswitch = 0.
+!               IF( v_ice(ji,2) < 0._wp ) rswitch = 1.
+!               ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  &
+!                  &           +      vmask(ji,2,1)   *  &
+!                  &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) &
+!                  &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) )
+!               ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1)
+!            END DO
+!         END IF
+!         !
+!         ! Treatment of corners
+!         IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south
+!         IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north
+!         IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south
+!         IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north
+!
+!         ! retrieve ice tracers
+!         jm = 1
+!         DO jl = 1, jpl
+!            a_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
+!            v_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
+!            v_s  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
+!            smv_i(i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
+!            oa_i (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
+!            DO jk = 1, nlay_s
+!               e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
+!            ENDDO
+!            DO jk = 1, nlay_i
+!               e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
+!            ENDDO
+!         ENDDO
+       
+         ! integrated values
+         vt_i (i1:i2,j1:j2) = SUM( v_i(i1:i2,j1:j2,:), dim=3 )
+         vt_s (i1:i2,j1:j2) = SUM( v_s(i1:i2,j1:j2,:), dim=3 )
+         at_i (i1:i2,j1:j2) = SUM( a_i(i1:i2,j1:j2,:), dim=3 )
+         et_s(i1:i2,j1:j2)  = SUM( SUM( e_s(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
+         et_i(i1:i2,j1:j2)  = SUM( SUM( e_i(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
+
       ENDIF
       
Index: /trunk/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90	(revision 7761)
@@ -52,16 +52,8 @@
       !!----------------------------------------------------------------------
       !
-!      IF( ( MOD( kt-nit000, Agrif_irhot() * Agrif_Parent(nn_fsbc) ) /=0 ) .AND. (kt /= 0) ) THEN
-!         PRINT *, 'clem NOT udpate, kt=',kt,Agrif_NbStepint()
-!      ELSE
-!         PRINT *, 'clem     UPDATE, kt=',kt,Agrif_NbStepint()
-!      ENDIF
-
-      !! clem: I think the update should take place each time the ocean sees the surface forcings (but maybe I am wrong and we should update every rhot time steps) 
-      IF( ( MOD( kt-nit000, Agrif_irhot() * Agrif_Parent(nn_fsbc) ) /=0 ) .AND. (kt /= 0) ) RETURN ! do not update if nb of child time steps differ from time refinement
-                                                                                                   ! i.e. update only at the parent time step
-      !! clem: this condition is clearly wrong if nn_fsbc/=1 (==> Agrif_NbStepint /= (Agrif_irhot()-1) all the time)
-      !!IF( ( Agrif_NbStepint() .NE. (Agrif_irhot()-1) ) .AND. (kt /= 0) )  RETURN
-      
+      !! clem: I think the update should take place each time the ocean sees the surface forcings
+      !!       (but maybe I am wrong and we should update every rhot time steps) 
+      IF( ( MOD( (kt-nit000)/nn_fsbc + 1, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) /=0 ) .AND. (kt /= 0) ) RETURN ! do not update if nb of child time steps differ from time refinement
+                                                                                                                           ! i.e. update only at the parent time step
       Agrif_UseSpecialValueInUpdate = .TRUE.
       Agrif_SpecialValueFineGrid = -9999.
@@ -98,5 +90,5 @@
       INTEGER  :: jk, jl, jm
       !!-----------------------------------------------------------------------
-      ! clem: it is ok not to multiply by e1 e2 since we conserve tracers here (cf ce qui est fait dans opa).
+      ! it is ok not to multiply by e1*e2 since we conserve tracers here (same as in the ocean).
       IF( before ) THEN
          jm = 1
@@ -135,4 +127,11 @@
          ENDDO
 
+         ! integrated values
+         vt_i (i1:i2,j1:j2) = SUM( v_i(i1:i2,j1:j2,:), dim=3 )
+         vt_s (i1:i2,j1:j2) = SUM( v_s(i1:i2,j1:j2,:), dim=3 )
+         at_i (i1:i2,j1:j2) = SUM( a_i(i1:i2,j1:j2,:), dim=3 )
+         et_s(i1:i2,j1:j2)  = SUM( SUM( e_s(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
+         et_i(i1:i2,j1:j2)  = SUM( SUM( e_i(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
+         
       ENDIF
       !
Index: /trunk/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 7761)
@@ -37,5 +37,5 @@
       jpim1   = jpi-1 
       jpjm1   = jpj-1 
-      jpkm1   = jpk-1                                         
+      jpkm1   = MAX( 1, jpk-1 )                                         
       jpij    = jpi*jpj 
       nperio  = 0
@@ -83,4 +83,7 @@
 # if defined key_top
    CALL Agrif_InitValues_cont_top
+# endif
+# if defined key_lim3
+   CALL Agrif_InitValues_cont_lim3
 # endif
    !
@@ -290,9 +293,9 @@
          ! check if vmask agree with parent along northern and southern boundaries:
          CALL Agrif_Bc_variable(vmsk_id,calledweight=1.,procname=interpvmsk)
-	 ! check if tmask and vertical scale factors agree with parent over first two coarse grid points:
+         ! check if tmask and vertical scale factors agree with parent over first two coarse grid points:
          CALL Agrif_Bc_variable(e3t_id,calledweight=1.,procname=interpe3t)
          !
          IF (lk_mpp) CALL mpp_sum( kindic_agr )
-         IF( kindic_agr /= 0 ) THEN             		
+         IF( kindic_agr /= 0 ) THEN
             CALL ctl_stop('Child Bathymetry is not correct near boundaries.')
          ELSE
@@ -457,5 +460,5 @@
 !   CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average)
 !   CALL Agrif_Set_Updatetype(sshn_id, update = Agrif_Update_Full_Weighting)
- 
+
    !
 END SUBROUTINE agrif_declare_var
@@ -579,5 +582,12 @@
    CALL agrif_declare_var_lim3
 
-   ! Controls (clem)
+   ! Controls
+
+   ! clem: For some reason, nn_fsbc(child)/=1 does not work properly (signal is largely degraded by the agrif zoom)
+   !          the run must satisfy CFL=Uice/(dx/dt) < 0.6/nn_fsbc(child)
+   !          therefore, if nn_fsbc(child)>1 one must reduce the time-step in proportion to nn_fsbc(child), which is not acceptable
+   !       If a solution is found, the following stop could be removed
+   IF( nn_fsbc > 1 )  CALL ctl_stop('nn_fsbc(child) must be set to 1 otherwise agrif and lim3 do not work properly')
+
    ! stop if rhot * nn_fsbc(parent) /= N * nn_fsbc(child) with N being integer
    IF( MOD( Agrif_irhot() * Agrif_Parent(nn_fsbc), nn_fsbc ) /= 0 )  THEN
@@ -591,5 +601,6 @@
    ! First Interpolations (using "after" ice subtime step => lim_nbstep=1)
    !----------------------------------------------------------------------
-   lim_nbstep = 1
+!!   lim_nbstep = 1
+   lim_nbstep = ( Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) ! clem: to have calledweight=1 in interp (otherwise the western border of the zoom is wrong)
    CALL agrif_interp_lim3('U') ! interpolation of ice velocities
    CALL agrif_interp_lim3('V') ! interpolation of ice velocities
@@ -618,4 +629,8 @@
    ! 1. Declaration of the type of variable which have to be interpolated (parent=>child)
    !       agrif_declare_variable(position,1st point index,--,--,dimensions,name)
+   !           ex.:  position=> 1,1 = not-centered (in i and j)
+   !                            2,2 =     centered (    -     )
+   !                 index   => 1,1 = one ghost line
+   !                            2,2 = two ghost lines
    !-------------------------------------------------------------------------------------
    CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpl*(5+nlay_s+nlay_i)/),tra_ice_id )
@@ -637,5 +652,5 @@
    ! 4. Set update type in case 2 ways (child=>parent) (normal & tangent to the grid cell for velocities)
    !--------------------------------------------------
-   CALL Agrif_Set_Updatetype(tra_ice_id, update = AGRIF_Update_Average) ! clem je comprends pas average/copy
+   CALL Agrif_Set_Updatetype(tra_ice_id, update = AGRIF_Update_Average)
    CALL Agrif_Set_Updatetype(u_ice_id  ,update1 = Agrif_Update_Copy   , update2 = Agrif_Update_Average)
    CALL Agrif_Set_Updatetype(v_ice_id  ,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy   )
Index: /trunk/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90	(revision 7761)
@@ -540,7 +540,5 @@
       !
       ! lfax contains the set of allowed factors.
-      data (ilfax(jl),jl=1,ntest) / 16384, 8192, 4096, 2048, 1024, 512, 256,  &
-         &                            128,   64,   32,   16,    8,   4,   2  /
-      !!----------------------------------------------------------------------
+      ilfax(:) = (/(2**jl,jl=ntest,1,-1)/)
 
       ! Clear the error flag and initialise output vars
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 7761)
@@ -40,9 +40,8 @@
    REAL(wp), PUBLIC ::   rtt      = 273.16_wp        !: triple point of temperature   [Kelvin]
    REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water [Kelvin]
+   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin]
 #if defined key_lim3
-   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin]
    REAL(wp), PUBLIC ::   rt0_ice  = 273.15_wp        !: melting point of ice          [Kelvin]
 #else
-   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin]
    REAL(wp), PUBLIC ::   rt0_ice  = 273.05_wp        !: melting point of ice          [Kelvin]
 #endif
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 7761)
@@ -1330,5 +1330,5 @@
       jpi   = size(fsp,1)
       jpj   = size(fsp,2)
-      jpkm1 = size(fsp,3) - 1
+      jpkm1 = MAX( 1, size(fsp,3) - 1 )
       !
       IF (polynomial_type == 1) THEN     ! Constrained Cubic Spline
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7761)
@@ -1590,12 +1590,21 @@
          zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
          zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * zicefr(:,:)
-               CALL iom_put( 'rain'         ,   frcv(jpr_rain)%z3(:,:,1)                                                         )  ! liquid precipitation 
+         IF( iom_use('precip') )          &
+            &  CALL iom_put( 'precip'       ,   frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1)                              )  ! total  precipitation
+         IF( iom_use('rain') )            &
+            &  CALL iom_put( 'rain'         ,   frcv(jpr_rain)%z3(:,:,1)                                                         )  ! liquid precipitation 
+         IF( iom_use('rain_ao_cea') )   &
+            &  CALL iom_put( 'rain_ao_cea'  , frcv(jpr_rain)%z3(:,:,1)* p_frld(:,:) * tmask(:,:,1)      )   ! liquid precipitation 
          IF( iom_use('hflx_rain_cea') )   &
-            &  CALL iom_put( 'hflx_rain_cea',   frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:)                                            )  ! heat flux from liq. precip. 
+            CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) * tmask(:,:,1))   ! heat flux from liq. precip. 
+         IF( iom_use('hflx_prec_cea') )   &
+            CALL iom_put( 'hflx_prec_cea', ztprecip * zcptn(:,:) * tmask(:,:,1) * p_frld(:,:) )   ! heat content flux from all precip  (cell avg)
+         IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') )   &
+            ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
          IF( iom_use('evap_ao_cea'  ) )   &
-            &  CALL iom_put( 'evap_ao_cea'  ,   frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)                )  ! ice-free oce evap (cell average)
+            CALL iom_put( 'evap_ao_cea'  , ztmp * tmask(:,:,1)                  )   ! ice-free oce evap (cell average)
          IF( iom_use('hflx_evap_cea') )   &
-            &  CALL iom_put( 'hflx_evap_cea', ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) ) * zcptn(:,:) )  ! heat flux from from evap (cell average)
-      CASE( 'oce and ice' )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
+            CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * zcptn(:,:) * tmask(:,:,1) )   ! heat flux from from evap (cell average)
+      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
          zemp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
          zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * zicefr(:,:)
@@ -1660,4 +1669,6 @@
       ! runoffs and calving (put in emp_tot)
       IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+      IF( iom_use('hflx_rnf_cea') )   &
+         CALL iom_put( 'hflx_rnf_cea' , rnf(:,:) * zcptn(:,:) )
       IF( srcv(jpr_cal)%laction ) THEN 
          zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 7761)
@@ -116,11 +116,4 @@
       IF( nn_timing == 1 )   CALL timing_start('sbc_ice_lim')
 
-      ! clem: it is important to initialize agrif_lim3 variables here and not in sbc_lim_init
-# if defined key_agrif
-      IF( kt == nit000 ) THEN
-         IF( .NOT. Agrif_Root() )   CALL Agrif_InitValues_cont_lim3
-      ENDIF
-# endif
-
       !-----------------------!
       ! --- Ice time step --- !
@@ -194,5 +187,5 @@
          IF( .NOT. Agrif_Root() )     CALL agrif_interp_lim3('T')
 #endif
-         IF( ln_limthd .AND. ln_bdy )     CALL bdy_ice_lim( kt )   ! -- bdy ice thermo 
+         IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice_lim( kt )   ! -- bdy ice thermo 
          ! previous lead fraction and ice volume for flux calculations
                                       CALL sbc_lim_bef                        
@@ -246,11 +239,13 @@
                                       CALL lim_var_agg( 2 )     ! necessary calls (at least for coupling)
                                       !
-# if defined key_agrif
-!!         IF( .NOT. Agrif_Root() )     CALL Agrif_ChildGrid_To_ParentGrid()  ! clem: should be called at the update frequency only (cf agrif_lim3_update)
-# endif
+!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
+!!       moreover it should only be called at the update frequency (as in agrif_lim3_update.F90)
+!!# if defined key_agrif
+!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ChildGrid_To_ParentGrid()
+!!# endif
                                       CALL lim_sbc_flx( kt )    ! -- Update surface ocean mass, heat and salt fluxes
-# if defined key_agrif
-!!         IF( .NOT. Agrif_Root() )     CALL Agrif_ParentGrid_To_ChildGrid()  ! clem: should be called at the update frequency only (cf agrif_lim3_update)
-# endif
+!!# if defined key_agrif
+!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid()
+!!# endif
          IF( ln_limdiahsb )           CALL lim_diahsb( kt )     ! -- Diagnostics and outputs 
          !
@@ -425,5 +420,5 @@
       !
       IF( lwp .AND. ln_bdy .AND. ln_limdiachk )  &
-      &   CALL ctl_warn('online conservation check activated but it does not work with BDY')
+         &   CALL ctl_warn('online conservation check activated but it does not work with BDY')
       !
       IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice  = ', rdt_ice
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 7761)
@@ -376,5 +376,5 @@
       jpim1 = jpi-1                                            ! inner domain indices
       jpjm1 = jpj-1                                            !   "           "
-      jpkm1 = jpk-1                                            !   "           "
+      jpkm1 = MAX( 1, jpk-1 )                                  !   "           "
       jpij  = jpi*jpj                                          !  jpi x j
 
Index: /trunk/NEMOGCM/NEMO/SAS_SRC/daymod.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/SAS_SRC/daymod.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/SAS_SRC/daymod.F90	(revision 7761)
@@ -80,8 +80,8 @@
       IF( MOD( rday , 2.  ) /= 0. )   CALL ctl_stop( 'the number of second of in a day must be an even number'    )
       IF( MOD( rdt  , 2.  ) /= 0. )   CALL ctl_stop( 'the time step (in second) must be an even number'           )
-      nsecd   = NINT(rday         )
-      nsecd05 = NINT(0.5 * rday   )
-      ndt     = NINT(      rdt    )
-      ndt05   = NINT(0.5 * rdt    )
+      nsecd   = NINT(rday       )
+      nsecd05 = NINT(0.5 * rday )
+      ndt     = NINT(      rdt  )
+      ndt05   = NINT(0.5 * rdt  )
 
       ! ==> clem: here we read the ocean restart for the date (only if it exists)
Index: /trunk/NEMOGCM/NEMO/SAS_SRC/diawri.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/SAS_SRC/diawri.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/SAS_SRC/diawri.F90	(revision 7761)
@@ -40,4 +40,6 @@
 #if defined key_lim2
    USE limwri_2 
+#elif defined key_lim3
+   USE limwri
 #endif
    USE lib_mpp         ! MPP library
@@ -80,5 +82,5 @@
    !!   Default option                                   NetCDF output file
    !!----------------------------------------------------------------------
-#if defined key_iomput
+# if defined key_iomput
    !!----------------------------------------------------------------------
    !!   'key_iomput'                                        use IOM library
@@ -95,8 +97,13 @@
       !! ** Method  :  use iom_put
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !!
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
       !!----------------------------------------------------------------------
       ! 
-      !! no relevant 2D arrays to write in iomput case
+      ! Output the initial state and forcings
+      IF( ninist == 1 ) THEN
+         CALL dia_wri_state( 'output.init', kt )
+         ninist = 0
+      ENDIF
       !
    END SUBROUTINE dia_wri
@@ -392,4 +399,6 @@
 #if defined key_lim2
       CALL lim_wri_state_2( kt, id_i, nh_i )
+#elif defined key_lim3
+      CALL lim_wri_state( kt, id_i, nh_i )
 #else
       CALL histend( id_i, snc4chunks=snc4set )
Index: /trunk/NEMOGCM/NEMO/SAS_SRC/nemogcm.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/SAS_SRC/nemogcm.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/SAS_SRC/nemogcm.F90	(revision 7761)
@@ -308,5 +308,5 @@
       jpim1 = jpi-1                                            ! inner domain indices
       jpjm1 = jpj-1                                            !   "           "
-      jpkm1 = jpk-1                                            !   "           "
+      jpkm1 = MAX( 1, jpk-1 )                                  !   "           "
       jpij  = jpi*jpj                                          !  jpi x j
 
@@ -370,5 +370,5 @@
       !           the environment of ocean BDY. Therefore bdy is called in both OPA and SAS modules. 
       !           This is not clean and should be changed in the future. 
-                            CALL     bdy_init
+                            CALL bdy_init
       ! ==>
                             CALL icb_init( rdt, nit000)   ! initialise icebergs instance
Index: /trunk/NEMOGCM/NEMO/SAS_SRC/sbcssm.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/SAS_SRC/sbcssm.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/SAS_SRC/sbcssm.F90	(revision 7761)
@@ -88,11 +88,19 @@
          ! 
          IF( ln_3d_uve ) THEN
-            IF( .NOT. ln_linssh )   e3t_m(:,:) = sf_ssm_3d(jf_e3t)%fnow(:,:,1) * tmask(:,:,1)    ! v-velocity 
-                                    ssu_m(:,:) = sf_ssm_3d(jf_usp)%fnow(:,:,1) * umask(:,:,1)    ! u-velocity
-                                    ssv_m(:,:) = sf_ssm_3d(jf_vsp)%fnow(:,:,1) * vmask(:,:,1)    ! v-velocity 
+            IF( .NOT. ln_linssh ) THEN
+               e3t_m(:,:) = sf_ssm_3d(jf_e3t)%fnow(:,:,1) * tmask(:,:,1) ! vertical scale factor 
+            ELSE
+               e3t_m(:,:) = e3t_0(:,:,1)                                 ! vertical scale factor
+            ENDIF
+            ssu_m(:,:) = sf_ssm_3d(jf_usp)%fnow(:,:,1) * umask(:,:,1)    ! u-velocity
+            ssv_m(:,:) = sf_ssm_3d(jf_vsp)%fnow(:,:,1) * vmask(:,:,1)    ! v-velocity 
          ELSE
-            IF( .NOT. ln_linssh )   e3t_m(:,:) = sf_ssm_2d(jf_e3t)%fnow(:,:,1) * tmask(:,:,1)    ! v-velocity 
-                                    ssu_m(:,:) = sf_ssm_2d(jf_usp)%fnow(:,:,1) * umask(:,:,1)    ! u-velocity
-                                    ssv_m(:,:) = sf_ssm_2d(jf_vsp)%fnow(:,:,1) * vmask(:,:,1)    ! v-velocity 
+            IF( .NOT. ln_linssh ) THEN
+               e3t_m(:,:) = sf_ssm_2d(jf_e3t)%fnow(:,:,1) * tmask(:,:,1) ! vertical scale factor 
+            ELSE
+               e3t_m(:,:) = e3t_0(:,:,1)                                 ! vertical scale factor
+            ENDIF
+            ssu_m(:,:) = sf_ssm_2d(jf_usp)%fnow(:,:,1) * umask(:,:,1)    ! u-velocity
+            ssv_m(:,:) = sf_ssm_2d(jf_vsp)%fnow(:,:,1) * vmask(:,:,1)    ! v-velocity 
          ENDIF
          !
@@ -111,7 +119,7 @@
          ssv_m(:,:) = 0._wp
          ssh_m(:,:) = 0._wp
-         e3t_m(:,:) = e3t_0(:,:,1) !clem: necessary at least for sas2D
-         frq_m(:,:) = 1._wp        !              - -
-         sshn (:,:) = 0._wp        !              - -
+         IF( .NOT. ln_linssh ) e3t_m(:,:) = e3t_0(:,:,1) !clem: necessary at least for sas2D
+         frq_m(:,:) = 1._wp                              !              - -
+         sshn (:,:) = 0._wp                              !              - -
       ENDIF
       
@@ -173,5 +181,5 @@
       NAMELIST/namsbc_sas/l_sasread, cn_dir, ln_3d_uve, ln_read_frq, sn_tem, sn_sal, sn_usp, sn_vsp, sn_ssh, sn_e3t, sn_frq
       !!----------------------------------------------------------------------
-      
+
       IF( ln_rstart .AND. nn_components == jp_iam_sas ) RETURN
       
@@ -306,5 +314,4 @@
  
       CALL sbc_ssm( nit000 )   ! need to define ss?_m arrays used in limistate
-      IF( .NOT. ln_read_frq )   frq_m(:,:) = 1.
       l_initdone = .TRUE.
       !
Index: /trunk/NEMOGCM/NEMO/SAS_SRC/step.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/SAS_SRC/step.F90	(revision 7760)
+++ /trunk/NEMOGCM/NEMO/SAS_SRC/step.F90	(revision 7761)
@@ -53,5 +53,5 @@
 
 #if defined key_agrif
-   SUBROUTINE stp( )
+   RECURSIVE SUBROUTINE stp( )
       INTEGER             ::   kstp   ! ocean time-step index
 #else
