Index: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4737)
+++ branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5837)
@@ -24,5 +24,6 @@
    !!             -   !  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
+   !!                 !  2012-07  (J. Simeon, G. Madec, C. Ethe) Online coarsening of outputs
+   !!            3.7  !  2014-04  (F. Roquet, G. Madec) New equations of state
    !!----------------------------------------------------------------------
 
@@ -31,4 +32,5 @@
    !!----------------------------------------------------------------------
    USE step_oce         ! time stepping definition modules
+   USE iom
 
    IMPLICIT NONE
@@ -39,7 +41,7 @@
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-#  include "zdfddm_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+!!gm   #  include "zdfddm_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -81,16 +83,17 @@
       IF ( kstp == (nit000 + 1) ) lk_agrif_fstep = .FALSE.
 # if defined key_iomput
-      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( "nemo" )
+      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
 # endif
 #endif
                              indic = 0           ! reset to no error condition
       IF( kstp == nit000 ) THEN
-                      CALL iom_init( "nemo" )      ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS)
-         IF( ln_crs ) CALL iom_init( "nemo_crs" )  ! initialize context for coarse grid
+         ! 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
 
       IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
-                             CALL iom_setkt( kstp - nit000 + 1, "nemo"     )   ! say to iom that we are at time step kstp
-      IF( ln_crs     )       CALL iom_setkt( kstp - nit000 + 1, "nemo_crs" )   ! say to iom that we are at time step kstp
+                             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
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -98,14 +101,23 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( lk_tide    )   CALL sbc_tide( kstp )
-      IF( lk_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
-
+      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
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update stochastic parameters and random T/S fluctuations
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+                        CALL sto_par( kstp )          ! Stochastic parameters
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Ocean physics update                (ua, va, tsa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL bn2( tsb, rn2b )        ! before Brunt-Vaisala frequency
-                         CALL bn2( tsn, rn2  )        ! now    Brunt-Vaisala frequency
+      !  THERMODYNAMICS
+                         CALL eos_rab( tsb, rab_b )       ! before local thermal/haline expension ratio at T-points
+                         CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points
+                         CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
+                         CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
       !
       !  VERTICAL PHYSICS
@@ -117,7 +129,7 @@
       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 * tmask(:,:,:)
-         avmu(:,:,:) = rn_avm0 * umask(:,:,:)
-         avmv(:,:,:) = rn_avm0 * vmask(:,:,:)
+         avt (:,:,:) = rn_avt0 * wmask (:,:,:)
+         avmu(:,:,:) = rn_avm0 * wumask(:,:,:)
+         avmv(:,:,:) = rn_avm0 * wvmask(:,:,:)
       ENDIF
       IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
@@ -140,7 +152,13 @@
       !
       IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
-                         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
+         IF(ln_sto_eos ) CALL sto_pts( tsn )          ! Random T/S fluctuations
+                         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)
+            &                                          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 )
@@ -152,8 +170,8 @@
       IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )      ! eddy induced velocity coefficient
 #endif
-#if defined key_traldf_c3d && key_traldf_smag
+#if defined key_traldf_c3d && defined key_traldf_smag
                           CALL ldf_tra_smag( kstp )      ! eddy induced velocity coefficient
 #  endif
-#if defined key_dynldf_c3d && key_dynldf_smag
+#if defined key_dynldf_c3d && defined key_dynldf_smag
                           CALL ldf_dyn_smag( kstp )      ! eddy induced velocity coefficient
 #  endif
@@ -170,7 +188,13 @@
           ! 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, 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
+            IF(ln_sto_eos ) CALL sto_pts( tsn )                             ! Random T/S fluctuations
+                            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)
+               &                                          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
@@ -201,15 +225,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( lk_diafwb  )   CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
-      IF( ln_diaptr  )   CALL dia_ptr( kstp )         ! Poleward TRansports 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
@@ -219,4 +241,5 @@
                          CALL trc_stp( kstp )         ! time-stepping
 #endif
+
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -236,4 +259,7 @@
       IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
                              CALL tra_ldf    ( kstp )       ! lateral mixing
+
+      IF( ln_diaptr      )   CALL dia_ptr                   ! Poleward adv/ldf TRansports diagnostics
+
 #if defined key_agrif
       IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge
@@ -244,13 +270,24 @@
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
+            IF( ln_sto_eos ) CALL sto_pts( tsn )                 ! Random T/S fluctuations
                              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
-
+            IF( ln_zps .AND. .NOT. ln_isfcav)                                &
+               &             CALL zps_hde    ( kstp, jpts, tsa, 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, tsa, gtsu, gtsv,  &    ! 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 last ocean level
       ELSE                                                  ! centered hpg  (eos then time stepping)
          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
+            IF( ln_sto_eos ) CALL sto_pts( tsn )    ! Random T/S fluctuations
+                             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)
+               &                                           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
          ENDIF
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
@@ -302,5 +339,4 @@
       IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
 
-      IF( lrst_oce .AND. ln_diahsb )   CALL dia_hsb_rst( kstp, 'WRITE' )
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Control and restarts
@@ -314,27 +350,17 @@
                  CALL iom_close( numror )     ! close input  ocean restart file
          IF(lwm) CALL FLUSH    ( numond )     ! flush output namelist oce
-         IF(lwm) CALL FLUSH    ( numoni )     ! flush output namelist ice    
+         IF( lwm.AND.numoni /= -1 ) CALL FLUSH    ( numoni )     ! flush output namelist ice
       ENDIF
       IF( lrst_oce         )   CALL rst_write    ( kstp )   ! write output ocean restart file
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Trends                              (ua, va, tsa used as workspace)
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( nstop == 0 ) THEN
-         IF( lk_trddyn     )   CALL trd_dwr( kstp )         ! trends: dynamics
-         IF( lk_trdtra     )   CALL trd_twr( kstp )         ! trends: active tracers
-         IF( lk_trdmld     )   CALL trd_mld( kstp )         ! trends: Mixed-layer
-         IF( lk_trdvor     )   CALL trd_vor( kstp )         ! trends: vorticity budget
-      ENDIF
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_cpl           )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
+      IF( lk_oasis         )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
       !
 #if defined key_iomput
       IF( kstp == nitend .OR. indic < 0 ) THEN 
-                      CALL iom_context_finalize( "nemo"     ) ! needed for XIOS+AGRIF
-         IF( ln_crs ) CALL iom_context_finalize( "nemo_crs" ) ! 
+                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
+         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 
       ENDIF
 #endif
