Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90	(revision 4490)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90	(revision 4491)
@@ -32,4 +32,5 @@
    LOGICAL , PUBLIC :: spongedoneT = .FALSE.   !: tracer   sponge layer indicator
    LOGICAL , PUBLIC :: spongedoneU = .FALSE.   !: dynamics sponge layer indicator
+   LOGICAL , PUBLIC :: lk_agrif_fstep = .TRUE. !: if true: first step
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE,  DIMENSION(:,:) ::   spe1ur , spe2vr , spbtr2   !: ???
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 4490)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 4491)
@@ -145,4 +145,21 @@
          END DO
       ELSE
+         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
+                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
+                              & + atfp * ( tabres(ji,jj,jk,jn) &
+                              &             - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
+                        ENDIF
+                     ENDDO
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
          DO jn = n1,n2
             DO jk=k1,k2
@@ -150,5 +167,5 @@
                   DO ji=i1,i2
                      IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 
-                         tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) * tmask(ji,jj,jk)
+                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) * tmask(ji,jj,jk)
                      END IF
                   END DO
@@ -179,5 +196,5 @@
                DO ji=i1,i2
                   tabres(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk)
-                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * fse3u(ji,jj,jk)
+                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * fse3u_n(ji,jj,jk)
                END DO
             END DO
@@ -188,7 +205,12 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  un(ji,jj,jk) = tabres(ji,jj,jk) / (e2u(ji,jj))
-                  un(ji,jj,jk) = un(ji,jj,jk) * umask(ji,jj,jk)
-                  un(ji,jj,jk) = un(ji,jj,jk) / fse3u(ji,jj,jk)
+                  tabres(ji,jj,jk) = tabres(ji,jj,jk) / e2u(ji,jj) / fse3u_n(ji,jj,jk)
+                  !
+                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
+                     ub(ji,jj,jk) = ub(ji,jj,jk) & 
+                       & + atfp * ( tabres(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
+                  ENDIF
+                  !
+                  un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk)
                END DO
             END DO
@@ -217,5 +239,5 @@
                DO ji=i1,i2
                   tabres(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk)
-                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * fse3v(ji,jj,jk)
+                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * fse3v_n(ji,jj,jk)
                END DO
             END DO
@@ -226,7 +248,12 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  vn(ji,jj,jk) = tabres(ji,jj,jk) / (e1v(ji,jj))
-                  vn(ji,jj,jk) = vn(ji,jj,jk) * vmask(ji,jj,jk)
-                  vn(ji,jj,jk) = vn(ji,jj,jk) / fse3v(ji,jj,jk)
+                  tabres(ji,jj,jk) = tabres(ji,jj,jk) / e1v(ji,jj) / fse3v_n(ji,jj,jk)
+                  !
+                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
+                     vb(ji,jj,jk) = vb(ji,jj,jk) & 
+                       & + atfp * ( tabres(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
+                  ENDIF
+                  !
+                  vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk)
                END DO
             END DO
@@ -248,5 +275,5 @@
       INTEGER :: ji, jj, jk
       REAL(wp) :: zrhoy
-      REAL(wp) :: zhinv
+      REAL(wp) :: zcorr
 
       IF (before) THEN
@@ -261,18 +288,39 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               IF(umask(ji,jj,1) .NE. 0.) THEN             
-                  spgu(ji,jj) = 0.e0
-                  DO jk=1,jpk
-                     spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk)
-                  END DO
-                  spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
-                  zhinv = (tabres(ji,jj)-spgu(ji,jj))/(hu(ji,jj)*e2u(ji,jj))
-                  Do jk=1,jpk              
-                     un(ji,jj,jk) = un(ji,jj,jk) + zhinv
-                     un(ji,jj,jk) = un(ji,jj,jk) * umask(ji,jj,jk)            
-                  END DO
-               ENDIF
+               tabres(ji,jj) =  tabres(ji,jj) * hur(ji,jj) / e2u(ji,jj)  
+               !    
+               ! Update "now" 3d velocities:
+               spgu(ji,jj) = 0.e0
+               DO jk=1,jpkm1
+                  spgu(ji,jj) = spgu(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk)
+               END DO
+               spgu(ji,jj) = spgu(ji,jj) * hur(ji,jj)
+               !
+               zcorr = tabres(ji,jj) - spgu(ji,jj)
+               DO jk=1,jpkm1              
+                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
+               END DO
+               !
                ! Update barotropic velocities:
-               un_b(ji,jj) = tabres(ji,jj) * hur(ji,jj) / e2u(ji,jj)
+#if defined key_dynspg_ts
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
+                  zcorr = tabres(ji,jj) - un_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) * 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) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)
+               END DO
+               spgu(ji,jj) = spgu(ji,jj) * hur_b(ji,jj)
+               !
+               zcorr = ub_b(ji,jj) - spgu(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
@@ -292,5 +340,5 @@
       INTEGER :: ji, jj, jk
       REAL(wp) :: zrhox
-      REAL(wp) :: zhinv
+      REAL(wp) :: zcorr
 
       IF (before) THEN
@@ -305,18 +353,39 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               IF(vmask(ji,jj,1) .NE. 0.) THEN             
-                  spgv(ji,jj) = 0.
-                  DO jk=1,jpk
-                     spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk)
-                  END DO
-                  spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
-                  zhinv = (tabres(ji,jj)-spgv(ji,jj))/(hv(ji,jj)*e1v(ji,jj))
-                  DO jk=1,jpk             
-                     vn(ji,jj,jk) = vn(ji,jj,jk) + zhinv
-                     vn(ji,jj,jk) = vn(ji,jj,jk) * vmask(ji,jj,jk)
-                  END DO
-               ENDIF
+               tabres(ji,jj) =  tabres(ji,jj) * hvr(ji,jj) / e1v(ji,jj)  
+               !    
+               ! Update "now" 3d velocities:
+               spgv(ji,jj) = 0.e0
+               DO jk=1,jpkm1
+                  spgv(ji,jj) = spgv(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk)
+               END DO
+               spgv(ji,jj) = spgv(ji,jj) * hvr(ji,jj)
+               !
+               zcorr = tabres(ji,jj) - spgv(ji,jj)
+               DO jk=1,jpkm1              
+                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
+               END DO
+               !
                ! Update barotropic velocities:
-               vn_b(ji,jj) = tabres(ji,jj) * hvr(ji,jj) / e1v(ji,jj)
+#if defined key_dynspg_ts
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
+                  zcorr = tabres(ji,jj) - vn_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) * 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) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)
+               END DO
+               spgv(ji,jj) = spgv(ji,jj) * hvr_b(ji,jj)
+               !
+               zcorr = vb_b(ji,jj) - spgv(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
@@ -344,4 +413,15 @@
          END DO
       ELSE
+
+#if ! defined key_dynspg_ts
+         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
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90	(revision 4490)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90	(revision 4491)
@@ -78,4 +78,21 @@
             ENDDO
          ELSE
+            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
+                              trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn) & 
+                                 & + atfp * ( tabres(ji,jj,jk,jn) &
+                                               - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
+                           ENDIF
+                        ENDDO
+                     ENDDO
+                  ENDDO
+               ENDDO
+            ENDIF
+
             DO jn = n1, n2
                DO jk = k1, k2
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4490)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4491)
@@ -79,4 +79,5 @@
 !      IF ( Agrif_Root() .and. lwp) Write(*,*) '---'
 !      IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp
+      IF ( kstp == (nit000 + 1) ) lk_agrif_fstep = .FALSE.
 # if defined key_iomput
       IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( "nemo" )
