Index: trunk/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5656)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5836)
@@ -24,6 +24,8 @@
    !!             -   !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
    !!            3.4  !  2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal
-   !!                 !  2012-07  (J. Simeon, G. Madec, C. Ethe) Online coarsening of outputs
-   !!            3.7  !  2014-04  (F. Roquet, G. Madec) New equations of state
+   !!            3.6  !  2012-07  (J. Simeon, G. Madec. C. Ethe)  Online coarsening of outputs
+   !!            3.6  !  2014-04  (F. Roquet, G. Madec) New equations of state
+   !!            3.7  !  2014-10  (G. Madec)  LDF simplication 
+   !!             -   !  2014-12  (G. Madec) remove KPP scheme
    !!----------------------------------------------------------------------
 
@@ -37,5 +39,5 @@
    PRIVATE
 
-   PUBLIC   stp   ! called by opa.F90
+   PUBLIC   stp   ! called by nemogcm.F90
 
    !! * Substitutions
@@ -43,5 +45,5 @@
 !!gm   #  include "zdfddm_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -69,5 +71,5 @@
       !!              -5- Compute the momentum trends
       !!              -6- Update the horizontal velocity
-      !!              -7- Compute the diagnostics variables (rd,N2, div,cur,w)
+      !!              -7- Compute the diagnostics variables (rd,N2, hdiv,w)
       !!              -8- Outputs and diagnostics
       !!----------------------------------------------------------------------
@@ -76,43 +78,40 @@
       INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt)
       !! ---------------------------------------------------------------------
-
 #if defined key_agrif
       kstp = nit000 + Agrif_Nb_Step()
-      IF ( lk_agrif_debug ) THEN
-         IF ( Agrif_Root() .and. lwp) Write(*,*) '---'
-         IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp, 'int tstep',Agrif_NbStepint()
-      ENDIF
-
-      IF ( kstp == (nit000 + 1) ) lk_agrif_fstep = .FALSE.
-
+      IF( lk_agrif_debug ) THEN
+         IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---'
+         IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
+      ENDIF
+      IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
 # if defined key_iomput
       IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
 # endif
 #endif
-                             indic = 0           ! reset to no error condition
-      IF( kstp == nit000 ) THEN
-         ! must be done after nemo_init for AGRIF+XIOS+OASIS
-                      CALL iom_init(      cxios_context          )  ! iom_put initialization
-         IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" )  ! initialize context for coarse grid
-      ENDIF
-
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! update I/O and calendar 
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                             indic = 0                ! reset to no error condition
+                             
+      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
+                             CALL iom_init(      cxios_context          )  ! for model grid (including passible AGRIF zoom)
+         IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
+      ENDIF
       IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
-                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell iom we are at time step kstp
-      IF( ln_crs     )       CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell iom we are at time step kstp
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Update data, open boundaries, surface boundary condition (including sea-ice)
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_tide    )   CALL sbc_tide( kstp )
-      IF( lk_bdy     )  THEN
-         IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib 
-                         CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
-      ENDIF
-                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
-                                                      ! clem: moved here for bdy ice purpose
+                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
+      IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( lk_tide    )   CALL sbc_tide( kstp )                   ! update tide potential
+      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                   ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 
+      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)
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update stochastic parameters and random T/S fluctuations
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-                        CALL sto_par( kstp )          ! Stochastic parameters
+                         CALL sto_par( kstp )          ! Stochastic parameters
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -124,4 +123,5 @@
                          CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
                          CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
+
       !
       !  VERTICAL PHYSICS
@@ -131,5 +131,4 @@
       IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz
       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)
          avt (:,:,:) = rn_avt0 * wmask (:,:,:)
@@ -137,6 +136,7 @@
          avmv(:,:,:) = rn_avm0 * wvmask(:,:,:)
       ENDIF
+
       IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
-         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
+         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
       ENDIF
       IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
@@ -144,6 +144,5 @@
       IF( lk_zdftmx  )   CALL zdf_tmx( kstp )         ! tidal vertical mixing
 
-      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
-         &               CALL zdf_ddm( kstp )         ! double diffusive mixing
+      IF( lk_zdfddm  )   CALL zdf_ddm( kstp )         ! double diffusive mixing
 
                          CALL zdf_mxl( kstp )         ! mixed layer depth
@@ -155,34 +154,32 @@
       !  LATERAL  PHYSICS
       !
-      IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
+      IF( l_ldfslp ) THEN                             ! slope of lateral mixing
+!!gm : why this here ????
          IF(ln_sto_eos ) CALL sto_pts( tsn )          ! Random T/S fluctuations
+!!gm
                          CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
+
          IF( ln_zps .AND. .NOT. ln_isfcav)                               &
             &            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
+
          IF( ln_zps .AND.       ln_isfcav)                               &
-            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
+            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi,   &    ! Partial steps for top cell (ISF)
             &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-            &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level
-         IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator
-                         CALL ldf_slp_grif( kstp )
-         ELSE
-                         CALL ldf_slp( kstp, rhd, rn2b )     ! before slope for Madec operator
+            &                                   grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level
+
+         IF( ln_traldf_triad ) THEN 
+                         CALL ldf_slp_triad( kstp )                       ! before slope for triad operator
+         ELSE     
+                         CALL ldf_slp     ( kstp, rhd, rn2b )             ! before slope for standard operator
          ENDIF
       ENDIF
-#if defined key_traldf_c2d
-      IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )      ! eddy induced velocity coefficient
-#endif
-#if defined key_traldf_c3d && defined key_traldf_smag
-                          CALL ldf_tra_smag( kstp )      ! eddy induced velocity coefficient
-#  endif
-#if defined key_dynldf_c3d && defined key_dynldf_smag
-                          CALL ldf_dyn_smag( kstp )      ! eddy induced velocity coefficient
-#  endif
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      !  Ocean dynamics : hdiv, rot, ssh, e3, wn
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur)
+      !                                               ! eddy diffusivity coeff. and/or eiv coeff.
+      IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp ) 
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  Ocean dynamics : hdiv, ssh, e3, wn
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor)
       IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
                          CALL wzv           ( kstp )  ! now cross-level velocity 
@@ -192,24 +189,25 @@
           ! Note that the computation of vertical velocity above, hence "after" sea level
           ! is necessary to compute momentum advection for the rhs of barotropic loop:
+!!gm : why also here ????
             IF(ln_sto_eos ) CALL sto_pts( tsn )                             ! Random T/S fluctuations
+!!gm
                             CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation
-            IF( ln_zps .AND. .NOT. ln_isfcav)                               &
-               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
-               &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level
-            IF( ln_zps .AND.       ln_isfcav)                               &
-               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
+                            
+            IF( ln_zps .AND. .NOT. ln_isfcav)   &                           ! Partial steps: bottom before horizontal gradient
+               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &  ! of t, s, rd at the last ocean level
+               &                                          rhd, gru , grv    )
+            IF( ln_zps .AND.       ln_isfcav)   &                           ! Partial steps: top & bottom before horizontal gradient
+               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,   & 
                &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-               &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
-
-                                  ua(:,:,:) = 0.e0             ! set dynamics trends to zero
-                                  va(:,:,:) = 0.e0
+               &                                               grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    )
+
+                                  ua(:,:,:) = 0._wp            ! set dynamics trends to zero
+                                  va(:,:,:) = 0._wp
           IF(  lk_asminc .AND. 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
@@ -221,5 +219,5 @@
                                   va_sv(:,:,:) = va(:,:,:)
 
-                                  CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case)
+                                  CALL div_hor( kstp )         ! Horizontal divergence  (2nd call in time-split case)
           IF( lk_vvl     )        CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
                                   CALL wzv           ( kstp )  ! now cross-level velocity 
@@ -229,13 +227,13 @@
       ! diagnostics and outputs             (ua, va, tsa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_floats  )      CALL flo_stp( kstp )         ! drifting Floats
-      IF( lk_diahth  )      CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
-      IF( .NOT. ln_cpl )    CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
-      IF( lk_diadct  )      CALL dia_dct( kstp )         ! Transports
-      IF( lk_diaar5  )      CALL dia_ar5( kstp )         ! ar5 diag
-      IF( lk_diaharm )      CALL dia_harm( kstp )        ! Tidal harmonic analysis
-                            CALL dia_wri( kstp )         ! ocean model: outputs
-      !
-      IF( ln_crs     )      CALL crs_fld( kstp )         ! ocean model: online field coarsening & output
+      IF( lk_floats  )   CALL flo_stp( kstp )         ! drifting Floats
+      IF( lk_diahth  )   CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
+      IF(.NOT.ln_cpl )   CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
+      IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports
+      IF( lk_diaar5  )   CALL dia_ar5( kstp )         ! ar5 diag
+      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
+                         CALL dia_wri( kstp )         ! ocean model: outputs
+      !
+      IF( ln_crs     )   CALL crs_fld( kstp )         ! ocean model: online field coarsening & output
 
 #if defined key_top
@@ -246,9 +244,8 @@
 #endif
 
-
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Active tracers                              (ua, va used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                             tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero
+                             tsa(:,:,:,:) = 0._wp           ! set tracer trends to zero
 
       IF(  lk_asminc .AND. ln_asmiau .AND. &
@@ -261,8 +258,9 @@
       IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends
                              CALL tra_adv    ( kstp )       ! horizontal & vertical advection
-      IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
                              CALL tra_ldf    ( kstp )       ! lateral mixing
 
+!!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
       IF( ln_diaptr      )   CALL dia_ptr                   ! Poleward adv/ldf TRansports diagnostics
+!!gm
 
 #if defined key_agrif
@@ -274,5 +272,7 @@
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
+!!gm : why again a call to sto_pts ???
             IF( ln_sto_eos ) CALL sto_pts( tsn )                 ! Random T/S fluctuations
+!!gm
                              CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation
             IF( ln_zps .AND. .NOT. ln_isfcav)                                &
@@ -280,18 +280,20 @@
                &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level
             IF( ln_zps .AND.       ln_isfcav)                                &
-               &             CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
+               &             CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv, gtui, gtvi,  &    ! Partial steps for top/bottom cells
                &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
+               &                                                grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    )
       ELSE                                                  ! centered hpg  (eos then time stepping)
          IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case
+!!gm : why again a call to sto_pts ???
             IF( ln_sto_eos ) CALL sto_pts( tsn )    ! Random T/S fluctuations
+!!gm
                              CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
          IF( ln_zps .AND. .NOT. ln_isfcav)                                   &
-               &             CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
-               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level
+               &             CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: bottom before horizontal gradient
+               &                                           rhd, gru , grv    )    ! of t, s, rd at the last ocean level
          IF( ln_zps .AND.       ln_isfcav)                                   & 
-               &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
+               &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,   &    ! Partial steps for top/bottom cells
                &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
+               &                                    grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
          ENDIF
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
@@ -307,24 +309,18 @@
                                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
+                               ua(:,:,:) = 0._wp            ! set dynamics trends to zero
+                               va(:,:,:) = 0._wp
 
         IF(  lk_asminc .AND. 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
@@ -340,4 +336,7 @@
       IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
       !
+
+!!gm : This does not only concern the dynamics ==>>> add a new title
+!!gm2: why ouput restart before AGRIF update?
       IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file
 
@@ -367,5 +366,6 @@
                  CALL iom_close( numror )     ! close input  ocean restart file
          IF(lwm) CALL FLUSH    ( numond )     ! flush output namelist oce
-         IF( lwm.AND.numoni /= -1 ) CALL FLUSH    ( numoni )     ! flush output namelist ice
+         IF(lwm.AND.numoni /= -1 )   &
+            &    CALL FLUSH    ( numoni )     ! flush output namelist ice (if exist)
       ENDIF
 
@@ -373,4 +373,5 @@
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!!gm why lk_oasis and not lk_cpl ????
       IF( lk_oasis         )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
       !
@@ -383,5 +384,4 @@
       !
       IF( nn_timing == 1 .AND.  kstp == nit000  )   CALL timing_reset
-      !     
       !
    END SUBROUTINE stp
