Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4291)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4292)
@@ -95,16 +95,50 @@
       ! Update data, open boundaries, surface boundary condition (including sea-ice)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( lk_tide    )   CALL sbc_tide( kstp )
+      IF( lk_obc     )   CALL obc_dta ( kstp )        ! update dynamic and tracer data at open boundaries
+      IF( lk_obc     )   CALL obc_rad ( kstp )        ! compute phase velocities at open boundaries
+      IF( lk_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
+
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
-
-      IF( lk_tide.AND.(kstp /= nit000 ))   CALL tide_init ( kstp )
-      IF( lk_tide    )   CALL sbc_tide( kstp )
-      IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries
-      IF( lk_obc     )   CALL obc_rad( kstp )         ! compute phase velocities at open boundaries
-      IF( lk_bdy     )   CALL bdy_dta( kstp, time_offset=+1 ) ! update dynamic and tracer data at open boundaries
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      !  Ocean dynamics : ssh, wn, hdiv, rot                                 !
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL ssh_wzv( kstp )         ! after ssh & vertical velocity
+                                                      ! clem: moved here for bdy ice purpose
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  Ocean dynamics : hdiv, rot, ssh, e3, wn
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL zdf_bfr( kstp )         ! bottom friction (if quadratic)
+                         CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur)
+      IF( lk_dynspg_ts ) THEN
+                                  CALL wzv           ( kstp )  ! now cross-level velocity 
+          ! In case the time splitting case, update almost all momentum trends here:
+          ! Note that the computation of vertical velocity above, hence "after" sea level
+          ! is necessary to compute momentum advection for the rhs of barotropic loop:
+                                  CALL eos    ( tsn, rhd, rhop )                 ! now in situ density for hpg computation
+          IF( ln_zps      )       CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &  ! zps: now hor. derivative
+                &                                          rhd, gru , grv  )     ! of t, s, rd at the last ocean level
+
+                                  ua(:,:,:) = 0.e0             ! set dynamics trends to zero
+                                  va(:,:,:) = 0.e0
+          IF(  ln_asmiau .AND. &
+             & ln_dyninc       )  CALL dyn_asm_inc  ( kstp )   ! apply dynamics assimilation increment
+          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! subtract Neptune velocities (simplified)
+          IF( lk_bdy           )  CALL bdy_dyn3d_dmp( kstp )   ! bdy damping trends
+                                  CALL dyn_adv      ( kstp )   ! advection (vector or flux form)
+                                  CALL dyn_vor      ( kstp )   ! vorticity term including Coriolis
+                                  CALL dyn_ldf      ( kstp )   ! lateral mixing
+          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified)
+#if defined key_agrif
+          IF(.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn        ! momentum sponge
+#endif
+                                  CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
+                                  CALL dyn_spg( kstp, indic )  ! surface pressure gradient
+
+                                  hdivb(:,:,:) = hdivn(:,:,:)  ! Store now divergence and rot temporarly, revert to these below 
+                                  rotb(:,:,:)  = rotn(:,:,:)     
+                                  ua_sv(:,:,:) = ua(:,:,:)     ! Save trends (barotropic trend has been fully updated)
+                                  va_sv(:,:,:) = va(:,:,:)
+
+                                  CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case)
+      ENDIF
+      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors
+                         CALL wzv           ( kstp )  ! now cross-level velocity (original)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -115,6 +149,4 @@
       !
       !  VERTICAL PHYSICS
-                         CALL zdf_bfr( kstp )         ! bottom friction
-
       !                                               ! Vertical eddy viscosity and diffusivity coefficients
       IF( lk_zdfric  )   CALL zdf_ric( kstp )            ! Richardson number dependent Kz
@@ -122,5 +154,5 @@
       IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
       IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz
-      IF( lk_zdfcst  )   THEN                            ! Constant Kz (reset avt, avm[uv] to the background value)
+      IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value)
          avt (:,:,:) = rn_avt0 * tmask(:,:,:)
          avmu(:,:,:) = rn_avm0 * umask(:,:,:)
@@ -146,5 +178,5 @@
       !
       IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
-                         CALL eos( tsb, rhd )                ! before in situ density
+                         CALL eos( tsb, rhd, gdept_0(:,:,:) )             ! before in situ density
          IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
             &                                      rhd, gru , grv  )      ! of t, s, rd at the last ocean level
@@ -193,7 +225,4 @@
                              tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero
 
-!write(numout,*) "MAV kt",kstp
-!write(numout,'(a5,3(1x,f21.18))') "INIn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(24,11)
-!write(numout,'(a5,3(1x,f21.18))') "INIa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)
       IF(  ln_asmiau .AND. &
          & ln_trainc     )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
@@ -205,36 +234,26 @@
       IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends
                              CALL tra_adv    ( kstp )       ! horizontal & vertical advection
-!write(numout,'(a5,3(1x,f21.18))') "ADVn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(24,11)
-!write(numout,'(a5,3(1x,f21.18))') "ADVa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)
       IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
                              CALL tra_ldf    ( kstp )       ! lateral mixing
-!write(numout,'(a5,3(1x,f21.18))') "LDFn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(24,11)
-!write(numout,'(a5,3(1x,f21.18))') "LDFa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)
 #if defined key_agrif
       IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge
 #endif
                              CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields
-!do jk=1,jpk
-!write(numout,'(a5,3(1x,f21.18))') "ZDFn:",tsn(5,10,jk,jp_tem),tsn(5,10,jk,jp_sal),tmask(5,10,jk)
-!write(numout,'(a5,3(1x,f21.18))') "ZDFa:",tsa(5,10,jk,jp_tem),tsa(5,10,jk,jp_sal),ssha(5,10)
-!end do
 
       IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg (time stepping then eos)
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
-                             CALL eos    ( tsa, rhd, rhop )      ! Time-filtered in situ density for hpg computation
+                             CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation
          IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &    ! zps: time filtered hor. derivative
             &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
 
       ELSE                                                  ! centered hpg  (eos then time stepping)
-                             CALL eos    ( tsn, rhd, rhop )      ! now in situ density for hpg computation
-         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! zps: now hor. derivative
+         IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case
+                                CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
+            IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! zps: now hor. derivative
             &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
-!write(numout,'(a5,3(1x,f21.18))') "ZPSn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(24,11)
-!write(numout,'(a5,3(1x,f21.18))') "ZPSa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)
+         ENDIF
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
-!write(numout,'(a5,3(1x,f21.18))') "NXTn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(25,11)
-!write(numout,'(a5,3(1x,f21.18))') "NXTa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)
       ENDIF
 
@@ -242,18 +261,31 @@
       ! Dynamics                                    (tsa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( lk_dynspg_ts   )  THEN
+                                                             ! revert to previously computed momentum tendencies
+                                                             ! (not using ua, va as temporary arrays during tracers' update could avoid that)
+                               ua(:,:,:) = ua_sv(:,:,:)
+                               va(:,:,:) = va_sv(:,:,:)
+                                                             ! Revert now divergence and rotational to previously computed ones 
+                                                             !(needed because of the time swap in div_cur, at the beginning of each time step)
+                               hdivn(:,:,:) = hdivb(:,:,:)
+                               rotn(:,:,:)  = rotb(:,:,:) 
+
+                               CALL dyn_bfr( kstp )         ! bottom friction
+                               CALL dyn_zdf( kstp )         ! vertical diffusion
+      ELSE
                                ua(:,:,:) = 0.e0             ! set dynamics trends to zero
                                va(:,:,:) = 0.e0
 
-      IF(  ln_asmiau .AND. &
-         & ln_dyninc       )   CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
-      IF( ln_bkgwri )          CALL asm_bkg_wri( kstp )     ! output background fields
-      IF( ln_neptsimp )        CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified)
-      IF( lk_bdy           )   CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends
+        IF(  ln_asmiau .AND. &
+           & ln_dyninc      )  CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
+        IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields
+        IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified)
+        IF( lk_bdy          )  CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends
                                CALL dyn_adv( kstp )         ! advection (vector or flux form)
                                CALL dyn_vor( kstp )         ! vorticity term including Coriolis
                                CALL dyn_ldf( kstp )         ! lateral mixing
-      IF( ln_neptsimp )        CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified)
-#if defined key_agrif
-      IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn        ! momemtum sponge
+        IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified)
+#if defined key_agrif
+        IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn        ! momemtum sponge
 #endif
                                CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
@@ -261,7 +293,9 @@
                                CALL dyn_zdf( kstp )         ! vertical diffusion
                                CALL dyn_spg( kstp, indic )  ! surface pressure gradient
+      ENDIF
                                CALL dyn_nxt( kstp )         ! lateral velocity at next time step
 
-                               CALL ssh_nxt( kstp )         ! sea surface height at next time step
+                               CALL ssh_swp( kstp )         ! swap of sea surface height
+      IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
 
       IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
