Index: branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 6019)
+++ branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 6020)
@@ -19,4 +19,5 @@
    !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
    !!            3.7  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
+   !!            3.7  !  2015-11  (J. Chanut) Free surface simplification
    !!-------------------------------------------------------------------------
   
@@ -28,5 +29,4 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE phycst          ! physical constants
-   USE dynspg_oce      ! type of surface pressure gradient
    USE dynadv          ! dynamics: vector invariant versus flux form
    USE domvvl          ! variable volume
@@ -68,16 +68,11 @@
       !!                  ***  ROUTINE dyn_nxt  ***
       !!                   
-      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary 
-      !!             condition on the after velocity, achieved the time stepping 
+      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary 
+      !!             condition on the after velocity, achieve the time stepping 
       !!             by applying the Asselin filter on now fields and swapping 
       !!             the fields.
       !!
-      !! ** Method  : * After velocity is compute using a leap-frog scheme:
-      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
-      !!             Note that with flux form advection and variable volume layer
-      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
-      !!             velocity.
-      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
-      !!             the time stepping has already been done in dynspg module
+      !! ** Method  : * Ensure after velocities transport matches time splitting
+      !!             estimate (ln_dynspg_ts=T)
       !!
       !!              * Apply lateral boundary conditions on after velocity 
@@ -101,8 +96,5 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   iku, ikv     ! local integers
-#if ! defined key_dynspg_flt
-      REAL(wp) ::   z2dt         ! temporary scalar
-#endif
-      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec      ! local scalars
+      REAL(wp) ::   zue3a, zue3n, zue3b, zuf           ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
       REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
@@ -112,6 +104,7 @@
       IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
       !
-      CALL wrk_alloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
-      IF( lk_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve )
+      IF( ln_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve )
+      IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
+      IF( l_trddyn     )   CALL wrk_alloc( jpi,jpj,jpk, zua, zva )
       !
       IF( kt == nit000 ) THEN
@@ -121,86 +114,48 @@
       ENDIF
 
-#if defined key_dynspg_flt
-      !
-      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
-      ! -------------
-
-      ! Update after velocity on domain lateral boundaries      (only local domain required)
-      ! --------------------------------------------------
-      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
-      CALL lbc_lnk( va, 'V', -1. ) 
-      !
-#else
-
-# if defined key_dynspg_exp
-      ! Next velocity :   Leap-frog time stepping
-      ! -------------
-      z2dt = 2. * rdt                                 ! Euler or leap-frog time step 
-      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
-      !
-      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
+      IF ( ln_dynspg_ts ) THEN
+         ! Ensure below that barotropic velocities match time splitting estimate
+         ! Compute actual transport and replace it with ts estimate at "after" time step
+         zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
+         zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
+         DO jk = 2, jpkm1
+            zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
+            zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
+         END DO
          DO jk = 1, jpkm1
-            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
-            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
-         END DO
-      ELSE                                            ! applied on thickness weighted velocity
-         DO jk = 1, jpkm1
-            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
-               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
-               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
-            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
-               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
-               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
-         END DO
-      ENDIF
-# endif
-
-# if defined key_dynspg_ts
-!!gm IF ( lk_dynspg_ts ) THEN ....
-      ! Ensure below that barotropic velocities match time splitting estimate
-      ! Compute actual transport and replace it with ts estimate at "after" time step
-      zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
-      zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
-      DO jk = 2, jpkm1
-         zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
-         zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
-      END DO
-      DO jk = 1, jpkm1
-         ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
-         va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
-      END DO
-
-      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
-         ! Remove advective velocity from "now velocities" 
-         ! prior to asselin filtering     
-         ! In the forward case, this is done below after asselin filtering   
-         ! so that asselin contribution is removed at the same time 
-         DO jk = 1, jpkm1
-            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
-            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
-         END DO  
-      ENDIF
-!!gm ENDIF
-# endif
+            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
+            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
+         END DO
+
+         IF ( .NOT.ln_bt_fw ) THEN
+            ! Remove advective velocity from "now velocities" 
+            ! prior to asselin filtering     
+            ! In the forward case, this is done below after asselin filtering   
+            ! so that asselin contribution is removed at the same time 
+            DO jk = 1, jpkm1
+               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
+               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
+            END DO  
+         ENDIF
+      ENDIF
 
       ! Update after velocity on domain lateral boundaries
       ! --------------------------------------------------      
-      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
-      CALL lbc_lnk( va, 'V', -1. ) 
-      !
-# if defined key_bdy
-      !                                !* BDY open boundaries
-      IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
-      IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
-
-!!$   Do we need a call to bdy_vol here??
-      !
-# endif
-      !
 # if defined key_agrif
       CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
 # endif
-#endif
-
+      !
+      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
+      CALL lbc_lnk( va, 'V', -1. ) 
+      !
+# if defined key_bdy
+      !                                !* BDY open boundaries
+      IF( lk_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt )
+      IF( lk_bdy .AND. ln_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
+
+!!$   Do we need a call to bdy_vol here??
+      !
+# endif
+      !
       IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
          z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step 
@@ -259,13 +214,25 @@
             ! (used as a now filtered scale factor until the swap)
             ! ----------------------------------------------------
-            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
+            IF (ln_dynspg_ts.AND.ln_bt_fw) THEN
                ! No asselin filtering on thicknesses if forward time splitting
-                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
+               fse3t_b(:,:,:) = fse3t_n(:,:,:)
             ELSE
                fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) )
                ! Add volume filter correction: compatibility with tracer advection scheme
                ! => time filter + conservation correction (only at the first level)
-               fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) &
-                              &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
+               IF ( nn_isf == 0) THEN   ! if no ice shelf melting
+                  fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) &
+                                 &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
+               ELSE                     ! if ice shelf melting
+                  DO jj = 1,jpj
+                     DO ji = 1,jpi
+                        jk = mikt(ji,jj)
+                        fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                       &
+                                          &                          * ( (emp_b(ji,jj)    - emp(ji,jj)   ) &
+                                          &                            - (rnf_b(ji,jj)    - rnf(ji,jj)   ) &
+                                          &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk)
+                     END DO
+                  END DO
+               END IF
             ENDIF
             !
@@ -324,5 +291,5 @@
          ENDIF
          !
-         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
+         IF (ln_dynspg_ts.AND.ln_bt_fw) THEN
             ! Revert "before" velocities to time split estimate
             ! Doing it here also means that asselin filter contribution is removed  
@@ -388,7 +355,8 @@
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
-      ! 
-      CALL wrk_dealloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
-      IF( lk_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve )
+      !
+      IF( ln_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve )
+      IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
+      IF( l_trddyn     )   CALL wrk_dealloc( jpi,jpj,jpk, zua, zva )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
