Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 6827)
@@ -21,5 +21,4 @@
    USE phycst          ! physical constants
    USE in_out_manager  ! I/O manager
-   USE sbc_oce         ! ocean surface boundary conditions
    USE lib_fortran,    ONLY: glob_sum, DDPDD
    USE lbclnk          ! lateral boundary condition - MPP exchanges
@@ -31,7 +30,4 @@
 
    PUBLIC dom_clo      ! routine called by domain module
-   PUBLIC sbc_clo      ! routine called by step module
-   PUBLIC clo_rnf      ! routine called by sbcrnf module
-   PUBLIC clo_ups      ! routine called in traadv_cen2(_jki) module
    PUBLIC clo_bat      ! routine called in domzgr module
 
@@ -185,247 +181,4 @@
 
 
-   SUBROUTINE sbc_clo( kt )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_clo  ***
-      !!                    
-      !! ** Purpose :   Special handling of closed seas
-      !!
-      !! ** Method  :   Water flux is forced to zero over closed sea
-      !!      Excess is shared between remaining ocean, or
-      !!      put as run-off in open ocean.
-      !!
-      !! ** Action  :   emp updated surface freshwater fluxes and associated heat content at kt
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean model time step
-      !
-      INTEGER             ::   ji, jj, jc, jn   ! dummy loop indices
-      REAL(wp), PARAMETER ::   rsmall = 1.e-20_wp    ! Closed sea correction epsilon
-      REAL(wp)            ::   zze2, ztmp, zcorr     ! 
-      REAL(wp)            ::   zcoef, zcoef1         ! 
-      COMPLEX(wp)         ::   ctmp 
-      REAL(wp), DIMENSION(jpncs) ::   zfwf   ! 1D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('sbc_clo')
-      !                                                   !------------------!
-      IF( kt == nit000 ) THEN                             !  Initialisation  !
-         !                                                !------------------!
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*)'sbc_clo : closed seas '
-         IF(lwp) WRITE(numout,*)'~~~~~~~'
-
-         surf(:) = 0.e0_wp
-         !
-         surf(jpncs+1) = glob_sum( e1e2t(:,:) )   ! surface of the global ocean
-         !
-         !                                        ! surface of closed seas 
-         IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation
-            DO jc = 1, jpncs
-               ctmp = CMPLX( 0.e0, 0.e0, wp )
-               DO jj = ncsj1(jc), ncsj2(jc)
-                  DO ji = ncsi1(jc), ncsi2(jc)
-                     ztmp = e1e2t(ji,jj) * tmask_i(ji,jj)
-                     CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
-                  END DO 
-               END DO 
-               IF( lk_mpp )   CALL mpp_sum( ctmp )
-               surf(jc) = REAL(ctmp,wp)
-            END DO
-         ELSE                                          ! Standard calculation           
-            DO jc = 1, jpncs
-               DO jj = ncsj1(jc), ncsj2(jc)
-                  DO ji = ncsi1(jc), ncsi2(jc)
-                     surf(jc) = surf(jc) + e1e2t(ji,jj) * tmask_i(ji,jj)      ! surface of closed seas
-                  END DO 
-               END DO 
-            END DO 
-            IF( lk_mpp )   CALL mpp_sum ( surf, jpncs )       ! mpp: sum over all the global domain
-         ENDIF
-
-         IF(lwp) WRITE(numout,*)'     Closed sea surfaces'
-         DO jc = 1, jpncs
-            IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc)
-         END DO
-
-         ! jpncs+1 : surface of sea, closed seas excluded
-         DO jc = 1, jpncs
-            surf(jpncs+1) = surf(jpncs+1) - surf(jc)
-         END DO           
-         !
-      ENDIF
-      !                                                   !--------------------!
-      !                                                   !  update emp        !
-      zfwf = 0.e0_wp                                      !--------------------!
-      IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation
-         DO jc = 1, jpncs
-            ctmp = CMPLX( 0.e0, 0.e0, wp )
-            DO jj = ncsj1(jc), ncsj2(jc)
-               DO ji = ncsi1(jc), ncsi2(jc)
-                  ztmp = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)
-                  CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
-               END DO  
-            END DO 
-            IF( lk_mpp )   CALL mpp_sum( ctmp )
-            zfwf(jc) = REAL(ctmp,wp)
-         END DO
-      ELSE                                          ! Standard calculation           
-         DO jc = 1, jpncs
-            DO jj = ncsj1(jc), ncsj2(jc)
-               DO ji = ncsi1(jc), ncsi2(jc)
-                  zfwf(jc) = zfwf(jc) + e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj) 
-               END DO  
-            END DO 
-         END DO
-         IF( lk_mpp )   CALL mpp_sum ( zfwf(:) , jpncs )       ! mpp: sum over all the global domain
-      ENDIF
-
-      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN      ! Black Sea case for ORCA_R2 configuration
-         zze2    = ( zfwf(3) + zfwf(4) ) * 0.5_wp
-         zfwf(3) = zze2
-         zfwf(4) = zze2
-      ENDIF
-
-      zcorr = 0._wp
-
-      DO jc = 1, jpncs
-         !
-         ! The following if avoids the redistribution of the round off
-         IF ( ABS(zfwf(jc) / surf(jpncs+1) ) > rsmall) THEN
-            !
-            IF( ncstt(jc) == 0 ) THEN           ! water/evap excess is shared by all open ocean
-               zcoef    = zfwf(jc) / surf(jpncs+1)
-               zcoef1   = rcp * zcoef
-               emp(:,:) = emp(:,:) + zcoef
-               qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
-               ! accumulate closed seas correction
-               zcorr    = zcorr    + zcoef
-               !
-            ELSEIF( ncstt(jc) == 1 ) THEN       ! Excess water in open sea, at outflow location, excess evap shared
-               IF ( zfwf(jc) <= 0.e0_wp ) THEN 
-                   DO jn = 1, ncsnr(jc)
-                     ji = mi0(ncsir(jc,jn))
-                     jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
-                     IF (      ji > 1 .AND. ji < jpi   &
-                         .AND. jj > 1 .AND. jj < jpj ) THEN 
-                         zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) )
-                         zcoef1     = rcp * zcoef
-                         emp(ji,jj) = emp(ji,jj) + zcoef
-                         qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)
-                     ENDIF 
-                   END DO 
-               ELSE 
-                   zcoef    = zfwf(jc) / surf(jpncs+1)
-                   zcoef1   = rcp * zcoef
-                   emp(:,:) = emp(:,:) + zcoef
-                   qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
-                   ! accumulate closed seas correction
-                   zcorr    = zcorr    + zcoef
-               ENDIF
-            ELSEIF( ncstt(jc) == 2 ) THEN       ! Excess e-p-r (either sign) goes to open ocean, at outflow location
-               DO jn = 1, ncsnr(jc)
-                  ji = mi0(ncsir(jc,jn))
-                  jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
-                  IF(      ji > 1 .AND. ji < jpi    &
-                     .AND. jj > 1 .AND. jj < jpj ) THEN 
-                     zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) *  e1e2t(ji,jj) )
-                     zcoef1     = rcp * zcoef
-                     emp(ji,jj) = emp(ji,jj) + zcoef
-                     qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)
-                  ENDIF 
-               END DO 
-            ENDIF 
-            !
-            DO jj = ncsj1(jc), ncsj2(jc)
-               DO ji = ncsi1(jc), ncsi2(jc)
-                  zcoef      = zfwf(jc) / surf(jc)
-                  zcoef1     = rcp * zcoef
-                  emp(ji,jj) = emp(ji,jj) - zcoef
-                  qns(ji,jj) = qns(ji,jj) + zcoef1 * sst_m(ji,jj)
-               END DO  
-            END DO 
-            !
-         END IF
-      END DO 
-
-      IF ( ABS(zcorr) > rsmall ) THEN      ! remove the global correction from the closed seas
-         DO jc = 1, jpncs                  ! only if it is large enough
-            DO jj = ncsj1(jc), ncsj2(jc)
-               DO ji = ncsi1(jc), ncsi2(jc)
-                  emp(ji,jj) = emp(ji,jj) - zcorr
-                  qns(ji,jj) = qns(ji,jj) + rcp * zcorr * sst_m(ji,jj)
-               END DO  
-             END DO 
-          END DO
-      ENDIF
-      !
-      emp (:,:) = emp (:,:) * tmask(:,:,1)
-      !
-      CALL lbc_lnk( emp , 'T', 1._wp )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('sbc_clo')
-      !
-   END SUBROUTINE sbc_clo
-
-
-   SUBROUTINE clo_rnf( p_rnfmsk )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_rnf  ***
-      !!                    
-      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
-      !!                to be the same as river mouth grid-points
-      !!
-      !! ** Method  :   set to 1 the runoff mask (mskrnf, see sbcrnf module)
-      !!                at the closed sea outflow grid-point.
-      !!
-      !! ** Action  :   update (p_)mskrnf (set 1 at closed sea outflow)
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array)
-      !
-      INTEGER  ::   jc, jn, ji, jj      ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      DO jc = 1, jpncs
-         IF( ncstt(jc) >= 1 ) THEN            ! runoff mask set to 1 at closed sea outflows
-             DO jn = 1, 4
-                DO jj =    mj0( ncsjr(jc,jn) ), mj1( ncsjr(jc,jn) )
-                   DO ji = mi0( ncsir(jc,jn) ), mi1( ncsir(jc,jn) )
-                      p_rnfmsk(ji,jj) = MAX( p_rnfmsk(ji,jj), 1.0_wp )
-                   END DO
-                END DO
-            END DO 
-         ENDIF 
-      END DO 
-      !
-   END SUBROUTINE clo_rnf
-
-   
-   SUBROUTINE clo_ups( p_upsmsk )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_rnf  ***
-      !!                    
-      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
-      !!                to be the same as river mouth grid-points
-      !!
-      !! ** Method  :   set to 0.5 the upstream mask (upsmsk, see traadv_cen2 
-      !!                module) over the closed seas.
-      !!
-      !! ** Action  :   update (p_)upsmsk (set 0.5 over closed seas)
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_upsmsk   ! upstream mask (upsmsk array)
-      !
-      INTEGER  ::   jc, ji, jj      ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      DO jc = 1, jpncs
-         DO jj = ncsj1(jc), ncsj2(jc)
-            DO ji = ncsi1(jc), ncsi2(jc)
-               p_upsmsk(ji,jj) = 0.5_wp         ! mixed upstream/centered scheme over closed seas
-            END DO 
-         END DO 
-       END DO 
-       !
-   END SUBROUTINE clo_ups
-   
-      
    SUBROUTINE clo_bat( pbat, kbat )
       !!---------------------------------------------------------------------
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90	(revision 6827)
@@ -33,7 +33,5 @@
    USE ioipsl  , ONLY :   ymds2ju   ! for calendar
    USE prtctl         ! Print control
-   USE trc_oce , ONLY : lk_offline ! offline flag
    USE timing         ! Timing
-   USE restart        ! restart
 
    IMPLICIT NONE
@@ -88,5 +86,4 @@
       ndt05   = NINT(0.5 * rdt  )
 
-      IF( .NOT. lk_offline ) CALL day_rst( nit000, 'READ' )
 
       ! set the calandar from ndastp (read in restart file and namelist)
@@ -281,6 +278,4 @@
       ENDIF
 
-      IF( .NOT. lk_offline ) CALL rst_opn( kt )               ! Open the restart file if needed and control lrst_oce
-      IF( lrst_oce         ) CALL day_rst( kt, 'WRITE' )      ! write day restart information
       !
       IF( nn_timing == 1 )  CALL timing_stop('day')
@@ -288,130 +283,4 @@
    END SUBROUTINE day
 
-
-   SUBROUTINE day_rst( kt, cdrw )
-      !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE ts_rst  ***
-      !!
-      !!  ** Purpose : Read or write calendar in restart file:
-      !!
-      !!  WRITE(READ) mode:
-      !!       kt        : number of time step since the begining of the experiment at the
-      !!                   end of the current(previous) run
-      !!       adatrj(0) : number of elapsed days since the begining of the experiment at the
-      !!                   end of the current(previous) run (REAL -> keep fractions of day)
-      !!       ndastp    : date at the end of the current(previous) run (coded as yyyymmdd integer)
-      !!
-      !!   According to namelist parameter nrstdt,
-      !!       nrstdt = 0  no control on the date (nit000 is  arbitrary).
-      !!       nrstdt = 1  we verify that nit000 is equal to the last
-      !!                   time step of previous run + 1.
-      !!       In both those options, the  exact duration of the experiment
-      !!       since the beginning (cumulated duration of all previous restart runs)
-      !!       is not stored in the restart and is assumed to be (nit000-1)*rdt.
-      !!       This is valid is the time step has remained constant.
-      !!
-      !!       nrstdt = 2  the duration of the experiment in days (adatrj)
-      !!                    has been stored in the restart file.
-      !!----------------------------------------------------------------------
-      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
-      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
-      !
-      REAL(wp) ::   zkt, zndastp, zdayfrac, ksecs, ktime
-      INTEGER  ::   ihour, iminute
-      !!----------------------------------------------------------------------
-
-      IF( TRIM(cdrw) == 'READ' ) THEN
-
-         IF( iom_varid( numror, 'kt', ldstop = .FALSE. ) > 0 ) THEN
-            ! Get Calendar informations
-            CALL iom_get( numror, 'kt', zkt )   ! last time-step of previous run
-            IF(lwp) THEN
-               WRITE(numout,*) ' *** Info read in restart : '
-               WRITE(numout,*) '   previous time-step                               : ', NINT( zkt )
-               WRITE(numout,*) ' *** restart option'
-               SELECT CASE ( nrstdt )
-               CASE ( 0 )   ;   WRITE(numout,*) ' nrstdt = 0 : no control of nit000'
-               CASE ( 1 )   ;   WRITE(numout,*) ' nrstdt = 1 : no control the date at nit000 (use ndate0 read in the namelist)'
-               CASE ( 2 )   ;   WRITE(numout,*) ' nrstdt = 2 : calendar parameters read in restart'
-               END SELECT
-               WRITE(numout,*)
-            ENDIF
-            ! Control of date
-            IF( nit000 - NINT( zkt ) /= 1 .AND. nrstdt /= 0 )                                         &
-                 &   CALL ctl_stop( ' ===>>>> : problem with nit000 for the restart',                 &
-                 &                  ' verify the restart file or rerun with nrstdt = 0 (namelist)' )
-            ! define ndastp and adatrj
-            IF ( nrstdt == 2 ) THEN
-               ! read the parameters corresponding to nit000 - 1 (last time step of previous run)
-               CALL iom_get( numror, 'ndastp', zndastp )
-               ndastp = NINT( zndastp )
-               CALL iom_get( numror, 'adatrj', adatrj  )
-	       CALL iom_get( numror, 'ntime', ktime )
-	       nn_time0=INT(ktime)
-               ! calculate start time in hours and minutes
-	       zdayfrac=adatrj-INT(adatrj)
-	       ksecs = NINT(zdayfrac*86400)	       ! Nearest second to catch rounding errors in adatrj	       
-	       ihour = INT(ksecs/3600)
-	       iminute = ksecs/60-ihour*60
- 	        
-               ! Add to nn_time0
-               nhour   =   nn_time0 / 100
-               nminute = ( nn_time0 - nhour * 100 )
-	       nminute=nminute+iminute
-	       
-	       IF( nminute >= 60 ) THEN
-	          nminute=nminute-60
-		  nhour=nhour+1
-	       ENDIF
-	       nhour=nhour+ihour
-	       IF( nhour >= 24 ) THEN
-		  nhour=nhour-24
-	          adatrj=adatrj+1
-	       ENDIF	       
-	       nn_time0 = nhour * 100 + nminute
-	       adatrj = INT(adatrj)                    ! adatrj set to integer as nn_time0 updated	       
-            ELSE
-               ! parameters corresponding to nit000 - 1 (as we start the step loop with a call to day)
-               ndastp = ndate0        ! ndate0 read in the namelist in dom_nam
-               nhour   =   nn_time0 / 100
-               nminute = ( nn_time0 - nhour * 100 )
-               IF( nhour*3600+nminute*60-ndt05 .lt. 0 )  ndastp=ndastp-1      ! Start hour is specified in the namelist (default 0)
-               adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday
-               ! note this is wrong if time step has changed during run
-            ENDIF
-         ELSE
-            ! parameters corresponding to nit000 - 1 (as we start the step loop with a call to day)
-            ndastp = ndate0           ! ndate0 read in the namelist in dom_nam
-            nhour   =   nn_time0 / 100
-	    nminute = ( nn_time0 - nhour * 100 )
-            IF( nhour*3600+nminute*60-ndt05 .lt. 0 )  ndastp=ndastp-1      ! Start hour is specified in the namelist (default 0)
-            adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday
-         ENDIF
-         IF( ABS(adatrj  - REAL(NINT(adatrj),wp)) < 0.1 / rday )   adatrj = REAL(NINT(adatrj),wp)   ! avoid truncation error
-         !
-         IF(lwp) THEN
-            WRITE(numout,*) ' *** Info used values : '
-            WRITE(numout,*) '   date ndastp                                      : ', ndastp
-            WRITE(numout,*) '   number of elapsed days since the begining of run : ', adatrj
-	    WRITE(numout,*) '   nn_time0                                         : ',nn_time0
-            WRITE(numout,*)
-         ENDIF
-         !
-      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
-         !
-         IF( kt == nitrst ) THEN
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'rst_write : write oce restart file  kt =', kt
-            IF(lwp) WRITE(numout,*) '~~~~~~~'
-         ENDIF
-         ! calendar control
-         CALL iom_rstput( kt, nitrst, numrow, 'kt'     , REAL( kt    , wp) )   ! time-step
-         CALL iom_rstput( kt, nitrst, numrow, 'ndastp' , REAL( ndastp, wp) )   ! date
-         CALL iom_rstput( kt, nitrst, numrow, 'adatrj' , adatrj            )   ! number of elapsed days since
-         !                                                                     ! the begining of the run [s]
-	 CALL iom_rstput( kt, nitrst, numrow, 'ntime'  , REAL( nn_time0, wp) ) ! time
-      ENDIF
-      !
-   END SUBROUTINE day_rst
 
    !!======================================================================
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 6827)
@@ -166,5 +166,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1e2f , r1_e1e2f                !: associated metrics at f-point
    !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ff                              !: coriolis factor                   [1/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ff_f, ff_t                      !: coriolis factor                   [1/s]
 
    !!----------------------------------------------------------------------
@@ -332,5 +332,5 @@
          &      e1e2v(jpi,jpj) , r1_e1e2v(jpi,jpj) , e1_e2v(jpi,jpj)                   ,     &
          &      e1e2f(jpi,jpj) , r1_e1e2f(jpi,jpj)                                     ,     &
-         &        ff (jpi,jpj)                                                         , STAT=ierr(3) )
+         &       ff_f(jpi,jpj) ,     ff_t(jpi,jpj)                                     , STAT=ierr(3) )
          !
       ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) ,     &
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6827)
@@ -24,5 +24,4 @@
    USE oce             ! ocean variables
    USE dom_oce         ! domain: ocean
-   USE sbc_oce         ! surface boundary condition: ocean
    USE phycst          ! physical constants
    USE closea          ! closed seas
@@ -33,8 +32,7 @@
    USE domwri          ! domain: write the meshmask file
    USE domvvl          ! variable volume
-   USE c1d             ! 1D vertical configuration
-   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
    !
    USE in_out_manager  ! I/O manager
+   USE iom             ! 
    USE wrk_nemo        ! Memory Allocation
    USE lib_mpp         ! distributed memory computing library
@@ -137,10 +135,5 @@
       ENDIF
       !
-      IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
-      !
-                             CALL dom_stp       ! time step
-      IF( nmsh /= 0 .AND. .NOT. ln_iscpl )                         CALL dom_wri      ! Create a domain file
-      IF( nmsh /= 0 .AND.       ln_iscpl .AND. .NOT. ln_rstart )   CALL dom_wri      ! Create a domain file
-      IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control
+      CALL cfg_write         ! create the configuration file
       !
       IF( nn_timing == 1 )   CALL timing_stop('dom_init')
@@ -465,4 +458,131 @@
    END SUBROUTINE dom_stiff
 
+
+   SUBROUTINE cfg_write
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE cfg_write  ***
+      !!                   
+      !! ** Purpose :   Create the "domain_cfg" file, a NetCDF file which 
+      !!              contains all the ocean domain informations required to 
+      !!              define an ocean configuration.
+      !!
+      !! ** Method  :   Write in a file all the arrays required to set up an
+      !!              ocean configuration.
+      !!
+      !! ** output file :   domain_cfg.nc : domain size, characteristics,
+      !horizontal mesh,
+      !!                              Coriolis parameter, depth and vertical
+      !scale factors
+      !!----------------------------------------------------------------------
+      INTEGER           ::   ji, jj, jk   ! dummy loop indices
+      INTEGER           ::   izco, izps, isco, icav
+      INTEGER           ::   inum     ! temprary units for 'domain_cfg.nc' file
+      CHARACTER(len=21) ::   clnam    ! filename (mesh and mask informations)
+      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! workspace
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'cfg_write : create the "domain_cfg.nc" file containing all required configuration information'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
+      !
+      !                       ! ============================= !
+      !                       !  create 'domain_cfg.nc' file  !
+      !                       ! ============================= !
+      !         
+      clnam = 'domain_cfg'  ! filename (configuration information)
+      CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib )
+      
+      !                             !==  global domain size  ==!
+      !
+      CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
+      CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
+      CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpk   , wp), ktype = jp_i4 )
+      !
+      !                             !==  domain characteristics  ==!
+      !
+      !                                   ! lateral boundary of the global
+      !                                   domain
+      CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
+      !
+      !                                   ! type of vertical coordinate
+      IF( ln_zco    ) THEN   ;   izco = 1   ;   ELSE   ;   izco = 0   ;   ENDIF
+      IF( ln_zps    ) THEN   ;   izps = 1   ;   ELSE   ;   izps = 0   ;   ENDIF
+      IF( ln_sco    ) THEN   ;   isco = 1   ;   ELSE   ;   isco = 0   ;   ENDIF
+      CALL iom_rstput( 0, 0, inum, 'ln_zco'   , REAL( izco, wp), ktype = jp_i4 )
+      CALL iom_rstput( 0, 0, inum, 'ln_zps'   , REAL( izps, wp), ktype = jp_i4 )
+      CALL iom_rstput( 0, 0, inum, 'ln_sco'   , REAL( isco, wp), ktype = jp_i4 )
+      !
+      !                                   ! ocean cavities under iceshelves
+      IF( ln_isfcav ) THEN   ;   icav = 1   ;   ELSE   ;   icav = 0   ;   ENDIF
+      CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
+      !
+      !                             !==  horizontal mesh  !
+      !
+      CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 )   ! latitude
+      CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
+      !                                
+      CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 )   ! longitude
+      CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
+      !                                
+      CALL iom_rstput( 0, 0, inum, 'e1t'  , e1t  , ktype = jp_r8 )   ! i-scale factors (e1.)
+      CALL iom_rstput( 0, 0, inum, 'e1u'  , e1u  , ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e1v'  , e1v  , ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e1f'  , e1f  , ktype = jp_r8 )
+      !
+      CALL iom_rstput( 0, 0, inum, 'e2t'  , e2t  , ktype = jp_r8 )   ! j-scale factors (e2.)
+      CALL iom_rstput( 0, 0, inum, 'e2u'  , e2u  , ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e2v'  , e2v  , ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e2f'  , e2f  , ktype = jp_r8 )
+      !
+      CALL iom_rstput( 0, 0, inum, 'ff_f' , ff_f , ktype = jp_r8 )   ! coriolis factor
+      CALL iom_rstput( 0, 0, inum, 'ff_t' , ff_t , ktype = jp_r8 )
+      !
+      !                             !==  vertical mesh - 3D mask  ==!
+      !                                                     
+      CALL iom_rstput( 0, 0, inum, 'gdept_1d', gdept_1d, ktype = jp_r8 )   !  reference 1D-coordinate
+      CALL iom_rstput( 0, 0, inum, 'gdepw_1d', gdepw_1d, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e3t_1d'  , e3t_1d  , ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e3w_1d'  , e3w_1d  , ktype = jp_r8 )
+      !                                                     
+      CALL iom_rstput( 0, 0, inum, 'gdept_0' , gdept_0 , ktype = jp_r8 )   !  depth (t- & w-points)
+      CALL iom_rstput( 0, 0, inum, 'gdepw_0' , gdepw_0 , ktype = jp_r8 )
+      !
+      CALL iom_rstput( 0, 0, inum, 'e3t_0'   , e3t_0   , ktype = jp_r8 )   !  vertical scale factors (e
+      CALL iom_rstput( 0, 0, inum, 'e3u_0'   , e3u_0   , ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e3v_0'   , e3v_0   , ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e3f_0'   , e3f_0   , ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e3w_0'   , e3w_0   , ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e3uw_0'  , e3uw_0  , ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e3vw_0'  , e3vw_0  , ktype = jp_r8 )
+      !                                         
+      !                             !==  ocean top and bottom level  ==!
+      !
+      CALL iom_rstput( 0, 0, inum, 'bottom level' , REAL( mbkt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points
+      CALL iom_rstput( 0, 0, inum, 'top    level' , REAL( mikt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points (ISF)
+      !
+      IF( ln_sco ) THEN             ! s-coordinate: store grid stiffness ratio (Not required anyway)
+         CALL dom_stiff
+         !SF  CALL dom_stiff( z2d ) !commented because at compilation error:  The
+                                   !number of actual arguments cannot be greater than the number of dummy
+                                   !arguments.   [DOM_STIFF]
+                                   !CALL dom_stiff( z2d )
+                                   ! --------------^ compilation aborted for
+                                   !/workgpfs/rech/gzi/rgzi011/commit-simplif2-TOOLS/NEMOGCM/CONFIG/TEST/BLD/ppsrc/nemo/domain.f90
+         CALL iom_rstput( 0, 0, inum, 'stiffness', z2d )        !    ! Max. grid stiffness ratio
+      ENDIF
+      !
+      !                                ! ============================
+      !                                !        close the files 
+      !                                ! ============================
+      CALL iom_close( inum )
+      !
+   END SUBROUTINE cfg_write
+
+
+
    !!======================================================================
 END MODULE domain
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90	(revision 6827)
@@ -16,6 +16,4 @@
    USE lib_mpp         ! distributed memory computing library
    USE timing          ! Timing
-   USE c1d             ! 1D configuration
-   USE domc1d          ! 1D configuration: column location
 
    IMPLICIT NONE
@@ -82,5 +80,4 @@
       !!----------------------------------------------------------------------
       !                              ! recalculate jpizoom/jpjzoom given lat/lon
-      IF( lk_c1d .AND. ln_c1d_locpt )  CALL dom_c1d( rn_lat1d, rn_lon1d )
       !
       !                        ! ============== !
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6827)
@@ -377,11 +377,13 @@
       CASE ( 0, 1, 4 )               ! mesh on the sphere
          !
-         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
+         ff_f(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
+         ff_t(:,:) = 2. * omega * SIN( rad * gphit(:,:) )     !    -        -       - at t-point
          !
       CASE ( 2 )                     ! f-plane at ppgphi0 
          !
-         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
-         !
-         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
+         ff_f(:,:) = 2. * omega * SIN( rad * ppgphi0 )
+         ff_t(:,:) = 2. * omega * SIN( rad * ppgphi0 )
+         !
+         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff_f(1,1)
          !
       CASE ( 3 )                     ! beta-plane
@@ -399,14 +401,15 @@
          zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
          !
-         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
+         ff_f(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
+         ff_t(:,:) = ( zf0  + zbeta * gphit(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
          !
          IF(lwp) THEN
             WRITE(numout,*) 
-            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
-            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
+            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff_f(nldi,nldj)
+            WRITE(numout,*) '          Coriolis parameter varies from ', ff_f(nldi,nldj),' to ', ff_f(nldi,nlej)
          ENDIF
          IF( lk_mpp ) THEN 
-            zminff=ff(nldi,nldj)
-            zmaxff=ff(nldi,nlej)
+            zminff=ff_f(nldi,nldj)
+            zmaxff=ff_f(nldi,nlej)
             CALL mpp_min( zminff )   ! min over the global domain
             CALL mpp_max( zmaxff )   ! max over the global domain
@@ -420,15 +423,16 @@
          zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
          !
-         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
+         ff_f(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
+         ff_t(:,:) = ( zf0 + zbeta * ABS( gphit(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
          !
          IF(lwp) THEN
             WRITE(numout,*) 
             WRITE(numout,*) '          Beta-plane and rotated domain : '
-            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
+            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff_f(nldi,nldj),' to ', ff_f(nldi,nlej)
          ENDIF
          !
          IF( lk_mpp ) THEN 
-            zminff=ff(nldi,nldj)
-            zmaxff=ff(nldi,nlej)
+            zminff=ff_f(nldi,nldj)
+            zmaxff=ff_f(nldi,nlej)
             CALL mpp_min( zminff )   ! min over the global domain
             CALL mpp_max( zmaxff )   ! max over the global domain
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 6827)
@@ -22,7 +22,4 @@
    USE phycst          ! physical constant
    USE dom_oce         ! ocean space and time domain
-   USE sbc_oce         ! ocean surface boundary condition
-   USE wet_dry         ! wetting and drying
-   USE restart         ! ocean restart
    !
    USE in_out_manager  ! I/O manager
@@ -37,7 +34,4 @@
 
    PUBLIC  dom_vvl_init       ! called by domain.F90
-   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
-   PUBLIC  dom_vvl_sf_swp     ! called by step.F90
-   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
 
    !                                                      !!* Namelist nam_vvl
@@ -133,5 +127,4 @@
       !
       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
-      CALL dom_vvl_rst( nit000, 'READ' )
       e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
       !
@@ -246,421 +239,4 @@
 
 
-   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
-      !!                   
-      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
-      !!                 tranxt and dynspg routines
-      !!
-      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
-      !!               - z_tilde_case: after scale factor increment = 
-      !!                                    high frequency part of horizontal divergence
-      !!                                  + retsoring towards the background grid
-      !!                                  + thickness difusion
-      !!                               Then repartition of ssh INCREMENT proportionnaly
-      !!                               to the "baroclinic" level thickness.
-      !!
-      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
-      !!               - tilde_e3t_a: after increment of vertical scale factor 
-      !!                              in z_tilde case
-      !!               - e3(t/u/v)_a
-      !!
-      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in )           ::   kt      ! time step
-      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
-      !
-      INTEGER                ::   ji, jj, jk            ! dummy loop indices
-      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
-      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
-      LOGICAL                ::   ll_do_bclinic         ! local logical
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze3t
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zht, z_scale, zwu, zwv, zhdiv
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
-      !
-      IF( nn_timing == 1 )   CALL timing_start('dom_vvl_sf_nxt')
-      !
-      CALL wrk_alloc( jpi,jpj,zht,   z_scale, zwu, zwv, zhdiv )
-      CALL wrk_alloc( jpi,jpj,jpk,   ze3t )
-
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
-      ENDIF
-
-      ll_do_bclinic = .TRUE.
-      IF( PRESENT(kcall) ) THEN
-         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
-      ENDIF
-
-      ! ******************************* !
-      ! After acale factors at t-points !
-      ! ******************************* !
-      !                                                ! --------------------------------------------- !
-      !                                                ! z_star coordinate and barotropic z-tilde part !
-      !                                                ! --------------------------------------------- !
-      !
-      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
-      DO jk = 1, jpkm1
-         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
-         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
-      END DO
-      !
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
-         !                                                            ! ------baroclinic part------ !
-         ! I - initialization
-         ! ==================
-
-         ! 1 - barotropic divergence
-         ! -------------------------
-         zhdiv(:,:) = 0._wp
-         zht(:,:)   = 0._wp
-         DO jk = 1, jpkm1
-            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
-            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
-
-         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
-         ! --------------------------------------------------
-         IF( ln_vvl_ztilde ) THEN
-            IF( kt > nit000 ) THEN
-               DO jk = 1, jpkm1
-                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
-                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
-               END DO
-            ENDIF
-         ENDIF
-
-         ! II - after z_tilde increments of vertical scale factors
-         ! =======================================================
-         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
-
-         ! 1 - High frequency divergence term
-         ! ----------------------------------
-         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
-            DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
-            END DO
-         ELSE                         ! layer case
-            DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
-            END DO
-         ENDIF
-
-         ! 2 - Restoring term (z-tilde case only)
-         ! ------------------
-         IF( ln_vvl_ztilde ) THEN
-            DO jk = 1, jpk
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
-            END DO
-         ENDIF
-
-         ! 3 - Thickness diffusion term
-         ! ----------------------------
-         zwu(:,:) = 0._wp
-         zwv(:,:) = 0._wp
-         DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
-                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
-                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
-                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
-                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
-                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         DO jj = 1, jpj          ! b - correction for last oceanic u-v points
-            DO ji = 1, jpi
-               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
-               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
-            END DO
-         END DO
-         DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
-                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
-                     &                                            ) * r1_e1e2t(ji,jj)
-               END DO
-            END DO
-         END DO
-         !                       ! d - thickness diffusion transport: boundary conditions
-         !                             (stored for tracer advction and continuity equation)
-         CALL lbc_lnk( un_td , 'U' , -1._wp)
-         CALL lbc_lnk( vn_td , 'V' , -1._wp)
-
-         ! 4 - Time stepping of baroclinic scale factors
-         ! ---------------------------------------------
-         ! Leapfrog time stepping
-         ! ~~~~~~~~~~~~~~~~~~~~~~
-         IF( neuler == 0 .AND. kt == nit000 ) THEN
-            z2dt =  rdt
-         ELSE
-            z2dt = 2.0_wp * rdt
-         ENDIF
-         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )
-         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
-
-         ! Maximum deformation control
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         ze3t(:,:,jpk) = 0._wp
-         DO jk = 1, jpkm1
-            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
-         END DO
-         z_tmax = MAXVAL( ze3t(:,:,:) )
-         IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain
-         z_tmin = MINVAL( ze3t(:,:,:) )
-         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
-         ! - ML - test: for the moment, stop simulation for too large e3_t variations
-         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
-            IF( lk_mpp ) THEN
-               CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
-               CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
-            ELSE
-               ijk_max = MAXLOC( ze3t(:,:,:) )
-               ijk_max(1) = ijk_max(1) + nimpp - 1
-               ijk_max(2) = ijk_max(2) + njmpp - 1
-               ijk_min = MINLOC( ze3t(:,:,:) )
-               ijk_min(1) = ijk_min(1) + nimpp - 1
-               ijk_min(2) = ijk_min(2) + njmpp - 1
-            ENDIF
-            IF (lwp) THEN
-               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
-               WRITE(numout, *) 'at i, j, k=', ijk_max
-               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
-               WRITE(numout, *) 'at i, j, k=', ijk_min            
-               CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
-            ENDIF
-         ENDIF
-         ! - ML - end test
-         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
-         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
-         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
-
-         !
-         ! "tilda" change in the after scale factor
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
-         END DO
-         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
-         ! ==================================================================================
-         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
-         ! - ML - baroclinicity error should be better treated in the future
-         !        i.e. locally and not spread over the water column.
-         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
-         zht(:,:) = 0.
-         DO jk = 1, jpkm1
-            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
-         DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
-         END DO
-
-      ENDIF
-
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
-      !                                           ! ---baroclinic part--------- !
-         DO jk = 1, jpkm1
-            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-      ENDIF
-
-      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
-         !
-         IF( lwp ) WRITE(numout, *) 'kt =', kt
-         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
-            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
-            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
-         END IF
-         !
-         zht(:,:) = 0.0_wp
-         DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
-         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax
-         !
-         zht(:,:) = 0.0_wp
-         DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
-         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax
-         !
-         zht(:,:) = 0.0_wp
-         DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
-         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
-         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
-         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
-         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
-      END IF
-
-      ! *********************************** !
-      ! After scale factors at u- v- points !
-      ! *********************************** !
-
-      CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
-      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
-
-      ! *********************************** !
-      ! After depths at u- v points         !
-      ! *********************************** !
-
-      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
-      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
-      DO jk = 2, jpkm1
-         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
-         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
-      END DO
-      !                                        ! Inverse of the local depth
-!!gm BUG ?  don't understand the use of umask_i here .....
-      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
-      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
-      !
-      CALL wrk_dealloc( jpi,jpj,       zht, z_scale, zwu, zwv, zhdiv )
-      CALL wrk_dealloc( jpi,jpj,jpk,   ze3t )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_nxt')
-      !
-   END SUBROUTINE dom_vvl_sf_nxt
-
-
-   SUBROUTINE dom_vvl_sf_swp( kt )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE dom_vvl_sf_swp  ***
-      !!                   
-      !! ** Purpose :  compute time filter and swap of scale factors 
-      !!               compute all depths and related variables for next time step
-      !!               write outputs and restart file
-      !!
-      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
-      !!               - reconstruct scale factor at other grid points (interpolate)
-      !!               - recompute depths and water height fields
-      !!
-      !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
-      !!               - Recompute:
-      !!                    e3(u/v)_b       
-      !!                    e3w_n           
-      !!                    e3(u/v)w_b      
-      !!                    e3(u/v)w_n      
-      !!                    gdept_n, gdepw_n  and gde3w_n
-      !!                    h(u/v) and h(u/v)r
-      !!
-      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
-      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! time step
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcoef        ! local scalar
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
-      !
-      IF( kt == nit000 )   THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
-      ENDIF
-      !
-      ! Time filter and swap of scale factors
-      ! =====================================
-      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-         IF( neuler == 0 .AND. kt == nit000 ) THEN
-            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
-         ELSE
-            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
-            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
-         ENDIF
-         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
-      ENDIF
-      gdept_b(:,:,:) = gdept_n(:,:,:)
-      gdepw_b(:,:,:) = gdepw_n(:,:,:)
-
-      e3t_n(:,:,:) = e3t_a(:,:,:)
-      e3u_n(:,:,:) = e3u_a(:,:,:)
-      e3v_n(:,:,:) = e3v_a(:,:,:)
-
-      ! Compute all missing vertical scale factor and depths
-      ! ====================================================
-      ! Horizontal scale factor interpolations
-      ! --------------------------------------
-      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
-      ! - JC - hu_b, hv_b, hur_b, hvr_b also
-      
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
-      
-      ! Vertical scale factor interpolations
-      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
-      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
-
-      ! t- and w- points depth (set the isf depth as it is in the initial step)
-      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-      gdepw_n(:,:,1) = 0.0_wp
-      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
-      DO jk = 2, jpk
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
-                                                                 ! 1 for jk = mikt
-               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
-                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
-               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
-            END DO
-         END DO
-      END DO
-
-      ! Local depth and Inverse of the local depth of the water
-      ! -------------------------------------------------------
-      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
-      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
-      !
-      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
-      DO jk = 2, jpkm1
-         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-      END DO
-
-      ! write restart file
-      ! ==================
-      IF( lrst_oce )   CALL dom_vvl_rst( kt, 'WRITE' )
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('dom_vvl_sf_swp')
-      !
-   END SUBROUTINE dom_vvl_sf_swp
-
-
    SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
       !!---------------------------------------------------------------------
@@ -684,9 +260,5 @@
       IF( nn_timing == 1 )   CALL timing_start('dom_vvl_interpol')
       !
-      IF(ln_wd) THEN
-        zlnwd = 1.0_wp
-      ELSE
-        zlnwd = 0.0_wp
-      END IF
+      zlnwd = 0.0_wp
       !
       SELECT CASE ( pout )    !==  type of interpolation  ==!
@@ -772,157 +344,4 @@
       !
    END SUBROUTINE dom_vvl_interpol
-
-
-   SUBROUTINE dom_vvl_rst( kt, cdrw )
-      !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE dom_vvl_rst  ***
-      !!                     
-      !! ** Purpose :   Read or write VVL file in restart file
-      !!
-      !! ** Method  :   use of IOM library
-      !!                if the restart does not contain vertical scale factors,
-      !!                they are set to the _0 values
-      !!                if the restart does not contain vertical scale factors increments (z_tilde),
-      !!                they are set to 0.
-      !!----------------------------------------------------------------------
-      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
-      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
-      !
-      INTEGER ::   ji, jj, jk
-      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_rst')
-      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
-         !                                   ! ===============
-         IF( ln_rstart ) THEN                   !* Read the restart file
-            CALL rst_read_open                  !  open the restart file if necessary
-            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    )
-            !
-            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
-            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
-            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
-            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
-            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
-            !                             ! --------- !
-            !                             ! all cases !
-            !                             ! --------- !
-            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) )
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) )
-               ! needed to restart if land processor not computed 
-               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
-               WHERE ( tmask(:,:,:) == 0.0_wp ) 
-                  e3t_n(:,:,:) = e3t_0(:,:,:)
-                  e3t_b(:,:,:) = e3t_0(:,:,:)
-               END WHERE
-               IF( neuler == 0 ) THEN
-                  e3t_b(:,:,:) = e3t_n(:,:,:)
-               ENDIF
-            ELSE IF( id1 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
-               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
-               IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) )
-               e3t_n(:,:,:) = e3t_b(:,:,:)
-               neuler = 0
-            ELSE IF( id2 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
-               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
-               IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) )
-               e3t_b(:,:,:) = e3t_n(:,:,:)
-               neuler = 0
-            ELSE
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
-               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
-               IF(lwp) write(numout,*) 'neuler is forced to 0'
-               DO jk = 1, jpk
-                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
-                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
-                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
-               END DO
-               e3t_b(:,:,:) = e3t_n(:,:,:)
-               neuler = 0
-            ENDIF
-            !                             ! ----------- !
-            IF( ln_vvl_zstar ) THEN       ! z_star case !
-               !                          ! ----------- !
-               IF( MIN( id3, id4 ) > 0 ) THEN
-                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
-               ENDIF
-               !                          ! ----------------------- !
-            ELSE                          ! z_tilde and layer cases !
-               !                          ! ----------------------- !
-               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
-                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
-                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
-               ELSE                            ! one at least array is missing
-                  tilde_e3t_b(:,:,:) = 0.0_wp
-                  tilde_e3t_n(:,:,:) = 0.0_wp
-               ENDIF
-               !                          ! ------------ !
-               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
-                  !                       ! ------------ !
-                  IF( id5 > 0 ) THEN  ! required array exists
-                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
-                  ELSE                ! array is missing
-                     hdiv_lf(:,:,:) = 0.0_wp
-                  ENDIF
-               ENDIF
-            ENDIF
-            !
-         ELSE                                   !* Initialize at "rest"
-            e3t_b(:,:,:) = e3t_0(:,:,:)
-            e3t_n(:,:,:) = e3t_0(:,:,:)
-            sshn(:,:) = 0.0_wp
-
-            IF( ln_wd ) THEN
-              DO jj = 1, jpj
-                DO ji = 1, jpi
-                  IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN
-                     e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 
-                     e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 
-                     e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 
-                     sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj)
-                     sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj)
-                     ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj)
-                  ENDIF
-                ENDDO
-              ENDDO
-            END IF
-
-            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
-               tilde_e3t_b(:,:,:) = 0.0_wp
-               tilde_e3t_n(:,:,:) = 0.0_wp
-               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
-            END IF
-         ENDIF
-         !
-      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
-         !                                   ! ===================
-         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
-         !                                           ! --------- !
-         !                                           ! all cases !
-         !                                           ! --------- !
-         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:) )
-         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) )
-         !                                           ! ----------------------- !
-         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
-            !                                        ! ----------------------- !
-            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
-            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
-         END IF
-         !                                           ! -------------!    
-         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
-            !                                        ! ------------ !
-            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
-         ENDIF
-         !
-      ENDIF
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
-      !
-   END SUBROUTINE dom_vvl_rst
 
 
@@ -990,5 +409,4 @@
       !
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
-      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
       !
       IF(lwp) THEN                   ! Print the choice
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 6827)
@@ -249,5 +249,6 @@
       CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 )
       
-      CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 )           !    ! coriolis factor
+      CALL iom_rstput( 0, 0, inum3, 'ff_f', ff_f, ktype = jp_r8 )           !    ! coriolis factor
+      CALL iom_rstput( 0, 0, inum3, 'ff_t', ff_t, ktype = jp_r8 )           !    ! coriolis factor
       
       ! note that mbkt is set to 1 over land ==> use surface tmask
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 6827)
@@ -37,7 +37,5 @@
    USE oce               ! ocean variables
    USE dom_oce           ! ocean domain
-   USE wet_dry           ! wetting and drying
    USE closea            ! closed seas
-   USE c1d               ! 1D vertical configuration
    !
    USE in_out_manager    ! I/O manager
@@ -147,5 +145,4 @@
                           CALL zgr_z            ! Reference z-coordinate system (always called)
                           CALL zgr_bat          ! Bathymetry fields (levels and meters)
-      IF( lk_c1d      )   CALL lbc_lnk( bathy , 'T', 1._wp )   ! 1D config.: same bathy value over the 3x3 domain
       IF( ln_zco      )   CALL zgr_zco          ! z-coordinate
       IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate
@@ -155,12 +152,7 @@
       ! -----------------------------------
       IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain
-      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points
+                          CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points
                           CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points
                           CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points
-      !
-      IF( lk_c1d ) THEN                         ! 1D config.: same mbathy value over the 3x3 domain
-         ibat = mbathy(2,2)
-         mbathy(:,:) = ibat
-      END IF
       !
       IF( nprint == 1 .AND. lwp )   THEN
@@ -1950,5 +1942,4 @@
       bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
 
-      IF( .NOT.ln_wd ) THEN                      
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -1956,5 +1947,4 @@
             END DO
          END DO
-      END IF
       !                                        ! =============================
       !                                        ! Define the envelop bathymetry   (hbatt)
@@ -1963,5 +1953,4 @@
       zenv(:,:) = bathy(:,:)
       !
-      IF( .NOT.ln_wd ) THEN    
       ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
          DO jj = 1, jpj
@@ -1986,5 +1975,4 @@
             END DO
          END DO
-      END IF
 
       ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
@@ -2080,10 +2068,5 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( .NOT.ln_wd ) THEN
            WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
-         ELSE
-           WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min
-           WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld
-         ENDIF
       ENDIF
       hbatu(:,:) = rn_sbot_min
@@ -2099,18 +2082,4 @@
       END DO
 
-      IF( ln_wd ) THEN               !avoid the zero depth on T- (U-,V-,F-) points
-        DO jj = 1, jpj
-          DO ji = 1, jpi
-            IF(ABS(hbatt(ji,jj)) < rn_wdmin1) &
-              & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1
-            IF(ABS(hbatu(ji,jj)) < rn_wdmin1) &
-              & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1
-            IF(ABS(hbatv(ji,jj)) < rn_wdmin1) &
-              & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1
-            IF(ABS(hbatf(ji,jj)) < rn_wdmin1) &
-              & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1
-          END DO
-        END DO
-      END IF
       ! 
       ! Apply lateral boundary condition
@@ -2146,10 +2115,8 @@
 
 !!bug:  key_helsinki a verifer
-      IF( .NOT.ln_wd ) THEN
         hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
         hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
         hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
         hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
-      END IF
 
       IF( nprint == 1 .AND. lwp )   THEN
@@ -2193,5 +2160,4 @@
       CALL lbc_lnk( e3vw_0, 'V', 1._wp )
       !
-      IF( .NOT.ln_wd ) THEN
         WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp
         WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp
@@ -2201,5 +2167,4 @@
         WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp
         WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp
-      END IF
 
 #if defined key_agrif
@@ -2234,13 +2199,4 @@
                IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
             END DO
-            IF( ln_wd ) THEN
-              IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN
-                mbathy(ji,jj) = 0
-              ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN
-                mbathy(ji,jj) = 2
-              ENDIF
-            ELSE
-              IF( scobot(ji,jj) == 0._wp )   mbathy(ji,jj) = 0
-            ENDIF
          END DO
       END DO
@@ -2425,33 +2381,18 @@
                    & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
             DO jk = 1, jpk
-               IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN
-                 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
-                 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
-               ELSE
                  z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
                         &              / ztmpu
                  z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
                         &              / ztmpu
-               END IF
-
-               IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN
-                 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
-                 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
-               ELSE
+
                  z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
                         &              / ztmpv
                  z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
                         &              / ztmpv
-               END IF
-
-               IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN
-                 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj  ,jk) + z_esigt3(ji+1,jj  ,jk)               &
-                        &                        + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
-               ELSE
+
                  z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
                         &            +   hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
                         &            +   hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
                         &            +   hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
-               END IF
 
                !
@@ -2582,33 +2523,17 @@
                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
            DO jk = 1, jpk
-              IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN
-                z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
-                z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
-              ELSE
                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
                        &              / ztmpu
                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
                        &              / ztmpu
-              END IF
-
-              IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN
-                z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
-                z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
-              ELSE
                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
                        &              / ztmpv
                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
                        &              / ztmpv
-              END IF
-
-              IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN
-                z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)                 &
-                       &                        + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
-              ELSE
+
                 z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
                        &              + hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
                        &              + hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
                        &              + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
-              END IF
 
              ! Code prior to wetting and drying option (for reference)
Index: anches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 6826)
+++ 	(revision )
@@ -1,265 +1,0 @@
-MODULE dtatsd
-   !!======================================================================
-   !!                     ***  MODULE  dtatsd  ***
-   !! Ocean data  :  read ocean Temperature & Salinity Data from gridded data
-   !!======================================================================
-   !! History :  OPA  ! 1991-03  ()  Original code
-   !!             -   ! 1992-07  (M. Imbard)
-   !!            8.0  ! 1999-10  (M.A. Foujols, M. Imbard)  NetCDF FORMAT 
-   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
-   !!            3.3  ! 2010-10  (C. Bricaud, S. Masson)  use of fldread
-   !!            3.4  ! 2010-11  (G. Madec, C. Ethe) Merge of dtatem and dtasal + suppression of CPP keys
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dta_tsd      : read and time interpolated ocean Temperature & Salinity Data
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE fldread         ! read input fields
-   USE in_out_manager  ! I/O manager
-   USE phycst          ! physical constants
-   USE lib_mpp         ! MPP library
-   USE wrk_nemo        ! Memory allocation
-   USE timing          ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   dta_tsd_init   ! called by opa.F90
-   PUBLIC   dta_tsd        ! called by istate.F90 and tradmp.90
-
-   LOGICAL , PUBLIC ::   ln_tsd_init      !: T & S data flag
-   LOGICAL , PUBLIC ::   ln_tsd_tradmp    !: internal damping toward input data flag
-
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read)
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$ 
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE dta_tsd_init( ld_tradmp )
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE dta_tsd_init  ***
-      !!                    
-      !! ** Purpose :   initialisation of T & S input data 
-      !! 
-      !! ** Method  : - Read namtsd namelist
-      !!              - allocates T & S data structure 
-      !!----------------------------------------------------------------------
-      LOGICAL, INTENT(in), OPTIONAL ::   ld_tradmp   ! force the initialization when tradp is used
-      !
-      INTEGER ::   ierr0, ierr1, ierr2, ierr3   ! temporary integers
-      !
-      CHARACTER(len=100)            ::   cn_dir          ! Root directory for location of ssr files
-      TYPE(FLD_N), DIMENSION( jpts) ::   slf_i           ! array of namelist informations on the fields to read
-      TYPE(FLD_N)                   ::   sn_tem, sn_sal
-      !!
-      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal
-      INTEGER  ::   ios
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dta_tsd_init')
-      !
-      !  Initialisation
-      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0
-      !
-      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist : 
-      READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtsd in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run
-      READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtsd in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namtsd )
-
-      IF( PRESENT( ld_tradmp ) )   ln_tsd_tradmp = .TRUE.     ! forces the initialization when tradmp is used
-      
-      IF(lwp) THEN                  ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'dta_tsd_init : Temperature & Salinity data '
-         WRITE(numout,*) '~~~~~~~~~~~~ '
-         WRITE(numout,*) '   Namelist namtsd'
-         WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init   = ', ln_tsd_init
-         WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_tradmp = ', ln_tsd_tradmp
-         WRITE(numout,*)
-         IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_tradmp ) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) '   T & S data not used'
-         ENDIF
-      ENDIF
-      !
-      IF( ln_rstart .AND. ln_tsd_init ) THEN
-         CALL ctl_warn( 'dta_tsd_init: ocean restart and T & S data intialisation, ',   &
-            &           'we keep the restart T & S values and set ln_tsd_init to FALSE' )
-         ln_tsd_init = .FALSE.
-      ENDIF
-      !
-      !                             ! allocate the arrays (if necessary)
-      IF(  ln_tsd_init .OR. ln_tsd_tradmp  ) THEN
-         !
-         ALLOCATE( sf_tsd(jpts), STAT=ierr0 )
-         IF( ierr0 > 0 ) THEN
-            CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' )   ;   RETURN
-         ENDIF
-         !
-                                ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 )
-         IF( sn_tem%ln_tint )   ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )
-                                ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 )
-         IF( sn_sal%ln_tint )   ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
-         !
-         IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN
-            CALL ctl_stop( 'dta_tsd : unable to allocate T & S data arrays' )   ;   RETURN
-         ENDIF
-         !                         ! fill sf_tsd with sn_tem & sn_sal and control print
-         slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal
-         CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd' )
-         !
-      ENDIF
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dta_tsd_init')
-      !
-   END SUBROUTINE dta_tsd_init
-
-
-   SUBROUTINE dta_tsd( kt, ptsd )
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE dta_tsd  ***
-      !!                    
-      !! ** Purpose :   provides T and S data at kt
-      !! 
-      !! ** Method  : - call fldread routine
-      !!              - ORCA_R2: add some hand made alteration to read data  
-      !!              - 'key_orca_lev10' interpolates on 10 times more levels
-      !!              - s- or mixed z-s coordinate: vertical interpolation on model mesh
-      !!              - ln_tsd_tradmp=F: deallocates the T-S data structure
-      !!                as T-S data are no are used
-      !!
-      !! ** Action  :   ptsd   T-S data on medl mesh and interpolated at time-step kt
-      !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   kt     ! ocean time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data
-      !
-      INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies
-      INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1   ! local integers
-      REAL(wp)::   zl, zi
-      REAL(wp), POINTER, DIMENSION(:) ::  ztp, zsp   ! 1D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dta_tsd')
-      !
-      CALL fld_read( kt, 1, sf_tsd )      !==   read T & S data at kt time step   ==!
-      !
-      !
-      !                                   !==   ORCA_R2 configuration and T & S damping   ==! 
-      IF( cp_cfg == "orca" .AND. jp_cfg == 2 .AND. ln_tsd_tradmp ) THEN    ! some hand made alterations
-         !
-         ij0 = 101   ;   ij1 = 109                       ! Reduced T & S in the Alboran Sea
-         ii0 = 141   ;   ii1 = 155
-         DO jj = mj0(ij0), mj1(ij1)
-            DO ji = mi0(ii0), mi1(ii1)
-               sf_tsd(jp_tem)%fnow(ji,jj,13:13) = sf_tsd(jp_tem)%fnow(ji,jj,13:13) - 0.20_wp
-               sf_tsd(jp_tem)%fnow(ji,jj,14:15) = sf_tsd(jp_tem)%fnow(ji,jj,14:15) - 0.35_wp
-               sf_tsd(jp_tem)%fnow(ji,jj,16:25) = sf_tsd(jp_tem)%fnow(ji,jj,16:25) - 0.40_wp
-               !
-               sf_tsd(jp_sal)%fnow(ji,jj,13:13) = sf_tsd(jp_sal)%fnow(ji,jj,13:13) - 0.15_wp
-               sf_tsd(jp_sal)%fnow(ji,jj,14:15) = sf_tsd(jp_sal)%fnow(ji,jj,14:15) - 0.25_wp
-               sf_tsd(jp_sal)%fnow(ji,jj,16:17) = sf_tsd(jp_sal)%fnow(ji,jj,16:17) - 0.30_wp
-               sf_tsd(jp_sal)%fnow(ji,jj,18:25) = sf_tsd(jp_sal)%fnow(ji,jj,18:25) - 0.35_wp
-            END DO
-         END DO
-         ij0 =  87   ;   ij1 =  96                          ! Reduced temperature in Red Sea
-         ii0 = 148   ;   ii1 = 160
-         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 ) = 7.0_wp
-         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp
-         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp
-      ENDIF
-      !
-      ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)    ! NO mask
-      ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) 
-      !
-      IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
-         !
-         CALL wrk_alloc( jpk, ztp, zsp )
-         !
-         IF( kt == nit000 .AND. lwp )THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'dta_tsd: interpolates T & S data onto the s- or mixed s-z-coordinate mesh'
-         ENDIF
-         !
-         DO jj = 1, jpj                         ! vertical interpolation of T & S
-            DO ji = 1, jpi
-               DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
-                  zl = gdept_0(ji,jj,jk)
-                  IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data
-                     ztp(jk) =  ptsd(ji,jj,1    ,jp_tem)
-                     zsp(jk) =  ptsd(ji,jj,1    ,jp_sal)
-                  ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data
-                     ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem)
-                     zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal)
-                  ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
-                     DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
-                        IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
-                           zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
-                           ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi 
-                           zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi
-                        ENDIF
-                     END DO
-                  ENDIF
-               END DO
-               DO jk = 1, jpkm1
-                  ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
-                  ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk)
-               END DO
-               ptsd(ji,jj,jpk,jp_tem) = 0._wp
-               ptsd(ji,jj,jpk,jp_sal) = 0._wp
-            END DO
-         END DO
-         ! 
-         CALL wrk_dealloc( jpk, ztp, zsp )
-         ! 
-      ELSE                                !==   z- or zps- coordinate   ==!
-         !                             
-         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:)    ! Mask
-         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:)
-         !
-         IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ik = mbkt(ji,jj) 
-                  IF( ik > 1 ) THEN
-                     zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
-                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
-                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
-                  ENDIF
-                  ik = mikt(ji,jj)
-                  IF( ik > 1 ) THEN
-                     zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 
-                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
-                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
-                  END IF
-               END DO
-            END DO
-         ENDIF
-         !
-      ENDIF
-      !
-      IF( .NOT.ln_tsd_tradmp ) THEN                   !==   deallocate T & S structure   ==! 
-         !                                              (data used only for initialisation)
-         IF(lwp) WRITE(numout,*) 'dta_tsd: deallocte T & S arrays as they are only use to initialize the run'
-                                        DEALLOCATE( sf_tsd(jp_tem)%fnow )     ! T arrays in the structure
-         IF( sf_tsd(jp_tem)%ln_tint )   DEALLOCATE( sf_tsd(jp_tem)%fdta )
-                                        DEALLOCATE( sf_tsd(jp_sal)%fnow )     ! S arrays in the structure
-         IF( sf_tsd(jp_sal)%ln_tint )   DEALLOCATE( sf_tsd(jp_sal)%fdta )
-                                        DEALLOCATE( sf_tsd              )     ! the structure itself
-      ENDIF
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dta_tsd')
-      !
-   END SUBROUTINE dta_tsd
-
-   !!======================================================================
-END MODULE dtatsd
Index: anches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 6826)
+++ 	(revision )
@@ -1,321 +1,0 @@
-MODULE iscplhsb
-   !!======================================================================
-   !!                       ***  MODULE  iscplhsb***
-   !! Ocean forcing: ice sheet/ocean coupling (conservation)
-   !!=====================================================================
-   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   iscpl_alloc    : variable allocation
-   !!   iscpl_hsb      : compute and store the input of heat/salt/volume 
-   !!                    into the system due to the coupling process
-   !!   iscpl_div      : correction of divergence to keep volume conservation
-   !!----------------------------------------------------------------------
-   USE dom_oce         ! ocean space and time domain
-   USE domwri          ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE sbc_oce         ! surface boundary condition variables
-   USE oce             ! global tra/dyn variable
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! MPP library
-   USE lib_fortran     ! MPP library
-   USE wrk_nemo        ! Memory allocation
-   USE lbclnk          !
-   USE domngb          !
-   USE iscplini
-
-   IMPLICIT NONE
-   PRIVATE
-   
-   PUBLIC   iscpl_div   
-   PUBLIC   iscpl_cons       
-   !! * Substitutions  
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE iscpl_cons(ptmask_b, psmask_b, pe3t_b, pts_flx, pvol_flx, prdt_iscpl)
-      !!---------------------------------------------------------------------- 
-      !!                   ***  ROUTINE iscpl_cons  ***
-      !! 
-      !! ** Purpose :   compute input into the system during the coupling step
-      !!                compute the correction term
-      !!                compute where the correction have to be applied
-      !! 
-      !! ** Method  :   compute tsn*e3t-tsb*e3tb and e3t-e3t_b
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b    !! mask before
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b      !! scale factor before
-      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b    !! mask before
-      REAL(wp), DIMENSION(:,:,:,:), INTENT(out) :: pts_flx     !! corrective flux to have tracer conservation
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(out) :: pvol_flx    !! corrective flux to have volume conservation
-      REAL(wp),                     INTENT(in ) :: prdt_iscpl  !! coupling period 
-      !!
-      INTEGER :: ji, jj, jk                                    !! loop index
-      INTEGER :: jip1, jim1, jjp1, jjm1
-      !!
-      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
-      REAL(wp):: r1_rdtiscpl
-      REAL(wp):: zjip1_ratio  , zjim1_ratio  , zjjp1_ratio  , zjjm1_ratio
-      !!
-      REAL(wp):: zde3t, zdtem, zdsal
-      REAL(wp), DIMENSION(:,:), POINTER :: zdssh
-      !!
-      REAL(wp), DIMENSION(:), ALLOCATABLE :: zlon, zlat
-      REAL(wp), DIMENSION(:), ALLOCATABLE :: zcorr_vol, zcorr_tem, zcorr_sal
-      INTEGER , DIMENSION(:), ALLOCATABLE :: ixpts, iypts, izpts, inpts
-      INTEGER :: jpts, npts
-
-      CALL wrk_alloc(jpi,jpj, zdssh )
-
-      ! get imbalance (volume heat and salt)
-      ! initialisation difference
-      zde3t = 0.0_wp; zdsal = 0.0_wp ; zdtem = 0.0_wp
-
-      ! initialisation correction term
-      pvol_flx(:,:,:  ) = 0.0_wp
-      pts_flx (:,:,:,:) = 0.0_wp
-      
-      r1_rdtiscpl = 1._wp / prdt_iscpl 
-
-      ! mask tsn and tsb 
-      tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)*ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)*tmask(:,:,:);
-      tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)*ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)*tmask(:,:,:);
-
-      !==============================================================================
-      ! diagnose the heat, salt and volume input and compute the correction variable
-      !==============================================================================
-
-      ! 
-      zdssh(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
-      IF (.NOT. ln_linssh ) zdssh = 0.0_wp ! already included in the levels by definition
-      
-      DO jk = 1,jpk-1
-         DO jj = 2,jpj-1
-            DO ji = fs_2,fs_jpim1
-               IF (tmask_h(ji,jj) == 1._wp) THEN
-
-                  ! volume differences
-                  zde3t = e3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
-
-                  ! heat diff
-                  zdtem = tsn(ji,jj,jk,jp_tem) * e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
-                        - tsb(ji,jj,jk,jp_tem) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
-                  ! salt diff
-                  zdsal = tsn(ji,jj,jk,jp_sal) * e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
-                        - tsb(ji,jj,jk,jp_sal) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
-               
-                  ! shh changes
-                  IF ( ptmask_b(ji,jj,jk) == 1._wp .OR. tmask(ji,jj,jk) == 1._wp ) THEN 
-                     zde3t = zde3t + zdssh(ji,jj) ! zdssh = 0 if vvl
-                     zdssh(ji,jj) = 0._wp
-                  END IF
-
-                  ! volume, heat and salt differences in each cell 
-                  pvol_flx(ji,jj,jk)       =   pvol_flx(ji,jj,jk)        + zde3t * r1_rdtiscpl
-                  pts_flx (ji,jj,jk,jp_sal)=   pts_flx (ji,jj,jk,jp_sal) + zdsal * r1_rdtiscpl 
-                  pts_flx (ji,jj,jk,jp_tem)=   pts_flx (ji,jj,jk,jp_tem) + zdtem * r1_rdtiscpl
-
-                  ! case where we close a cell: check if the neighbour cells are wet 
-                  IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
-
-                     jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ;
-
-                     zsum =   e1e2t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e1e2t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) &
-                       &    + e1e2t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e1e2t(jip1,jj  ) * tmask(jip1,jj  ,jk)
-
-                     IF ( zsum /= 0._wp ) THEN
-                        zjip1_ratio   = e1e2t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum
-                        zjim1_ratio   = e1e2t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum
-                        zjjp1_ratio   = e1e2t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum
-                        zjjm1_ratio   = e1e2t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) / zsum
-
-                        pvol_flx(ji  ,jjp1,jk       ) = pvol_flx(ji  ,jjp1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjp1_ratio
-                        pvol_flx(ji  ,jjm1,jk       ) = pvol_flx(ji  ,jjm1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjm1_ratio
-                        pvol_flx(jip1,jj  ,jk       ) = pvol_flx(jip1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjip1_ratio
-                        pvol_flx(jim1,jj  ,jk       ) = pvol_flx(jim1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjim1_ratio
-                        pts_flx (ji  ,jjp1,jk,jp_sal) = pts_flx (ji  ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio
-                        pts_flx (ji  ,jjm1,jk,jp_sal) = pts_flx (ji  ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio
-                        pts_flx (jip1,jj  ,jk,jp_sal) = pts_flx (jip1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio
-                        pts_flx (jim1,jj  ,jk,jp_sal) = pts_flx (jim1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio
-                        pts_flx (ji  ,jjp1,jk,jp_tem) = pts_flx (ji  ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio
-                        pts_flx (ji  ,jjm1,jk,jp_tem) = pts_flx (ji  ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio
-                        pts_flx (jip1,jj  ,jk,jp_tem) = pts_flx (jip1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio
-                        pts_flx (jim1,jj  ,jk,jp_tem) = pts_flx (jim1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio
-
-                        ! set to 0 the cell we distributed over neigbourg cells
-                        pvol_flx(ji,jj,jk       ) = 0._wp
-                        pts_flx (ji,jj,jk,jp_sal) = 0._wp
-                        pts_flx (ji,jj,jk,jp_tem) = 0._wp
-
-                     ELSE IF (zsum == 0._wp ) THEN
-                        ! case where we close a cell and no adjacent cell open
-                        ! check if the cell beneath is wet
-                        IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
-                           pvol_flx(ji,jj,jk+1)       =  pvol_flx(ji,jj,jk+1)        + pvol_flx(ji,jj,jk)
-                           pts_flx (ji,jj,jk+1,jp_sal)=  pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal)
-                           pts_flx (ji,jj,jk+1,jp_tem)=  pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem)
-
-                           ! set to 0 the cell we distributed over neigbourg cells
-                           pvol_flx(ji,jj,jk       ) = 0._wp
-                           pts_flx (ji,jj,jk,jp_sal) = 0._wp
-                           pts_flx (ji,jj,jk,jp_tem) = 0._wp
-                        ELSE
-                        ! case no adjacent cell on the horizontal and on the vertical
-                           IF ( lwp ) THEN   ! JMM : cAution this warning may occur on any mpp subdomain but numout is only
-                                             ! open for narea== 1 (lwp=T)
-                           WRITE(numout,*) 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
-                           WRITE(numout,*) '                     ',mig(ji),' ',mjg(jj),' ',jk
-                           WRITE(numout,*) '                     ',ji,' ',jj,' ',jk,' ',narea
-                           WRITE(numout,*) ' we are now looking for the closest wet cell on the horizontal '
-                           ENDIF
-                        ! We deal with these points later.
-                        END IF
-                     END IF
-                  END IF
-               END IF
-            END DO
-         END DO
-      END DO
-
-      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.)
-      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
-      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
-
-      ! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point 
-      ! allocation and initialisation of the list of problematic point
-      ALLOCATE(inpts(jpnij))
-      inpts(:)=0
-
-      ! fill narea location with the number of problematic point
-      DO jk = 1,jpk-1
-         DO jj = 2,jpj-1
-            DO ji = fs_2,fs_jpim1
-               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
-                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
-                  inpts(narea) = inpts(narea) + 1 
-               END IF
-            END DO
-         END DO
-      END DO
-
-      ! build array of total problematic point on each cpu (share to each cpu)
-      CALL mpp_max(inpts,jpnij) 
-
-      ! size of the new variable
-      npts  = SUM(inpts)    
-      
-      ! allocation of the coordinates, correction, index vector for the problematic points 
-      ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
-      ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20_wp ; zlat(:) = -1.0e20_wp
-      zcorr_vol(:) = -1.0e20_wp
-      zcorr_sal(:) = -1.0e20_wp
-      zcorr_tem(:) = -1.0e20_wp
-
-      ! fill new variable
-      jpts = SUM(inpts(1:narea-1))
-      DO jk = 1,jpk-1
-         DO jj = 2,jpj-1
-            DO ji = fs_2,fs_jpim1
-               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
-                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
-                  jpts = jpts + 1  ! positioning in the inpts vector for the area narea
-                  ixpts(jpts) = ji           ; iypts(jpts) = jj ; izpts(jpts) = jk
-                  zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj)
-                  zcorr_vol(jpts) = pvol_flx(ji,jj,jk)
-                  zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal)
-                  zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem)
-
-                  ! set flx to 0 (safer)
-                  pvol_flx(ji,jj,jk       ) = 0.0_wp
-                  pts_flx (ji,jj,jk,jp_sal) = 0.0_wp
-                  pts_flx (ji,jj,jk,jp_tem) = 0.0_wp
-               END IF
-            END DO
-         END DO
-      END DO
-
-      ! build array of total problematic point on each cpu (share to each cpu)
-      ! point coordinates
-      CALL mpp_max(zlat ,npts)
-      CALL mpp_max(zlon ,npts)
-      CALL mpp_max(izpts,npts)
-
-      ! correction values 
-      CALL mpp_max(zcorr_vol,npts)
-      CALL mpp_max(zcorr_sal,npts)
-      CALL mpp_max(zcorr_tem,npts)
-
-      ! put correction term in the closest cell          
-      DO jpts = 1,npts
-         CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
-         DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
-            DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
-               jk = izpts(jpts)
-
-               IF (tmask_h(ji,jj) == 1._wp) THEN
-                  ! correct the vol_flx in the closest cell
-                  pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts)
-                  pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts)
-                  pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts)
-
-                  ! set correction to 0
-                  zcorr_vol(jpts) = 0.0_wp
-                  zcorr_sal(jpts) = 0.0_wp
-                  zcorr_tem(jpts) = 0.0_wp
-               END IF
-            END DO
-         END DO
-      END DO
-
-      ! deallocate variables 
-      DEALLOCATE(inpts)
-      DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
-    
-      ! add contribution store on the hallo (lbclnk remove one of the contribution)
-      pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:)
-      pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
-      pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
-
-      ! compute sum over the halo and set it to 0.
-      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1._wp)
-      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1._wp)
-      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1._wp)
-
-      ! deallocate variables
-      CALL wrk_dealloc(jpi,jpj, zdssh ) 
-
-   END SUBROUTINE iscpl_cons
-
-   SUBROUTINE iscpl_div( phdivn )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE iscpl_div  ***
-      !!
-      !! ** Purpose :   update the horizontal divergenc
-      !!
-      !! ** Method  :
-      !!                CAUTION : iscpl is positive (inflow) and expressed in m/s
-      !!
-      !! ** Action  :   phdivn   increase by the iscpl correction term
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
-      !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / e3t_n(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-   END SUBROUTINE iscpl_div
-
-END MODULE iscplhsb
Index: anches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90	(revision 6826)
+++ 	(revision )
@@ -1,85 +1,0 @@
-MODULE iscplini
-   !!======================================================================
-   !!                       ***  MODULE  sbciscpl***
-   !! Ocean forcing:  river runoff
-   !!=====================================================================
-   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   iscpl_init     : initialisation routine (namelist)
-   !!   iscpl_alloc    : allocation of correction variables
-   !!----------------------------------------------------------------------
-   USE dom_oce         ! ocean space and time domain
-   USE oce             ! global tra/dyn variable
-   USE lib_mpp         ! MPP library
-   USE lib_fortran     ! MPP library
-   USE in_out_manager  ! I/O manager
-
-   IMPLICIT NONE
-   PRIVATE
-   
-   PUBLIC   iscpl_init      
-   PUBLIC   iscpl_alloc 
-   !!                                                      !!* namsbc_iscpl namelist *
-   LOGICAL , PUBLIC                                        ::   ln_hsb
-   INTEGER , PUBLIC                                        ::   nn_fiscpl, nstp_iscpl
-   INTEGER , PUBLIC                                        ::   nn_drown
-   REAL(wp), PUBLIC                                        ::   rdt_iscpl
-   !!                                                      !!* namsbc_iscpl namelist *
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) ::   hdiv_iscpl
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   htsc_iscpl
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   INTEGER FUNCTION iscpl_alloc()
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE sbc_iscpl_alloc  ***
-      !!----------------------------------------------------------------------
-      ALLOCATE( htsc_iscpl(jpi,jpj,jpk,jpts) , hdiv_iscpl(jpi,jpj,jpk) , STAT=iscpl_alloc )
-         !
-      IF( lk_mpp          )   CALL mpp_sum ( iscpl_alloc )
-      IF( iscpl_alloc > 0 )   CALL ctl_warn('iscpl_alloc: allocation of arrays failed')
-   END FUNCTION iscpl_alloc
-
-   SUBROUTINE iscpl_init()
-      INTEGER ::   ios           ! Local integer output status for namelist read
-      NAMELIST/namsbc_iscpl/nn_fiscpl,ln_hsb,nn_drown
-      !!----------------------------------------------------------------------
-      !                                   ! ============
-      !                                   !   Namelist
-      !                                   ! ============
-      !
-      nn_fiscpl = 0
-      ln_hsb    = .FALSE.
-      REWIND( numnam_ref )              ! Namelist namsbc_iscpl in reference namelist : Ice sheet coupling
-      READ  ( numnam_ref, namsbc_iscpl, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_iscpl in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namsbc_iscpl in configuration namelist : Ice Sheet coupling
-      READ  ( numnam_cfg, namsbc_iscpl, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_iscpl in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namsbc_iscpl )
-      !
-      nstp_iscpl=MIN(nn_fiscpl, nitend-nit000+1) ! the coupling period have to be less or egal than the total number of time step
-      rdt_iscpl = nstp_iscpl * rn_rdt
-      !
-      IF (lwp) THEN
-         WRITE(numout,*) 'iscpl_rst:'
-         WRITE(numout,*) '~~~~~~~~~'
-         WRITE(numout,*) ' coupling     flag (ln_iscpl )            = ', ln_iscpl
-         WRITE(numout,*) ' conservation flag (ln_hsb   )            = ', ln_hsb
-         WRITE(numout,*) ' nb of stp for cons (rn_fiscpl)           = ', nstp_iscpl
-         IF (nstp_iscpl .NE. nn_fiscpl) WRITE(numout,*) 'W A R N I N G: nb of stp for cons has been modified &
-            &                                           (larger than run length)'
-         WRITE(numout,*) ' coupling time step                       = ', rdt_iscpl
-         WRITE(numout,*) ' number of call of the extrapolation loop = ', nn_drown
-      END IF
-
-   END SUBROUTINE iscpl_init
-
-END MODULE iscplini
Index: anches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 6826)
+++ 	(revision )
@@ -1,441 +1,0 @@
-MODULE iscplrst
-   !!======================================================================
-   !!                       ***  MODULE  iscplrst***
-   !! Ocean forcing: update the restart file in case of ice sheet/ocean coupling
-   !!=====================================================================
-   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   iscpl_stp          : step management
-   !!   iscpl_rst_interpol : restart interpolation in case of coupling with ice sheet
-   !!----------------------------------------------------------------------
-   USE dom_oce         ! ocean space and time domain
-   USE domwri          ! ocean space and time domain
-   USE domvvl, ONLY : dom_vvl_interpol
-   USE phycst          ! physical constants
-   USE sbc_oce         ! surface boundary condition variables
-   USE oce             ! global tra/dyn variable
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O module
-   USE lib_mpp         ! MPP library
-   USE lib_fortran     ! MPP library
-   USE wrk_nemo        ! Memory allocation
-   USE lbclnk          ! communication
-   USE iscplini        ! ice sheet coupling: initialisation
-   USE iscplhsb        ! ice sheet coupling: conservation
-
-   IMPLICIT NONE
-   PRIVATE
-   
-   PUBLIC   iscpl_stp          ! step management 
-   !!
-   !! * Substitutions  
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE iscpl_stp
-      !!---------------------------------------------------------------------- 
-      !!                   ***  ROUTINE iscpl_stp  ***
-      !! 
-      !! ** Purpose : compute initialisation
-      !!              compute extrapolation of restart variable un, vn, tsn, sshn (wetting/drying)   
-      !!              compute correction term if needed
-      !! 
-      !!----------------------------------------------------------------------
-      INTEGER  ::   inum0
-      REAL(wp), DIMENSION(:,:  ), POINTER :: zsmask_b
-      REAL(wp), DIMENSION(:,:,:), POINTER :: ztmask_b, zumask_b, zvmask_b
-      REAL(wp), DIMENSION(:,:,:), POINTER :: ze3t_b  , ze3u_b  , ze3v_b  
-      REAL(wp), DIMENSION(:,:,:), POINTER :: zdepw_b
-      CHARACTER(20) :: cfile
-      !!----------------------------------------------------------------------
-
-      CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! mask before
-      CALL wrk_alloc(jpi,jpj,jpk, ze3t_b  , ze3u_b  , ze3v_b  ) ! e3   before
-      CALL wrk_alloc(jpi,jpj,jpk, zdepw_b )
-      CALL wrk_alloc(jpi,jpj,     zsmask_b                    )
-
-
-      !! get restart variable
-      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b   ) ! need to extrapolate T/S
-      CALL iom_get( numror, jpdom_autoglo, 'umask'  , zumask_b   ) ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'vmask'  , zvmask_b   ) ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'smask'  , zsmask_b   ) ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'e3t_n'  , ze3t_b(:,:,:) )  ! need to compute temperature correction
-      CALL iom_get( numror, jpdom_autoglo, 'e3u_n'  , ze3u_b(:,:,:) )  ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'e3v_n'  , ze3v_b(:,:,:) )  ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:) ) ! need to interpol vertical profile (vvl)
-
-      !! read namelist
-      CALL iscpl_init()
-
-      !!  ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
-      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
-
-      !! compute correction if conservation needed
-      IF ( ln_hsb ) THEN
-         IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
-         CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl)
-      END IF
-         
-      !! print mesh/mask
-      IF( nmsh /= 0 .AND. ln_iscpl )   CALL dom_wri      ! Create a domain file
-
-      IF ( ln_hsb ) THEN
-         cfile='correction'
-         cfile = TRIM( cfile )
-         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
-         CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
-         CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
-         CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
-         CALL iom_close ( inum0 )
-      END IF
-
-      CALL wrk_dealloc(jpi,jpj,jpk, ztmask_b,zumask_b,zvmask_b )  
-      CALL wrk_dealloc(jpi,jpj,jpk, ze3t_b  ,ze3u_b  ,ze3v_b   )  
-      CALL wrk_dealloc(jpi,jpj,jpk, zdepw_b                    )
-      CALL wrk_dealloc(jpi,jpj,     zsmask_b                   )
-
-      !! next step is an euler time step
-      neuler = 0
-
-      !! set _b and _n variables equal
-      tsb (:,:,:,:) = tsn (:,:,:,:)
-      ub  (:,:,:  ) = un  (:,:,:  )
-      vb  (:,:,:  ) = vn  (:,:,:  )
-      sshb(:,:    ) = sshn(:,:)
-
-      !! set _b and _n vertical scale factor equal
-      e3t_b (:,:,:) = e3t_n (:,:,:)
-      e3u_b (:,:,:) = e3u_n (:,:,:)
-      e3v_b (:,:,:) = e3v_n (:,:,:)
-
-      e3uw_b(:,:,:)  = e3uw_n(:,:,:)
-      e3vw_b(:,:,:)  = e3vw_n(:,:,:)
-      gdept_b(:,:,:)  = gdept_n(:,:,:)
-      gdepw_b(:,:,:) = gdepw_n(:,:,:)
-      hu_b (:,:) = hu_n(:,:)
-      hv_b (:,:) = hv_n(:,:)
-      r1_hu_b(:,:) = r1_hu_n(:,:)
-      r1_hv_b(:,:) = r1_hv_n(:,:)
-      !
-   END SUBROUTINE iscpl_stp
-   
-   SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
-      !!---------------------------------------------------------------------- 
-      !!                   ***  ROUTINE iscpl_rst_interpol  ***
-      !! 
-      !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves 
-      !!                compute 2d fields of heat, salt and volume correction
-      !! 
-      !! ** Method  :   tn, sn : extrapolation from neigbourg cells
-      !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
-      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
-      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
-      !!
-      INTEGER :: ji, jj, jk, iz          !! loop index
-      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
-      !!
-      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
-      REAL(wp):: zdz, zdzm1, zdzp1
-      !!
-      REAL(wp), DIMENSION(:,:    ), POINTER :: zdmask , zdsmask, zvcorr, zucorr, zde3t
-      REAL(wp), DIMENSION(:,:    ), POINTER :: zbub   , zbvb   , zbun  , zbvn
-      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  , zssh1, zhu1, zhv1
-      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1
-      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp
-      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zwmaskn, zwmaskb, ztmp3d
-      REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
-      !!----------------------------------------------------------------------
-
-      !! allocate variables
-      CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   )
-      CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d        ) 
-      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb                       ) 
-      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       ) 
-      CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
-      CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         ) 
-      CALL wrk_alloc(jpi,jpj,       zssh0  , zssh1, zhu1, zhv1             ) 
-
-      !! mask value to be sure
-      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
-      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
-      
-      ! compute wmask
-      zwmaskn(:,:,1) = tmask   (:,:,1)
-      zwmaskb(:,:,1) = ptmask_b(:,:,1)
-      DO jk = 2,jpk
-         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
-         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
-      END DO
-           
-      ! compute new ssh if we open a full water column (average of the closest neigbourgs)  
-      sshb (:,:)=sshn(:,:)
-      zssh0(:,:)=sshn(:,:)
-      zsmask0(:,:) = psmask_b(:,:)
-      zsmask1(:,:) = psmask_b(:,:)
-      DO iz = 1,10    ! need to be tuned (configuration dependent) (OK for ISOMIP+)
-         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
-         DO jj = 2,jpj-1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               jip1=ji+1; jim1=ji-1;
-               jjp1=jj+1; jjm1=jj-1;
-               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
-               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
-                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
-                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
-                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
-                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
-                  zsmask1(ji,jj)=1._wp
-               END IF
-            END DO
-         END DO
-         CALL lbc_lnk(sshn,'T',1._wp)
-         CALL lbc_lnk(zsmask1,'T',1._wp)
-         zssh0   = sshn
-         zsmask0 = zsmask1
-      END DO
-      sshn(:,:) = sshn(:,:) * ssmask(:,:)
-
-!=============================================================================
-!PM: Is this needed since introduction of VVL by default?
-      IF (.NOT.ln_linssh) THEN
-      ! Reconstruction of all vertical scale factors at now time steps
-      ! =============================================================================
-      ! Horizontal scale factor interpolations
-      ! --------------------------------------
-         DO jk = 1,jpk
-            DO jj=1,jpj
-               DO ji=1,jpi
-                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
-                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) / ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-         CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
-
-      ! Vertical scale factor interpolations
-      ! ------------------------------------
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
-         CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-         CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
-
-      ! t- and w- points depth
-      ! ----------------------
-         gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-         gdepw_n(:,:,1) = 0.0_wp
-         gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               DO jk = 2,mikt(ji,jj)-1
-                  gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
-                  gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
-               END DO
-               IF (mikt(ji,jj) > 1) THEN
-                  jk = mikt(ji,jj)
-                  gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
-                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj)
-               END IF
-               DO jk = mikt(ji,jj)+1, jpk
-                  gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn (ji,jj)
-               END DO
-            END DO
-         END DO
-
-      ! t-, u- and v- water column thickness
-      ! ------------------------------------
-         ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp
-         DO jk = 1, jpkm1
-            hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-            hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
-            ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-         !                                        ! Inverse of the local depth
-         r1_hu_n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
-         r1_hv_n(:,:) = 1._wp / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
-
-      END IF
-
-!=============================================================================
-! compute velocity
-! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
-      ub(:,:,:)=un(:,:,:)
-      vb(:,:,:)=vn(:,:,:)
-      DO jk = 1,jpk
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               un(ji,jj,jk) = ub(ji,jj,jk)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u_n(ji,jj,jk)*umask(ji,jj,jk);
-               vn(ji,jj,jk) = vb(ji,jj,jk)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v_n(ji,jj,jk)*vmask(ji,jj,jk);
-            END DO
-         END DO
-      END DO
-
-! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
-! compute barotropic velocity now and after 
-      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
-      zbub(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
-      zbvb(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:); 
-      zbun(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:); 
-      zbvn(:,:)   = SUM(ztrp,DIM=3)
-
-      ! new water column
-      zhu1=0.0_wp ;
-      zhv1=0.0_wp ;
-      DO jk  = 1,jpk
-        zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-        zhv1(:,:) = zhv1(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
-      END DO
-      
-      ! compute correction      
-      zucorr = 0._wp
-      zvcorr = 0._wp
-      DO jj = 1,jpj
-         DO ji = 1,jpi
-            IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
-               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
-            END IF
-            IF (zbvn(ji,jj) /= zbvb(ji,jj) .AND. zhv1(ji,jj) /= 0._wp ) THEN
-               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/zhv1(ji,jj)
-            END IF
-         END DO
-      END DO 
-      
-      ! update velocity
-      DO jk  = 1,jpk
-         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
-         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
-      END DO
-
-!=============================================================================
-      ! compute temp and salt
-      ! compute new tn and sn if we open a new cell
-      tsb (:,:,:,:) = tsn(:,:,:,:)
-      zts0(:,:,:,:) = tsn(:,:,:,:)
-      ztmask1(:,:,:) = ptmask_b(:,:,:)
-      ztmask0(:,:,:) = ptmask_b(:,:,:)
-      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
-          DO jk = 1,jpk-1
-              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
-              DO jj = 2,jpj-1
-                 DO ji = fs_2,fs_jpim1
-                      jip1=ji+1; jim1=ji-1;
-                      jjp1=jj+1; jjm1=jj-1;
-                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
-                      IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
-                      !! horizontal basic extrapolation
-                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
-                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
-                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
-                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
-                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
-                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
-                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
-                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
-                         ztmask1(ji,jj,jk)=1
-                      ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
-                      !! vertical extrapolation if horizontal extrapolation failed
-                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
-                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
-                         IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
-                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
-                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
-                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
-                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
-                            ztmask1(ji,jj,jk)=1._wp
-                         END IF
-                      END IF
-                  END DO
-              END DO
-          END DO
-          
-          CALL lbc_lnk(tsn(:,:,:,1),'T',1._wp)
-          CALL lbc_lnk(tsn(:,:,:,2),'T',1._wp)
-          CALL lbc_lnk(ztmask1,     'T',1._wp)
-
-          ! update
-          zts0(:,:,:,:) = tsn(:,:,:,:)
-          ztmask0 = ztmask1
-
-      END DO
-
-      ! mask new tsn field
-      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
-      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
-
-      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
-      !PM: Is this IF needed since change to VVL by default
-      IF (.NOT.ln_linssh) THEN
-         DO jk = 2,jpk-1
-            DO jj = 1,jpj
-               DO ji = 1,jpi
-                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND. (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN
-                     !compute weight
-                     zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
-                     zdz   =           gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
-                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw_n(ji,jj,jk  ))
-                     IF (zdz .LT. 0._wp) THEN 
-                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
-                     END IF
-                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
-                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
-                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/e3t_n(ji,jj,jk)
-                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
-                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
-                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/e3t_n(ji,jj,jk)
-                  END IF
-               END DO
-            END DO
-         END DO               
-      END IF
-
-      ! closed pool
-      ! -----------------------------------------------------------------------------------------
-      ! case we open a cell but no neigbour cells available to get an estimate of T and S
-      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
-         tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
-         tmask(:,:,:) = 0._wp    ! set mask to 0 to run
-         umask(:,:,:) = 0._wp
-         vmask(:,:,:) = 0._wp
-      END WHERE
-      
-      ! set mbkt and mikt to 1 in thiese location
-      WHERE (SUM(tmask,dim=3) == 0)
-         mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
-         mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
-      END WHERE
-      ! -------------------------------------------------------------------------------------------
-      ! compute new tn and sn if we close cell 
-      ! nothing to do
-      ! 
-      ! deallocation tmp arrays
-      CALL wrk_dealloc(jpi,jpj,jpk,2, zts0                                   )
-      CALL wrk_dealloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp                ) 
-      CALL wrk_dealloc(jpi,jpj,jpk,   zwmaskn, zwmaskb , ztmp3d              ) 
-      CALL wrk_dealloc(jpi,jpj,       zsmask0, zsmask1                       ) 
-      CALL wrk_dealloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
-      CALL wrk_dealloc(jpi,jpj,       zbub   , zbvb    , zbun  , zbvn        ) 
-      CALL wrk_dealloc(jpi,jpj,       zssh0  , zssh1  , zhu1 , zhv1          ) 
-
-   END SUBROUTINE iscpl_rst_interpol
-
-END MODULE iscplrst
Index: anches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 6826)
+++ 	(revision )
@@ -1,506 +1,0 @@
-MODULE istate
-   !!======================================================================
-   !!                     ***  MODULE  istate  ***
-   !! Ocean state   :  initial state setting
-   !!=====================================================================
-   !! History :  OPA  !  1989-12  (P. Andrich)  Original code
-   !!            5.0  !  1991-11  (G. Madec)  rewritting
-   !!            6.0  !  1996-01  (G. Madec)  terrain following coordinates
-   !!            8.0  !  2001-09  (M. Levy, M. Ben Jelloul)  istate_eel
-   !!            8.0  !  2001-09  (M. Levy, M. Ben Jelloul)  istate_uvg
-   !!   NEMO     1.0  !  2003-08  (G. Madec, C. Talandier)  F90: Free form, modules + EEL R5
-   !!             -   !  2004-05  (A. Koch-Larrouy)  istate_gyre 
-   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
-   !!            3.3  !  2010-10  (C. Ethe) merge TRC-TRA
-   !!            3.4  !  2011-04  (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   istate_init   : initial state setting
-   !!   istate_tem    : analytical profile for initial Temperature
-   !!   istate_sal    : analytical profile for initial Salinity
-   !!   istate_eel    : initial state setting of EEL R5 configuration
-   !!   istate_gyre   : initial state setting of GYRE configuration
-   !!   istate_uvg    : initial velocity in geostropic balance
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers 
-   USE dom_oce         ! ocean space and time domain 
-   USE c1d             ! 1D vertical configuration
-   USE daymod          ! calendar
-   USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine)
-   USE ldftra          ! lateral physics: ocean active tracers
-   USE zdf_oce         ! ocean vertical physics
-   USE phycst          ! physical constants
-   USE dtatsd          ! data temperature and salinity   (dta_tsd routine)
-   USE dtauvd          ! data: U & V current             (dta_uvd routine)
-   USE domvvl          ! varying vertical mesh
-   USE iscplrst        ! ice sheet coupling
-   !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O library
-   USE lib_mpp         ! MPP library
-   USE restart         ! restart
-   USE wrk_nemo        ! Memory allocation
-   USE timing          ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   istate_init   ! routine called by step.F90
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE istate_init
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE istate_init  ***
-      !! 
-      !! ** Purpose :   Initialization of the dynamics and tracer fields.
-      !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )   CALL timing_start('istate_init')
-      !
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
-
-                     CALL dta_tsd_init        ! Initialisation of T & S input data
-      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
-
-      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
-      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
-      tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk
-      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
-
-      IF( ln_rstart ) THEN                    ! Restart from a file
-         !                                    ! -------------------
-         CALL rst_read                        ! Read the restart file
-         IF (ln_iscpl)       CALL iscpl_stp   ! extraloate restart to wet and dry
-         CALL day_init                        ! model calendar (using both namelist and restart infos)
-      ELSE
-         !                                    ! Start from rest
-         !                                    ! ---------------
-         numror = 0                              ! define numror = 0 -> no restart file to read
-         neuler = 0                              ! Set time-step indicator at nit000 (euler forward)
-         CALL day_init                           ! model calendar (using both namelist and restart infos)
-         !                                       ! Initialization of ocean to zero
-         !   before fields      !       now fields     
-         sshb (:,:)   = 0._wp   ;   sshn (:,:)   = 0._wp
-         ub   (:,:,:) = 0._wp   ;   un   (:,:,:) = 0._wp
-         vb   (:,:,:) = 0._wp   ;   vn   (:,:,:) = 0._wp  
-                                    hdivn(:,:,:) = 0._wp
-         !
-         IF( cp_cfg == 'eel' ) THEN
-            CALL istate_eel                      ! EEL   configuration : start from pre-defined U,V T-S fields
-         ELSEIF( cp_cfg == 'gyre' ) THEN         
-            CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields
-         ELSE                                    ! Initial T-S, U-V fields read in files
-            IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000
-               CALL dta_tsd( nit000, tsb )  
-               tsn(:,:,:,:) = tsb(:,:,:,:)
-               !
-            ELSE                                 ! Initial T-S fields defined analytically
-               CALL istate_t_s
-            ENDIF
-            IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
-               CALL wrk_alloc( jpi,jpj,jpk,2,   zuvd )
-               CALL dta_uvd( nit000, zuvd )
-               ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
-               vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
-               CALL wrk_dealloc( jpi,jpj,jpk,2,   zuvd )
-            ENDIF
-         ENDIF
-         !
-!!gm This is to be changed !!!!
-         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
-         IF( .NOT.ln_linssh ) THEN
-            DO jk = 1, jpk
-               e3t_b(:,:,jk) = e3t_n(:,:,jk)
-            END DO
-         ENDIF
-!!gm 
-         ! 
-      ENDIF
-      ! 
-      ! Initialize "now" and "before" barotropic velocities:
-      ! Do it whatever the free surface method, these arrays being eventually used
-      !
-      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
-      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
-      !
-!!gm  the use of umsak & vmask is not necessary belox as un, vn, ub, vb are always masked
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
-               vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
-               !
-               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
-               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-      un_b(:,:) = un_b(:,:) * r1_hu_n(:,:)
-      vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:)
-      !
-      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
-      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('istate_init')
-      !
-   END SUBROUTINE istate_init
-
-
-   SUBROUTINE istate_t_s
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE istate_t_s  ***
-      !!   
-      !! ** Purpose :   Intialization of the temperature field with an 
-      !!      analytical profile or a file (i.e. in EEL configuration)
-      !!
-      !! ** Method  : - temperature: use Philander analytic profile
-      !!              - salinity   : use to a constant value 35.5
-      !!
-      !! References :  Philander ???
-      !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zsal = 35.50_wp
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'istate_t_s : Philander s initial temperature profile'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~   and constant salinity (',zsal,' psu)'
-      !
-      DO jk = 1, jpk
-         tsn(:,:,jk,jp_tem) = (  ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((gdept_n(:,:,jk)-80.)/30.) )   &
-            &                + 10. * ( 5000. - gdept_n(:,:,jk) ) /5000.)  ) * tmask(:,:,jk)
-         tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
-      END DO
-      tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
-      tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
-      !
-   END SUBROUTINE istate_t_s
-
-
-   SUBROUTINE istate_eel
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE istate_eel  ***
-      !! 
-      !! ** Purpose :   Initialization of the dynamics and tracers for EEL R5
-      !!      configuration (channel with or without a topographic bump)
-      !!
-      !! ** Method  : - set temprature field
-      !!              - set salinity field
-      !!              - set velocity field including horizontal divergence
-      !!                and relative vorticity fields
-      !!----------------------------------------------------------------------
-      USE divhor     ! hor. divergence      (div_hor routine)
-      USE iom
-      !
-      INTEGER  ::   inum              ! temporary logical unit
-      INTEGER  ::   ji, jj, jk        ! dummy loop indices
-      INTEGER  ::   ijloc
-      REAL(wp) ::   zh1, zh2, zslope, zcst, zfcor   ! temporary scalars
-      REAL(wp) ::   zt1  = 15._wp                   ! surface temperature value (EEL R5)
-      REAL(wp) ::   zt2  =  5._wp                   ! bottom  temperature value (EEL R5)
-      REAL(wp) ::   zsal = 35.0_wp                  ! constant salinity (EEL R2, R5 and R6)
-      REAL(wp) ::   zueel = 0.1_wp                  ! constant uniform zonal velocity (EEL R5)
-      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zssh  ! initial ssh over the global domain
-      !!----------------------------------------------------------------------
-      !
-      SELECT CASE ( jp_cfg ) 
-         !                                              ! ====================
-         CASE ( 5 )                                     ! EEL R5 configuration
-            !                                           ! ====================
-            !
-            ! set temperature field with a linear profile
-            ! -------------------------------------------
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: linear temperature profile'
-            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
-            !
-            zh1 = gdept_1d(  1  )
-            zh2 = gdept_1d(jpkm1)
-            !
-            zslope = ( zt1 - zt2 ) / ( zh1 - zh2 )
-            zcst   = ( zt1 * ( zh1 - zh2) - ( zt1 - zt2 ) * zh1 ) / ( zh1 - zh2 )
-            !
-            DO jk = 1, jpk
-               tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - gdept_n(:,:,jk) / 1000 ) ) * tmask(:,:,jk)
-               tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
-            END DO
-            !
-            ! set salinity field to a constant value
-            ! --------------------------------------
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
-            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
-            !
-            tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
-            tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
-            !
-            ! set the dynamics: U,V, hdiv (and ssh if necessary)
-            ! ----------------
-            ! Start EEL5 configuration with barotropic geostrophic velocities 
-            ! according the sshb and sshn SSH imposed.
-            ! we assume a uniform grid (hence the use of e1t(1,1) for delta_y)
-            ! we use the Coriolis frequency at mid-channel.   
-            ub(:,:,:) = zueel * umask(:,:,:)
-            un(:,:,:) = ub(:,:,:)
-            ijloc = mj0(INT(jpjglo-1)/2)
-            zfcor = ff(1,ijloc)
-            !
-            DO jj = 1, jpjglo
-               zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)*zueel*e1t(1,1)*zfcor/grav 
-            END DO
-            !
-            IF(lwp) THEN
-               WRITE(numout,*) ' Uniform zonal velocity for EEL R5:',zueel
-               WRITE(numout,*) ' Geostrophic SSH profile as a function of y:'
-               WRITE(numout,'(12(1x,f6.2))') zssh(1,:)
-            ENDIF
-            !
-            DO jj = 1, nlcj
-               DO ji = 1, nlci
-                  sshb(ji,jj) = zssh( mig(ji) , mjg(jj) ) * tmask(ji,jj,1)
-               END DO
-            END DO
-            sshb(nlci+1:jpi,      :   ) = 0.e0      ! set to zero extra mpp columns
-            sshb(      :   ,nlcj+1:jpj) = 0.e0      ! set to zero extra mpp rows
-            !
-            sshn(:,:) = sshb(:,:)                   ! set now ssh to the before value
-            !
-            IF( nn_rstssh /= 0 ) THEN  
-               nn_rstssh = 0                        ! hand-made initilization of ssh 
-               CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' )
-            ENDIF
-            !
-!!gm  Check  here call to div_hor should not be necessary
-!!gm         div_hor call runoffs  not sure they are defined at that level
-            CALL div_hor( nit000 )                  ! horizontal divergence and relative vorticity (curl)
-            ! N.B. the vertical velocity will be computed from the horizontal divergence field
-            ! in istate by a call to wzv routine
-
-
-            !                                     ! ==========================
-         CASE ( 2 , 6 )                           ! EEL R2 or R6 configuration
-            !                                     ! ==========================
-            !
-            ! set temperature field with a NetCDF file
-            ! ----------------------------------------
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R2 or R6: read initial temperature in a NetCDF file'
-            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
-            !
-            CALL iom_open ( 'eel.initemp', inum )
-            CALL iom_get ( inum, jpdom_data, 'initemp', tsb(:,:,:,jp_tem) ) ! read before temprature (tb)
-            CALL iom_close( inum )
-            !
-            tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem)                            ! set nox temperature to tb
-            !
-            ! set salinity field to a constant value
-            ! --------------------------------------
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
-            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
-            !
-            tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
-            tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
-            !
-            !                                    ! ===========================
-         CASE DEFAULT                            ! NONE existing configuration
-            !                                    ! ===========================
-            WRITE(ctmp1,*) 'EEL with a ', jp_cfg,' km resolution is not coded'
-            CALL ctl_stop( ctmp1 )
-            !
-      END SELECT
-      !
-   END SUBROUTINE istate_eel
-
-
-   SUBROUTINE istate_gyre
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE istate_gyre  ***
-      !! 
-      !! ** Purpose :   Initialization of the dynamics and tracers for GYRE
-      !!      configuration (double gyre with rotated domain)
-      !!
-      !! ** Method  : - set temprature field
-      !!              - set salinity   field
-      !!----------------------------------------------------------------------
-      INTEGER :: ji, jj, jk  ! dummy loop indices
-      INTEGER            ::   inum          ! temporary logical unit
-      INTEGER, PARAMETER ::   ntsinit = 0   ! (0/1) (analytical/input data files) T&S initialization
-      !!----------------------------------------------------------------------
-      !
-      SELECT CASE ( ntsinit)
-      !
-      CASE ( 0 )                  ! analytical T/S profil deduced from LEVITUS
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS '
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-         !
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  tsn(ji,jj,jk,jp_tem) = (  16. - 12. * TANH( (gdept_n(ji,jj,jk) - 400) / 700 )         )   &
-                       &           * (-TANH( (500-gdept_n(ji,jj,jk)) / 150 ) + 1) / 2               &
-                       &       + (      15. * ( 1. - TANH( (gdept_n(ji,jj,jk)-50.) / 1500.) )       &
-                       &                - 1.4 * TANH((gdept_n(ji,jj,jk)-100.) / 100.)               &    
-                       &                + 7.  * (1500. - gdept_n(ji,jj,jk)) / 1500.             )   & 
-                       &           * (-TANH( (gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2
-                  tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
-                  tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem)
-
-                  tsn(ji,jj,jk,jp_sal) =  (  36.25 - 1.13 * TANH( (gdept_n(ji,jj,jk) - 305) / 460 )  )  &
-                     &              * (-TANH((500 - gdept_n(ji,jj,jk)) / 150) + 1) / 2          &
-                     &          + (  35.55 + 1.25 * (5000. - gdept_n(ji,jj,jk)) / 5000.         &
-                     &                - 1.62 * TANH( (gdept_n(ji,jj,jk) - 60.  ) / 650. )       &
-                     &                + 0.2  * TANH( (gdept_n(ji,jj,jk) - 35.  ) / 100. )       &
-                     &                + 0.2  * TANH( (gdept_n(ji,jj,jk) - 1000.) / 5000.)    )  &
-                     &              * (-TANH((gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2 
-                  tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
-                  tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal)
-               END DO
-            END DO
-         END DO
-         !
-      CASE ( 1 )                  ! T/S data fields read in dta_tem.nc/data_sal.nc files
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'istate_gyre : initial T and S read from dta_tem.nc/data_sal.nc files'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-         IF(lwp) WRITE(numout,*) '              NetCDF FORMAT'
-
-         ! Read temperature field
-         ! ----------------------
-         CALL iom_open ( 'data_tem', inum )
-         CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) ) 
-         CALL iom_close( inum )
-
-         tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:) 
-         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
-
-         ! Read salinity field
-         ! -------------------
-         CALL iom_open ( 'data_sal', inum )
-         CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) ) 
-         CALL iom_close( inum )
-
-         tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 
-         tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
-         !
-      END SELECT
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '              Initial temperature and salinity profiles:'
-         WRITE(numout, "(9x,' level   gdept_1d   temperature   salinity   ')" )
-         WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_1d(jk), tsn(2,2,jk,jp_tem), tsn(2,2,jk,jp_sal), jk = 1, jpk )
-      ENDIF
-      !
-   END SUBROUTINE istate_gyre
-
-
-   SUBROUTINE istate_uvg
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE istate_uvg  ***
-      !!
-      !! ** Purpose :   Compute the geostrophic velocities from (tn,sn) fields
-      !!
-      !! ** Method  :   Using the hydrostatic hypothesis the now hydrostatic 
-      !!      pressure is computed by integrating the in-situ density from the
-      !!      surface to the bottom.
-      !!                 p=integral [ rau*g dz ]
-      !!----------------------------------------------------------------------
-      USE divhor          ! hor. divergence                       (div_hor routine)
-      USE lbclnk          ! ocean lateral boundary condition (or mpp link)
-      !
-      INTEGER ::   ji, jj, jk        ! dummy loop indices
-      REAL(wp) ::   zmsv, zphv, zmsu, zphu, zalfg     ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zprn)
-      !
-      IF(lwp) WRITE(numout,*) 
-      IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
-
-      ! Compute the now hydrostatic pressure
-      ! ------------------------------------
-
-      zalfg = 0.5 * grav * rau0
-      
-      zprn(:,:,1) = zalfg * e3w_n(:,:,1) * ( 1 + rhd(:,:,1) )       ! Surface value
-
-      DO jk = 2, jpkm1                                              ! Vertical integration from the surface
-         zprn(:,:,jk) = zprn(:,:,jk-1)   &
-            &         + zalfg * e3w_n(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
-      END DO  
-
-      ! Compute geostrophic balance
-      ! ---------------------------
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vertor opt.
-               zmsv = 1. / MAX(  umask(ji-1,jj+1,jk) + umask(ji  ,jj+1,jk)   &
-                               + umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) , 1.  )
-               zphv = ( zprn(ji  ,jj+1,jk) - zprn(ji-1,jj+1,jk) ) * umask(ji-1,jj+1,jk) / e1u(ji-1,jj+1)   &
-                    + ( zprn(ji+1,jj+1,jk) - zprn(ji  ,jj+1,jk) ) * umask(ji  ,jj+1,jk) / e1u(ji  ,jj+1)   &
-                    + ( zprn(ji  ,jj  ,jk) - zprn(ji-1,jj  ,jk) ) * umask(ji-1,jj  ,jk) / e1u(ji-1,jj  )   &
-                    + ( zprn(ji+1,jj  ,jk) - zprn(ji  ,jj  ,jk) ) * umask(ji  ,jj  ,jk) / e1u(ji  ,jj  )
-               zphv = 1. / rau0 * zphv * zmsv * vmask(ji,jj,jk)
-
-               zmsu = 1. / MAX(  vmask(ji+1,jj  ,jk) + vmask(ji  ,jj  ,jk)   &
-                               + vmask(ji+1,jj-1,jk) + vmask(ji  ,jj-1,jk) , 1.  )
-               zphu = ( zprn(ji+1,jj+1,jk) - zprn(ji+1,jj  ,jk) ) * vmask(ji+1,jj  ,jk) / e2v(ji+1,jj  )   &
-                    + ( zprn(ji  ,jj+1,jk) - zprn(ji  ,jj  ,jk) ) * vmask(ji  ,jj  ,jk) / e2v(ji  ,jj  )   &
-                    + ( zprn(ji+1,jj  ,jk) - zprn(ji+1,jj-1,jk) ) * vmask(ji+1,jj-1,jk) / e2v(ji+1,jj-1)   &
-                    + ( zprn(ji  ,jj  ,jk) - zprn(ji  ,jj-1,jk) ) * vmask(ji  ,jj-1,jk) / e2v(ji  ,jj-1)
-               zphu = 1. / rau0 * zphu * zmsu * umask(ji,jj,jk)
-
-               ! Compute the geostrophic velocities
-               un(ji,jj,jk) = -2. * zphu / ( ff(ji,jj) + ff(ji  ,jj-1) )
-               vn(ji,jj,jk) =  2. * zphv / ( ff(ji,jj) + ff(ji-1,jj  ) )
-            END DO
-         END DO
-      END DO
-
-      IF(lwp) WRITE(numout,*) '         we force to zero bottom velocity'
-
-      ! Susbtract the bottom velocity (level jpk-1 for flat bottom case)
-      ! to have a zero bottom velocity
-
-      DO jk = 1, jpkm1
-         un(:,:,jk) = ( un(:,:,jk) - un(:,:,jpkm1) ) * umask(:,:,jk)
-         vn(:,:,jk) = ( vn(:,:,jk) - vn(:,:,jpkm1) ) * vmask(:,:,jk)
-      END DO
-
-      CALL lbc_lnk( un, 'U', -1. )
-      CALL lbc_lnk( vn, 'V', -1. )
-      
-      ub(:,:,:) = un(:,:,:)
-      vb(:,:,:) = vn(:,:,:)
-      
-      !
-!!gm  Check  here call to div_hor should not be necessary
-!!gm         div_hor call runoffs  not sure they are defined at that level
-      CALL div_hor( nit000 )            ! now horizontal divergence
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zprn)
-      !
-   END SUBROUTINE istate_uvg
-
-   !!=====================================================================
-END MODULE istate
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 6827)
@@ -20,27 +20,14 @@
    !!--------------------------------------------------------------------
    USE dom_oce         ! ocean space and time domain
-   USE c1d             ! 1D vertical configuration
-   USE flo_oce         ! floats module declarations
    USE lbclnk          ! lateal boundary condition / mpp exchanges
    USE iom_def         ! iom variables definitions
    USE iom_nf90        ! NetCDF format with native NetCDF library
    USE in_out_manager  ! I/O manager
-   USE lib_mpp           ! MPP library
-#if defined key_iomput
-   USE sbc_oce, ONLY :   nn_fsbc         ! ocean space and time domain
-   USE trc_oce, ONLY :   nn_dttrc        !  !: frequency of step on passive tracers
-   USE icb_oce, ONLY :   nclasses, class_num       !  !: iceberg classes
-#if defined key_lim3
-   USE ice    , ONLY :   jpl
-#elif defined key_lim2
-   USE par_ice_2
-#endif
+   USE lib_mpp         ! MPP library
    USE domngb          ! ocean space and time domain
    USE phycst          ! physical constants
-   USE dianam          ! build name of file
+ !SF  USE dianam          ! build name of file
    USE xios
-# endif
    USE ioipsl, ONLY :  ju2ymds    ! for calendar
-   USE crs             ! Grid coarsening
 
    IMPLICIT NONE
@@ -139,27 +126,4 @@
       ENDIF
 
-      IF( TRIM(cdname) == TRIM(cxios_context)//"_crs" ) THEN  
-         CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
-         !
-         CALL set_grid( "T", glamt_crs, gphit_crs ) 
-         CALL set_grid( "U", glamu_crs, gphiu_crs ) 
-         CALL set_grid( "V", glamv_crs, gphiv_crs ) 
-         CALL set_grid( "W", glamt_crs, gphit_crs ) 
-         CALL set_grid_znl( gphit_crs )
-          !
-         CALL dom_grid_glo   ! Return to parent grid domain
-         !
-         IF( ln_cfmeta ) THEN   ! Add additional grid metadata
-            CALL iom_set_domain_attr("grid_T", area = e1e2t_crs(nldi:nlei, nldj:nlej))
-            CALL iom_set_domain_attr("grid_U", area = e1u_crs(nldi:nlei, nldj:nlej) * e2u_crs(nldi:nlei, nldj:nlej))
-            CALL iom_set_domain_attr("grid_V", area = e1v_crs(nldi:nlei, nldj:nlej) * e2v_crs(nldi:nlei, nldj:nlej))
-            CALL iom_set_domain_attr("grid_W", area = e1e2t_crs(nldi:nlei, nldj:nlej))
-            CALL set_grid_bounds( "T", glamf_crs, gphif_crs, glamt_crs, gphit_crs )
-            CALL set_grid_bounds( "U", glamv_crs, gphiv_crs, glamu_crs, gphiu_crs )
-            CALL set_grid_bounds( "V", glamu_crs, gphiu_crs, glamv_crs, gphiv_crs )
-            CALL set_grid_bounds( "W", glamf_crs, gphif_crs, glamt_crs, gphit_crs )
-         ENDIF
-      ENDIF
-
       ! vertical grid definition
       CALL iom_set_axis_attr( "deptht", gdept_1d )
@@ -180,11 +144,4 @@
       CALL iom_set_axis_attr( "depthw", bounds=z_bnds )
 
-# if defined key_floats
-      CALL iom_set_axis_attr( "nfloat", (/ (REAL(ji,wp), ji=1,nfloat) /) )
-# endif
-#if defined key_lim3 || defined key_lim2
-      CALL iom_set_axis_attr( "ncatice", (/ (REAL(ji,wp), ji=1,jpl) /) )
-#endif
-      CALL iom_set_axis_attr( "icbcla", class_num )
       CALL iom_set_axis_attr( "iax_20C", (/ REAL(20,wp) /) )
       CALL iom_set_axis_attr( "iax_28C", (/ REAL(28,wp) /) )
@@ -882,8 +839,4 @@
             ENDIF
             
-            ! C1D case : always call lbc_lnk to replicate the central value over the whole 3X3 domain
-            IF( lk_c1d .AND. PRESENT(pv_r2d) )   CALL lbc_lnk( pv_r2d,'Z',1. )
-            IF( lk_c1d .AND. PRESENT(pv_r3d) )   CALL lbc_lnk( pv_r3d,'Z',1. )
-    
             !--- Apply scale_factor and offset
             zscf = iom_file(kiomid)%scf(idvar)      ! scale factor
@@ -1459,8 +1412,4 @@
       ! frequency of the call of iom_put (attribut: freq_op)
       WRITE(cl1,'(i1)')        1   ;   CALL iom_set_field_attr('field_definition', freq_op = cl1//'ts', freq_offset='0ts')
-      WRITE(cl1,'(i1)')  nn_fsbc   ;   CALL iom_set_field_attr('SBC'             , freq_op = cl1//'ts', freq_offset='0ts')
-      WRITE(cl1,'(i1)')  nn_fsbc   ;   CALL iom_set_field_attr('SBC_scalar'      , freq_op = cl1//'ts', freq_offset='0ts')
-      WRITE(cl1,'(i1)') nn_dttrc   ;   CALL iom_set_field_attr('ptrc_T'          , freq_op = cl1//'ts', freq_offset='0ts')
-      WRITE(cl1,'(i1)') nn_dttrc   ;   CALL iom_set_field_attr('diad_T'          , freq_op = cl1//'ts', freq_offset='0ts')
        
       ! output file names (attribut: name)
Index: anches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 6826)
+++ 	(revision )
@@ -1,271 +1,0 @@
-MODULE restart
-   !!======================================================================
-   !!                     ***  MODULE  restart  ***
-   !! Ocean restart :  write the ocean restart file
-   !!======================================================================
-   !! History :  OPA  !  1999-11  (M. Imbard)  Original code
-   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form
-   !!            2.0  !  2006-07  (S. Masson)  use IOM for restart
-   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
-   !!            - -  !  2010-10  (C. Ethe, G. Madec) TRC-TRA merge (T-S in 4D)
-   !!            3.7  !  2014-01  (G. Madec) suppression of curl and hdiv from the restart
-   !!             -   !  2014-12  (G. Madec) remove KPP scheme
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   rst_opn    : open the ocean restart file
-   !!   rst_write  : write the ocean restart file
-   !!   rst_read   : read the ocean restart file
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers 
-   USE dom_oce         ! ocean space and time domain
-   USE sbc_ice         ! only lk_lim3 
-   USE phycst          ! physical constants
-   USE eosbn2          ! equation of state            (eos bn2 routine)
-   USE trdmxl_oce      ! ocean active mixed layer tracers trends variables
-   !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O module
-   USE diurnal_bulk
-   
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   rst_opn         ! routine called by step module
-   PUBLIC   rst_write       ! routine called by step module
-   PUBLIC   rst_read        ! routine called by istate module
-   PUBLIC   rst_read_open   ! routine called in rst_read and (possibly) in dom_vvl_init
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE rst_opn( kt )
-      !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE rst_opn  ***
-      !!                     
-      !! ** Purpose : + initialization (should be read in the namelist) of nitrst 
-      !!              + open the restart when we are one time step before nitrst
-      !!                   - restart header is defined when kt = nitrst-1
-      !!                   - restart data  are written when kt = nitrst
-      !!              + define lrst_oce to .TRUE. when we need to define or write the restart
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt     ! ocean time-step
-      !!
-      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character
-      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
-      CHARACTER(lc)       ::   clpath   ! full path to ocean output restart file
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN   ! default definitions
-         lrst_oce = .FALSE.   
-         IF( ln_rst_list ) THEN
-            nrst_lst = 1
-            nitrst = nstocklist( nrst_lst )
-         ELSE
-            nitrst = nitend
-         ENDIF
-      ENDIF
-
-      ! frequency-based restart dumping (nn_stock)
-      IF( .NOT. ln_rst_list .AND. MOD( kt - 1, nstock ) == 0 ) THEN   
-         ! we use kt - 1 and not kt - nit000 to keep the same periodicity from the beginning of the experiment
-         nitrst = kt + nstock - 1                  ! define the next value of nitrst for restart writing
-         IF( nitrst > nitend )   nitrst = nitend   ! make sure we write a restart at the end of the run
-      ENDIF
-      ! to get better performances with NetCDF format:
-      ! we open and define the ocean restart file one time step before writing the data (-> at nitrst - 1)
-      ! except if we write ocean restart files every time step or if an ocean restart file was writen at nitend - 1
-      IF( kt == nitrst - 1 .OR. nstock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN
-         IF( nitrst <= nitend .AND. nitrst > 0 ) THEN 
-            ! beware of the format used to write kt (default is i8.8, that should be large enough...)
-            IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
-            ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst
-            ENDIF
-            ! create the file
-            clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_ocerst_out)
-            clpath = TRIM(cn_ocerst_outdir)
-            IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
-            IF(lwp) THEN
-               WRITE(numout,*)
-               SELECT CASE ( jprstlib )
-               CASE DEFAULT         ;   WRITE(numout,*)                            &
-                   '             open ocean restart NetCDF file: ',TRIM(clpath)//clname
-               END SELECT
-               IF ( snc4set%luse )      WRITE(numout,*) '             opened for NetCDF4 chunking and compression'
-               IF( kt == nitrst - 1 ) THEN   ;   WRITE(numout,*) '             kt = nitrst - 1 = ', kt
-               ELSE                          ;   WRITE(numout,*) '             kt = '             , kt
-               ENDIF
-            ENDIF
-            !
-            CALL iom_open( TRIM(clpath)//TRIM(clname), numrow, ldwrt = .TRUE., kiolib = jprstlib )
-            lrst_oce = .TRUE.
-         ENDIF
-      ENDIF
-      !
-   END SUBROUTINE rst_opn
-
-
-   SUBROUTINE rst_write( kt )
-      !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE rstwrite  ***
-      !!                     
-      !! ** Purpose :   Write restart fields in the format corresponding to jprstlib
-      !!
-      !! ** Method  :   Write in numrow when kt == nitrst in NetCDF
-      !!              file, save fields which are necessary for restart
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step
-      !!----------------------------------------------------------------------
-
-                     CALL iom_rstput( kt, nitrst, numrow, 'rdt'    , rdt       )   ! dynamics and tracer time step
-
-      IF ( .NOT. ln_diurnal_only ) THEN
-                     CALL iom_rstput( kt, nitrst, numrow, 'ub'     , ub        )     ! before fields
-                     CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vb        )
-                     CALL iom_rstput( kt, nitrst, numrow, 'tb'     , tsb(:,:,:,jp_tem) )
-                     CALL iom_rstput( kt, nitrst, numrow, 'sb'     , tsb(:,:,:,jp_sal) )
-                     CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb      )
-                     !
-                     CALL iom_rstput( kt, nitrst, numrow, 'un'     , un        )     ! now fields
-                     CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vn        )
-                     CALL iom_rstput( kt, nitrst, numrow, 'tn'     , tsn(:,:,:,jp_tem) )
-                     CALL iom_rstput( kt, nitrst, numrow, 'sn'     , tsn(:,:,:,jp_sal) )
-                     CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , sshn      )
-                     CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop      )
-
-                  ! extra variable needed for the ice sheet coupling
-                  IF ( ln_iscpl ) THEN 
-                     CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask     ) ! need to extrapolate T/S
-                     CALL iom_rstput( kt, nitrst, numrow, 'umask'  , umask     ) ! need to correct barotropic velocity
-                     CALL iom_rstput( kt, nitrst, numrow, 'vmask'  , vmask     ) ! need to correct barotropic velocity
-                     CALL iom_rstput( kt, nitrst, numrow, 'smask'  , ssmask    ) ! need to correct barotropic velocity
-                     CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) )   ! need to compute temperature correction
-                     CALL iom_rstput( kt, nitrst, numrow, 'e3u_n', e3u_n(:,:,:) )   ! need to compute bt conservation
-                     CALL iom_rstput( kt, nitrst, numrow, 'e3v_n', e3v_n(:,:,:) )   ! need to compute bt conservation
-                     CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw_n(:,:,:) ) ! need to compute extrapolation if vvl
-                  END IF
-      ENDIF
-      
-      IF (ln_diurnal) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst    )  
-
-      IF( kt == nitrst ) THEN
-         CALL iom_close( numrow )     ! close the restart file (only at last time step)
-!!gm         IF( .NOT. lk_trdmld )   lrst_oce = .FALSE.
-!!gm  not sure what to do here   ===>>>  ask to Sebastian
-         lrst_oce = .FALSE.
-            IF( ln_rst_list ) THEN
-               nrst_lst = MIN(nrst_lst + 1, SIZE(nstocklist,1))
-               nitrst = nstocklist( nrst_lst )
-            ENDIF
-            lrst_oce = .FALSE.
-      ENDIF
-      !
-   END SUBROUTINE rst_write
-
-
-   SUBROUTINE rst_read_open
-      !!---------------------------------------------------------------------- 
-      !!                   ***  ROUTINE rst_read_open  ***
-      !! 
-      !! ** Purpose :   Open read files for restart (format fixed by jprstlib )
-      !! 
-      !! ** Method  :   Use a non-zero, positive value of numror to assess whether or not
-      !!                the file has already been opened
-      !!----------------------------------------------------------------------
-      INTEGER        ::   jlibalt = jprstlib
-      LOGICAL        ::   llok
-      CHARACTER(lc)  ::   clpath   ! full path to ocean output restart file
-      !!----------------------------------------------------------------------
-      !
-      IF( numror <= 0 ) THEN
-         IF(lwp) THEN                                             ! Contol prints
-            WRITE(numout,*)
-            SELECT CASE ( jprstlib )
-            CASE ( jpnf90    )   ;   WRITE(numout,*) 'rst_read : read oce NetCDF restart file'
-            END SELECT
-            IF ( snc4set%luse )      WRITE(numout,*) 'rst_read : configured with NetCDF4 support'
-            WRITE(numout,*) '~~~~~~~~'
-         ENDIF
-
-         clpath = TRIM(cn_ocerst_indir)
-         IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
-         CALL iom_open( TRIM(clpath)//cn_ocerst_in, numror, kiolib = jlibalt )
-      ENDIF
-   END SUBROUTINE rst_read_open
-
-
-   SUBROUTINE rst_read
-      !!---------------------------------------------------------------------- 
-      !!                   ***  ROUTINE rst_read  ***
-      !! 
-      !! ** Purpose :   Read files for restart (format fixed by jprstlib )
-      !! 
-      !! ** Method  :   Read in restart.nc file fields which are necessary for restart
-      !!----------------------------------------------------------------------
-      REAL(wp) ::   zrdt
-      INTEGER  ::   jk
-      !!----------------------------------------------------------------------
-
-      CALL rst_read_open           ! open restart for reading (if not already opened)
-
-      ! Check dynamics and tracer time-step consistency and force Euler restart if changed
-      IF( iom_varid( numror, 'rdt', ldstop = .FALSE. ) > 0 )   THEN
-         CALL iom_get( numror, 'rdt', zrdt )
-         IF( zrdt /= rdt )   neuler = 0
-      ENDIF
-
-      ! Diurnal DSST 
-      IF( ln_diurnal ) CALL iom_get( numror, jpdom_autoglo, 'Dsst' , x_dsst  ) 
-      IF ( ln_diurnal_only ) THEN 
-         IF(lwp) WRITE( numout, * ) &
-         &   "rst_read:- ln_diurnal_only set, setting rhop=rau0" 
-         rhop = rau0
-         CALL iom_get( numror, jpdom_autoglo, 'tn'     , tsn(:,:,1,jp_tem) ) 
-         RETURN 
-      ENDIF  
-      
-      IF( iom_varid( numror, 'ub', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'ub'     , ub      )   ! before fields
-         CALL iom_get( numror, jpdom_autoglo, 'vb'     , vb      )
-         CALL iom_get( numror, jpdom_autoglo, 'tb'     , tsb(:,:,:,jp_tem) )
-         CALL iom_get( numror, jpdom_autoglo, 'sb'     , tsb(:,:,:,jp_sal) )
-         CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb    )
-      ELSE
-         neuler = 0
-      ENDIF
-      !
-      CALL iom_get( numror, jpdom_autoglo, 'un'     , un      )   ! now    fields
-      CALL iom_get( numror, jpdom_autoglo, 'vn'     , vn      )
-      CALL iom_get( numror, jpdom_autoglo, 'tn'     , tsn(:,:,:,jp_tem) )
-      CALL iom_get( numror, jpdom_autoglo, 'sn'     , tsn(:,:,:,jp_sal) )
-      CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    )
-      IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'rhop'   , rhop    )   ! now    potential density
-      ELSE
-         CALL eos( tsn, rhd, rhop, gdept_n(:,:,:) )   
-      ENDIF
-      !
-      IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0)
-         tsb  (:,:,:,:) = tsn  (:,:,:,:)                             ! all before fields set to now values
-         ub   (:,:,:)   = un   (:,:,:)
-         vb   (:,:,:)   = vn   (:,:,:)
-         sshb (:,:)     = sshn (:,:)
-         !
-         IF( .NOT.ln_linssh ) THEN
-            DO jk = 1, jpk
-               e3t_b(:,:,jk) = e3t_n(:,:,jk)
-            END DO
-         ENDIF
-         !
-      ENDIF
-      !
-   END SUBROUTINE rst_read
-
-   !!=====================================================================
-END MODULE restart
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 6827)
@@ -39,6 +39,4 @@
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
-   USE stopar         ! Stochastic T/S fluctuations
-   USE stopts         ! Stochastic T/S fluctuations
    !
    USE in_out_manager ! I/O manager
@@ -335,68 +333,4 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         ! Stochastic equation of state
-         IF ( ln_sto_eos ) THEN
-            ALLOCATE(zn0_sto(1:2*nn_sto_eos))
-            ALLOCATE(zn_sto(1:2*nn_sto_eos))
-            ALLOCATE(zsign(1:2*nn_sto_eos))
-            DO jsmp = 1, 2*nn_sto_eos, 2
-              zsign(jsmp)   = 1._wp
-              zsign(jsmp+1) = -1._wp
-            END DO
-            !
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     !
-                     ! compute density (2*nn_sto_eos) times:
-                     ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)
-                     ! (2) for t-dt, s-ds (with the opposite fluctuation)
-                     DO jsmp = 1, nn_sto_eos*2
-                        jdof   = (jsmp + 1) / 2
-                        zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth
-                        zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature
-                        zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp)
-                        zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity
-                        ztm    = tmask(ji,jj,jk)                                         ! tmask
-                        !
-                        zn3 = EOS013*zt   &
-                           &   + EOS103*zs+EOS003
-                           !
-                        zn2 = (EOS022*zt   &
-                           &   + EOS112*zs+EOS012)*zt   &
-                           &   + (EOS202*zs+EOS102)*zs+EOS002
-                           !
-                        zn1 = (((EOS041*zt   &
-                           &   + EOS131*zs+EOS031)*zt   &
-                           &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
-                           &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
-                           &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
-                           !
-                        zn0_sto(jsmp) = (((((EOS060*zt   &
-                           &   + EOS150*zs+EOS050)*zt   &
-                           &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
-                           &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
-                           &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
-                           &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
-                           &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
-                           !
-                        zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp)
-                     END DO
-                     !
-                     ! compute stochastic density as the mean of the (2*nn_sto_eos) densities
-                     prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp
-                     DO jsmp = 1, nn_sto_eos*2
-                        prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface
-                        !
-                        prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked)
-                     END DO
-                     prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
-                     prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos
-                  END DO
-               END DO
-            END DO
-            DEALLOCATE(zn0_sto,zn_sto,zsign)
-         ! Non-stochastic equation of state
-         ELSE
             DO jk = 1, jpkm1
                DO jj = 1, jpj
@@ -437,5 +371,4 @@
                END DO
             END DO
-         ENDIF
          
       CASE( np_seos )                !==  simplified EOS  ==!
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 6827)
@@ -48,49 +48,11 @@
    USE mppini         ! shared/distributed memory setting (mpp_init routine)
    USE domain         ! domain initialization             (dom_init routine)
-#if defined key_nemocice_decomp
-   USE ice_domain_size, only: nx_global, ny_global
-#endif
-   USE tideini        ! tidal components initialization   (tide_ini routine)
-   USE bdyini         ! open boundary cond. setting       (bdy_init routine)
-   USE bdydta         ! open boundary cond. setting   (bdy_dta_init routine)
-   USE bdytides       ! open boundary cond. setting   (bdytide_init routine)
-   USE sbctide, ONLY  : lk_tide
-   USE istate         ! initial state setting          (istate_init routine)
-   USE ldfdyn         ! lateral viscosity setting      (ldfdyn_init routine)
-   USE ldftra         ! lateral diffusivity setting    (ldftra_init routine)
-   USE zdfini         ! vertical physics setting          (zdf_init routine)
    USE phycst         ! physical constant                  (par_cst routine)
-   USE trdini         ! dyn/tra trends initialization     (trd_init routine)
-   USE asminc         ! assimilation increments     
-   USE asmbkg         ! writing out state trajectory
-   USE diaptr         ! poleward transports           (dia_ptr_init routine)
-   USE diadct         ! sections transports           (dia_dct_init routine)
-   USE diaobs         ! Observation diagnostics       (dia_obs_init routine)
-   USE diacfl         ! CFL diagnostics               (dia_cfl_init routine)
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
-   USE step           ! NEMO time-stepping                 (stp     routine)
-   USE icbini         ! handle bergs, initialisation
-   USE icbstp         ! handle bergs, calving, themodynamics and transport
-   USE cpl_oasis3     ! OASIS3 coupling
-   USE c1d            ! 1D configuration
-   USE step_c1d       ! Time stepping loop for the 1D configuration
-   USE dyndmp         ! Momentum damping
-   USE stopar         ! Stochastic param.: ???
-   USE stopts         ! Stochastic param.: ???
-#if defined key_top
-   USE trcini         ! passive tracer initialisation
-#endif
    USE lib_mpp        ! distributed memory computing
-   USE diurnal_bulk    ! diurnal bulk SST 
-   USE step_diu        ! diurnal bulk SST timestepping (called from here if run offline)
 #if defined key_iomput
    USE xios           ! xIOserver
 #endif
-   USE crsini         ! initialise grid coarsening utility
    USE lbcnfd , ONLY  : isendto, nsndto, nfsloop, nfeloop   ! Setup of north fold exchanges 
-   USE sbc_oce, ONLY  : lk_oasis
-   USE diatmb          ! Top,middle,bottom output
-   USE dia25h          ! 25h mean output
-   USE wet_dry         ! Wetting and drying setting   (wad_init routine)
 
    IMPLICIT NONE
@@ -137,10 +99,4 @@
       CALL Agrif_Declare_Var_dom   ! AGRIF: set the meshes for DOM
       CALL Agrif_Declare_Var       !  "      "   "   "      "  DYN/TRA 
-# if defined key_top
-      CALL Agrif_Declare_Var_top   !  "      "   "   "      "  TOP
-# endif
-# if defined key_lim2
-      CALL Agrif_Declare_Var_lim2  !  "      "   "   "      "  LIM
-# endif
 #endif
       ! check that all process are still there... If some process have an error,
@@ -151,49 +107,6 @@
 
       !                            !-----------------------!
-      !                            !==   time stepping   ==!
-      !                            !-----------------------!
-      istp = nit000
-#if defined key_c1d
-         DO WHILE ( istp <= nitend .AND. nstop == 0 )
-            CALL stp_c1d( istp )
-            istp = istp + 1
-         END DO
-#else
-          IF( lk_asminc ) THEN
-             IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 )    ! Output background fields
-             IF( ln_asmdin ) THEN                        ! Direct initialization
-                IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1 )    ! Tracers
-                IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics
-                IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH
-             ENDIF
-          ENDIF
-
-#if defined key_agrif
-          CALL Agrif_Regrid()
-#endif
-
-         DO WHILE ( istp <= nitend .AND. nstop == 0 )
-#if defined key_agrif
-            CALL stp                         ! AGRIF: time stepping
-#else
-            IF ( .NOT. ln_diurnal_only ) THEN 
-               CALL stp( istp )                 ! standard time stepping 
-            ELSE 
-               CALL stp_diurnal( istp )        ! time step only the diurnal SST 
-            ENDIF 
-#endif
-            istp = istp + 1
-            IF( lk_mpp )   CALL mpp_max( nstop )
-         END DO
-#endif
-
-      IF( ln_diaobs   )   CALL dia_obs_wri
-      !
-      IF( ln_icebergs )   CALL icb_end( nitend )
-
-      !                            !------------------------!
       !                            !==  finalize the run  ==!
       !                            !------------------------!
-      IF(lwp) WRITE(numout,cform_aaa)   ! Flag AAAAAAA
       !
       IF( nstop /= 0 .AND. lwp ) THEN   ! error print
@@ -204,6 +117,5 @@
 #if defined key_agrif
       IF( .NOT. Agrif_Root() ) THEN
-                         CALL Agrif_ParentGrid_To_ChildGrid()
-         IF( ln_diaobs ) CALL dia_obs_wri
+                                CALL Agrif_ParentGrid_To_ChildGrid()
          IF( nn_timing == 1 )   CALL timing_finalize
                                 CALL Agrif_ChildGrid_To_ParentGrid()
@@ -216,11 +128,6 @@
 #if defined key_iomput
       CALL xios_finalize                  ! end mpp communications with xios
-      IF( lk_oasis )   CALL cpl_finalize  ! end coupling and mpp communications with OASIS
 #else
-      IF( lk_oasis ) THEN 
-         CALL cpl_finalize              ! end coupling and mpp communications with OASIS
-      ELSE
-         IF( lk_mpp )   CALL mppstop    ! end mpp communications
-      ENDIF
+      IF( lk_mpp )   CALL mppstop    ! end mpp communications
 #endif
       !
@@ -294,25 +201,12 @@
 #if defined key_iomput
       IF( Agrif_Root() ) THEN
-         IF( lk_oasis ) THEN
-            CALL cpl_init( "oceanx", ilocal_comm )                     ! nemo local communicator given by oasis
-            CALL xios_initialize( "not used",local_comm=ilocal_comm )    ! send nemo communicator to xios
-         ELSE
             CALL  xios_initialize( "for_xios_mpi_id",return_comm=ilocal_comm )    ! nemo local communicator given by xios
-         ENDIF
       ENDIF
       ! Nodes selection (control print return in cltxt)
       narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )
 #else
-      IF( lk_oasis ) THEN
-         IF( Agrif_Root() ) THEN
-            CALL cpl_init( "oceanx", ilocal_comm )                      ! nemo local communicator given by oasis
-         ENDIF
-         ! Nodes selection (control print return in cltxt)
-         narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )
-      ELSE
          ilocal_comm = 0
          ! Nodes selection (control print return in cltxt)
          narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop )
-      ENDIF
 #endif
       narea = narea + 1                                     ! mynode return the rank of proc (0 --> jpnij -1 )
@@ -402,106 +296,7 @@
                             CALL     phy_cst    ! Physical constants
                             CALL     eos_init   ! Equation of state
-      IF( lk_c1d        )   CALL     c1d_init   ! 1D column configuration
-                            CALL     wad_init   ! Wetting and drying options
                             CALL     dom_cfg    ! Domain configuration
                             CALL     dom_init   ! Domain
-      IF( ln_crs        )   CALL     crs_init   ! coarsened grid: domain initialization 
-      IF( ln_nnogather )    CALL nemo_northcomms! northfold neighbour lists (must be done after the masks are defined)
       IF( ln_ctl        )   CALL prt_ctl_init   ! Print control
-      
-      CALL diurnal_sst_bulk_init            ! diurnal sst
-      IF ( ln_diurnal ) CALL diurnal_sst_coolskin_init   ! cool skin   
-      
-      ! IF ln_diurnal_only, then we only want a subset of the initialisation routines
-      IF ( ln_diurnal_only ) THEN
-         CALL  istate_init   ! ocean initial state (Dynamics and tracers)
-         CALL     sbc_init   ! Forcings : surface module
-         CALL tra_qsr_init   ! penetrative solar radiation qsr
-         IF( ln_diaobs     ) THEN                  ! Observation & model comparison
-            CALL dia_obs_init            ! Initialize observational data
-            CALL dia_obs( nit000 - 1 )   ! Observation operator for restart
-         ENDIF     
-         !                                     ! Assimilation increments
-         IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments
-                 
-         IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler
-         RETURN
-      ENDIF
-      
-                            CALL  istate_init   ! ocean initial state (Dynamics and tracers)
-
-      !                                      ! external forcing 
-!!gm to be added : creation and call of sbc_apr_init
-      IF( lk_tide       )   CALL    tide_init   ! tidal harmonics
-                            CALL     sbc_init   ! surface boundary conditions (including sea-ice)
-!!gm ==>> bdy_init should call bdy_dta_init and bdytide_init  NOT in nemogcm !!!
-      IF( lk_bdy        )   CALL     bdy_init   ! Open boundaries initialisation
-      IF( lk_bdy        )   CALL bdy_dta_init   ! Open boundaries initialisation of external data arrays
-      IF( lk_bdy .AND. lk_tide )   &
-         &                  CALL bdytide_init   ! Open boundaries initialisation of tidal harmonic forcing
-         
-      !                                      ! Ocean physics
-      !                                         ! Vertical physics
-                            CALL     zdf_init      ! namelist read
-                            CALL zdf_bfr_init      ! bottom friction
-      IF( lk_zdfric     )   CALL zdf_ric_init      ! Richardson number dependent Kz
-      IF( lk_zdftke     )   CALL zdf_tke_init      ! TKE closure scheme
-      IF( lk_zdfgls     )   CALL zdf_gls_init      ! GLS closure scheme
-      IF( lk_zdftmx     )   CALL zdf_tmx_init      ! tidal vertical mixing
-      IF( lk_zdfddm     )   CALL zdf_ddm_init      ! double diffusive mixing
-         
-      !                                         ! Lateral physics
-                            CALL ldf_tra_init      ! Lateral ocean tracer physics
-                            CALL ldf_eiv_init      ! eddy induced velocity param.
-                            CALL ldf_dyn_init      ! Lateral ocean momentum physics
-
-      !                                         ! Active tracers
-                            CALL tra_qsr_init      ! penetrative solar radiation qsr
-                            CALL tra_bbc_init      ! bottom heat flux
-      IF( lk_trabbl     )   CALL tra_bbl_init      ! advective (and/or diffusive) bottom boundary layer scheme
-                            CALL tra_dmp_init      ! internal tracer damping
-                            CALL tra_adv_init      ! horizontal & vertical advection
-                            CALL tra_ldf_init      ! lateral mixing
-                            CALL tra_zdf_init      ! vertical mixing and after tracer fields
-
-      !                                         ! Dynamics
-      IF( lk_c1d        )   CALL dyn_dmp_init      ! internal momentum damping
-                            CALL dyn_adv_init      ! advection (vector or flux form)
-                            CALL dyn_vor_init      ! vorticity term including Coriolis
-                            CALL dyn_ldf_init      ! lateral mixing
-                            CALL dyn_hpg_init      ! horizontal gradient of Hydrostatic pressure
-                            CALL dyn_zdf_init      ! vertical diffusion
-                            CALL dyn_spg_init      ! surface pressure gradient
-
-#if defined key_top
-      !                                      ! Passive tracers
-                            CALL     trc_init
-#endif
-      IF( l_ldfslp      )   CALL ldf_slp_init   ! slope of lateral mixing
-
-      !                                      ! Icebergs
-                            CALL icb_init( rdt, nit000)   ! initialise icebergs instance
-
-      !                                      ! Misc. options
-                            CALL sto_par_init   ! Stochastic parametrization
-      IF( ln_sto_eos     )  CALL sto_pts_init   ! RRandom T/S fluctuations
-     
-      !                                      ! Diagnostics
-      IF( lk_floats     )   CALL     flo_init   ! drifting Floats
-                            CALL dia_cfl_init   ! Initialise CFL diagnostics
-      IF( lk_diaar5     )   CALL dia_ar5_init   ! ar5 diag
-                            CALL dia_ptr_init   ! Poleward TRansports initialization
-      IF( lk_diadct     )   CALL dia_dct_init   ! Sections tranports
-                            CALL dia_hsb_init   ! heat content, salt content and volume budgets
-                            CALL     trd_init   ! Mixed-layer/Vorticity/Integral constraints trends
-                            CALL dia_obs_init            ! Initialize observational data
-      IF( ln_diaobs     )   CALL dia_obs( nit000 - 1 )   ! Observation operator for restart
-
-      !                                         ! Assimilation increments
-      IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments
-      IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler
-                            CALL dia_tmb_init  ! TMB outputs
-                            CALL dia_25h_init  ! 25h mean  outputs
-
       !
    END SUBROUTINE nemo_init
@@ -600,12 +395,4 @@
       ENDIF
       !
-      IF( nbench == 1 ) THEN              ! Benchmark
-         SELECT CASE ( cp_cfg )
-         CASE ( 'gyre' )   ;   CALL ctl_warn( ' The Benchmark is activated ' )
-         CASE DEFAULT      ;   CALL ctl_stop( ' The Benchmark is based on the GYRE configuration:',   &
-            &                                 ' cp_cfg = "gyre" in namelist &namcfg or set nbench = 0' )
-         END SELECT
-      ENDIF
-      !
       IF( 1_wp /= SIGN(1._wp,-0._wp)  )   CALL ctl_stop( 'nemo_ctl: The intrinsec SIGN function follows ',  &
          &                                               'f2003 standard. '                              ,  &
@@ -653,13 +440,5 @@
       !! ** Method  :
       !!----------------------------------------------------------------------
-      USE diawri    , ONLY: dia_wri_alloc
       USE dom_oce   , ONLY: dom_oce_alloc
-      USE trc_oce   , ONLY: trc_oce_alloc
-#if defined key_diadct 
-      USE diadct    , ONLY: diadct_alloc 
-#endif 
-#if defined key_bdy
-      USE bdy_oce   , ONLY: bdy_oce_alloc
-#endif
       !
       INTEGER :: ierr
@@ -667,16 +446,5 @@
       !
       ierr =        oce_alloc       ()          ! ocean
-      ierr = ierr + dia_wri_alloc   ()
       ierr = ierr + dom_oce_alloc   ()          ! ocean domain
-      ierr = ierr + zdf_oce_alloc   ()          ! ocean vertical physics
-      !
-      ierr = ierr + trc_oce_alloc   ()          ! shared TRC / TRA arrays
-      !
-#if defined key_diadct 
-      ierr = ierr + diadct_alloc    ()          ! 
-#endif 
-#if defined key_bdy
-      ierr = ierr + bdy_oce_alloc   ()          ! bdy masks (incl. initialization)
-#endif
       !
       IF( lk_mpp    )   CALL mpp_sum( ierr )
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 6826)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 6827)
@@ -9,96 +9,10 @@
    USE oce              ! ocean dynamics and tracers variables
    USE dom_oce          ! ocean space and time domain variables
-   USE zdf_oce          ! ocean vertical physics variables
 
    USE daymod           ! calendar                         (day     routine)
 
-   USE sbc_oce          ! surface boundary condition: ocean
-   USE sbcmod           ! surface boundary condition       (sbc     routine)
-   USE sbcrnf           ! surface boundary condition: runoff variables
-   USE sbccpl           ! surface boundary condition: coupled formulation (call send at end of step)
-   USE sbcapr           ! surface boundary condition: atmospheric pressure
-   USE sbctide          ! Tide initialisation
-
-   USE traqsr           ! solar radiation penetration      (tra_qsr routine)
-   USE trasbc           ! surface boundary condition       (tra_sbc routine)
-   USE trabbc           ! bottom boundary condition        (tra_bbc routine)
-   USE trabbl           ! bottom boundary layer            (tra_bbl routine)
-   USE tradmp           ! internal damping                 (tra_dmp routine)
-   USE traadv           ! advection scheme control     (tra_adv_ctl routine)
-   USE traldf           ! lateral mixing                   (tra_ldf routine)
-   USE trazdf           ! vertical mixing                  (tra_zdf routine)
-   USE tranxt           ! time-stepping                    (tra_nxt routine)
-   USE tranpc           ! non-penetrative convection       (tra_npc routine)
-
    USE eosbn2           ! equation of state                (eos_bn2 routine)
 
-   USE divhor           ! horizontal divergence            (div_hor routine)
-   USE dynadv           ! advection                        (dyn_adv routine)
-   USE dynbfr           ! Bottom friction terms            (dyn_bfr routine)
-   USE dynvor           ! vorticity term                   (dyn_vor routine)
-   USE dynhpg           ! hydrostatic pressure grad.       (dyn_hpg routine)
-   USE dynldf           ! lateral momentum diffusion       (dyn_ldf routine)
-   USE dynzdf           ! vertical diffusion               (dyn_zdf routine)
-   USE dynspg           ! surface pressure gradient        (dyn_spg routine)
-
-   USE dynnxt           ! time-stepping                    (dyn_nxt routine)
-
-   USE stopar           ! Stochastic parametrization       (sto_par routine)
-   USE stopts 
-
-   USE bdy_par          ! for lk_bdy
-   USE bdy_oce          ! for dmp logical
-   USE bdydta           ! open boundary condition data     (bdy_dta routine)
-   USE bdytra           ! bdy cond. for tracers            (bdy_tra routine)
-   USE bdydyn3d         ! bdy cond. for baroclinic vel.  (bdy_dyn3d routine)
-
-   USE sshwzv           ! vertical velocity and ssh        (ssh_nxt routine)
-   !                                                       (ssh_swp routine)
-   !                                                       (wzv     routine)
-   USE domvvl           ! variable vertical scale factors  (dom_vvl_sf_nxt routine)
-   !                                                       (dom_vvl_sf_swp routine)
-
-   USE ldfslp           ! iso-neutral slopes               (ldf_slp routine)
-   USE ldfdyn           ! lateral eddy viscosity coef.     (ldf_dyn routine)
-   USE ldftra           ! lateral eddy diffusive coef.     (ldf_tra routine)
-
-   USE zdftmx           ! tide-induced vertical mixing     (zdf_tmx routine)
-   USE zdfbfr           ! bottom friction                  (zdf_bfr routine)
-   USE zdftke           ! TKE vertical mixing              (zdf_tke routine)
-   USE zdfgls           ! GLS vertical mixing              (zdf_gls routine)
-   USE zdfddm           ! double diffusion mixing          (zdf_ddm routine)
-   USE zdfevd           ! enhanced vertical diffusion      (zdf_evd routine)
-   USE zdfric           ! Richardson vertical mixing       (zdf_ric routine)
-   USE zdfmxl           ! Mixed-layer depth                (zdf_mxl routine)
-
-   USE step_diu        ! Time stepping for diurnal sst
-   USE diurnal_bulk    ! diurnal SST bulk routines  (diurnal_sst_takaya routine) 
-   USE cool_skin       ! diurnal cool skin correction (diurnal_sst_coolskin routine)   
-   USE sbc_oce         ! surface fluxes  
-   
-   USE zpshde           ! partial step: hor. derivative     (zps_hde routine)
-
-   USE diawri           ! Standard run outputs             (dia_wri routine)
-   USE diaptr           ! poleward transports              (dia_ptr routine)
-   USE diadct           ! sections transports              (dia_dct routine)
-   USE diaar5           ! AR5 diagnosics                   (dia_ar5 routine)
-   USE diahth           ! thermocline depth                (dia_hth routine)
-   USE diafwb           ! freshwater budget                (dia_fwb routine)
-   USE diahsb           ! heat, salt and volume budgets    (dia_hsb routine)
-   USE diaharm
-   USE diacfl
-   USE flo_oce          ! floats variables
-   USE floats           ! floats computation               (flo_stp routine)
-
-   USE crsfld           ! Standard output on coarse grid   (crs_fld routine)
-
-   USE asminc           ! assimilation increments      (tra_asm_inc routine)
-   !                                                   (dyn_asm_inc routine)
-   USE asmbkg
-   USE stpctl           ! time stepping control            (stp_ctl routine)
-   USE restart          ! ocean restart                    (rst_wri routine)
    USE prtctl           ! Print control                    (prt_ctl routine)
-
-   USE diaobs           ! Observation operator
 
    USE in_out_manager   ! I/O manager
@@ -110,11 +24,4 @@
    USE xios
 #endif
-#if defined key_agrif
-   USE agrif_opa_sponge ! Momemtum and tracers sponges
-   USE agrif_opa_update ! Update (2-way nesting)
-#endif
-#if defined key_top
-   USE trcstp           ! passive tracer time-stepping      (trc_stp routine)
-#endif
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.7 , NEMO Consortium (2014)
