Index: trunk/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 508)
@@ -5,4 +5,8 @@
    !!                 (please no more than 2 lines)
    !!=====================================================================
+   !! History :  9.0  !  03-09  (C. Talandir, G. Madec)  Original code
+   !!            9.0  !  06-01  (A. Biastoch)  Allow sub-basins computation
+   !!----------------------------------------------------------------------
+
    !!----------------------------------------------------------------------
    !!   dia_ptr      : Poleward Transport Diagnostics module
@@ -14,9 +18,4 @@
    !!                : flux array; Generic interface: ptr_vj_3d, ptr_vj_2d
    !!----------------------------------------------------------------------
-   !! History :
-   !!   9.0  !  03-09  (C. Talandir, G. Madec)  Original code
-   !!   9.0  !  06-01  (A. Biastoch)  Allow sub-basins computation
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce           ! ocean dynamics and active tracers
    USE dom_oce       ! ocean space and time domain
@@ -26,5 +25,6 @@
    USE dianam
    USE phycst
-   USE ioipsl          ! NetCDF IPSL library
+   USE iom
+   USE ioipsl         
    USE daymod
 
@@ -36,39 +36,34 @@
    END INTERFACE
 
-   !! *  Routine accessibility
-   PUBLIC dia_ptr_init   ! call in opa module
-   PUBLIC dia_ptr        ! call in step module
-   PUBLIC ptr_vj         ! call by tra_ldf & tra_adv routines
-   PUBLIC ptr_vjk        ! call by tra_ldf & tra_adv routines
-
-   !! * Share Module variables
-   LOGICAL, PUBLIC ::       & !!! ** init namelist (namptr) **
-      ln_diaptr = .FALSE.,  &  !: Poleward transport flag (T) or not (F)
-      ln_subbas = .FALSE.      !: Atlantic/Pacific/Indian basins calculation
-   INTEGER, PUBLIC ::       & !!: ** ptr namelist (namptr) **
-      nf_ptr = 15              !: frequency of ptr computation
-   REAL(wp), PUBLIC, DIMENSION(jpj) ::   &   !!: poleward transport
-      pht_adv, pst_adv,     &  !: heat and salt: advection
-      pht_ove, pst_ove,     &  !: heat and salt: overturning
-      pht_ldf, pst_ldf,     &  !: heat and salt: lateral diffusion
-#if defined key_diaeiv
-      pht_eiv, pst_eiv,     &  !: heat and salt: bolus advection
-#endif
-      ht_atl,ht_ind,ht_pac, &  !: heat
-      st_atl,st_ind,st_pac     !: salt
-   REAL(wp),DIMENSION(jpi,jpj) ::   &
-      abasin,pbasin,ibasin     !: return function value
+   PUBLIC   dia_ptr_init   ! call in opa module
+   PUBLIC   dia_ptr        ! call in step module
+   PUBLIC   ptr_vj         ! call by tra_ldf & tra_adv routines
+   PUBLIC   ptr_vjk        ! call by tra_ldf & tra_adv routines
+
+   !!! ** init namelist (namptr)
+   LOGICAL , PUBLIC                 ::   ln_diaptr = .FALSE.   !: Poleward transport flag (T) or not (F)
+   LOGICAL , PUBLIC                 ::   ln_subbas = .FALSE.   !: Atlantic/Pacific/Indian basins calculation
+   INTEGER , PUBLIC                 ::   nf_ptr = 15           !: frequency of ptr computation
+
+   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_adv, pst_adv      !: heat and salt poleward transport: advection
+   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_ove, pst_ove      !: heat and salt poleward transport: overturning
+   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_ldf, pst_ldf      !: heat and salt poleward transport: lateral diffusion
+#if defined key_diaeiv
+   REAL(wp), PUBLIC, DIMENSION(jpj) ::   pht_eiv, pst_eiv      !: heat and salt poleward transport: bolus advection
+#endif
+   REAL(wp), PUBLIC, DIMENSION(jpj) ::   ht_atl,ht_ind,ht_pac  !: heat
+   REAL(wp), PUBLIC, DIMENSION(jpj) ::   st_atl,st_ind,st_pac  !: salt
+
      
 
-   !! Module variables
-   REAL(wp), DIMENSION(jpj,jpk) ::   &  
-      tn_jk  , sn_jk  ,  &  !: "zonal" mean temperature and salinity
-      v_msf_atl       ,  &  !: "meridional" Stream-Function
-      v_msf_glo       ,  &  !: "meridional" Stream-Function
-      v_msf_ipc       ,  &  !: "meridional" Stream-Function
-#if defined key_diaeiv
-      v_msf_eiv       ,  &  !: bolus "meridional" Stream-Function
-#endif
-      surf_jk_r             !: inverse of the ocean "zonal" section surface
+   REAL(wp), DIMENSION(jpj,jpk) ::   tn_jk  , sn_jk  ,  &  !: "zonal" mean temperature and salinity
+      &                              v_msf_atl       ,  &  !: "meridional" Stream-Function
+      &                              v_msf_glo       ,  &  !: "meridional" Stream-Function
+      &                              v_msf_ipc       ,  &  !: "meridional" Stream-Function
+      &                              surf_jk_r             !: inverse of the ocean "zonal" section surface
+#if defined key_diaeiv
+   REAL(wp), DIMENSION(jpj,jpk) ::   v_msf_eiv                  !: bolus "meridional" Stream-Function
+#endif
+   REAL(wp), DIMENSION(jpi,jpj) ::   abasin, pbasin, ibasin     !: return function value
 
    !! * Substitutions
@@ -78,5 +73,5 @@
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -94,17 +89,12 @@
       !!
       !! ** Action  : - p_fval: i-k-mean poleward flux of pva
-      !!
-      !!----------------------------------------------------------------------
-      !! * arguments
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   &
-         pva                         ! mask flux array at V-point
-
-      !! * local declarations
-      INTEGER  ::   ji, jj, jk        ! dummy loop arguments
-      INTEGER  ::   ijpj             ! ???
-      REAL(wp),DIMENSION(jpj) ::   &
-         p_fval                       ! function value
+      !!----------------------------------------------------------------------
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pva   ! mask flux array at V-point
+      !!
+      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
+      INTEGER                  ::   ijpj         ! ???
+      REAL(wp), DIMENSION(jpj) ::   p_fval       ! function value
       !!--------------------------------------------------------------------
-
+      !
       ijpj = jpj
       p_fval(:) = 0.e0
@@ -116,9 +106,8 @@
          END DO
       END DO
-
+      !
       IF( lk_mpp )   CALL mpp_sum( p_fval, ijpj )     !!bug  I presume
-
+      !
    END FUNCTION ptr_vj_3d
-
 
 
@@ -134,17 +123,12 @@
       !!
       !! ** Action  : - p_fval: i-k-mean poleward flux of pva
-      !!
-      !!----------------------------------------------------------------------
-      !! * arguments
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   &
-         pva                         ! mask flux array at V-point
-
-      !! * local declarations
-      INTEGER  ::   ji,jj             ! dummy loop arguments
-      INTEGER  ::   ijpj             ! ???
-      REAL(wp),DIMENSION(jpj) ::   &
-         p_fval                       ! function value
+      !!----------------------------------------------------------------------
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pva   ! mask flux array at V-point
+      !!
+      INTEGER                  ::   ji,jj    ! dummy loop arguments
+      INTEGER                  ::   ijpj     ! ???
+      REAL(wp), DIMENSION(jpj) ::   p_fval   ! function value
       !!--------------------------------------------------------------------
-      
+      ! 
       ijpj = jpj
       p_fval(:) = 0.e0
@@ -154,9 +138,8 @@
          END DO
       END DO
-
+      !
       IF( lk_mpp )   CALL mpp_sum( p_fval, ijpj )     !!bug  I presume
- 
-    END FUNCTION ptr_vj_2d
-
+      ! 
+   END FUNCTION ptr_vj_2d
 
 
@@ -171,22 +154,16 @@
       !!
       !! ** Action  : - p_fval: i-k-mean poleward flux of pva
-      !!
-      !!----------------------------------------------------------------------
-      !! * arguments
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   &
-         pva                         ! mask flux array at V-point
-
-      !! * local declarations
-      INTEGER  ::   ji, jj, jk        ! dummy loop arguments
-      INTEGER, DIMENSION (1) :: ish
-      INTEGER, DIMENSION (2) :: ish2
-      REAL(wp),DIMENSION(jpj*jpk) ::   &
-         zwork                        ! temporary vector for mpp_sum
-      REAL(wp),DIMENSION(jpj,jpk) ::   &
-         p_fval                       ! return function value
+      !!----------------------------------------------------------------------
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pva   ! mask flux array at V-point
+      !!
+      INTEGER                      ::   ji, jj, jk   ! dummy loop arguments
+      INTEGER , DIMENSION (1)      ::   ish
+      INTEGER , DIMENSION (2)      ::   ish2
+      REAL(wp), DIMENSION(jpj*jpk) ::   zwork        ! temporary vector for mpp_sum
+      REAL(wp), DIMENSION(jpj,jpk) ::   p_fval       ! return function value
       !!--------------------------------------------------------------------
- 
+      ! 
       p_fval(:,:) = 0.e0
-
+      !
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
@@ -197,5 +174,5 @@
          END DO
       END DO
-
+      !
       IF(lk_mpp) THEN
          ish(1) = jpj*jpk  ;  ish2(1) = jpj  ;  ish2(2) = jpk
@@ -204,6 +181,7 @@
          p_fval(:,:)= RESHAPE( zwork, ish2 )
       END IF
-
+      !
    END FUNCTION ptr_vjk
+
 
    FUNCTION ptr_vtjk( pva )   RESULT ( p_fval )
@@ -218,20 +196,14 @@
       !!
       !! ** Action  : - p_fval: i-k-mean poleward flux of pva
-      !!
-      !!----------------------------------------------------------------------
-      !! * arguments
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   &
-         pva                         ! mask flux array at V-point
- 
-      !! * local declarations
-      INTEGER  ::   ji, jj, jk        ! dummy loop arguments
-      INTEGER, DIMENSION (1) :: ish
-      INTEGER, DIMENSION (2) :: ish2
-      REAL(wp),DIMENSION(jpj*jpk) ::   &
-         zwork                        ! temporary vector for mpp_sum
-      REAL(wp),DIMENSION(jpj,jpk) ::   &
-         p_fval                       ! return function value
+      !!----------------------------------------------------------------------
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pva   ! mask flux array at V-point
+      !!
+      INTEGER                     ::   ji, jj, jk   ! dummy loop arguments
+      INTEGER, DIMENSION (1)      ::   ish
+      INTEGER, DIMENSION (2)      ::   ish2
+      REAL(wp),DIMENSION(jpj*jpk) ::   zwork        ! temporary vector for mpp_sum
+      REAL(wp),DIMENSION(jpj,jpk) ::   p_fval       ! return function value
       !!-------------------------------------------------------------------- 
-
+      !
       p_fval(:,:) = 0.e0
       DO jk = 1, jpkm1
@@ -251,5 +223,5 @@
          p_fval(:,:)= RESHAPE(zwork,ish2)
       END IF
-
+      !
    END FUNCTION ptr_vtjk
 
@@ -259,16 +231,10 @@
       !!                  ***  ROUTINE dia_ptr  ***
       !!----------------------------------------------------------------------
-      !! * Moudules used
-      USE ioipsl
-
-      !! * Argument
       INTEGER, INTENT(in) ::   kt   ! ocean time step index
-
-      !! * Local variables
-      INTEGER ::   jk,jj,ji               ! dummy loop
-      REAL(wp) ::    &
-         zsverdrup,  &              ! conversion from m3/s to Sverdrup
-         zpwatt,     &              ! conversion from W    to PW
-         zggram                     ! conversion from g    to Pg
+      !!
+      INTEGER  ::   jk, jj, ji               ! dummy loop
+      REAL(wp) ::   zsverdrup,  &              ! conversion from m3/s to Sverdrup
+         &          zpwatt,     &              ! conversion from W    to PW
+         &          zggram                     ! conversion from g    to Pg
 
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &
@@ -277,13 +243,5 @@
          vs_atl, vs_pac, vs_ind,           &
          zv_eiv
-      CHARACTER (len=32) ::   &
-         clnam = 'subbasins.nc'                
-      INTEGER ::  itime,inum,ipi,ipj,ipk       ! temporary integer
-      INTEGER, DIMENSION (1) ::   istep
-      REAL(wp) ::    zdate0,zsecond,zdt        ! temporary scalars
-      REAL(wp), DIMENSION(jpidta,jpjdta) ::   &
-         zlamt, zphit, zdta             ! temporary workspace (NetCDF read)
-      REAL(wp), DIMENSION(jpk) ::   &
-         zdept                          ! temporary workspace (NetCDF read)
+      INTEGER ::  inum       ! temporary logical unit
       !!----------------------------------------------------------------------
 
@@ -293,10 +251,4 @@
          zpwatt    = 1.e-15
          zggram    = 1.e-6
-         ipi       = jpidta
-         ipj       = jpjdta
-         ipk       = 1
-         itime     = 1
-         zsecond   = 0.e0
-         zdate0    = 0.e0
    
 # if defined key_diaeiv
@@ -315,50 +267,10 @@
          IF( ln_subbas ) THEN              ! Basins computation
 
-            IF( kt == nit000 ) THEN                ! load basin mask
-               itime = 1
-               ipi   = jpidta
-               ipj   = jpjdta
-               ipk   = 1
-               zdt   = 0.e0
-               istep = 0
-               clnam = 'subbasins.nc'
-
-               CALL flinopen(clnam,1,jpidta,1,jpjdta,.FALSE.,ipi,ipj, &
-                  &          ipk,zlamt,zphit,zdept,itime,istep,zdate0,zdt,inum)
-
-               ! get basins:
-               abasin (:,:) = 0.e0
-               pbasin (:,:) = 0.e0
-               ibasin (:,:) = 0.e0
-
-               ! Atlantic basin
-               CALL flinget(inum,'atlmsk',jpidta,jpjdta,1,itime,1,   &
-                  &         0,1,jpidta,1,jpjdta,zdta(:,:))
-               DO jj = 1, nlcj                                 ! interior values
-                  DO ji = 1, nlci
-                     abasin (ji,jj) = zdta( mig(ji), mjg(jj) )
-                  END DO
-               END DO
-
-               ! Pacific basin
-               CALL flinget(inum,'pacmsk',jpidta,jpjdta,1,itime,1,   &
-                  &         0,1,jpidta,1,jpjdta,zdta(:,:))
-               DO jj = 1, nlcj                                 ! interior values
-                  DO ji = 1, nlci
-                     pbasin (ji,jj) = zdta( mig(ji), mjg(jj) )
-                  END DO
-               END DO
-
-               ! Indian basin
-               CALL flinget(inum,'indmsk',jpidta,jpjdta,1,itime,1,   &
-                  &         0,1,jpidta,1,jpjdta,zdta(:,:))
-               DO jj = 1, nlcj                                 ! interior values
-                  DO ji = 1, nlci
-                     ibasin (ji,jj) = zdta( mig(ji), mjg(jj) )
-                  END DO
-               END DO
-
-               CALL flinclo(inum)
-
+            IF( kt == nit000 ) THEN                ! load sub-basin mask
+               CALL iom_open( 'subbasins', inum )
+               CALL iom_get( inum, jpdom_data, 'atlmsk', abasin )      ! Atlantic basin
+               CALL iom_get( inum, jpdom_data, 'pacmsk', pbasin )      ! Pacific basin
+               CALL iom_get( inum, jpdom_data, 'indmsk', ibasin )      ! Indian basin
+               CALL iom_close( inum )
             ENDIF
 
@@ -396,12 +308,12 @@
 #endif
          IF( ln_subbas ) THEN
-            v_msf_atl(:,:) = ptr_vjk( v_atl(:,:,:) ) 
-            v_msf_ipc(:,:) = ptr_vjk( v_ipc(:,:,:) ) 
-            ht_atl(:) = SUM(ptr_vjk( vt_atl(:,:,:)),2 )
-            ht_pac(:) = SUM(ptr_vjk( vt_pac(:,:,:)),2 )
-            ht_ind(:) = SUM(ptr_vjk( vt_ind(:,:,:)),2 )
-            st_atl(:) = SUM(ptr_vjk( vs_atl(:,:,:)),2 )
-            st_pac(:) = SUM(ptr_vjk( vs_pac(:,:,:)),2 )
-            st_ind(:) = SUM(ptr_vjk( vs_ind(:,:,:)),2 )
+            v_msf_atl(:,:) = ptr_vjk( v_atl (:,:,:) ) 
+            v_msf_ipc(:,:) = ptr_vjk( v_ipc (:,:,:) ) 
+            ht_atl(:) = SUM( ptr_vjk( vt_atl(:,:,:)), 2 )
+            ht_pac(:) = SUM( ptr_vjk( vt_pac(:,:,:)), 2 )
+            ht_ind(:) = SUM( ptr_vjk( vt_ind(:,:,:)), 2 )
+            st_atl(:) = SUM( ptr_vjk( vs_atl(:,:,:)), 2 )
+            st_pac(:) = SUM( ptr_vjk( vs_pac(:,:,:)), 2 )
+            st_ind(:) = SUM( ptr_vjk( vs_ind(:,:,:)), 2 )
          ENDIF
 
@@ -466,5 +378,5 @@
       ! Close the file
       IF( kt == nitend ) CALL histclo( numptr )
-
+      !
    END SUBROUTINE dia_ptr
 
@@ -475,5 +387,4 @@
       !!                   
       !! ** Purpose :   Initialization, namelist read
-      !!
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_1         ! temporary workspace
@@ -485,5 +396,4 @@
       REWIND ( numnam )
       READ   ( numnam, namptr )
-
 
       ! Control print
@@ -513,19 +423,12 @@
       !!
       !! ** Method  :   NetCDF file
-      !!
-      !!----------------------------------------------------------------------
-      !! * Arguments
+      !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-
-      !! * Save variables   
+      !!
       INTEGER, SAVE ::   nhoridz, ndepidzt, ndepidzw, ndex(1)
 
-      !! * Local variables
-      CHARACTER (len=40) ::   &
-         clhstnam, clop             ! temporary names
-      INTEGER ::   iline, it, ji    !
-      REAL(wp) ::   &
-         zsto, zout, zdt, zmax, &   ! temporary scalars
-         zjulian
+      CHARACTER (len=40)       ::   clhstnam, clop                   ! temporary names
+      INTEGER                  ::   iline, it, ji                    !
+      REAL(wp)                 ::   zsto, zout, zdt, zmax, zjulian   ! temporary scalars
       REAL(wp), DIMENSION(jpj) ::   zphi, zfoo
       !!----------------------------------------------------------------------
@@ -720,5 +623,5 @@
  
       ENDIF
-
+      !
    END SUBROUTINE dia_ptr_wri
 
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 508)
@@ -4,14 +4,23 @@
    !! Ocean dynamics:  surface pressure gradient trend
    !!======================================================================
-#if ( defined key_dynspg_flt && ! defined key_mpp_omp )   ||   defined key_esopa
+   !! History    8.0  !  98-05  (G. Roullet)  free surface
+   !!                 !  98-10  (G. Madec, M. Imbard)  release 8.2
+   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module
+   !!            " "  !  02-11  (C. Talandier, A-M Treguier) Open boundaries
+   !!            9.0  !  04-08  (C. Talandier) New trends organization
+   !!            " "  !  05-11  (V. Garnier) Surface pressure gradient organization
+   !!            " "  !  06-07  (S. Masson)  distributed restart using iom
+   !!----------------------------------------------------------------------
+#if ( defined key_dynspg_flt  && ! defined key_mpp_omp )  ||   defined key_esopa  
    !!----------------------------------------------------------------------
    !!   'key_dynspg_flt'                              filtered free surface
    !!   NOT 'key_mpp_omp'                          k-j-i loop (vector opt.)
+   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   dyn_spg_flt  : update the momentum trend with the surface pressure
    !!                  gradient in the filtered free surface case with
    !!                  vector optimization
-   !!----------------------------------------------------------------------
-   !! * Modules used
+   !!   flt_rst      : read/write the time-splitting restart fields in the ocean restart file
+   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers 
    USE dom_oce         ! ocean space and time domain 
@@ -21,4 +30,5 @@
    USE flxrnf          ! ocean runoffs
    USE sol_oce         ! ocean elliptic solver
+   USE solver          ! solver initialization
    USE solpcg          ! preconditionned conjugate gradient solver
    USE solsor          ! Successive Over-relaxation solver
@@ -32,12 +42,13 @@
    USE cla_dynspg      ! cross land advection
    USE prtctl          ! Print control
-   USE in_out_manager  ! I/O manager
    USE solmat          ! matrix construction for elliptic solvers
    USE agrif_opa_interp
+   USE in_out_manager  ! I/O manager
+   USE iom
+   USE restart         ! only for lrst_oce
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Accessibility
    PUBLIC dyn_spg_flt  ! routine called by step.F90
 
@@ -48,5 +59,5 @@
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    !!----------------------------------------------------------------------
 
@@ -96,28 +107,15 @@
       !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
       !!
-      !! References :
-      !!      Roullet and Madec 1999, JGR.
-      !!
-      !! History :
-      !!        !  98-05 (G. Roullet)  Original code
-      !!        !  98-10 (G. Madec, M. Imbard)  release 8.2
-      !!   8.5  !  02-08 (G. Madec)  F90: Free form and module
-      !!        !  02-11 (C. Talandier, A-M Treguier) Open boundaries
-      !!   9.0  !  04-08 (C. Talandier) New trends organization
-      !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
+      !! References : Roullet and Madec 1999, JGR.
       !!---------------------------------------------------------------------
-      !! * Arguments
       INTEGER, INTENT( in )  ::   kt         ! ocean time-step index
-      INTEGER, INTENT( out ) ::   kindic     ! solver convergence flag
-                                             ! if the solver doesn t converge
-                                             ! the flag is < 0
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk               ! dummy loop indices
-      REAL(wp) ::                         & 
-         z2dt, z2dtg, zraur, znugdt,      &  ! temporary scalars
-         znurau, zssha, zgcb, zbtd,       &  !   "          "
-         ztdgu, ztdgv                        !   "          "
+      INTEGER, INTENT( out ) ::   kindic     ! solver convergence flag (<0 if not converge)
+      !!                                   
+      INTEGER  ::   ji, jj, jk                          ! dummy loop indices
+      REAL(wp) ::   z2dt, z2dtg, zraur, znugdt,      &  ! temporary scalars
+         &          znurau, zssha, zgcb, zbtd,       &  !   "          "
+         &          ztdgu, ztdgv                        !   "          "
       !!----------------------------------------------------------------------
-
+      !
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
@@ -128,4 +126,9 @@
          spgu(:,:) = 0.e0                     ! surface pressure gradient (i-direction)
          spgv(:,:) = 0.e0                     ! surface pressure gradient (j-direction)
+         CALL solver_init( nit000 )           ! Elliptic solver initialisation
+
+         ! read filtered free surface arrays in restart file
+         CALL flt_rst( nit000, 'READ' )       ! read or initialize the following fields:
+         !                                    ! gcx, gcxb, sshb, sshn
       ENDIF
 
@@ -168,13 +171,9 @@
 
 #if defined key_obc
-      ! Update velocities on each open boundary with the radiation algorithm
-      CALL obc_dyn( kt )
-      ! Correction of the barotropic componant velocity to control the volume of the system
-      CALL obc_vol( kt )
+      CALL obc_dyn( kt )      ! Update velocities on each open boundary with the radiation algorithm
+      CALL obc_vol( kt )      ! Correction of the barotropic componant velocity to control the volume of the system
 #endif
 #if defined key_agrif
-      ! Update velocities on each coarse/fine interfaces
-
-      CALL Agrif_dyn( kt )
+      CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces 
 #endif
 #if defined key_orca_r2
@@ -243,12 +242,12 @@
       IF( .NOT. AGRIF_ROOT() ) THEN
          ! add contribution of gradient of after barotropic transport divergence 
-         IF( (nbondi == -1) .OR. (nbondi == 2) ) gcb(3,:) = gcb(3,:) &
-            &            -znugdt * z2dt*laplacu(2,:)*gcdprc(3,:)*hu(2,:)*e2u(2,:)
-         IF( (nbondi ==  1) .OR. (nbondi == 2) )  gcb(nlci-2,:) = gcb(nlci-2,:) &
-            &           +znugdt * z2dt*laplacu(nlci-2,:)*gcdprc(nlci-2,:)*hu(nlci-2,:)*e2u(nlci-2,:)
-         IF( (nbondj == -1) .OR. (nbondj == 2) ) gcb(:,3) = gcb(:,3) &
-            &           -znugdt * z2dt*laplacv(:,2)*gcdprc(:,3)*hv(:,2)*e1v(:,2)
-         IF( (nbondj ==  1) .OR. (nbondj == 2) )  gcb(:,nlcj-2) = gcb(:,nlcj-2) &
-            &           +znugdt * z2dt*laplacv(:,nlcj-2)*gcdprc(:,nlcj-2)*hv(:,nlcj-2)*e1v(:,nlcj-2)
+         IF( nbondi == -1 .OR. nbondi == 2 )   gcb(3     ,:) =   &
+            &    gcb(3     ,:) - znugdt * z2dt * laplacu(2     ,:) * gcdprc(3     ,:) * hu(2     ,:) * e2u(2     ,:)
+         IF( nbondi ==  1 .OR. nbondi == 2 )   gcb(nlci-2,:) =   &
+            &    gcb(nlci-2,:) + znugdt * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:)
+         IF( nbondj == -1 .OR. nbondj == 2 )   gcb(:     ,3) =   &
+            &    gcb(:,3     ) - znugdt * z2dt * laplacv(:,2     ) * gcdprc(:,3     ) * hv(:,2     ) * e1v(:,2     )
+         IF( nbondj ==  1 .OR. nbondj == 2 )   gcb(:,nlcj-2) =   &
+            &    gcb(:,nlcj-2) + znugdt * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2)
       ENDIF
 #endif
@@ -263,5 +262,5 @@
       epsr = eps * eps * rnorme
       ncut = 0
-      ! if rnorme is 0, the solution is 0, the solver isn't called
+      ! if rnorme is 0, the solution is 0, the solver is not called
       IF( rnorme == 0.e0 ) THEN
          gcx(:,:) = 0.e0
@@ -313,8 +312,8 @@
       IF( .NOT. Agrif_Root() ) THEN
          ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
-         IF( (nbondi == -1) .OR. (nbondi == 2) ) spgu(2,:) = znugdt * z2dt * laplacu(2,:) * umask(2,:,1)
-         IF( (nbondi ==  1) .OR. (nbondi == 2) ) spgu(nlci-2,:) = znugdt * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
-         IF( (nbondj == -1) .OR. (nbondj == 2) ) spgv(:,2) = znugdt * z2dt * laplacv(:,2) * vmask(:,2,1)
-         IF( (nbondj ==  1) .OR. (nbondj == 2) ) spgv(:,nlcj-2) = znugdt * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
+         IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2     ,:) = znugdt * z2dt * laplacu(2     ,:) * umask(2     ,:,1)
+         IF( nbondi ==  1 .OR. nbondi == 2 ) spgu(nlci-2,:) = znugdt * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
+         IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2     ) = znugdt * z2dt * laplacv(:,2     ) * vmask(:     ,2,1)
+         IF( nbondj ==  1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = znugdt * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
       ENDIF
 #endif      
@@ -323,5 +322,5 @@
       !     ( c a u t i o n : (ua,va) here are the after velocity not the
       !                       trend, the leap-frog time stepping will not
-      !                       be done in dynnxt.F routine)
+      !                       be done in dynnxt.F90 routine)
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
@@ -332,5 +331,4 @@
          END DO
       END DO
-
 
       ! Sea surface elevation time stepping
@@ -358,8 +356,52 @@
       ENDIF
 
-      !                       ! print sum trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg  - ssh: ', mask1=tmask )
-
+      ! write filtered free surface arrays in restart file
+      ! --------------------------------------------------
+      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
+
+      ! print sum trends (used for debugging)
+      IF(ln_ctl)     CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg  - ssh: ', mask1=tmask )
+      !
    END SUBROUTINE dyn_spg_flt
+
+
+   SUBROUTINE flt_rst( kt, cdrw )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE ts_rst  ***
+     !!
+     !! ** Purpose : Read or write filtered free surface arrays in restart file
+     !!----------------------------------------------------------------------
+     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+     !!----------------------------------------------------------------------
+
+     IF( TRIM(cdrw) == 'READ' ) THEN
+        IF( iom_varid( numror, 'gcx' ) > 0 ) THEN
+! Caution : extra-hallow
+! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
+           CALL iom_get( numror, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) )
+           CALL iom_get( numror, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) )
+           CALL iom_get( numror, jpdom_local, 'sshb', sshb(:,:)         )
+           CALL iom_get( numror, jpdom_local, 'sshn', sshn(:,:)         )
+           IF( neuler == 0 ) THEN
+              sshb(:,:) = sshn(:,:)
+              gcxb(:,:) = gcx (:,:)
+           ENDIF
+        ELSE
+           gcx (:,:) = 0.e0
+           gcxb(:,:) = 0.e0
+           sshb(:,:) = 0.e0
+           sshn(:,:) = 0.e0
+        ENDIF
+     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
+! Caution : extra-hallow
+! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
+        CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) )
+        CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
+        CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:)         )
+        CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:)         )
+     ENDIF
+     !
+   END SUBROUTINE flt_rst
 
 #else
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90	(revision 508)
@@ -25,4 +25,5 @@
    USE solfet          ! FETI solver
    USE solsor_e        ! Successive Over-relaxation solver with MPP optimization
+   USE solver
    USE obc_oce         ! Lateral open boundary condition
    USE obcdyn          ! ocean open boundary condition (obc_dyn routines)
@@ -35,4 +36,6 @@
    USE solmat          ! matrix construction for elliptic solvers
    USE agrif_opa_interp
+   USE restart         ! only for lrst_oce
+   USE iom
 
    IMPLICIT NONE
@@ -112,7 +115,14 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~  (free surface constant volume, autotasking case)'
 
-         ! set to zero free surface specific arrays 
-         spgu(:,:) = 0.e0      ! surface pressure gradient (i-direction)
-         spgv(:,:) = 0.e0      ! surface pressure gradient (j-direction)
+         ! set to zero free surface specific arrays
+         spgu(:,:) = 0.e0                     ! surface pressure gradient (i-direction)
+         
+         spgv(:,:) = 0.e0                     ! surface pressure gradient (j-direction)
+         
+         CALL solver_init( nit000 )           ! Elliptic solver initialisation
+
+         ! read filtered free surface arrays in restart file
+         CALL flt_rst( nit000, 'READ' )       ! read or initialize the following fields:
+         !                                    ! gcx, gcxb, sshb, sshn
       ENDIF
 
@@ -354,4 +364,8 @@
       CALL lbc_lnk( sshn, 'T', 1. )
 
+      ! write filtered free surface arrays in restart file
+      ! --------------------------------------------------
+      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
+
       IF(ln_ctl) THEN         ! print sum trends (used for debugging)
          CALL prt_ctl( tab3d_1=ua  , clinfo1=' spg  - Ua : ', mask1=umask,   &
@@ -362,4 +376,43 @@
    END SUBROUTINE dyn_spg_flt_jki
 
+   SUBROUTINE flt_rst( kt, cdrw )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE ts_rst  ***
+     !!
+     !! ** Purpose : Read or write filtered free surface arrays in restart file
+     !!----------------------------------------------------------------------
+     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+     !!----------------------------------------------------------------------
+
+     IF( TRIM(cdrw) == 'READ' ) THEN
+        IF( iom_varid( numror, 'gcx' ) > 0 ) THEN
+! Caution : extra-hallow
+! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
+           CALL iom_get( numror, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) )
+           CALL iom_get( numror, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) )
+           CALL iom_get( numror, jpdom_local, 'sshb', sshb(:,:)         )
+           CALL iom_get( numror, jpdom_local, 'sshn', sshn(:,:)         )
+           IF( neuler == 0 ) THEN
+              sshb(:,:) = sshn(:,:)
+              gcxb(:,:) = gcx (:,:)
+           ENDIF
+        ELSE
+           gcx (:,:) = 0.e0
+           gcxb(:,:) = 0.e0
+           sshb(:,:) = 0.e0
+           sshn(:,:) = 0.e0
+        ENDIF
+     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
+! Caution : extra-hallow
+! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
+        CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) )
+        CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
+        CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:)         )
+        CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:)         )
+     ENDIF
+     !
+   END SUBROUTINE flt_rst
+
 #else
    !!----------------------------------------------------------------------
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_rl.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_rl.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_rl.F90	(revision 508)
@@ -4,4 +4,17 @@
    !! Ocean dynamics:  surface pressure gradient trend
    !!======================================================================
+   !! History :  7.0  !  96-05  (G. Madec, M. Imbard, M. Guyon)  rewitting in 1
+   !!                           routine, without pointers, and with the same matrix
+   !!                           for sor and pcg, mpp exchange, and symmetric conditions
+   !!            " "  !  96-07  (A. Weaver)  Euler forward step
+   !!            " "  !  96-11  (A. Weaver)  correction to preconditioning:
+   !!            8.0  !  98-02  (M. Guyon)  FETI method
+   !!            " "  !  97-09  (J.-M. Molines)  Open boundaries
+   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module
+   !!                 !  02-11  (C. Talandier, A-M Treguier) Open boundaries
+   !!            9.0  !  04-08  (C. Talandier)  New trends organization
+   !!            " "  !  05-11  (V. Garnier) Surface pressure gradient organization
+   !!            " "  !  06-07  (S. Masson)  distributed restart using iom
+   !!---------------------------------------------------------------------
 #if   defined key_dynspg_rl   ||   defined key_esopa
    !!----------------------------------------------------------------------
@@ -10,6 +23,6 @@
    !!   dyn_spg_rl   : update the momentum trend with the surface pressure
    !!                  for the rigid-lid case.
-   !!----------------------------------------------------------------------
-   !! * Modules used
+   !!   rl_rst       : read/write the rigid-lid restart fields in the ocean restart file
+   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -19,4 +32,5 @@
    USE zdf_oce         ! ocean vertical physics
    USE sol_oce         ! ocean elliptic solver
+   USE solver          ! solver initialization
    USE solpcg          ! preconditionned conjugate gradient solver
    USE solsor          ! Successive Over-relaxation solver
@@ -28,9 +42,10 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
+   USE iom
+   USE restart         ! only for lrst_oce
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Accessibility
    PUBLIC dyn_spg_rl   ! called by step.F90
 
@@ -42,5 +57,5 @@
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    !!----------------------------------------------------------------------
 
@@ -76,26 +91,8 @@
       !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
       !!
-      !! References :
-      !!      Madec et al. 1988, ocean modelling, issue 78, 1-6.
-      !!
-      !! History :
-      !!        !  96-05  (G. Madec, M. Imbard, M. Guyon)  rewitting in 1
-      !!                  routine, without pointers, and with the same matrix
-      !!                  for sor and pcg, mpp exchange, and symmetric conditions
-      !!        !  96-07  (A. Weaver)  Euler forward step
-      !!        !  96-11  (A. Weaver)  correction to preconditioning:
-      !!        !  98-02  (M. Guyon)  FETI method
-      !!        !  98-05  (G. Roullet)  free surface
-      !!        !  97-09  (J.-M. Molines)  Open boundaries
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
-      !!        !  02-11  (C. Talandier, A-M Treguier) Open boundaries
-      !!   9.0  !  04-08  (C. Talandier)  New trends organization
-      !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
+      !! References : Madec et al. 1988, ocean modelling, issue 78, 1-6.
       !!---------------------------------------------------------------------
-      !! * Arguments
       INTEGER, INTENT( in  ) ::   kt       ! ocean time-step index
-      INTEGER, INTENT( out ) ::   kindic   ! solver flag, take a negative value
-      !                                    ! when the solver doesnot converge
-      !! * Local declarations
+      INTEGER, INTENT( out ) ::   kindic   ! solver flag (<0 when the solver does not converge)
       INTEGER ::   ji, jj, jk    ! dummy loop indices
       REAL(wp) ::  zbsfa, zgcx, z2dt
@@ -114,32 +111,23 @@
 
          ! set to zero rigid-lid specific arrays
-         spgu(:,:) = 0.e0      ! surface pressure gradient (i-direction) 
-         spgv(:,:) = 0.e0      ! surface pressure gradient (j-direction)
-      ENDIF
-
-      ! 0. Initializations:
-      ! -------------------
-# if defined key_obc
-      ! space index on boundary arrays
-      ib = 1
-      ibm = 2
-      ibm2 = 3
-      ! time index on boundary arrays
-      it = 1
-      itm = 2
-      itm2 = 3
-# endif
-
+         spgu(:,:) = 0.e0                   ! surface pressure gradient (i-direction) 
+         spgv(:,:) = 0.e0                   ! surface pressure gradient (j-direction)
+
+         CALL solver_init( nit000 )         ! Elliptic solver initialisation
+
+         ! read rigid lid arrays in restart file
+         CALL rl_rst( nit000, 'READ' )      ! read or initialize the following fields:
+         !                                  ! gcx, gcxb, bsfb, bsfn, bsfd
+      ENDIF
+
+      !  Vertically averaged momentum trend
+      ! ------------------------------------
       !                                                ! ===============
       DO jj = 2, jpjm1                                 !  Vertical slab
          !                                             ! ===============
-
-         ! 1. Vertically averaged momentum trend
-         ! -------------------------------------
-         ! initialization to zero
-         spgu(:,jj) = 0.
+         
+         spgu(:,jj) = 0.                          ! initialization to zero
          spgv(:,jj) = 0.
-         ! vertical sum
-         DO jk = 1, jpk
+         DO jk = 1, jpk                           ! vertical sum
             DO ji = 2, jpim1
                spgu(ji,jj) = spgu(ji,jj) + ua(ji,jj,jk) * fse3u(ji,jj,jk) * umask(ji,jj,jk)
@@ -147,6 +135,5 @@
             END DO 
          END DO 
-         ! divide by the depth
-         spgu(:,jj) = spgu(:,jj) * hur(:,jj)
+         spgu(:,jj) = spgu(:,jj) * hur(:,jj)      ! divide by the depth 
          spgv(:,jj) = spgv(:,jj) * hvr(:,jj)
 
@@ -155,27 +142,20 @@
       !                                                ! ===============
 
-      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-
       ! Boundary conditions on (spgu,spgv)
       CALL lbc_lnk( spgu, 'U', -1. )
       CALL lbc_lnk( spgv, 'V', -1. )
 
-      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-
-      ! 2. Barotropic streamfunction trend (bsfd)
+      !  Barotropic streamfunction trend (bsfd)
       ! ----------------------------------
-
-      ! Curl of the vertically averaged velocity
-      DO jj = 2, jpjm1
+      DO jj = 2, jpjm1                            ! Curl of the vertically averaged velocity 
          DO ji = 2, jpim1
             gcb(ji,jj) = -gcdprc(ji,jj)   &
-                       *(  ( e2v(ji+1,jj  )*spgv(ji+1,jj  ) - e2v(ji,jj)*spgv(ji,jj) )   &
-                          -( e1u(ji  ,jj+1)*spgu(ji  ,jj+1) - e1u(ji,jj)*spgu(ji,jj) )   ) 
+               &       *(  ( e2v(ji+1,jj  )*spgv(ji+1,jj  ) - e2v(ji,jj)*spgv(ji,jj) )   &
+               &          -( e1u(ji  ,jj+1)*spgu(ji  ,jj+1) - e1u(ji,jj)*spgu(ji,jj) )   ) 
          END DO
       END DO
 
 # if defined key_obc
-      ! Open boundary contribution
-      DO jj = 2, jpjm1
+      DO jj = 2, jpjm1                            ! Open boundary contribution 
          DO ji = 2, jpim1
            gcb(ji,jj) = gcb(ji,jj) - gcdprc(ji,jj) * gcbob(ji,jj)
@@ -198,6 +178,4 @@
       ! applied the lateral boundary conditions
       IF( nsolv == 4)   CALL lbc_lnk_e( gcb, c_solver_pt, 1. )   
-
-      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
 
       ! Relative precision (computation on one processor)
@@ -234,16 +212,11 @@
       ENDIF
 
-
       ! bsf trend update
       ! ----------------
-
       bsfd(1:nlci,1:nlcj) = gcx(1:nlci,1:nlcj)
-
       
       ! update bsf trend with islands trend
       ! -----------------------------------
-
       IF( lk_isl )   CALL isl_dyn_spg         ! update bsfd
-
 
 # if defined key_obc
@@ -337,5 +310,4 @@
       ! 4. Barotropic stream function and array swap
       ! --------------------------------------------
-
       ! Leap-frog time scheme, time filter and array swap
       IF( neuler == 0 .AND. kt == nit000 ) THEN
@@ -362,4 +334,11 @@
 
 # if defined key_obc
+      ib   = 1      ! space index on boundary arrays
+      ibm  = 2
+      ibm2 = 3
+      it   = 1      ! time index on boundary arrays
+      itm  = 2
+      itm2 = 3
+
       ! Swap of boundary arrays
       IF( lp_obc_east ) THEN
@@ -499,12 +478,8 @@
       ENDIF
 # endif
-      !
-      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
       
       !  add the surface pressure trend to the general trend
       ! -----------------------------------------------------
-      
       DO jj = 2, jpjm1
-
          ! update the surface pressure gradient with the barotropic trend
          DO ji = 2, jpim1
@@ -519,8 +494,53 @@
             END DO
          END DO
-
-      END DO
+      END DO
+
+      ! write rigid lid arrays in restart file
+      ! --------------------------------------
+      IF( lrst_oce ) CALL rl_rst( kt, 'WRITE' )
 
    END SUBROUTINE dyn_spg_rl
+
+
+   SUBROUTINE rl_rst( kt, cdrw )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE rl_rst  ***
+     !!
+     !! ** Purpose : Read or write rigid-lid arrays in ocean restart file
+     !!----------------------------------------------------------------------
+     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+     !!----------------------------------------------------------------------
+     !
+     IF( TRIM(cdrw) == 'READ' ) THEN
+        IF( iom_varid( numror, 'gcx' ) > 0 ) THEN
+     ! Caution : extra-hallow
+     ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
+           CALL iom_get( numror, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) )
+           CALL iom_get( numror, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) )
+           CALL iom_get( numror, jpdom_local, 'bsfb', bsfb(:,:)         )
+           CALL iom_get( numror, jpdom_local, 'bsfn', bsfn(:,:)         )
+           CALL iom_get( numror, jpdom_local, 'bsfd', bsfd(:,:)         )
+           IF( neuler == 0 ) THEN
+              gcxb(:,:) = gcx (:,:)
+              bsfb(:,:) = bsfn(:,:)
+           ENDIF
+        ELSE
+           gcx (:,:) = 0.e0
+           gcxb(:,:) = 0.e0
+           bsfb(:,:) = 0.e0
+           bsfn(:,:) = 0.e0
+           bsfd(:,:) = 0.e0
+        ENDIF
+     ELSEIF(  TRIM(cdrw) == 'WRITE' ) THEN
+        ! Caution : extra-hallow, gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
+        CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
+        CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
+        CALL iom_rstput( kt, nitrst, numrow, 'bsfb', bsfb(:,:)         )
+        CALL iom_rstput( kt, nitrst, numrow, 'bsfn', bsfn(:,:)         )
+        CALL iom_rstput( kt, nitrst, numrow, 'bsfd', bsfd(:,:)         )
+     ENDIF
+     !
+   END SUBROUTINE rl_rst
 
 #else
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 508)
@@ -4,4 +4,8 @@
    !! Ocean dynamics:  surface pressure gradient trend
    !!======================================================================
+   !! History :   9.0  !  04-12  (L. Bessieres, G. Madec)  Original code
+   !!             " "  !  05-11  (V. Garnier, G. Madec)  optimization
+   !!             9.0  !  06-08  (S. Masson)  distributed restart using iom
+   !!---------------------------------------------------------------------
 #if ( defined key_dynspg_ts && ! defined key_mpp_omp ) ||   defined key_esopa
    !!----------------------------------------------------------------------
@@ -9,6 +13,8 @@
    !!   NOT 'key_mpp_omp'                          k-j-i loop (vector opt.)
    !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
    !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
    !!                 splitting scheme and add to the general trend 
+   !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file
    !!----------------------------------------------------------------------
    !! * Modules used
@@ -27,10 +33,15 @@
    USE dynspg_oce      ! surface pressure gradient variables
    USE in_out_manager  ! I/O manager
+   USE iom
+   USE restart         ! only for lrst_oce
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Accessibility
    PUBLIC dyn_spg_ts  ! routine called by step.F90
+
+   REAL(wp), DIMENSION(jpi,jpj) ::  ftnw, ftne,   &  ! triad of coriolis parameter
+      &                             ftsw, ftse       ! (only used with een vorticity scheme)
+
 
    !! * Substitutions
@@ -74,12 +85,6 @@
       !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
       !!
-      !! References :
-      !!   Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
-      !!
-      !! History :
-      !!   9.0  !  04-12  (L. Bessieres, G. Madec)  Original code
-      !!        !  05-11  (V. Garnier, G. Madec)  optimization
+      !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
       !!---------------------------------------------------------------------
-      !! * Arguments
       INTEGER, INTENT( in )  ::   kt           ! ocean time-step index
 
@@ -97,6 +102,4 @@
          zsshb_e, zub_e, zvb_e,             &  !     "        "
          zun_e, zvn_e                          !     "        "
-      REAL(wp), DIMENSION(jpi,jpj),SAVE ::  &
-         ztnw, ztne, ztsw, ztse
       !!----------------------------------------------------------------------
 
@@ -109,5 +112,5 @@
 
       IF( kt == nit000 ) THEN
-
+         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
@@ -115,25 +118,7 @@
          IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', FLOOR( 2*rdt/rdtbt )
 
-         IF( .NOT. ln_rstart ) THEN
-            ! initialize barotropic specific arrays
-            sshb_b(:,:) = sshb(:,:)
-            sshn_b(:,:) = sshn(:,:)
-            un_b(:,:)   = 0.e0
-            vn_b(:,:)   = 0.e0
-            ! vertical sum
-            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
-               DO jk = 1, jpkm1
-                  DO ji = 1, jpij
-                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
-                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
-                  END DO
-               END DO
-            ELSE                             ! No  vector opt.
-               DO jk = 1, jpkm1
-                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
-                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
-               END DO
-            ENDIF
-         ENDIF
+         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields:
+         !                               ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b
+
          ssha_e(:,:) = sshn(:,:)
          ua_e(:,:)   = un_b(:,:)
@@ -141,17 +126,17 @@
 
          IF( ln_dynvor_een ) THEN
-            ztne(1,:) = 0.e0   ;   ztnw(1,:) = 0.e0   ;   ztse(1,:) = 0.e0   ;   ztsw(1,:) = 0.e0
+            ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0
             DO jj = 2, jpj
                DO ji = fs_2, jpi   ! vector opt.
-                  ztne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
-                  ztnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
-                  ztse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3.
-                  ztsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3.
+                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
+                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
+                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3.
+                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3.
                END DO
             END DO
          ENDIF
-
+         !
       ENDIF
-    
+
       ! Local constant initialization
       ! --------------------------------
@@ -216,5 +201,5 @@
             END DO
          END DO
-
+         !
       ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
          DO jj = 2, jpjm1
@@ -228,5 +213,5 @@
             END DO
          END DO
-
+         !
       ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
          zfac25 = 0.25
@@ -241,5 +226,5 @@
             END DO
          END DO
-
+         !
       ENDIF
 
@@ -300,5 +285,4 @@
       DO jit = 1, icycle                                   !  sub-time-step loop  !
          !                                                 ! ==================== !
-
          z2dt_e = 2. * rdtbt
          IF ( jit == 1 )   z2dt_e = rdtbt
@@ -360,5 +344,5 @@
                END DO
             END DO
-
+            !
          ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
             DO jj = 2, jpjm1
@@ -379,5 +363,5 @@
                END DO
             END DO
-
+            !
          ELSEIF ( ln_dynvor_een ) THEN                    ! energy and enstrophy conserving scheme
             zfac25 = 0.25
@@ -397,4 +381,5 @@
                END DO
             END DO
+            ! 
          ENDIF
 
@@ -504,9 +489,72 @@
       END DO
 
-      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
-         CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
+      ! write filtered free surface arrays in restart file
+      ! --------------------------------------------------
+      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
+
+      ! print sum trends (used for debugging)
+      IF( ln_ctl )     CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
+      !
+   END SUBROUTINE dyn_spg_ts
+
+
+   SUBROUTINE ts_rst( kt, cdrw )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE ts_rst  ***
+      !!
+      !! ** Purpose : Read or write time-splitting arrays in restart file
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+      !
+      INTEGER ::  ji, jk        ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN
+         IF( iom_varid( numror, 'sshn' ) > 0 ) THEN
+            CALL iom_get( numror, jpdom_local, 'sshb'  , sshb(:,:)   )
+            CALL iom_get( numror, jpdom_local, 'sshn'  , sshn(:,:)   )
+            IF( neuler == 0 ) sshb(:,:) = sshn(:,:)
+         ELSE
+            sshb(:,:) = 0.e0
+            sshn(:,:) = 0.e0
+         ENDIF
+         IF( iom_varid( numror, 'sshn_b' ) > 0 ) THEN
+            CALL iom_get( numror, jpdom_local, 'sshb_b', sshb_b(:,:) )   ! free surface issued
+            CALL iom_get( numror, jpdom_local, 'sshn_b', sshn_b(:,:) )   ! from time-splitting loop
+            CALL iom_get( numror, jpdom_local, 'un_b'  , un_b  (:,:) )   ! horizontal transports issued
+            CALL iom_get( numror, jpdom_local, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
+            IF( neuler == 0 ) sshb_b(:,:) = sshn_b(:,:)
+         ELSE
+            sshb_b(:,:) = sshb(:,:)
+            sshn_b(:,:) = sshn(:,:)
+            un_b  (:,:) = 0.e0
+            vn_b  (:,:) = 0.e0
+            ! vertical sum
+            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
+               DO jk = 1, jpkm1
+                  DO ji = 1, jpij
+                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
+                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
+                  END DO
+               END DO
+            ELSE                             ! No  vector opt.
+               DO jk = 1, jpkm1
+                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
+                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
+               END DO
+            ENDIF
+         ENDIF
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
+         CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb  (:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn  (:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) )   ! free surface issued
+         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) )   ! from barotropic loop
+         CALL iom_rstput( kt, nitrst, numrow, 'un_b'  , un_b  (:,:) )   ! horizontal transports issued
+         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
       ENDIF
-      
-   END SUBROUTINE dyn_spg_ts
+      !
+   END SUBROUTINE ts_rst
+
 #else
    !!----------------------------------------------------------------------
Index: trunk/NEMO/OPA_SRC/SOL/solmat.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SOL/solmat.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/SOL/solmat.F90	(revision 508)
@@ -4,4 +4,16 @@
    !! solver       : construction of the matrix 
    !!======================================================================
+   !! History :   1.0  !  88-04  (G. Madec)  Original code
+   !!                  !  93-03  (M. Guyon)  symetrical conditions
+   !!                  !  93-06  (M. Guyon)  suppress pointers
+   !!                  !  96-05  (G. Madec)  merge sor and pcg formulations
+   !!                  !  96-11  (A. Weaver)  correction to preconditioning
+   !!                  !  98-02  (M. Guyon)  FETI method
+   !!             8.5  !  02-08  (G. Madec)  F90: Free form
+   !!                  !  02-11  (C. Talandier, A-M. Treguier) Free surface & Open boundaries
+   !!             9.0  !  05-09  (R. Benshila)  add sol_exd for extra outer halo
+   !!             9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
+   !!             9.0  !  06-07  (S. Masson)  distributed restart using iom
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -29,5 +41,5 @@
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    !!----------------------------------------------------------------------
 
@@ -55,16 +67,4 @@
       !!              - gcdmat : preconditioning matrix (diagonal elements)
       !!              - gcdprc : inverse of the preconditioning matrix
-      !!
-      !! History :
-      !!   1.0  !  88-04  (G. Madec)  Original code
-      !!        !  91-11  (G. Madec)
-      !!        !  93-03  (M. Guyon)  symetrical conditions
-      !!        !  93-06  (M. Guyon)  suppress pointers
-      !!        !  96-05  (G. Madec)  merge sor and pcg formulations
-      !!        !  96-11  (A. Weaver)  correction to preconditioning
-      !!        !  98-02  (M. Guyon)  FETI method
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form
-      !!        !  02-11  (C. Talandier, A-M. Treguier) Free surface & Open boundaries
-      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
       !!----------------------------------------------------------------------
       !! * Arguments
Index: trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 508)
@@ -5,67 +5,75 @@
    !!                 turbulent closure parameterization
    !!=====================================================================
+   !! History :   6.0  !  91-03  (b. blanke)  Original code
+   !!             7.0  !  91-11  (G. Madec)   bug fix
+   !!             7.1  !  92-10  (G. Madec)   new mixing length and eav
+   !!             7.2  !  93-03  (M. Guyon)   symetrical conditions
+   !!             7.3  !  94-08  (G. Madec, M. Imbard)   npdl flag
+   !!             7.5  !  96-01  (G. Madec)   s-coordinates
+   !!             8.0  !  97-07  (G. Madec)   lbc
+   !!             8.1  !  99-01  (E. Stretta) new option for the mixing length
+   !!             8.5  !  02-06  (G. Madec) add zdf_tke_init routine
+   !!             8.5  !  02-08  (G. Madec)  ri_c and Free form, F90
+   !!             9.0  !  04-10  (C. Ethe )  1D configuration
+   !!             9.0  !  02-08  (G. Madec)  autotasking optimization
+   !!             9.0  !  06-07  (S. Masson)  distributed restart using iom
+   !!----------------------------------------------------------------------
 #if defined key_zdftke   ||   defined key_esopa
    !!----------------------------------------------------------------------
-   !!   'key_zdftke'                                             TKE scheme
+   !!   'key_zdftke'                                   TKE vertical physics
+   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   zdf_tke      : update momentum and tracer Kz from a tke scheme
    !!   zdf_tke_init : initialization, namelist read, and parameters control
+   !!   tke_rst      : read/write tke restart in ocean restart file
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and active tracers 
    USE dom_oce         ! ocean space and time domain
    USE zdf_oce         ! ocean vertical physics
-   USE in_out_manager  ! I/O manager
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE phycst          ! physical constants
    USE taumod          ! surface stress
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
+   USE in_out_manager  ! I/O manager
+   USE iom
+   USE restart         ! only for lrst_oce
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC zdf_tke        ! routine called in step module
-   PUBLIC zdf_tke_init   ! routine called in zdftke_jki module
-
-   !! * Share Module variables
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .TRUE.    !: TKE vertical mixing flag
-   LOGICAL, PUBLIC ::         & !!: ** tke namelist (namtke) **
-     ln_rstke = .FALSE.          !: =T restart with tke from a run without tke with 
-     !                           !  a none zero initial value for en
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
-      en                         !: now turbulent kinetic energy
-
-   INTEGER, PUBLIC ::         & !!! ** tke namelist (namtke) **
-      nitke = 50 ,            &  ! number of restart iterative loops
-      nmxl  =  2 ,            &  ! = 0/1/2/3 flag for the type of mixing length used
-      npdl  =  1 ,            &  ! = 0/1/2 flag on prandtl number on vert. eddy coeff.
-      nave  =  1 ,            &  ! = 0/1 flag for horizontal average on avt, avmu, avmv
-      navb  =  0                 ! = 0/1 flag for constant or profile background avt
-   REAL(wp), PUBLIC ::        & !!! ** tke namlist (namtke) **
-      ediff = 0.1_wp       ,  &  ! coeff. for vertical eddy coef.; avt=ediff*mxl*sqrt(e)
-      ediss = 0.7_wp       ,  &  ! coef. of the Kolmogoroff dissipation 
-      ebb   = 3.75_wp      ,  &  ! coef. of the surface input of tke
-      efave = 1._wp        ,  &  ! coef. for the tke vert. diff. coeff.; avtke=efave*avm
-      emin  = 0.7071e-6_wp ,  &  ! minimum value of tke (m2/s2)
-      emin0 = 1.e-4_wp     ,  &  ! surface minimum value of tke (m2/s2)
-      ri_c  = 2._wp / 9._wp      ! critic Richardson number
-   REAL(wp), PUBLIC ::        &
-      eboost                     ! multiplicative coeff of the shear product.
-
-   !! caution vectopt_memory change the solution (last digit of the solver stat)
+   PUBLIC   zdf_tke        ! routine called in step module
+   PUBLIC   zdf_tke_init   ! routine also called in zdftke_jki module
+   PUBLIC   tke_rst        ! routine also called in zdftke_jki module
+
+   LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag
+   REAL(wp), PUBLIC                         ::   eboost              !: multiplicative coeff of the shear product.
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   en                  !: now turbulent kinetic energy
 # if defined key_vectopt_memory
-   REAL(wp), DIMENSION(jpi,jpj,jpk), PUBLIC ::   &
-      etmean,    &  ! coefficient used for horizontal smoothing
-      eumean,    &  ! at t-, u- and v-points
-      evmean        !
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   etmean              !: coefficient used for horizontal smoothing
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   eumean, evmean      !: at t-, u- and v-points
 # endif
 
+   !! * Namelist (namtke)
+   LOGICAL , PUBLIC ::   ln_rstke = .FALSE.         !: =T restart with tke from a run without tke with 
+     !                                              !  a none zero initial value for en
+   INTEGER , PUBLIC ::   nitke = 50 ,            &  !: number of restart iterative loops
+      &                  nmxl  =  2 ,            &  !: = 0/1/2/3 flag for the type of mixing length used
+      &                  npdl  =  1 ,            &  !: = 0/1/2 flag on prandtl number on vert. eddy coeff.
+      &                  nave  =  1 ,            &  !: = 0/1 flag for horizontal average on avt, avmu, avmv
+      &                  navb  =  0                 !: = 0/1 flag for constant or profile background avt
+   REAL(wp), PUBLIC ::   ediff = 0.1_wp       ,  &  !: coeff. for vertical eddy coef.; avt=ediff*mxl*sqrt(e)
+      &                  ediss = 0.7_wp       ,  &  !: coef. of the Kolmogoroff dissipation 
+      &                  ebb   = 3.75_wp      ,  &  !: coef. of the surface input of tke
+      &                  efave = 1._wp        ,  &  !: coef. for the tke vert. diff. coeff.; avtke=efave*avm
+      &                  emin  = 0.7071e-6_wp ,  &  !: minimum value of tke (m2/s2)
+      &                  emin0 = 1.e-4_wp     ,  &  !: surface minimum value of tke (m2/s2)
+      &                  ri_c  = 2._wp / 9._wp      !: critic Richardson number
+   NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0,   &
+      &             ri_c, nitke, nmxl, npdl, nave, navb
+
 # if defined key_cfg_1d
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   
-      e_dis,    &   ! dissipation turbulent lengh scale
-      e_mix,    &   ! mixing turbulent lengh scale
-      e_pdl,    &   ! prandl number
-      e_ric         ! local Richardson number
+   !                                                                   ! 1D cfg only
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix,      &  ! dissipation and mixing turbulent lengh scales
+      &                                          e_pdl, e_ric          ! prandl and local Richardson numbers
 #endif
 
@@ -74,10 +82,12 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
+   !! $Header$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
 
-   SUBROUTINE zdf_tke ( kt )
+   SUBROUTINE zdf_tke( kt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE zdf_tke  ***
@@ -136,45 +146,27 @@
       !!                update avt, avmu, avmv (before vertical eddy coef.)
       !!
-      !! References :
-      !!      Gaspar et al., jgr, 95, 1990,
-      !!      Blanke and Delecluse, jpo, 1991
-      !! History :
-      !!   6.0  !  91-03  (b. blanke)  Original code
-      !!   7.0  !  91-11  (G. Madec)   bug fix
-      !!   7.1  !  92-10  (G. Madec)   new mixing length and eav
-      !!   7.2  !  93-03  (M. Guyon)   symetrical conditions
-      !!   7.3  !  94-08  (G. Madec, M. Imbard)   npdl flag
-      !!   7.5  !  96-01  (G. Madec)   s-coordinates
-      !!   8.0  !  97-07  (G. Madec)   lbc
-      !!   8.1  !  99-01  (E. Stretta) new option for the mixing length
-      !!   8.5  !  02-08  (G. Madec)  ri_c and Free form, F90
-      !!   9.0  !  04-10  (C. Ethe )  1D configuration
+      !! References : Gaspar et al., jgr, 95, 1990,
+      !!              Blanke and Delecluse, jpo, 1991
       !!----------------------------------------------------------------------
-      !! * Modules used
       USE oce     , zwd   => ua,  &  ! use ua as workspace
          &          zmxlm => ta,  &  ! use ta as workspace
          &          zmxld => sa      ! use sa as workspace
-
-      !! * arguments
-      INTEGER, INTENT( in  ) ::   kt ! ocean time step
-
-      !! * local declarations
-      INTEGER ::   ji, jj, jk        ! dummy loop arguments
-      REAL(wp) ::   &
-         zmlmin, zbbrau,          &  ! temporary scalars
-         zfact1, zfact2, zfact3,  &  !
-         zrn2, zesurf,            &  !
-         ztx2, zty2, zav,         &  !
-         zcoef, zcof, zsh2,       &  !
-         zdku, zdkv, zpdl, zri,   &  !
-         zsqen, zesh2,            &  !
-         zemxl, zemlm, zemlp
+      !
+      INTEGER, INTENT(in) ::   kt ! ocean time step
+      !
+      INTEGER  ::   ji, jj, jk                  ! dummy loop arguments
+      REAL(wp) ::   zmlmin, zbbrau,          &  ! temporary scalars
+         &          zfact1, zfact2, zfact3,  &  !
+         &          zrn2, zesurf,            &  !
+         &          ztx2, zty2, zav,         &  !
+         &          zcoef, zcof, zsh2,       &  !
+         &          zdku, zdkv, zpdl, zri,   &  !
+         &          zsqen, zesh2,            &  !
+         &          zemxl, zemlm, zemlp
       !!--------------------------------------------------------------------
 
-      ! Initialization (first time-step only)
-      ! --------------
-      IF( kt == nit000  )   CALL zdf_tke_init
-
-      ! Local constant initialization
+      IF( kt == nit000  )   CALL zdf_tke_init      ! Initialization (first time-step only)
+
+      !                                            ! Local constant initialization
       zmlmin = 1.e-8
       zbbrau =  .5 * ebb / rau0
@@ -183,9 +175,7 @@
       zfact3 = 0.5 * rdt * ediss
 
-
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! I.  Mixing length
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
 
       ! Buoyancy length scale: l=sqrt(2*e/n**2)
@@ -204,5 +194,4 @@
          END DO
       END DO
-
 
       ! Physical limits for the mixing length
@@ -291,9 +280,7 @@
 # endif
 
-
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! II  Tubulent kinetic energy time stepping
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
 
       ! 1. Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt
@@ -475,5 +462,4 @@
       CALL lbc_lnk( en , 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
 
-
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! III.  Before vertical eddy vicosity and diffusivity coefficients
@@ -601,4 +587,8 @@
       ! ------------------------------=====
       CALL lbc_lnk( avt, 'W', 1. )
+
+      ! write en in restart file
+      ! ------------------------
+      IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' )
 
       IF(ln_ctl) THEN
@@ -624,21 +614,14 @@
       !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
       !!
-      !! history :
-      !!  8.5  ! 02-06 (G. Madec) original code
       !!----------------------------------------------------------------------
-      !! * Module used
       USE dynzdf_exp
       USE trazdf_exp
-
-      !! * local declarations
-      !! caution vectopt_memory change the solution (last digit of the solver stat)
+      !
 # if defined key_vectopt_memory
-      INTEGER ::   ji, jj, jk, jit   ! dummy loop indices
+      ! caution vectopt_memory change the solution (last digit of the solver stat)
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
 # else
-      INTEGER ::           jk, jit   ! dummy loop indices
+      INTEGER ::           jk   ! dummy loop indices
 # endif
-
-      NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0,   &
-         ri_c, nitke, nmxl, npdl, nave, navb
       !!----------------------------------------------------------------------
 
@@ -681,5 +664,5 @@
       ! Check nmxl and npdl values
       IF( nmxl < 0 .OR. nmxl > 3 ) CALL ctl_stop( '          bad flag: nmxl is < 0 or > 3 ' )
-      IF ( npdl < 0 .OR. npdl > 1 ) CALL ctl_stop( '          bad flag: npdl is < 0 or > 1 ' )
+      IF( npdl < 0 .OR. npdl > 1 ) CALL ctl_stop( '          bad flag: npdl is < 0 or > 1 ' )
 
       ! Horizontal average : initialization of weighting arrays 
@@ -691,6 +674,6 @@
          IF(lwp) WRITE(numout,*) '          no horizontal average on avt, avmu, avmv'
          IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
-!! caution vectopt_memory change the solution (last digit of the solver stat)
 # if defined key_vectopt_memory
+         ! caution vectopt_memory change the solution (last digit of the solver stat)
          ! weighting mean arrays etmean, eumean and evmean
          !           ( 1  1 )                                          ( 1 )
@@ -720,6 +703,6 @@
       CASE ( 1 )                ! horizontal average 
          IF(lwp) WRITE(numout,*) '          horizontal average on avt, avmu, avmv'
-!! caution vectopt_memory change the solution (last digit of the solver stat)
 # if defined key_vectopt_memory
+         ! caution vectopt_memory change the solution (last digit of the solver stat)
          ! weighting mean arrays etmean, eumean and evmean
          !           ( 1  1 )              ( 1/2  1/2 )             ( 1/2  1  1/2 )
@@ -790,21 +773,48 @@
 
 
-      ! Initialization of turbulent kinetic energy ( en )
+      ! read or initialize turbulent kinetic energy ( en )
       ! -------------------------------------------------
-      IF( ln_rstart ) THEN
-         ! no en field in the restart file, en set by iterative loop
-         IF( ln_rstke ) THEN
-            en (:,:,:) = emin * tmask(:,:,:)
-            DO jit = 2, nitke+1
-               CALL zdf_tke( jit )
-            END DO
-         ENDIF
-         ! otherwise en is already read in dtrlec called by inidtr
-      ELSE
-         ! no restart: en set to emin
-         en(:,:,:) = emin * tmask(:,:,:)
-      ENDIF
-
+      CALL tke_rst( nit000, 'READ' )
+      !
    END SUBROUTINE zdf_tke_init
+
+
+   SUBROUTINE tke_rst( kt, cdrw )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE ts_rst  ***
+     !!                     
+     !! ** Purpose : Read or write filtered free surface arrays in restart file
+     !!
+     !! ** Method  : 
+     !!
+     !!----------------------------------------------------------------------
+     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+     !
+     INTEGER ::   jit   ! dummy loop indices
+     !!----------------------------------------------------------------------
+     !
+     IF( TRIM(cdrw) == 'READ' ) THEN
+        IF( ln_rstart ) THEN
+           IF( iom_varid( numror, 'en' ) > 0 .AND. .NOT.(ln_rstke) ) THEN 
+              CALL iom_get( numror, jpdom_local, 'en', en )
+           ELSE
+              IF(lwp .AND. iom_varid(numror,'en') > 0 ) WRITE(numout,*) ' ===>>>> : previous run without tke scheme'
+              IF(lwp .AND. ln_rstke ) WRITE(numout,*) ' ===>>>> : We do not use en from the restart file'
+              IF(lwp) WRITE(numout,*) ' ===>>>> : en set by iterative loop'
+              IF(lwp) WRITE(numout,*) ' =======             ========='
+              en (:,:,:) = emin * tmask(:,:,:)
+              DO jit = 2, nitke+1
+                 CALL zdf_tke( jit )
+              END DO
+           ENDIF
+        ELSE
+           en(:,:,:) = emin * tmask(:,:,:)      ! no restart: en set to emin
+        ENDIF
+     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
+        CALL iom_rstput( kt, nitrst, numrow, 'en', en )
+     ENDIF
+     !
+   END SUBROUTINE tke_rst
 
 #else
@@ -812,5 +822,5 @@
    !!   Dummy module :                                        NO TKE scheme
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
+   PUBLIC, LOGICAL, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
 CONTAINS
    SUBROUTINE zdf_tke( kt )          ! Empty routine
Index: trunk/NEMO/OPA_SRC/ZDF/zdftke_jki.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftke_jki.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftke_jki.F90	(revision 508)
@@ -5,4 +5,7 @@
    !!                 turbulent closure parameterization
    !!=====================================================================
+   !! History :
+   !!   9.0  !  02-08  (G. Madec)  autotasking optimization
+   !!----------------------------------------------------------------------
 #if defined key_zdftke   ||   defined key_esopa
    !!----------------------------------------------------------------------
@@ -23,4 +26,5 @@
    USE taumod          ! surface stress
    USE prtctl          ! Print control
+   USE restart         ! only for lrst_oce
 
    IMPLICIT NONE
@@ -99,6 +103,4 @@
       !!      Gaspar et al., jgr, 95, 1990,
       !!      Blanke and Delecluse, jpo, 1991
-      !! History :
-      !!   9.0  !  02-08  (G. Madec)  autotasking optimization
       !!----------------------------------------------------------------------
       !! * Modules used
@@ -518,4 +520,8 @@
       CALL lbc_lnk( avt, 'W', 1. )
 
+      ! write en in restart file
+      ! ------------------------
+      IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' )
+
       IF(ln_ctl) THEN
          CALL prt_ctl(tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt , clinfo2=' t: ', ovlap=1, kdim=jpk)
Index: trunk/NEMO/OPA_SRC/in_out_manager.F90
===================================================================
--- trunk/NEMO/OPA_SRC/in_out_manager.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/in_out_manager.F90	(revision 508)
@@ -1,65 +1,66 @@
 MODULE in_out_manager   
+   !!======================================================================
+   !!                       ***  MODULE  in_out_manager  ***
+   !! Ocean physics:  vertical mixing coefficient compute from the tke 
+   !!                 turbulent closure parameterization
+   !!=====================================================================
+   !! History :   8.5  !  02-06  (G. Madec)  original code
+   !!             9.0  !  06-07  (S. Masson)  iom, add ctl_stop, ctl_warn
+   !!----------------------------------------------------------------------
 
-   USE lib_print         ! formated print library
+   !!----------------------------------------------------------------------
+   !!   ctl_stop   : update momentum and tracer Kz from a tke scheme
+   !!   ctl_warn   : initialization, namelist read, and parameters control
+   !!----------------------------------------------------------------------
    USE par_kind
    USE par_oce
+   USE lib_print         ! formated print library
 
    PUBLIC
 
    !!----------------------------------------------------------------------
-   !! namelist parameters
-   !! -------------------------------------
-   ! namrun:  parameters of the run
-   CHARACTER (len=16) ::    &   !:
-      cexper = "exp0"           !: experiment name used for output filename
-   
-   LOGICAL ::   &              !!: * namelist namrun *
-      ln_rstart = .FALSE. ,  &  !: start from (F) rest or (T) a restart file
-      ln_ctl    = .FALSE.       !: run control for debugging
-   
-   INTEGER ::                & !!: * namelist namrun *
-      no     = 0        ,    &  !: job number
-      nrstdt = 0        ,    &  !: control of the time step (0, 1 or 2)
-      nit000 = 1        ,    &  !: index of the first time step
-      nitend = 10       ,    &  !: index of the last time step
-      ndate0 = 961115   ,    &  !: initial calendar date aammjj
-      nleapy = 0        ,    &  !: Leap year calendar flag (0/1 or 30)
-      ninist = 0        ,    &  !: initial state output flag (0/1)
-      nbench = 0                !: benchmark parameter (0/1)
+   !!                   namrun namelist parameters
+   !!----------------------------------------------------------------------
+   CHARACTER (len=16) ::   cexper    = "exp0"        !: experiment name used for output filename
+   LOGICAL            ::   ln_rstart = .FALSE. ,  &  !: start from (F) rest or (T) a restart file
+      &                    ln_ctl    = .FALSE.       !: run control for debugging
+   INTEGER            ::   no     = 0        ,    &  !: job number
+      &                    nrstdt = 0        ,    &  !: control of the time step (0, 1 or 2)
+      &                    nit000 = 1        ,    &  !: index of the first time step
+      &                    nitend = 10       ,    &  !: index of the last time step
+      &                    ndate0 = 961115   ,    &  !: initial calendar date aammjj
+      &                    nleapy = 0        ,    &  !: Leap year calendar flag (0/1 or 30)
+      &                    ninist = 0        ,    &  !: initial state output flag (0/1)
+      &                    nbench = 0                !: benchmark parameter (0/1)
    
    !!----------------------------------------------------------------------
-   !! output monitoring
-   !! -----------------------------------
-
-   INTEGER ::                &  !:
-      nstock =   10 ,        &  !: restart file frequency
-      nprint =    0 ,        &  !: level of print (0 no print)
-      nwrite =   10 ,        &  !: restart file frequency
-      nictls =    0 ,        &  !: Start i indice for the SUM control
-      nictle =    0 ,        &  !: End   i indice for the SUM control
-      njctls =    0 ,        &  !: Start j indice for the SUM control
-      njctle =    0 ,        &  !: End   j indice for the SUM control
-      isplt  =    1 ,        &  !: number of processors following i
-      jsplt  =    1 ,        &  !: number of processors following j
-      ijsplt =    1             !: nb of local domain = nb of processors
+   !!                    output monitoring
+   !!----------------------------------------------------------------------
+   INTEGER ::   nstock =   10 ,        &  !: restart file frequency
+      &         nprint =    0 ,        &  !: level of print (0 no print)
+      &         nwrite =   10 ,        &  !: restart file frequency
+      &         nictls =    0 ,        &  !: Start i indice for the SUM control
+      &         nictle =    0 ,        &  !: End   i indice for the SUM control
+      &         njctls =    0 ,        &  !: Start j indice for the SUM control
+      &         njctle =    0 ,        &  !: End   j indice for the SUM control
+      &         isplt  =    1 ,        &  !: number of processors following i
+      &         jsplt  =    1 ,        &  !: number of processors following j
+      &         ijsplt =    1             !: nb of local domain = nb of processors
 
    !!----------------------------------------------------------------------
-   !! logical units
-   !! ------------------------------
-   INTEGER ::                &  !:
-      numstp     =  1 ,      &  !: logical unit for time step
-      numout     =  2 ,      &  !: logical unit for output print
-      numnam     =  3 ,      &  !: logical unit for namelist
-      numnam_ice =  4 ,      &  !: logical unit for ice namelist
-      numevo_ice = 17 ,      &  !: logical unit for ice variables (temp. evolution)
-      numice_dmp = 18 ,      &  !: logical unit for ice variables (damping)
-      numsol     = 25 ,      &  !: logical unit for solver statistics
-      numwri     = 40 ,      &  !: logical unit for output write
-      numisp     = 41 ,      &  !: logical unit for island statistics
-      numgap     = 45 ,      &  !: logical unit for differences diagnostic
-      numbol     = 67 ,      &  !: logical unit for "bol" diagnostics
-      numptr     = 68 ,      &  !: logical unit for Poleward TRansports
-      numflo     = 69           !: logical unit for drifting floats
-      !                         !: * coupled units
+   !!                        logical units
+   !!----------------------------------------------------------------------
+   INTEGER ::   numstp     =  1 ,      &  !: logical unit for time step
+      &         numout     =  2 ,      &  !: logical unit for output print
+      &         numnam     =  3 ,      &  !: logical unit for namelist
+      &         numnam_ice =  4 ,      &  !: logical unit for ice namelist
+      &         numevo_ice = 17 ,      &  !: logical unit for ice variables (temp. evolution)
+      &         numsol     = 25 ,      &  !: logical unit for solver statistics
+      &         numwri     = 40 ,      &  !: logical unit for output write
+      &         numisp     = 41 ,      &  !: logical unit for island statistics
+      &         numgap     = 45 ,      &  !: logical unit for differences diagnostic
+      &         numbol     = 67 ,      &  !: logical unit for "bol" diagnostics
+      &         numptr     = 68 ,      &  !: logical unit for Poleward TRansports
+      &         numflo     = 69           !: logical unit for drifting floats
 
    !!----------------------------------------------------------------------
@@ -67,27 +68,18 @@
    !!----------------------------------------------------------------------
    
-   INTEGER ::                &  !:
-      nstop = 0 ,            &  !: e r r o r  flag (=number of reason for a
-      !                         !                   prematurely stop the run)
-      nwarn = 0                 !: w a r n i n g  flag (=number of warning
-      !                         !                       found during the run)
-
-   
-   CHARACTER(LEN=100) :: ctmp1, ctmp2, ctmp3    ! temporary character
-   CHARACTER (len=64) ::        &                                                    !:
-      cform_err="(/,' ===>>> : E R R O R',     /,'         ===========',/)"    ,   & !:
-      cform_war="(/,' ===>>> : W A R N I N G', /,'         ===============',/)"      !:
-   LOGICAL ::   &               !:
-      lwp                ,   &  !: boolean : true on the 1st processor only
-      lsp_area = .TRUE.         !: to make a control print over a specific area
+   INTEGER            ::   nstop = 0 ,           &  !: error flag (=number of reason for a premature stop run)
+      &                    nwarn = 0                !: warning flag (=number of warning found during the run)
+   CHARACTER(LEN=100) ::   ctmp1, ctmp2, ctmp3      !: temporary character
+   CHARACTER (len=64) ::   cform_err="(/,' ===>>> : E R R O R',     /,'         ===========',/)"    ,   &  !:
+      &                    cform_war="(/,' ===>>> : W A R N I N G', /,'         ===============',/)"       !:
+   LOGICAL            ::   lwp               ,   &  !: boolean : true on the 1st processor only
+      &                    lsp_area = .TRUE.        !: to make a control print over a specific area
    !!----------------------------------------------------------------------
    !!  OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
-
 CONTAINS
-
 
    SUBROUTINE ctl_stop( cd1, cd2, cd3, cd4, cd5,   &
@@ -96,10 +88,10 @@
       !!                  ***  ROUTINE  stop_opa  ***
       !!
-      !! ** Purpose : ???
-      !!
+      !! ** Purpose : ??? blah blah....
       !!-----------------------------------------------------------------------
-      CHARACTER(len=*),INTENT(in),OPTIONAL ::  cd1, cd2, cd3, cd4, cd5, cd6, cd7, cd8, cd9, cd10
+      CHARACTER(len=*), INTENT(in), OPTIONAL ::  cd1, cd2, cd3, cd4, cd5,   &
+         &                                       cd6, cd7, cd8, cd9, cd10
       !!-----------------------------------------------------------------------
-      
+      !
       nstop = nstop + 1 
       IF(lwp) THEN
@@ -117,5 +109,5 @@
       ENDIF
       CALL FLUSH(numout)
-
+      !
    END SUBROUTINE ctl_stop
 
@@ -124,12 +116,12 @@
       &                 cd6, cd7, cd8, cd9, cd10 )
       !!-----------------------------------------------------------------------
-      !!                  ***  ROUTINE  stop_opa  ***
+      !!                  ***  ROUTINE  stop_warn  ***
       !!
-      !! ** Purpose : ???
-      !!
+      !! ** Purpose : ???  blah blah....
       !!-----------------------------------------------------------------------
-      CHARACTER(len=*),INTENT(in),OPTIONAL ::  cd1, cd2, cd3, cd4, cd5, cd6, cd7, cd8, cd9, cd10
+      CHARACTER(len=*), INTENT(in), OPTIONAL ::  cd1, cd2, cd3, cd4, cd5,   &
+         &                                       cd6, cd7, cd8, cd9, cd10
       !!-----------------------------------------------------------------------
-      
+      ! 
       nwarn = nwarn + 1 
       IF(lwp) THEN
@@ -147,6 +139,7 @@
       ENDIF
       CALL FLUSH(numout)
-
+      !
    END SUBROUTINE ctl_warn
 
+   !!=====================================================================
 END MODULE in_out_manager
Index: trunk/NEMO/OPA_SRC/iom.F90
===================================================================
--- trunk/NEMO/OPA_SRC/iom.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/iom.F90	(revision 508)
@@ -2,68 +2,62 @@
    !!=====================================================================
    !!                    ***  MODULE  iom ***
-   !!
    !! Input/Output manager :  Library to read input files
-   !!
-   !! Ongoing work : This code is here to help discussions about I/O
-   !!                library in the NEMO system
    !!====================================================================
+   !! History :  9.0  ! 05 12  (J. Belier) Original code
+   !!            9.0  ! 06 02  (S. Masson) Adaptation to NEMO
+   !!--------------------------------------------------------------------
+   !!gm  caution add !DIR nec: improved performance to be checked as well as no result changes
+
    !!--------------------------------------------------------------------
    !!   iom_open       : open a file read only
    !!   iom_close      : close a file or all files opened by iom
-   !!   iom_get        : read a field : interface to several routines
+   !!   iom_get        : read a field (interfaced to several routines)
+   !!   iom_gettime    : read the time axis cdvar in the file               !!gm : never call ??????
    !!   iom_varid      : get the id of a variable in a file
-   !!   iom_get_gblatt : ???
+   !!   iom_rstput     : write a field in a restart file (interfaced to several routines)
    !!--------------------------------------------------------------------
-   !! History :  9.0  ! 05 12  (J. Belier) Original code
-   !!            9.0  ! 06 02  (S. Masson) Adaptation to NEMO
-   !!--------------------------------------------------------------------
-   !! * Modules used
    USE in_out_manager  ! I/O manager
    USE dom_oce         ! ocean space and time domain
-   USE lbclnk          ! ???
-   USE ioipsl          ! ???
+   USE lbclnk          ! lateal boundary condition / mpp exchanges
+   USE ioipsl          ! IOIPSL library
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC iom_open, iom_close, iom_get, iom_varid, iom_get_gblatt
-
-   !! * Interfaces
+   PUBLIC iom_open, iom_close, iom_get, iom_varid, iom_rstput, iom_gettime
+
    INTERFACE iom_get
-      MODULE PROCEDURE iom_get_r_1d, iom_get_r_2d, iom_get_r_3d
+      MODULE PROCEDURE iom_get_r_0d, iom_get_r_1d, iom_get_r_2d, iom_get_r_3d
    END INTERFACE
-
-   !! * Share module variables
-   INTEGER, PARAMETER, PUBLIC ::        &  !:
-      jpdom_data                  = 1,  &  !: ( 1  :jpidta, 1  :jpjdta)
-      jpdom_global                = 2,  &  !: ( 1  :jpiglo, 1  :jpjglo)
-      jpdom_local                 = 3,  &  !: One of the 3 following cases
-      jpdom_local_full            = 4,  &  !: ( 1  :jpi   , 1  :jpi   )
-      jpdom_local_noextra         = 5,  &  !: ( 1  :nlci  , 1  :nlcj  )
-      jpdom_local_noovlap         = 6,  &  !: (nldi:nlei  ,nldj:nlej  )
-      jpdom_unknown               = 7      !: No dimension checking
-
-   !! * Module variables
-   INTEGER, PARAMETER ::    &
-      jpmax_vars    = 50,   &  ! maximum number of variables in one file
-      jpmax_dims    =  5,   &  ! maximum number of dimensions for one variable
-      jpmax_digits  =  5       ! maximum number of digits in the file name to reference the cpu number
- 
+   INTERFACE iom_rstput
+      MODULE PROCEDURE iom_rstput_0d, iom_rstput_1d, iom_rstput_2d, iom_rstput_3d
+   END INTERFACE
+
+   INTEGER, PARAMETER, PUBLIC ::   jpdom_data          = 1   !: ( 1  :jpidta, 1  :jpjdta)
+   INTEGER, PARAMETER, PUBLIC ::   jpdom_global        = 2   !: ( 1  :jpiglo, 1  :jpjglo)
+   INTEGER, PARAMETER, PUBLIC ::   jpdom_local         = 3   !: One of the 3 following cases
+   INTEGER, PARAMETER, PUBLIC ::   jpdom_local_full    = 4   !: ( 1  :jpi   , 1  :jpi   )
+   INTEGER, PARAMETER, PUBLIC ::   jpdom_local_noextra = 5   !: ( 1  :nlci  , 1  :nlcj  )
+   INTEGER, PARAMETER, PUBLIC ::   jpdom_local_noovlap = 6   !: (nldi:nlei  ,nldj:nlej  )
+   INTEGER, PARAMETER, PUBLIC ::   jpdom_unknown       = 7   !: No dimension checking
+
+   INTEGER, PARAMETER ::   jpmax_vars   = 60,   &  ! maximum number of variables in one file
+      &                    jpmax_dims   =  4,   &  ! maximum number of dimensions for one variable
+      &                    jpmax_digits =  5       ! maximum number of digits for the cpu number in the file name
 !$AGRIF_DO_NOT_TREAT
-   INTEGER :: iom_init = 0
-
-   TYPE :: flio_file
+   INTEGER ::   iom_init = 0
+   TYPE    ::   flio_file
       CHARACTER(LEN=240)                        ::   name     ! name of the file
-      INTEGER                                   ::   iopen    ! 1/0 is the file is open/not open
+      INTEGER                                   ::   iopen    ! 1(0) if the file is open(not open)
       INTEGER                                   ::   nvars    ! number of identified varibles in the file
       INTEGER                                   ::   iduld    ! id of the unlimited dimension
       CHARACTER(LEN=16), DIMENSION(jpmax_vars)  ::   cn_var   ! names of the variables
       INTEGER, DIMENSION(jpmax_vars)            ::   ndims    ! number of dimensions of the variables
-      LOGICAL, DIMENSION(jpmax_vars)            ::   luld     ! variable including unlimited dimension
-      INTEGER, DIMENSION(jpmax_dims,jpmax_vars) ::   dimsz    ! size of the dimensions of the variables
+      LOGICAL, DIMENSION(jpmax_vars)            ::   luld     ! variable using the unlimited dimension
+      INTEGER, DIMENSION(jpmax_dims,jpmax_vars) ::   dimsz    ! size of variables dimensions 
       REAL(kind=wp), DIMENSION(jpmax_vars)      ::   scf      ! scale_factor of the variables
       REAL(kind=wp), DIMENSION(jpmax_vars)      ::   ofs      ! add_offset of the variables
    END TYPE flio_file
-   TYPE(flio_file), DIMENSION(flio_max_files)   :: iom_file ! array containing the info for all opened files
+   TYPE(flio_file), DIMENSION(flio_max_files)   ::   iom_file ! array containing the info for all opened files
 !$AGRIF_END_DO_NOT_TREAT
 
@@ -76,72 +70,143 @@
 CONTAINS
 
-   SUBROUTINE iom_open( cdname, knumfl, ldimg )
+   SUBROUTINE iom_open( cdname, knumfl, ldwrt, kdom, ldimg )
       !!---------------------------------------------------------------------
       !!                   ***  SUBROUTINE  iom_open  ***
       !!
       !! ** Purpose :  open an input file read only (return 0 if not found)
-      !!
-      !! ** Method :
-      !!
       !!---------------------------------------------------------------------
-      CHARACTER(len=*), INTENT(in )  ::   cdname   ! File name
-      INTEGER, INTENT(out)           ::   knumfl   ! Identifier of the opened file
-      LOGICAL, INTENT(in ), OPTIONAL ::   ldimg    ! flg to specify that we use dimg format
-
-      CHARACTER(LEN=100) :: clname   ! the name of the file based on cdname [[+clcpu]+clcpu]
-      CHARACTER(LEN=10) :: clsuffix  ! ".nc" or ".dimg"
-      CHARACTER(LEN=10) :: clcpu     ! the cpu number (max jpmax_digits digits)
-      LOGICAL :: llok                ! check the existence 
-      INTEGER :: icnt                ! counter for digits in clcpu (max = jpmax_digits)
+      CHARACTER(len=*), INTENT(in   )           ::   cdname   ! File name
+      INTEGER         , INTENT(  out)           ::   knumfl   ! Identifier of the opened file
+      LOGICAL         , INTENT(in   ), OPTIONAL ::   ldwrt    ! read or write the file?
+      INTEGER         , INTENT(in   ), OPTIONAL ::   kdom     ! Type of domain to be written
+      LOGICAL         , INTENT(in   ), OPTIONAL ::   ldimg    ! use dimg format?
+
+      CHARACTER(LEN=100)    ::   clname    ! the name of the file based on cdname [[+clcpu]+clcpu]
+      CHARACTER(LEN=100)    ::   cltmpn    ! tempory name to store clname (in writting mode)
+      CHARACTER(LEN=10)     ::   clsuffix  ! ".nc" or ".dimg"
+      CHARACTER(LEN=10)     ::   clcpu     ! the cpu number (max jpmax_digits digits)
+      LOGICAL               ::   llok      ! check the existence 
+      LOGICAL               ::   llwrt     ! 
+      INTEGER               ::   icnt      ! counter for digits in clcpu (max = jpmax_digits)
+      INTEGER               ::   iln, ils  ! lengths of character
+      INTEGER               ::   idom      ! type of domain
+      INTEGER               ::   ifliodom  ! model domain identifier (see flio_dom_set)
+      INTEGER, DIMENSION(2) ::   iszl      ! local number of points for x,y dimensions
+      INTEGER, DIMENSION(2) ::   ifst      ! position of first local point for x,y dimensions
+      INTEGER, DIMENSION(2) ::   ilst      ! position of last local point for x,y dimensions
+      INTEGER, DIMENSION(2) ::   ihst      ! start halo size for x,y dimensions
+      INTEGER, DIMENSION(2) ::   ihnd      ! end halo size for x,y dimensions
       !---------------------------------------------------------------------
-
-      ! find the file
+      ! if iom_open is called for the first time: initialize iom_file(:)%iopen to 0
+      ! (could be done when defining iom_file in f95 but not in f90)
+      IF( iom_init == 0 ) THEN
+         iom_file(:)%iopen = 0
+         iom_init = 1
+      ENDIF
+      ! do we read or write the file?
+      ! =============
+      IF( PRESENT(ldwrt) ) THEN   ;   llwrt = ldwrt
+      ELSE                        ;   llwrt = .FALSE.
+      ENDIF
+      ! create the file name by added, if needed, TRIM(Agrif_CFixed()) and TRIM(clsuffix)
       ! =============
       clname   = trim(cdname)
 #if defined key_agrif
       if ( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
-#endif                
+#endif    
+      ! which suffix should we use?
       clsuffix = ".nc"
-      IF( PRESENT(ldimg) ) THEN
-         IF ( ldimg ) clsuffix = ".dimg"
-      ENDIF
-      !
+      IF( PRESENT(ldimg) ) THEN   ;   IF( ldimg )   clsuffix = ".dimg"   ;   ENDIF
+      ! Add the suffix if needed
+      iln = LEN_TRIM(clname)
+      ils = LEN_TRIM(clsuffix)
+      IF( iln <= ils) clname = clname(1:iln)//TRIM(clsuffix)
+      IF( clname(iln-ils+1:iln) /= TRIM(clsuffix) )   clname = clname(1:iln)//TRIM(clsuffix)
+      cltmpn = clname   ! store this name
+      ! try to find if the file to be opened already exist
       INQUIRE( FILE = clname, EXIST = llok )
-      IF( .NOT.llok ) THEN                         ! try to complete the name with the suffix only
-         clname = TRIM(cdname)//TRIM(clsuffix)
-#if defined key_agrif
-         if ( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
-#endif                
-         INQUIRE( FILE = clname, EXIST = llok )
-         IF( .NOT.llok ) THEN                      ! try to complete the name with both cpu number and suffix
-            WRITE(clcpu,*) narea-1
-            clcpu  = trim(adjustl(clcpu))
-            clname = trim(cdname)//"_"
-#if defined key_agrif
-            if ( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
-#endif                
-            icnt = 0
-            INQUIRE( FILE = trim(clname)//trim(clcpu)//trim(clsuffix), EXIST = llok ) 
-            DO WHILE( .NOT.llok .AND. icnt < jpmax_digits )   ! we try fifferent formats for the cpu number by adding 0
-               clname = trim(clname)//"0"
-               INQUIRE( FILE = trim(clname)//trim(clcpu)//trim(clsuffix), EXIST = llok )
-               icnt = icnt + 1
-            END DO
-            IF( .NOT.llok ) THEN                   ! no way to find the files...
-               CALL ctl_stop( 'iom_open: file '//trim(clname)//'... not found' )
+      IF( .NOT.llok ) THEN
+         ! we try to add the cpu number to the name
+         WRITE(clcpu,*) narea-1
+         clcpu  = TRIM(ADJUSTL(clcpu))
+         iln = INDEX(clname,TRIM(clsuffix))
+         clname = clname(1:iln-1)//'_'//TRIM(clcpu)//TRIM(clsuffix)
+         icnt = 0
+         INQUIRE( FILE = clname, EXIST = llok ) 
+         ! we try different formats for the cpu number by adding 0
+         DO WHILE( .NOT.llok .AND. icnt < jpmax_digits )
+            clcpu  = "0"//trim(clcpu)
+            clname = clname(1:iln-1)//'_'//TRIM(clcpu)//TRIM(clsuffix)
+            INQUIRE( FILE = clname, EXIST = llok )
+            icnt = icnt + 1
+         END DO
+      ENDIF
+      !
+      IF( llok ) THEN      ! Open the file
+         !                 ! =============
+         IF( llwrt ) THEN 
+            IF(lwp) WRITE(numout,*) '          iom_open ~~~  open existing file: '//TRIM(clname)//' in WRITE mode'
+            CALL flioopfd( TRIM(clname), knumfl, "WRITE" )
+         ELSE
+            IF(lwp) WRITE(numout,*) '          iom_open ~~~  open existing file: '//TRIM(clname)//' in READ mode'
+            CALL flioopfd( TRIM(clname), knumfl )
+         ENDIF
+      ELSE                 ! no way to find the file...
+         !                 ! =======================
+         IF( llwrt ) THEN 
+            ! file opened in write mode
+            ! the file does not exist, we must create it...
+            ! =============
+            llok = .TRUE.
+            ! on which domain must the file be written??
+            ! check the domain definition
+            idom = jpdom_local_noovlap   ! default definition
+            IF( PRESENT(kdom) )   idom = kdom
+            ! create the domain informations
+            ! =============
+            SELECT CASE (idom)
+            CASE (jpdom_local_full)
+               iszl = (/ jpi             , jpj              /)
+               ifst = (/ nimpp           , njmpp            /)
+               ilst = (/ nimpp + jpi - 1 , njmpp + jpj - 1  /)
+               ihst = (/ nldi - 1        , nldj - 1         /)
+               ihnd = (/ jpi - nlei      , jpj - nlej       /)
+            CASE (jpdom_local_noextra)
+               iszl = (/ nlci            , nlcj             /)
+               ifst = (/ nimpp           , njmpp            /)
+               ilst = (/ nimpp + nlci - 1, njmpp + nlcj - 1 /)
+               ihst = (/ nldi - 1        , nldj - 1         /)
+               ihnd = (/ nlci - nlei     , nlcj - nlej      /)
+            CASE (jpdom_local_noovlap)
+               iszl = (/ nlei - nldi + 1 , nlej - nldj + 1  /)
+               ifst = (/ nimpp + nldi - 1, njmpp + nldj - 1 /)
+               ilst = (/ nimpp + nlei - 1, njmpp + nlej - 1 /)
+               ihst = (/ 0               , 0                /)
+               ihnd = (/ 0               , 0                /)
+            CASE DEFAULT
+               llok = .FALSE.
+               CALL ctl_stop( 'iom_open: wrong value of kdom, only jpdom_local* cases are accepted' )
+            END SELECT
+            IF( llok ) THEN
+               CALL flio_dom_set( jpnij, narea-1, (/1, 2/), (/jpiglo, jpjglo/)   &
+                    &                  , iszl, ifst, ilst, ihst, ihnd, 'BOX', ifliodom )        
+               ! create the file
+               ! =============
+               ! Note that fliocrfd may change the value of clname (add the cpu number...)
+               clname = cltmpn   ! get back the file name without the cpu number in it
+               IF(lwp) WRITE(numout,*) '          iom_open ~~~  create new file: '//trim(clname)//' in WRITE mode'
+               CALL fliocrfd( clname, (/'x'    , 'y'    , 'z', 't'/)   &
+                    &               , (/iszl(1), iszl(2), jpk, -1 /)   &
+                    &               , knumfl, ifliodom )
             ENDIF
-            clname = trim(clname)//trim(clcpu)//trim(clsuffix)
+         ELSE
+            ! the file is open for read-only, it must exist...
+            iln = INDEX( cltmpn,TRIM(clsuffix) )
+            CALL ctl_stop( 'iom_open: file '//cltmpn(1:iln-1)//'* not found' )
          ENDIF
       ENDIF
-
-      ! Open the file
+      ! start to fill the information of opened files
       ! =============
       IF( llok ) THEN                         
-         IF (lwp) WRITE(numout,*) 'iom_open ~~~  open file: '//trim(clname)
-         CALL flioopfd( trim(clname), knumfl )
-         IF( iom_init == 0 ) THEN
-            iom_file(:)%iopen = 0
-            iom_init = 1
-         ENDIF
          iom_file(knumfl)%iopen      = 1
          iom_file(knumfl)%name       = TRIM(clname)
@@ -152,8 +217,9 @@
          ! does the file contain time axis (that must be unlimitted) ?
          CALL flioinqf( knumfl, id_uld = iom_file(knumfl)%iduld )
+         IF(lwp) WRITE(numout,*) '                   ---> OK'
       ELSE
-         knumfl = 0
-      ENDIF
-  
+         knumfl = 0      ! return error flag
+      ENDIF
+      !
    END SUBROUTINE iom_open
 
@@ -164,16 +230,12 @@
       !!
       !! ** Purpose : close an input file, or all files opened by iom
-      !!
-      !! ** Method :
-      !!
       !!--------------------------------------------------------------------
-      INTEGER,INTENT(in), OPTIONAL ::   knumfl   ! Identifier of the file to be closed
-      !                                          ! If this argument is not present,
-      !                                          ! all the files opened by iom are closed.
-
-      INTEGER           ::   jf            ! dummy loop indices
-      INTEGER           ::   i_s, i_e      ! temporary integer
+      INTEGER, INTENT(in), OPTIONAL ::   knumfl   ! Identifier of the file to be closed
+      !                                           ! No argument : all the files opened by iom are closed
+
+      INTEGER ::   jf         ! dummy loop indices
+      INTEGER ::   i_s, i_e   ! temporary integer
       !---------------------------------------------------------------------
-
+      !
       IF( PRESENT(knumfl) ) THEN
          i_s = knumfl
@@ -183,8 +245,10 @@
          i_e = flio_max_files
       ENDIF
-      IF ( i_s > 0 ) THEN
+      
+      IF( i_s > 0 ) THEN
          DO jf = i_s, i_e
             IF( iom_file(jf)%iopen > 0 ) THEN
                CALL flioclo( jf )
+               IF(lwp) WRITE(numout,*) '          iom_close, close file: '//TRIM(iom_file(knumfl)%name)//' ok'
                iom_file(jf)%iopen      = 0
                iom_file(jf)%name       = 'NONE'
@@ -200,142 +264,118 @@
          END DO
       ENDIF
-          
+      !    
    END SUBROUTINE iom_close
  
 
    !!----------------------------------------------------------------------
-   !!                   INTERFACE iom_u_getv
+   !!                   INTERFACE iom_get_123d
    !!----------------------------------------------------------------------
-   SUBROUTINE iom_get_r_1d( knumfl, kdom , cdvar , pvar  ,   &
-      &                             ktime, kstart, kcount )
-      INTEGER               , INTENT(in )           ::   knumfl    ! Identifier of the file
-      INTEGER               , INTENT(in )           ::   kdom      ! Type of domain to be read
-      CHARACTER(len=*)      , INTENT(in )           ::   cdvar     ! Name of the variable
-      REAL(wp), DIMENSION(:), INTENT(out)           ::   pvar      ! read field
-      INTEGER               , INTENT(in ) ,OPTIONAL ::   ktime     ! record number
-      INTEGER , DIMENSION(:), INTENT(in ), OPTIONAL ::   kstart    ! start position of the reading in each axis 
-      INTEGER , DIMENSION(:), INTENT(in ), OPTIONAL ::   kcount    ! number of points to be read in each axis
-!
-      CHARACTER(LEN=100) :: clinfo                    ! info character
-!
-      clinfo = 'iom_get_r_1d, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
-      IF( PRESENT(kstart) ) THEN
-         IF ( SIZE(kstart) /= 1 ) CALL ctl_stop( trim(clinfo), 'kstart must be a 1 element vector' )
-      ENDIF
-      IF( PRESENT(kcount) ) THEN
-         IF ( SIZE(kcount) /= 1 ) CALL ctl_stop( trim(clinfo), 'kcount must be a 1 element vector' )
-      ENDIF
-      IF ( knumfl > 0 ) CALL iom_u_getv( knumfl, kdom       , cdvar        , pv_r1d=pvar,   &
-           &                                     ktime=ktime, kstart=kstart, kcount=kcount )
+   SUBROUTINE iom_get_r_0d( knumfl, cdvar, pvar )
+      INTEGER         , INTENT(in   )                 ::   knumfl    ! Identifier of the file
+      CHARACTER(len=*), INTENT(in   )                 ::   cdvar     ! Name of the variable
+      REAL(wp)        , INTENT(  out)                 ::   pvar      ! read field
+      !
+      IF( knumfl > 0 .AND. iom_varid( knumfl, cdvar ) > 0 )   CALL fliogetv( knumfl, cdvar, pvar )
+   END SUBROUTINE iom_get_r_0d
+
+   SUBROUTINE iom_get_r_1d( knumfl, kdom, cdvar, pvar, ktime, kstart, kcount )
+      INTEGER         , INTENT(in   )                         ::   knumfl    ! Identifier of the file
+      INTEGER         , INTENT(in   )                         ::   kdom      ! Type of domain to be read
+      CHARACTER(len=*), INTENT(in   )                         ::   cdvar     ! Name of the variable
+      REAL(wp)        , INTENT(  out), DIMENSION(:)           ::   pvar      ! read field
+      INTEGER         , INTENT(in   )              , OPTIONAL ::   ktime     ! record number
+      INTEGER         , INTENT(in   ), DIMENSION(1), OPTIONAL ::   kstart    ! start axis position of the reading 
+      INTEGER         , INTENT(in   ), DIMENSION(1), OPTIONAL ::   kcount    ! number of points in each axis
+      !
+      IF( knumfl > 0 ) CALL iom_get_123d( knumfl, kdom       , cdvar        , pv_r1d=pvar,   &
+         &                                        ktime=ktime, kstart=kstart, kcount=kcount )
    END SUBROUTINE iom_get_r_1d
-   SUBROUTINE iom_get_r_2d( knumfl, kdom , cdvar , pvar  ,   &
-      &                             ktime, kstart, kcount )
-      INTEGER,INTENT(in) :: knumfl
-      INTEGER,INTENT(in) :: kdom
-      CHARACTER(len=*),INTENT(in) :: cdvar
-      REAL(wp),INTENT(out),DIMENSION(:,:) :: pvar
-      INTEGER,INTENT(in),OPTIONAL :: ktime
-      INTEGER,DIMENSION(:),INTENT(in),OPTIONAL :: kstart
-      INTEGER,DIMENSION(:),INTENT(in),OPTIONAL :: kcount
-!
-      CHARACTER(LEN=100) :: clinfo                    ! info character
-!
-      clinfo = 'iom_get_r_2d, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
-      IF( PRESENT(kstart) ) THEN
-         IF ( size(kstart) /= 2 ) CALL ctl_stop(trim(clinfo), 'kstart must be a 2 element vector')
-      ENDIF
-      IF( PRESENT(kcount) ) THEN
-         IF ( size(kcount) /= 2 ) CALL ctl_stop(trim(clinfo), 'kcount must be a 2 element vector')
-      ENDIF
-      IF ( knumfl > 0 ) CALL iom_u_getv( knumfl, kdom       , cdvar        , pv_r2d=pvar,   &
-         &                                       ktime=ktime, kstart=kstart, kcount=kcount )
+
+   SUBROUTINE iom_get_r_2d( knumfl, kdom, cdvar, pvar, ktime, kstart, kcount )
+      INTEGER         , INTENT(in   )                           ::   knumfl    ! Identifier of the file
+      INTEGER         , INTENT(in   )                           ::   kdom      ! Type of domain to be read
+      CHARACTER(len=*), INTENT(in   )                           ::   cdvar     ! Name of the variable
+      REAL(wp)        , INTENT(  out), DIMENSION(:,:)           ::   pvar      ! read field
+      INTEGER         , INTENT(in   )                , OPTIONAL ::   ktime     ! record number
+      INTEGER         , INTENT(in   ), DIMENSION(2)  , OPTIONAL ::   kstart    ! start axis position of the reading 
+      INTEGER         , INTENT(in   ), DIMENSION(2)  , OPTIONAL ::   kcount    ! number of points in each axis
+      !
+      IF( knumfl > 0 ) CALL iom_get_123d( knumfl, kdom       , cdvar        , pv_r2d=pvar,   &
+         &                                        ktime=ktime, kstart=kstart, kcount=kcount )
    END SUBROUTINE iom_get_r_2d
-   SUBROUTINE iom_get_r_3d( knumfl, kdom , cdvar , pvar  ,   &
-      &                             ktime, kstart, kcount )
-      INTEGER,INTENT(in) :: knumfl
-      INTEGER,INTENT(in) :: kdom
-      CHARACTER(len=*),INTENT(in) :: cdvar
-      REAL(wp),INTENT(out),DIMENSION(:,:,:) :: pvar
-      INTEGER,INTENT(in),OPTIONAL :: ktime
-      INTEGER,DIMENSION(:),INTENT(in),OPTIONAL :: kstart
-      INTEGER,DIMENSION(:),INTENT(in),OPTIONAL :: kcount
-!
-      CHARACTER(LEN=100) :: clinfo                    ! info character
-!
-      clinfo = 'iom_get_r_3d, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
-      IF ( PRESENT(kstart) ) THEN
-         IF ( size(kstart) /= 3 ) CALL ctl_stop(trim(clinfo), 'kstart must be a 3 element vector')
-      ENDIF
-      IF ( PRESENT(kcount) ) THEN
-         IF ( size(kcount) /= 3 ) CALL ctl_stop(trim(clinfo), 'kcount must be a 3 element vector')
-      ENDIF
-      IF ( knumfl > 0 ) CALL iom_u_getv( knumfl, kdom       , cdvar        , pv_r3d=pvar,   &
-	     &                                   ktime=ktime, kstart=kstart, kcount=kcount )
+
+   SUBROUTINE iom_get_r_3d( knumfl, kdom, cdvar, pvar, ktime, kstart, kcount )
+      INTEGER         , INTENT(in   )                             ::   knumfl    ! Identifier of the file
+      INTEGER         , INTENT(in   )                             ::   kdom      ! Type of domain to be read
+      CHARACTER(len=*), INTENT(in   )                             ::   cdvar     ! Name of the variable
+      REAL(wp)        , INTENT(  out), DIMENSION(:,:,:)           ::   pvar      ! read field
+      INTEGER         , INTENT(in   )                  , OPTIONAL ::   ktime     ! record number
+      INTEGER         , INTENT(in   ), DIMENSION(3)    , OPTIONAL ::   kstart    ! start axis position of the reading 
+      INTEGER         , INTENT(in   ), DIMENSION(3)    , OPTIONAL ::   kcount    ! number of points in each axis
+      !
+      IF( knumfl > 0 ) CALL iom_get_123d( knumfl, kdom       , cdvar        , pv_r3d=pvar,   &
+         &                                        ktime=ktime, kstart=kstart, kcount=kcount )
    END SUBROUTINE iom_get_r_3d
    !!----------------------------------------------------------------------
 
-
-   SUBROUTINE iom_u_getv( knumfl, kdom  , cdvar , &
-        &                   pv_r1d, pv_r2d, pv_r3d, &
-        &                   ktime , kstart, kcount )
+   SUBROUTINE iom_get_123d( knumfl, kdom  , cdvar ,   &
+        &                   pv_r1d, pv_r2d, pv_r3d,   &
+        &                   ktime , kstart, kcount  )
      !!-----------------------------------------------------------------------
-     !!                  ***  ROUTINE  iom_u_getv  ***
+     !!                  ***  ROUTINE  iom_get_123d  ***
      !!
      !! ** Purpose : read a 1D/2D/3D variable
      !!
-     !! ** Method : read ONE time step at each CALL
-     !!
+     !! ** Method : read ONE record at each CALL
      !!-----------------------------------------------------------------------
-     INTEGER,                      INTENT(in )           ::   knumfl     ! Identifier of the file
-     INTEGER,                      INTENT(in )           ::   kdom       ! Type of domain to be read
-     CHARACTER(len=*),             INTENT(in )           ::   cdvar      ! Name of the variable
-     REAL(wp), DIMENSION(:)      , INTENT(out), OPTIONAL ::   pv_r1d     ! read field (1D case)
-     REAL(wp), DIMENSION(:,:)    , INTENT(out), OPTIONAL ::   pv_r2d     ! read field (2D case)
-     REAL(wp), DIMENSION(:,:,:)  , INTENT(out), OPTIONAL ::   pv_r3d     ! read field (3D case)
-     INTEGER                     , INTENT(in ), OPTIONAL ::   ktime      ! record number
-     INTEGER , DIMENSION(:)      , INTENT(in ), OPTIONAL ::   kstart     ! start position of the reading in each axis 
-     INTEGER , DIMENSION(:)      , INTENT(in ), OPTIONAL ::   kcount     ! number of points to be read in each axis
-
-     INTEGER                        :: jl            ! loop on number of dimension 
-     INTEGER                        :: idom,    &    ! type of domain
-          &                            idvar,   &    ! id of the variable
-          &                            inbdim,  &    ! number of dimensions of the variable
-          &                            idmspc,  &    ! number of spatial dimensions 
-          &                            itime,   &    ! record number
-          &                            istop         ! temporary value of nstop
-     INTEGER, DIMENSION(jpmax_dims) :: istart,  &    ! starting point to read for each axis
-          &                            icnt,    &    ! number of value to read along each axis 
-          &                            idimsz        ! size of the dimensions of the variable
-     REAL(wp)                       :: zscf, zofs    ! sacle_factor and add_offset
-     INTEGER                        ::  itmp         ! temporary integer
-     CHARACTER(LEN=100)             :: clinfo        ! info character
+     INTEGER                    , INTENT(in   )           ::   knumfl     ! Identifier of the file
+     INTEGER                    , INTENT(in   )           ::   kdom       ! Type of domain to be read
+     CHARACTER(len=*)           , INTENT(in   )           ::   cdvar      ! Name of the variable
+     REAL(wp), DIMENSION(:)     , INTENT(  out), OPTIONAL ::   pv_r1d     ! read field (1D case)
+     REAL(wp), DIMENSION(:,:)   , INTENT(  out), OPTIONAL ::   pv_r2d     ! read field (2D case)
+     REAL(wp), DIMENSION(:,:,:) , INTENT(  out), OPTIONAL ::   pv_r3d     ! read field (3D case)
+     INTEGER                    , INTENT(in   ), OPTIONAL ::   ktime      ! record number
+     INTEGER , DIMENSION(:)     , INTENT(in   ), OPTIONAL ::   kstart     ! start position of the reading in each axis 
+     INTEGER , DIMENSION(:)     , INTENT(in   ), OPTIONAL ::   kcount     ! number of points to be read in each axis
+     !
+     INTEGER                        ::   jl          ! loop on number of dimension 
+     INTEGER                        ::   idom,    &  ! type of domain
+          &                              idvar,   &  ! id of the variable
+          &                              inbdim,  &  ! number of dimensions of the variable
+          &                              idmspc,  &  ! number of spatial dimensions 
+          &                              itime,   &  ! record number
+          &                              istop       ! temporary value of nstop
+     INTEGER, DIMENSION(jpmax_dims) ::   istart,  &  ! starting point to read for each axis
+          &                              icnt,    &  ! number of value to read along each axis 
+          &                              idimsz      ! size of the dimensions of the variable
+     REAL(wp)                       ::   zscf, zofs  ! sacle_factor and add_offset
+     INTEGER                        ::   itmp        ! temporary integer
+     CHARACTER(LEN=100)             ::   clinfo      ! info character
      !---------------------------------------------------------------------
-     clinfo = 'iom_u_getv, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
+     !
+     clinfo = '          iom_get_123d, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
      ! local definition of the domain ?
      idom = kdom
      ! check kcount and kstart optionals parameters...
-     IF( PRESENT(kcount) .AND. (.NOT. PRESENT(kstart)) ) &
+     IF( PRESENT(kcount) .AND. (.NOT. PRESENT(kstart)) )   &
           CALL ctl_stop( trim(clinfo), 'kcount present needs kstart present' )
-     IF( PRESENT(kstart) .AND. (.NOT. PRESENT(kcount)) ) &
+     IF( PRESENT(kstart) .AND. (.NOT. PRESENT(kcount)) )   &
           CALL ctl_stop( trim(clinfo), 'kstart present needs kcount present' )
-     IF( PRESENT(kstart) .AND. idom /= jpdom_unknown ) &
+     IF( PRESENT(kstart) .AND. idom /= jpdom_unknown   )   &
           CALL ctl_stop( trim(clinfo), 'kstart present needs kdom = jpdom_unknown' )
 
      ! Search for the variable in the data base (eventually actualize data)
-     !-
      istop = nstop
      idvar = iom_varid( knumfl, cdvar )
      !
-     IF ( idvar > 0 ) THEN
+     IF( idvar > 0 ) THEN
         ! to write iom_file(knumfl)%dimsz in a shorter way !
         idimsz(:) = iom_file(knumfl)%dimsz(:, idvar) 
-        inbdim = iom_file(knumfl)%ndims(idvar)! number of dimensions in the file
-        idmspc = inbdim ! number of spatial dimensions in the file
-        IF( iom_file(knumfl)%luld(idvar) ) idmspc = inbdim - 1
-        IF( idmspc > 3 ) CALL ctl_stop(trim(clinfo), 'the file has more than 3 spatial dimensions',   &
-             &                    'this case is not coded...') 
-        ! Identify the domain in case of jpdom_local definition
-        !-
-        IF( idom == jpdom_local ) THEN
+        inbdim = iom_file(knumfl)%ndims(idvar)            ! number of dimensions in the file
+        idmspc = inbdim                                   ! number of spatial dimensions in the file
+        IF( iom_file(knumfl)%luld(idvar) )   idmspc = inbdim - 1
+        IF( idmspc > 3 )   CALL ctl_stop(trim(clinfo),   &
+           &                    'the file has more than 3 spatial dimensions this case is not coded...' ) 
+        IF( idom == jpdom_local ) THEN        ! Identify the domain in case of jpdom_local definition
            IF( idimsz(1) == jpi .AND. idimsz(2) == jpj ) THEN 
               idom = jpdom_local_full
@@ -348,10 +388,10 @@
            ENDIF
         ENDIF
-
+        !
         ! definition of istart and icnt
-        !-
+        !
         ! initializations
         istart(:) = 1
-        icnt(:) = 1
+        icnt  (:) = 1
         itime = 1
         IF( PRESENT(ktime) ) itime = ktime
@@ -383,5 +423,5 @@
            CASE (2)
               ! data is 2d array (+ maybe a temporal dimension)
-              IF ( PRESENT(kstart) ) THEN
+              IF( PRESENT(kstart) ) THEN
                  istart(1:3) = (/ kstart(1:2), itime /)
                  icnt(1:2) = kcount(1:2)
@@ -404,9 +444,9 @@
               ENDIF
            CASE DEFAULT
-              IF ( itime == 1 .AND. idimsz(3) == 1 .AND. idmspc == 3 ) THEN
-                 CALL ctl_warn( trim(clinfo), '2D array but 3 spatial dimensions for the data...', &
-                      &         'As the size of the z dimension is 1 and as we try to read the first reacord, ',     &
-                      &         'we accept this case even if there is a possible mix-up between z and time dimension...')           
-                 IF ( PRESENT(kstart) ) THEN
+              IF( itime == 1 .AND. idimsz(3) == 1 .AND. idmspc == 3 ) THEN
+                 CALL ctl_warn( trim(clinfo), '2D array but 3 spatial dimensions for the data...',                &
+                      &         'As the size of the z dimension is 1 and as we try to read the first record, ',   &
+                      &         'we accept this case even if there is a possible mix-up between z and time dimension' )           
+                 IF( PRESENT(kstart) ) THEN
                     istart(1:2) = kstart(1:2)
                     icnt(1:2) = kcount(1:2)
@@ -428,18 +468,18 @@
                  ENDIF
               ELSE
-                 CALL ctl_stop( trim(clinfo), 'to keep iom lisibility, when reading a 2D array,', &
-                      &         'we do not accept data with more than 2 spatial dimension',     &
-                      &         'Use ncwa -a to suppress the unnecessary dimensions')           
+                 CALL ctl_stop( trim(clinfo), 'to keep iom lisibility, when reading a 2D array,',           &
+                      &                       'we do not accept data with more than 2 spatial dimension',   &
+                      &                       'Use ncwa -a to suppress the unnecessary dimensions' )
               ENDIF
            END SELECT
         ELSEIF( PRESENT(pv_r3d) ) THEN
            SELECT CASE (idmspc)
-           CASE (1)
-              CALL ctl_stop(trim(clinfo), 'the file has only 1 spatial dimension',   &
-                   &        'it is impossible to read a 3d array from this file...')
-           CASE (2)
-              CALL ctl_stop(trim(clinfo), 'the file has only 2 spatial dimension',   &
-                   &        'it is impossible to read a 3d array from this file...')
-           CASE (3)
+           CASE( 1 )
+              CALL ctl_stop( trim(clinfo), 'the file has only 1 spatial dimension',            &
+                   &                       'it is impossible to read a 3d array from this file...' )
+           CASE( 2 )
+              CALL ctl_stop( trim(clinfo), 'the file has only 2 spatial dimension',            &
+                   &                       'it is impossible to read a 3d array from this file...' )
+           CASE( 3 )
               ! data is 3d array (+ maybe a temporal dimension)
               IF( PRESENT(kstart) ) THEN
@@ -469,7 +509,7 @@
               ENDIF
            CASE DEFAULT
-              CALL ctl_stop( trim(clinfo), 'to keep iom lisibility, when reading a 3D array,', &
-                   &         'we do not accept data with more than 3 spatial dimension',     &
-                   &         'Use ncwa -a to suppress the unnecessary dimensions')           
+              CALL ctl_stop( trim(clinfo), 'to keep iom lisibility, when reading a 3D array,',   &
+                   &         'we do not accept data with more than 3 spatial dimension',         &
+                   &         'Use ncwa -a to suppress the unnecessary dimensions' )           
            END SELECT
         ENDIF
@@ -491,5 +531,5 @@
            itmp = size(pv_r1d)
            WRITE(ctmp1,*) 'size(pv_r1d): ', itmp, ' /= icnt(1): ', icnt(1)
-           IF( itmp /= icnt(1) ) CALL ctl_stop( trim(clinfo), ctmp1 )
+           IF( itmp /= icnt(1) )   CALL ctl_stop( trim(clinfo), ctmp1 )
         ELSEIF( PRESENT(pv_r2d) ) THEN
            DO jl = 1, 2
@@ -501,5 +541,5 @@
                  WRITE(ctmp1,*) 'size(pv_r2d(nldi:nlei,nldj:nlej), ',jl,'): ',itmp,' /= icnt(',jl,'): ',icnt(jl)
               ENDIF
-              IF( itmp /= icnt(jl) ) CALL ctl_stop( trim(clinfo), ctmp1 )
+              IF( itmp /= icnt(jl) )   CALL ctl_stop( trim(clinfo), ctmp1 )
            END DO
         ELSEIF( PRESENT(pv_r3d) ) THEN
@@ -512,5 +552,5 @@
                  WRITE(ctmp1,*) 'size(pv_r3d(nldi:nlei,nldj:nlej), ',jl,'): ',itmp,' /= icnt(',jl,'): ',icnt(jl)
               ENDIF
-              IF( itmp /= icnt(jl) ) CALL ctl_stop( trim(clinfo), ctmp1 )
+              IF( itmp /= icnt(jl) )   CALL ctl_stop( trim(clinfo), ctmp1 )
            END DO
         ENDIF
@@ -520,6 +560,4 @@
      !-     
      IF( istop == nstop) THEN ! no additional errors until this point...
-        !
-        istop = nstop
         !
         zscf = iom_file(knumfl)%scf(idvar)      ! scale factor
@@ -529,12 +567,14 @@
            CALL fliogetv( knumfl, cdvar, pv_r1d(:), start=istart(1:inbdim), count=icnt(1:inbdim) )
            !--- Apply scale_factor and offset
-           IF( zscf /= 1. ) pv_r1d(:) = pv_r1d(:) * zscf 
-           IF( zofs /= 0. ) pv_r1d(:) = pv_r1d(:) + zofs
+           IF( zscf /= 1. )   pv_r1d(:) = pv_r1d(:) * zscf 
+           IF( zofs /= 0. )   pv_r1d(:) = pv_r1d(:) + zofs
         ELSEIF( PRESENT(pv_r2d) ) THEN
            IF( idom /= jpdom_unknown ) THEN
               CALL fliogetv( knumfl, cdvar, pv_r2d(nldi:nlei,nldj:nlej), start=istart(1:inbdim), count=icnt(1:inbdim) )
               !--- Apply scale_factor and offset
-              IF (zscf /= 1.) pv_r2d(nldi:nlei, nldj:nlej) = pv_r2d(nldi:nlei,nldj:nlej) * zscf
-              IF (zofs /= 0.) pv_r2d(nldi:nlei, nldj:nlej) = pv_r2d(nldi:nlei,nldj:nlej) + zofs
+!CDIR NOUNROLL
+              IF( zscf /= 1.)   pv_r2d(nldi:nlei, nldj:nlej) = pv_r2d(nldi:nlei,nldj:nlej) * zscf
+!CDIR NOUNROLL
+              IF( zofs /= 0.)   pv_r2d(nldi:nlei, nldj:nlej) = pv_r2d(nldi:nlei,nldj:nlej) + zofs
               !--- Fill the overlap areas and extra hallows (mpp)
               CALL lbc_lnk (pv_r2d,'Z',-999.,'no0')
@@ -542,13 +582,18 @@
               CALL fliogetv( knumfl, cdvar, pv_r2d(:,:), start=istart(1:inbdim), count=icnt(1:inbdim) )
               !--- Apply scale_factor and offset
-              IF (zscf /= 1.) pv_r2d(:,:) = pv_r2d(:,:) * zscf
-              IF (zofs /= 0.) pv_r2d(:,:) = pv_r2d(:,:) + zofs
+!CDIR COLLAPSE
+              IF( zscf /= 1.)   pv_r2d(:,:) = pv_r2d(:,:) * zscf
+!CDIR COLLAPSE
+              IF( zofs /= 0.)   pv_r2d(:,:) = pv_r2d(:,:) + zofs
            ENDIF
         ELSEIF( PRESENT(pv_r3d) ) THEN
            IF( idom /= jpdom_unknown ) THEN
-              CALL fliogetv( knumfl, cdvar, pv_r3d(nldi:nlei,nldj:nlej,:), start=istart(1:inbdim), count=icnt(1:inbdim) )
+              CALL fliogetv( knumfl, cdvar, pv_r3d(nldi:nlei,nldj:nlej,:), start=istart(1:inbdim),   &
+                 &                                                         count=icnt  (1:inbdim) )
               !--- Apply scale_factor and offset
-              IF( zscf /= 1. ) pv_r3d(nldi:nlei,nldj:nlej,:) = pv_r3d(nldi:nlei,nldj:nlej,:) * zscf
-              IF( zofs /= 0. ) pv_r3d(nldi:nlei,nldj:nlej,:) = pv_r3d(nldi:nlei,nldj:nlej,:) + zofs
+!CDIR NOUNROLL
+              IF( zscf /= 1. )   pv_r3d(nldi:nlei,nldj:nlej,:) = pv_r3d(nldi:nlei,nldj:nlej,:) * zscf
+!CDIR NOUNROLL
+              IF( zofs /= 0. )   pv_r3d(nldi:nlei,nldj:nlej,:) = pv_r3d(nldi:nlei,nldj:nlej,:) + zofs
               !--- Fill the overlap areas and extra hallows (mpp)
               IF( icnt(3) == jpk )   CALL lbc_lnk( pv_r3d,'Z',-999.,'no0' ) ! this if could be removed with the new lbc_lnk ...
@@ -556,186 +601,317 @@
               CALL fliogetv( knumfl, cdvar, pv_r3d(:,:,:), start=istart(1:inbdim), count=icnt(1:inbdim) )
               !--- Apply scale_factor and offset
-              IF (zscf /= 1.)   pv_r3d(:,:,:) = pv_r3d(:,:,:) * zscf
-              IF (zofs /= 0.)   pv_r3d(:,:,:) = pv_r3d(:,:,:) + zofs
+!CDIR COLLAPSE
+              IF( zscf /= 1.)   pv_r3d(:,:,:) = pv_r3d(:,:,:) * zscf
+!CDIR COLLAPSE
+              IF( zofs /= 0.)   pv_r3d(:,:,:) = pv_r3d(:,:,:) + zofs
            ENDIF
         ENDIF
         !
-        IF( istop == nstop .AND. lwp )  &
-             &  WRITE(numout,*) ' read '//trim(cdvar)//' in '//trim(iom_file(knumfl)%name)//' ok'
+        IF( istop == nstop .AND. lwp )   &
+           WRITE(numout,*) '           read '//trim(cdvar)//' in '//trim(iom_file(knumfl)%name)//' ok'
      ENDIF
      !
-   END SUBROUTINE iom_u_getv
+   END SUBROUTINE iom_get_123d
 
    
    SUBROUTINE iom_gettime( knumfl, cdvar, ptime )
-     !!--------------------------------------------------------------------
-     !!                   ***  SUBROUTINE  iom_close  ***
-     !!
-     !! ** Purpose : read the time axis cdvar in the file 
-     !!
-     !! ** Method :
-     !!
-     !!--------------------------------------------------------------------
-     INTEGER               , INTENT(in)  ::   knumfl   ! Identifier of the file to be closed
-     CHARACTER(len=*)      , INTENT(in)  ::   cdvar    ! time axis name
-     REAL(wp), DIMENSION(:), INTENT(out) ::   ptime    ! the time axis
-
-     INTEGER           :: idvar    ! id of the variable
-     CHARACTER(LEN=100) :: clinfo                    ! info character
-     !---------------------------------------------------------------------
-     clinfo = 'iom_gettime, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
-     idvar = iom_varid( knumfl, cdvar )
-     !
-     ptime(:) = 0. ! default definition
-     IF ( idvar > 0 ) THEN
-        IF ( iom_file(knumfl)%ndims(idvar) == 1 ) THEN
-           IF ( iom_file(knumfl)%luld(idvar) ) THEN
-              IF ( iom_file(knumfl)%dimsz(1,idvar) == size(ptime) ) THEN
-                 CALL fliogetv( knumfl, cdvar, ptime(:), start=(/ 1 /), &
-                      &                                   count=(/ iom_file(knumfl)%dimsz(1,idvar) /) )
-              ELSE
-                 WRITE(ctmp1,*) 'error with the size of ptime ',size(ptime),iom_file(knumfl)%dimsz(1,idvar)
-                 CALL ctl_stop( trim(clinfo), trim(ctmp1) )
-              ENDIF
-           ELSE
-              CALL ctl_stop( trim(clinfo), 'variable dimension is not unlimited... use iom_get' )
-           ENDIF
-        ELSE
-           CALL ctl_stop( trim(clinfo), 'the variable has more than 1 dimension' )
-        ENDIF
-     ELSE
-        CALL ctl_stop( trim(clinfo), 'variable not found in '//iom_file(knumfl)%name )
-     ENDIF
-
+      !!--------------------------------------------------------------------
+      !!                   ***  SUBROUTINE iom_gettime  ***
+      !!
+      !! ** Purpose : read the time axis cdvar in the file 
+      !!--------------------------------------------------------------------
+      INTEGER               , INTENT(in   ) ::   knumfl   ! file Identifier
+      CHARACTER(len=*)      , INTENT(in   ) ::   cdvar    ! time axis name
+      REAL(wp), DIMENSION(:), INTENT(  out) ::   ptime    ! the time axis
+      !
+      INTEGER            ::   idvar    ! id of the variable
+      CHARACTER(LEN=100) ::   clinfo   ! info character
+      !---------------------------------------------------------------------
+      !
+      clinfo = 'iom_gettime, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
+      idvar = iom_varid( knumfl, cdvar )
+      !
+      ptime(:) = 0. ! default definition
+      IF( idvar > 0 ) THEN
+         IF( iom_file(knumfl)%ndims(idvar) == 1 ) THEN
+            IF( iom_file(knumfl)%luld(idvar) ) THEN
+               IF( iom_file(knumfl)%dimsz(1,idvar) == size(ptime) ) THEN
+                  CALL fliogetv( knumfl, cdvar, ptime(:), start=(/ 1 /),   &
+                     &                                    count=(/ iom_file(knumfl)%dimsz(1,idvar) /) )
+               ELSE
+                  WRITE(ctmp1,*) 'error with the size of ptime ',size(ptime),iom_file(knumfl)%dimsz(1,idvar)
+                  CALL ctl_stop( trim(clinfo), trim(ctmp1) )
+               ENDIF
+            ELSE
+               CALL ctl_stop( trim(clinfo), 'variable dimension is not unlimited... use iom_get' )
+            ENDIF
+         ELSE
+            CALL ctl_stop( trim(clinfo), 'the variable has more than 1 dimension' )
+         ENDIF
+      ELSE
+         CALL ctl_stop( trim(clinfo), 'variable not found in '//iom_file(knumfl)%name )
+      ENDIF
+      !
    END SUBROUTINE iom_gettime
     
     
-    
    FUNCTION iom_varid ( knumfl, cdvar, kdimsz )  
-     !!-----------------------------------------------------------------------
-     !!                  ***  FUNCTION  iom_varid  ***
-     !!
-     !! ** Purpose : get the id of a variable in a file (return 0 if not found)
-     !!
-     !! ** Method : ???
-     !!
-     !!-----------------------------------------------------------------------
-     INTEGER              , INTENT(in)            ::   knumfl   ! file Identifier
-     CHARACTER(len=*)     , INTENT(in)            ::   cdvar    ! name of the variable
-     INTEGER, DIMENSION(:), INTENT(out), OPTIONAL ::   kdimsz   ! size of the dimensions
-     !
-     INTEGER                        :: iom_varid, idvar, i_nvd, ji
-     INTEGER, DIMENSION(jpmax_dims) :: idimid
-     LOGICAL                        :: ll_fnd
-     CHARACTER(LEN=100)             :: clinfo                   ! info character
-     !!-----------------------------------------------------------------------
-     clinfo = 'iom_varid, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
-     iom_varid = 0                         ! default definition
-     IF ( PRESENT(kdimsz) ) kdimsz(:) = 0  ! default definition
-     !
-     IF ( knumfl > 0 ) THEN
-        IF( iom_file(knumfl)%iopen == 0 ) THEN 
-           CALL ctl_stop( trim(clinfo), 'the file is not open' )
-        ELSE
-           !
-           ll_fnd  = .FALSE.
-           idvar = 0
-           !
-           DO WHILE ( .NOT.ll_fnd .AND. idvar < iom_file(knumfl)%nvars )
-              idvar = idvar + 1
-              ll_fnd  = ( TRIM(cdvar) == TRIM(iom_file(knumfl)%cn_var(idvar)) )
-           END DO
-           !
-           IF( .NOT.ll_fnd ) THEN
-              idvar = idvar + 1
-              IF( idvar <= jpmax_vars ) THEN
-                 CALL flioinqv( knumfl, cdvar, ll_fnd, nb_dims = i_nvd )
-                 IF( ll_fnd ) THEN
-                    IF( i_nvd <= jpmax_dims ) THEN
-                       iom_file(knumfl)%nvars           = idvar
-                       iom_file(knumfl)%cn_var(idvar) = trim(cdvar)
-                       iom_file(knumfl)%ndims(idvar)  = i_nvd
-                       CALL flioinqv( knumfl, cdvar, ll_fnd, &
-                            &         len_dims = iom_file(knumfl)%dimsz(1:i_nvd,idvar), &
-                            &         id_dims  = idimid(1:i_nvd) )
-                       DO ji = 1, i_nvd
-                          IF ( idimid(ji) == iom_file(knumfl)%iduld ) iom_file(knumfl)%luld(idvar) = .TRUE.
-                       END DO
-                       !----------
-                       !---------- Deal with scale_factor and offset
-                       CALL flioinqa( knumfl, cdvar, 'scale_factor', ll_fnd )
-                       IF (ll_fnd) THEN
-                          CALL fliogeta( knumfl, cdvar, 'scale_factor', iom_file(knumfl)%scf(idvar) )
-                       ELSE
-                          iom_file(knumfl)%scf(idvar) = 1.
-                       END IF
-                       CALL flioinqa( knumfl, cdvar, 'offset', ll_fnd )
-                       IF( ll_fnd ) THEN
-                          CALL fliogeta( knumfl, cdvar, 'offset', iom_file(knumfl)%ofs(idvar) )
-                       ELSE
-                          iom_file(knumfl)%ofs(idvar) = 0.
-                       END IF
-                       !
-                       iom_varid = idvar
-                       IF ( PRESENT(kdimsz) ) THEN 
-                          IF ( i_nvd == size(kdimsz) ) THEN
-                             kdimsz(:) = iom_file(knumfl)%dimsz(1:i_nvd,idvar)
-                          ELSE
-                             WRITE(ctmp1,*) i_nvd, size(kdimsz)
-                             CALL ctl_stop( trim(clinfo), 'error in kdimsz size'//trim(ctmp1) )
-                          ENDIF
-                       ENDIF
-                    ELSE
-                       CALL ctl_stop( trim(clinfo), 'Too many dimensions in the file '//iom_file(knumfl)%name,   &
-                            &           'increase the parameter jpmax_vars')
-                    ENDIF
-                 ELSE  
-                    CALL ctl_warn( trim(clinfo), 'Variable '//trim(cdvar)// &
-                         & ' is not found in the file '//trim(iom_file(knumfl)%name) )
-                 ENDIF
-              ELSE
-                 CALL ctl_stop( trim(clinfo), 'Too many variables in the file '//iom_file(knumfl)%name,   &
-                      &           'increase the parameter jpmax_vars')
-              ENDIF
-           ELSE
-              iom_varid = idvar
-              IF ( PRESENT(kdimsz) ) THEN 
-                 i_nvd = iom_file(knumfl)%ndims(idvar)
-                 IF ( i_nvd == size(kdimsz) ) THEN
-                    kdimsz(:) = iom_file(knumfl)%dimsz(1:i_nvd,idvar)
-                 ELSE
-                    WRITE(ctmp1,*) i_nvd, size(kdimsz)
-                    CALL ctl_stop( trim(clinfo), 'error in kdimsz size'//trim(ctmp1) )
-                 ENDIF
-              ENDIF
-           ENDIF
-        ENDIF
-     ENDIF
-
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  iom_varid  ***
+      !!
+      !! ** Purpose : get the id of a variable in a file (return 0 if not found)
+      !!-----------------------------------------------------------------------
+      INTEGER              , INTENT(in   )           ::   knumfl   ! file Identifier
+      CHARACTER(len=*)     , INTENT(in   )           ::   cdvar    ! name of the variable
+      INTEGER, DIMENSION(:), INTENT(  out), OPTIONAL ::   kdimsz   ! size of the dimensions
+      !
+      INTEGER                        ::   ji                       ! dummy loop index
+      INTEGER                        ::   iom_varid, idvar, i_nvd
+      INTEGER, DIMENSION(jpmax_dims) ::   idimid
+      LOGICAL                        ::   ll_fnd
+      CHARACTER(LEN=100)             ::   clinfo                   ! info character
+      !!-----------------------------------------------------------------------
+      clinfo = 'iom_varid, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar)
+      iom_varid = 0                         ! default definition
+      IF( PRESENT(kdimsz) ) kdimsz(:) = 0   ! default definition
+      !
+      IF( knumfl > 0 ) THEN
+         IF( iom_file(knumfl)%iopen == 0 ) THEN 
+            CALL ctl_stop( trim(clinfo), 'the file is not open' )
+         ELSE
+            ll_fnd  = .FALSE.
+            idvar = 0
+            !
+            DO WHILE ( .NOT.ll_fnd .AND. idvar < iom_file(knumfl)%nvars )
+               idvar = idvar + 1
+               ll_fnd  = ( TRIM(cdvar) == TRIM(iom_file(knumfl)%cn_var(idvar)) )
+            END DO
+            !
+            IF( .NOT.ll_fnd ) THEN
+               idvar = idvar + 1
+               IF( idvar <= jpmax_vars ) THEN
+                  CALL flioinqv( knumfl, cdvar, ll_fnd, nb_dims = i_nvd )
+                  IF( ll_fnd ) THEN
+                     IF( i_nvd <= jpmax_dims ) THEN
+                        iom_file(knumfl)%nvars           = idvar
+                        iom_file(knumfl)%cn_var(idvar) = trim(cdvar)
+                        iom_file(knumfl)%ndims(idvar)  = i_nvd
+                        CALL flioinqv( knumfl, cdvar, ll_fnd,   &
+                           &           len_dims = iom_file(knumfl)%dimsz(1:i_nvd,idvar), &
+                           &           id_dims  = idimid(1:i_nvd) )
+                        DO ji = 1, i_nvd
+                           IF( idimid(ji) == iom_file(knumfl)%iduld ) iom_file(knumfl)%luld(idvar) = .TRUE.
+                        END DO
+                        !----------
+                        !---------- Deal with scale_factor and offset
+                        CALL flioinqa( knumfl, cdvar, 'scale_factor', ll_fnd )
+                        IF( ll_fnd) THEN
+                           CALL fliogeta( knumfl, cdvar, 'scale_factor', iom_file(knumfl)%scf(idvar) )
+                        ELSE
+                           iom_file(knumfl)%scf(idvar) = 1.
+                        END IF
+                        CALL flioinqa( knumfl, cdvar, 'offset', ll_fnd )
+                        IF( ll_fnd ) THEN
+                           CALL fliogeta( knumfl, cdvar, 'offset', iom_file(knumfl)%ofs(idvar) )
+                        ELSE
+                           iom_file(knumfl)%ofs(idvar) = 0.
+                        END IF
+                        !
+                        iom_varid = idvar
+                        IF( PRESENT(kdimsz) ) THEN 
+                           IF( i_nvd == size(kdimsz) ) THEN
+                              kdimsz(:) = iom_file(knumfl)%dimsz(1:i_nvd,idvar)
+                           ELSE
+                              WRITE(ctmp1,*) i_nvd, size(kdimsz)
+                              CALL ctl_stop( trim(clinfo), 'error in kdimsz size'//trim(ctmp1) )
+                           ENDIF
+                        ENDIF
+                     ELSE
+                        CALL ctl_stop( trim(clinfo), 'Too many dimensions in the file '//iom_file(knumfl)%name,   &
+                           &                        'increase the parameter jpmax_vars')
+                     ENDIF
+!!$                  ELSE  
+!!$                     CALL ctl_warn( trim(clinfo), 'Variable '//trim(cdvar)// &
+!!$                        &                         ' is not found in the file '//trim(iom_file(knumfl)%name) )
+                  ENDIF
+               ELSE
+                  CALL ctl_stop( trim(clinfo), 'Too many variables in the file '//iom_file(knumfl)%name,   &
+                     &                         'increase the parameter jpmax_vars')
+               ENDIF
+            ELSE
+               iom_varid = idvar
+               IF( PRESENT(kdimsz) ) THEN 
+                  i_nvd = iom_file(knumfl)%ndims(idvar)
+                  IF( i_nvd == size(kdimsz) ) THEN
+                     kdimsz(:) = iom_file(knumfl)%dimsz(1:i_nvd,idvar)
+                  ELSE
+                     WRITE(ctmp1,*) i_nvd, size(kdimsz)
+                     CALL ctl_stop( trim(clinfo), 'error in kdimsz size'//trim(ctmp1) )
+                  ENDIF
+               ENDIF
+            ENDIF
+         ENDIF
+      ENDIF
+      !
    END FUNCTION iom_varid
 
-  
-   FUNCTION iom_get_gblatt( knumfl, kinfonum )
-      !!-----------------------------------------------------------------------
-      !!                  ***  FUNCTION  iom_get_gblatt  ***
+   !!----------------------------------------------------------------------
+   !!                   INTERFACE iom_rstput
+   !!----------------------------------------------------------------------
+   SUBROUTINE iom_rstput_0d( kt, kwrite, knumfl, cdvar, pvar )
+      INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
+      INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
+      INTEGER         , INTENT(in)                         ::   knumfl   ! Identifier of the file 
+      CHARACTER(len=*), INTENT(in)                         ::   cdvar    ! time axis name
+      REAL(wp)        , INTENT(in)                         ::   pvar     ! read field
+      IF( knumfl > 0 )   CALL iom_rstput_0123d( kt, kwrite, knumfl, cdvar, pv_r0d = pvar )
+   END SUBROUTINE iom_rstput_0d
+
+   SUBROUTINE iom_rstput_1d( kt, kwrite, knumfl, cdvar, pvar )
+      INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
+      INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
+      INTEGER         , INTENT(in)                         ::   knumfl   ! Identifier of the file 
+      CHARACTER(len=*), INTENT(in)                         ::   cdvar    ! time axis name
+      REAL(wp)        , INTENT(in), DIMENSION(        jpk) ::   pvar     ! read field
+      IF( knumfl > 0 )    CALL iom_rstput_0123d( kt, kwrite, knumfl, cdvar, pv_r1d = pvar )
+   END SUBROUTINE iom_rstput_1d
+
+   SUBROUTINE iom_rstput_2d( kt, kwrite, knumfl, cdvar, pvar )
+      INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
+      INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
+      INTEGER         , INTENT(in)                         ::   knumfl   ! Identifier of the file 
+      CHARACTER(len=*), INTENT(in)                         ::   cdvar    ! time axis name
+      REAL(wp)        , INTENT(in), DIMENSION(jpi,jpj    ) ::   pvar     ! read field
+      IF( knumfl > 0 )   CALL iom_rstput_0123d( kt, kwrite, knumfl, cdvar, pv_r2d = pvar )
+   END SUBROUTINE iom_rstput_2d
+
+   SUBROUTINE iom_rstput_3d( kt, kwrite, knumfl, cdvar, pvar )
+      INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
+      INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
+      INTEGER         , INTENT(in)                         ::   knumfl   ! Identifier of the file 
+      CHARACTER(len=*), INTENT(in)                         ::   cdvar    ! time axis name
+      REAL(wp)        , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvar     ! read field
+      IF( knumfl > 0 )   CALL iom_rstput_0123d( kt, kwrite, knumfl, cdvar, pv_r3d = pvar )
+   END SUBROUTINE iom_rstput_3d
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE iom_rstput_0123d( kt, kwrite, knumfl, cdvar ,   &
+      &                             pv_r0d, pv_r1d, pv_r2d, pv_r3d )
+      !!--------------------------------------------------------------------
+      !!                   ***  SUBROUTINE  iom_rstput  ***
       !!
-      !! ** Purpose : ???
-      !!
-      !! ** Method : ???
-      !!
-      !!-----------------------------------------------------------------------
-      INTEGER,INTENT(in) :: knumfl
-      INTEGER,          intent(in) :: kinfonum
-      !
-      CHARACTER(LEN=10) :: clinfo
-      REAL(wp) :: iom_get_gblatt
-      !!-----------------------------------------------------------------------
-
-      WRITE(clinfo,*) kinfonum
-      clinfo = 'info'//trim(adjustl(clinfo))
-      CALL fliogeta (knumfl, "?", clinfo, iom_get_gblatt)
-
-   END FUNCTION iom_get_gblatt
-
+      !! ** Purpose : read the time axis cdvar in the file 
+      !!--------------------------------------------------------------------
+      INTEGER                   , INTENT(in)           ::   kt       ! ocean time-step
+      INTEGER                   , INTENT(in)           ::   kwrite   ! writing time-step
+      INTEGER                   , INTENT(in)           ::   knumfl   ! Identifier of the file 
+      CHARACTER(len=*)          , INTENT(in)           ::   cdvar    ! time axis name
+      REAL(wp)                  , INTENT(in), OPTIONAL ::   pv_r0d   ! read field
+      REAL(wp), DIMENSION(:)    , INTENT(in), OPTIONAL ::   pv_r1d   ! read field
+      REAL(wp), DIMENSION(:,:)  , INTENT(in), OPTIONAL ::   pv_r2d   ! read field
+      REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL ::   pv_r3d   ! read field
+      !
+      INTEGER               :: idims, idvar        
+      INTEGER               :: ix1, ix2, iy1, iy2       
+      INTEGER, DIMENSION(4) :: idimsz, idimid     
+      CHARACTER(LEN=100)    :: clinfo        ! info character
+      !---------------------------------------------------------------------
+      !
+      clinfo = '          iom_rstput_0123d, file: '//TRIM(iom_file(knumfl)%name)//', var: '//TRIM(cdvar)
+
+      ! define dimension variables if it is not already done
+      ! ==========================
+      IF( iom_file(knumfl)%nvars == 0 ) THEN
+         ! define the dimension variables if it is not already done
+         CALL fliodefv( knumfl,'nav_lon', (/1,2/), v_t=flio_r4   , axis='X',   &
+            &                  long_name="Longitude", units="degrees_east" )
+         CALL fliodefv( knumfl,'nav_lat', (/1,2/), v_t=flio_r4   , axis='Y',   &
+            &                  long_name="Latitude", units="degrees_north" )
+         CALL fliodefv( knumfl,'nav_lev', (/3/)  , v_t=flio_i4   , axis='Z',   &
+            &                  long_name="Model levels",units="model_levels")
+         CALL fliodefv( knumfl,'time_counter', (/4/), v_t=flio_r4, axis='T',   &
+            &                  long_name="Time axis", units='seconds since 0001-01-01 00:00:00' )
+         ! update informations structure related the dimension variable we just added...
+         iom_file(knumfl)%nvars       = 4
+         iom_file(knumfl)%luld(1:4)   = (/ .FALSE., .FALSE., .FALSE., .TRUE. /)
+         iom_file(knumfl)%cn_var(1:3) = (/ 'nav_lon', 'nav_lat', 'nav_lev' /)
+         iom_file(knumfl)%cn_var(4)   = 'time_counter'
+         iom_file(knumfl)%ndims(1:4)  = (/ 2, 2, 1, 1 /)  
+         IF(lwp) WRITE(numout,*) TRIM(clinfo)//' define dimension variables done'
+      ENDIF
+
+      ! define the data if it is not already done
+      ! ===============
+      idvar = iom_varid( knumfl, cdvar )
+      IF( idvar <= 0 ) THEN
+         ! variable definition
+         IF(     PRESENT(pv_r0d) ) THEN   ;   idims = 0
+         ELSEIF( PRESENT(pv_r1d) ) THEN   ;   idims = 2   ;   idimid(1:idims) = (/    3,4/)
+         ELSEIF( PRESENT(pv_r2d) ) THEN   ;   idims = 3   ;   idimid(1:idims) = (/1,2  ,4/)
+         ELSEIF( PRESENT(pv_r3d) ) THEN   ;   idims = 4   ;   idimid(1:idims) = (/1,2,3,4/)
+         ENDIF
+         IF( PRESENT(pv_r0d) ) THEN   ;   CALL fliodefv (knumfl, cdvar                 , v_t = flio_r8)
+         ELSE                         ;   CALL fliodefv (knumfl, cdvar, idimid(1:idims), v_t = flio_r8)
+         ENDIF
+         ! update informations structure related the new variable we want to add...
+         idvar                          = iom_file(knumfl)%nvars + 1
+         iom_file(knumfl)%nvars         = idvar
+         iom_file(knumfl)%cn_var(idvar) = TRIM(cdvar)
+         iom_file(knumfl)%scf(idvar)    = 1.
+         iom_file(knumfl)%ofs(idvar)    = 0.
+         iom_file(knumfl)%ndims(idvar)  = idims
+         IF( .NOT. PRESENT(pv_r0d) ) THEN
+            iom_file(knumfl)%luld(idvar) = .TRUE.
+            CALL flioinqf( knumfl, ln_dim = idimsz )
+            iom_file(knumfl)%dimsz(1:idims-1,idvar) = idimsz(idimid(1:idims-1))
+         ENDIF
+         IF(lwp) WRITE(numout,*) TRIM(clinfo)//' defined ok'
+      ENDIF
+
+      ! time step kwrite : write the variable
+      IF( kt == kwrite ) THEN
+         ! on what kind of domain must the data be written?
+         IF( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) THEN
+            idimsz(1:2) = iom_file(knumfl)%dimsz(1:2,idvar)
+            IF(     idimsz(1) == (nlei - nldi + 1) .AND. idimsz(2) == (nlej - nldj + 1) ) THEN
+               ix1 = nldi   ;   ix2 = nlei   ;   iy1 = nldj   ;   iy2 = nlej
+            ELSEIF( idimsz(1) == nlci              .AND. idimsz(2) == nlcj              ) THEN
+               ix1 = 1      ;   ix2 = nlci   ;   iy1 = 1      ;   iy2 = nlcj
+            ELSEIF( idimsz(1) == jpi               .AND. idimsz(2) == jpj               ) THEN
+               ix1 = 1      ;   ix2 = jpi    ;   iy1 = 1      ;   iy2 = jpj
+            ELSE 
+               CALL ctl_stop( 'iom_rstput_0123d: should have been an impossible case...' )
+            ENDIF
+
+            ! write dimension variables if it is not already done
+            ! =============
+            IF( iom_file(knumfl)%dimsz(1, 1) == 0 ) THEN
+               CALL flioputv( knumfl, 'nav_lon'     , glamt(ix1:ix2, iy1:iy2) )
+               CALL flioputv( knumfl, 'nav_lat'     , gphit(ix1:ix2, iy1:iy2) )
+               CALL flioputv( knumfl, 'nav_lev'     , gdept_0 )
+               CALL flioputv( knumfl, 'time_counter', kt )   ! +++ WRONG VALUE: to be improved but not really useful...
+               ! update informations structure related the dimension variable
+               iom_file(knumfl)%dimsz(1:2, 1) = idimsz(1:2)
+               iom_file(knumfl)%dimsz(1:2, 2) = idimsz(1:2)
+               iom_file(knumfl)%dimsz(1, 3:4) = (/idimsz(3), 1/)
+               IF(lwp) WRITE(numout,*) TRIM(clinfo)//' write dimension variables done'
+            ENDIF
+         ENDIF
+
+         ! write the data
+         ! =============
+         IF(     PRESENT(pv_r0d) ) THEN   ;   CALL flioputv( knumfl, cdvar, pv_r0d                      )
+         ELSEIF( PRESENT(pv_r1d) ) THEN   ;   CALL flioputv( knumfl, cdvar, pv_r1d(                  :) )
+         ELSEIF( PRESENT(pv_r2d) ) THEN   ;   CALL flioputv( knumfl, cdvar, pv_r2d(ix1:ix2, iy1:iy2   ) )
+         ELSEIF( PRESENT(pv_r3d) ) THEN   ;   CALL flioputv( knumfl, cdvar, pv_r3d(ix1:ix2, iy1:iy2, :) )
+         ENDIF
+         ! add 1 to the size of the temporal dimension (not really useful...)
+         IF( iom_file(knumfl)%luld(idvar) )   iom_file(knumfl)%dimsz(iom_file(knumfl)%ndims(idvar), idvar)    &
+            &                               = iom_file(knumfl)%dimsz(iom_file(knumfl)%ndims(idvar), idvar) + 1
+         IF(lwp) WRITE(numout,*) TRIM(clinfo)//' written ok'
+      ENDIF
+      !     
+   END SUBROUTINE iom_rstput_0123d
+   
    !!======================================================================
 END MODULE iom
Index: trunk/NEMO/OPA_SRC/istate.F90
===================================================================
--- trunk/NEMO/OPA_SRC/istate.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/istate.F90	(revision 508)
@@ -4,4 +4,14 @@
    !! Ocean state   :  initial state setting
    !!=====================================================================
+   !! History :   4.0  !  89-12  (P. Andrich)  Original code
+   !!             5.0  !  91-11  (G. Madec)  rewritting
+   !!             6.0  !  96-01  (G. Madec)  terrain following coordinates
+   !!             8.0  !  01-09  (M. Levy, M. Ben Jelloul)  istate_eel
+   !!             8.0  !  01-09  (M. Levy, M. Ben Jelloul)  istate_uvg
+   !!             9.0  !  03-08  (G. Madec)  F90: Free form, modules
+   !!             9.0  !  03-09  (G. Madec, C. Talandier)  add EEL R5
+   !!             9.0  !  04-05  (A. Koch-Larrouy)  istate_gyre 
+   !!             9.0  !  06-07  (S. Masson)  distributed restart using iom
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -13,5 +23,4 @@
    !!   istate_uvg    : initial velocity in geostropic balance
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and active tracers 
    USE dom_oce         ! ocean space and time domain 
@@ -19,5 +28,4 @@
    USE ldftra_oce      ! ocean active tracers: lateral physics
    USE zdf_oce         ! ocean vertical physics
-   USE in_out_manager  ! I/O manager
    USE phycst          ! physical constants
    USE wzvmod          ! verctical velocity               (wzv     routine)
@@ -26,10 +34,11 @@
    USE restart         ! ocean restart                   (rst_read routine)
    USE solisl          ! ???
-
+   USE in_out_manager  ! I/O manager
+   USE iom
+   
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC istate_init   ! routine called by step.F90
+   PUBLIC   istate_init   ! routine called by step.F90
 
    !! * Substitutions
@@ -37,7 +46,7 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
    !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -48,46 +57,14 @@
       !!                   ***  ROUTINE istate_init  ***
       !! 
-      !! ** Purpose :   Initialization of the dynamics and tracers.
-      !!
-      !! ** Method  :
-      !!
-      !! History :
-      !!   4.0  !  91-03  ()  Original code
-      !!        !  91-11  (G. Madec)
-      !!   9.0  !  03-09  (G. Madec)  F90: Free form, modules, orthogonality
-      !!----------------------------------------------------------------------
-      USE iom
-      !! * Local declarations
-      !CT INTEGER ::   inum
-      !!----------------------------------------------------------------------
-
-
-      ! Initialization to zero
-      ! ----------------------
-
-      !     before fields       !       now fields        !      after fields       !
-      ;   ub   (:,:,:) = 0.e0   ;   un   (:,:,:) = 0.e0   ;   ua   (:,:,:) = 0.e0
-      ;   vb   (:,:,:) = 0.e0   ;   vn   (:,:,:) = 0.e0   ;   va   (:,:,:) = 0.e0
-      ;                         ;   wn   (:,:,:) = 0.e0   ;
-      ;   rotb (:,:,:) = 0.e0   ;   rotn (:,:,:) = 0.e0   ;
-      ;   hdivb(:,:,:) = 0.e0   ;   hdivn(:,:,:) = 0.e0   ;
-
-      ;   tb   (:,:,:) = 0.e0   ;   tn   (:,:,:) = 0.e0   ;   ta   (:,:,:) = 0.e0
-      ;   sb   (:,:,:) = 0.e0   ;   sn   (:,:,:) = 0.e0   ;   sa   (:,:,:) = 0.e0
+      !! ** Purpose :   Initialization of the dynamics and tracer fields.
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
 
       rhd  (:,:,:) = 0.e0
       rhop (:,:,:) = 0.e0
       rn2  (:,:,:) = 0.e0 
-
-#if defined key_dynspg_rl
-      ! rigid-lid formulation
-      bsfb(:,:) = 0.e0      ! before barotropic stream-function
-      bsfn(:,:) = 0.e0      ! now    barotropic stream-function
-      bsfd(:,:) = 0.e0      ! barotropic stream-function trend
-#endif
-      ! free surface formulation
-      sshb(:,:) = 0.e0      ! before sea-surface height
-      sshn(:,:) = 0.e0      ! now    sea-surface height
-
 
       IF( ln_rstart ) THEN                    ! Restart from a file
@@ -100,4 +77,11 @@
          neuler = 0                              ! Set time-step indicator at nit000 (euler forward)
          adatrj = 0._wp
+         !                                       ! Initialization of ocean to zero
+         !     before fields       !       now fields          
+         ;   ub   (:,:,:) = 0.e0   ;   un   (:,:,:) = 0.e0   
+         ;   vb   (:,:,:) = 0.e0   ;   vn   (:,:,:) = 0.e0    
+         ;   rotb (:,:,:) = 0.e0   ;   rotn (:,:,:) = 0.e0  
+         ;   hdivb(:,:,:) = 0.e0   ;   hdivn(:,:,:) = 0.e0  
+         !
          IF( cp_cfg == 'eel' ) THEN
             CALL istate_eel                      ! EEL   configuration : start from pre-defined
@@ -107,11 +91,5 @@
             !                                    !                       and salinity fields 
          ELSE
-         !                                       ! Other configurations: Initial temperature and salinity fields
-
-         !CT CALL iom_open ('ssh', inum) 
-         !CT CALL iom_get( inum, jpdom_local, 'sshb', sshb )     ! free surface formulation (ssh)
-         !CT sshn(:,:) = sshb(:,:)
-         !CT CALL iom_close (inum)
-
+            !                                    ! Other configurations: Initial temperature and salinity fields
 #if defined key_dtatem
             CALL dta_tem( nit000 )                  ! read 3D temperature data
@@ -120,5 +98,5 @@
 #else
             IF(lwp) WRITE(numout,*)                 ! analytical temperature profile
-            IF(lwp) WRITE(numout,*)' Temperature initialization using an analytic profile'
+            IF(lwp) WRITE(numout,*)'             Temperature initialization using an analytic profile'
             CALL istate_tem
 #endif
@@ -130,5 +108,5 @@
             ! No salinity data
             IF(lwp)WRITE(numout,*)                  ! analytical salinity profile
-            IF(lwp)WRITE(numout,*)' Salinity initialisation using a constant value'
+            IF(lwp)WRITE(numout,*)'             Salinity initialisation using a constant value'
             CALL istate_sal
 #endif
@@ -139,5 +117,5 @@
       !                                       ! -----------------
       CALL wzv( nit000 )                         ! from horizontal divergence
-
+      !
    END SUBROUTINE istate_init
 
@@ -153,14 +131,8 @@
       !!
       !! References :  Philander ???
-      !!
-      !! History :
-      !!   4.0  !  89-12  (P. Andrich)  Original code
-      !!   6.0  !  96-01  (G. Madec)  terrain following coordinates
-      !!   9.0  !  02-09  (G. Madec)  F90: Free form
-      !!----------------------------------------------------------------------
-      !! * Local declarations
+      !!----------------------------------------------------------------------
       INTEGER :: ji, jj, jk
       !!----------------------------------------------------------------------
-
+      !
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'istate_tem : initial temperature profile'
@@ -181,5 +153,5 @@
          &                 1     , jpi   , 5     , 1     , jpk   ,   &
          &                 1     , 1.    , numout                  )
-
+      !
    END SUBROUTINE istate_tem
 
@@ -194,10 +166,5 @@
       !!              
       !! ** Action  :   Initialize sn and sb
-      !!
-      !! History :
-      !!   4.0  !  89-12  (P. Andrich)  Original code
-      !!   8.5  !  02-09  (G. Madec)  F90: Free form
-      !!----------------------------------------------------------------------
-      !! * Local declarations
+      !!----------------------------------------------------------------------
       REAL(wp) ::   zsal = 35.50_wp
       !!----------------------------------------------------------------------
@@ -224,30 +191,19 @@
       !!              - set velocity field including horizontal divergence
       !!                and relative vorticity fields
-      !!
-      !! History :
-      !!   8.0  !  01-09  (M. Levy, M. Ben Jelloul)  read file for EEL 2 & 6
-      !!   9.0  !  03-09  (G. Madec, C. Talandier)  EEL 5
-      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-      !!----------------------------------------------------------------------
-      !! * Modules used
+      !!----------------------------------------------------------------------
       USE eosbn2     ! eq. of state, Brunt Vaisala frequency (eos     routine)
       USE divcur     ! hor. divergence & rel. vorticity      (div_cur routine)
       USE iom
  
-      !! * Local declarations
       INTEGER  ::   inum              ! temporary logical unit
       INTEGER  ::   ji, jj, jk        ! dummy loop indices
       INTEGER  ::   ijloc
-      REAL(wp) ::   &
-         zh1, zh2, zslope, zcst,   &  ! temporary scalars
-         zfcor
-      REAL(wp) ::   &
-         zt1  = 12._wp,            &  ! surface temperature value (EEL R5)
-         zt2  =  2._wp,            &  ! bottom  temperature value (EEL R5)
-         zsal = 35.5_wp,           &  ! constant salinity (EEL R2, R5 and R6)
-         zueel = 0.1_wp               ! constant uniform zonal velocity (EEL R5)
+      REAL(wp) ::   zh1, zh2, zslope, zcst, zfcor   ! temporary scalars
+      REAL(wp) ::   zt1  = 12._wp,               &  ! surface temperature value (EEL R5)
+         &          zt2  =  2._wp,               &  ! bottom  temperature value (EEL R5)
+         &          zsal = 35.5_wp,              &  ! constant salinity (EEL R2, R5 and R6)
+         &          zueel = 0.1_wp                  ! constant uniform zonal velocity (EEL R5)
 # if ! defined key_dynspg_rl
-      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   &
-         zssh                         ! initial ssh over the global domain
+      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zssh  ! initial ssh over the global domain
 # endif
       !!----------------------------------------------------------------------
@@ -389,17 +345,8 @@
       !! ** Method  : - set temprature field
       !!              - set salinity field
-      !!
-      !! ** History :                                     
-      !!      9.0  !  04-05  (A. Koch-Larrouy)  Original code 
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE iom
-
-      !! * Local variables
-      INTEGER  ::   inum              ! temporary logical unit
-      INTEGER, PARAMETER ::   &
-         ntsinit = 0         ! (0/1) (analytical/input data files) T&S initialization
-
+      !!----------------------------------------------------------------------
       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
       !!----------------------------------------------------------------------
 
@@ -469,7 +416,5 @@
       ENDIF
 
-     
    END SUBROUTINE istate_gyre
-
 
 
@@ -484,11 +429,5 @@
       !!      surface to the bottom.
       !!                 p=integral [ rau*g dz ]
-      !!
-      !! History :
-      !!   8.1  !  01-09  (M. Levy, M. Ben Jelloul)  Original code
-      !!   8.5  !  02-09  (G. Madec)  F90: Free form
-      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-      !!----------------------------------------------------------------------
-      !! * Modules used
+      !!----------------------------------------------------------------------
       USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine)
       USE dynspg          ! surface pressure gradient             (dyn_spg routine)
@@ -496,12 +435,8 @@
       USE lbclnk          ! ocean lateral boundary condition (or mpp link)
 
-      !! * Local declarations
       INTEGER ::   ji, jj, jk        ! dummy loop indices
       INTEGER ::   indic             ! ???
-      REAL(wp) ::   &
-         zmsv, zphv, zmsu, zphu,  &  ! temporary scalars
-         zalfg
-      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   &
-         zprn                        ! workspace
+      REAL(wp) ::   zmsv, zphv, zmsu, zphu, zalfg     ! temporary scalars
+      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zprn     ! workspace
       !!----------------------------------------------------------------------
 
@@ -514,9 +449,8 @@
 
       zalfg = 0.5 * grav * rau0
-      ! Surface value
-      zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) )
-
-      ! Vertical integration from the surface
-      DO jk = 2, jpkm1
+      
+      zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) )       ! Surface value
+
+      DO jk = 2, jpkm1                                              ! Vertical integration from the surface
          zprn(:,:,jk) = zprn(:,:,jk-1)   &
             &         + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
@@ -525,5 +459,4 @@
       ! Compute geostrophic balance
       ! ---------------------------
-
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
@@ -571,15 +504,9 @@
       ! after initializing u and v, we need to calculate the initial streamfunction bsf.
       ! Otherwise, only the trend will be computed and the model will blow up (inconsistency).
-      
       ! to do that, we call dyn_spg with a special trick:
-      ! we fill ua and va with the velocities divided by dt,
-      ! and the streamfunction will be brought to the right
-      ! value assuming the velocities have been set up in
-      ! one time step.
-      ! we then set bsfd to zero (first guess for next step
-      ! is d(psi)/dt = 0.)
-
-      !  sets up s false trend to calculate the barotropic
-      !  streamfunction.
+      ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the
+      ! right value assuming the velocities have been set up in one time step.
+      ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.)
+      !  sets up s false trend to calculate the barotropic streamfunction.
 
       ua(:,:,:) = ub(:,:,:) / rdt
@@ -612,5 +539,5 @@
       hdivb(:,:,:) = hdivn(:,:,:)       ! set the before to the now value
       rotb (:,:,:) = rotn (:,:,:)       ! set the before to the now value
-
+      !
    END SUBROUTINE istate_uvg
 
Index: trunk/NEMO/OPA_SRC/opa.F90
===================================================================
--- trunk/NEMO/OPA_SRC/opa.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/opa.F90	(revision 508)
@@ -48,5 +48,4 @@
    USE obc_par         ! open boundary cond. parameters
    USE obcini          ! open boundary cond. initialization (obc_ini routine)
-   USE solver          ! solver initialization          (solver_init routine)
    USE istate          ! initial state setting          (istate_init routine)
    USE eosbn2          ! equation of state            (eos bn2 routine)
@@ -69,5 +68,4 @@
 
    USE step            ! OPA time-stepping                  (stp     routine)
-   USE dynspg_oce      ! Control choice of surface pressure gradient schemes
    USE prtctl          ! Print control                 (prt_ctl_init routine)
    USE ini1d           ! re-initialization of u-v mask for the 1D configuration
@@ -245,8 +243,4 @@
       CALL istate_init                      ! ocean initial state (Dynamics and tracers)
 
-      IF( lk_dynspg_flt .OR. lk_dynspg_rl ) THEN
-         CALL solver_init( nit000 )         ! Elliptic solver
-      ENDIF
-
 !!add
                        CALL eos( tb, sb, rhd, rhop )        ! before potential and in situ densities
Index: trunk/NEMO/OPA_SRC/par_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/par_oce.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/par_oce.F90	(revision 508)
@@ -26,7 +26,7 @@
 
    INTEGER, PUBLIC, PARAMETER ::    &  !:
-      jpni   = 1,                   &  !: number of processors following i 
-      jpnj   = 1,                   &  !: number of processors following j
-      jpnij  = 1,                   &  !: nb of local domain = nb of processors 
+      jpni   = 2,                   &  !: number of processors following i 
+      jpnj   = 2,                   &  !: number of processors following j
+      jpnij  = 4,                   &  !: nb of local domain = nb of processors 
       !                                !  ( <= jpni x jpnj )
       jpr2di = 0,                   &  !: number of columns for extra outer halo 
Index: trunk/NEMO/OPA_SRC/restart.F90
===================================================================
--- trunk/NEMO/OPA_SRC/restart.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/restart.F90	(revision 508)
@@ -3,43 +3,88 @@
    !!                     ***  MODULE  restart  ***
    !! Ocean restart :  write the ocean restart file
-   !!=====================================================================
-
-   !!----------------------------------------------------------------------
-   !!   rst_write  : write of the restart file
-   !!   rst_read   : read the restart file
-   !!----------------------------------------------------------------------
-   !! * Modules used
+   !!======================================================================
+   !! History :        !  99-11  (M. Imbard)  Original code
+   !!             8.5  !  02-08  (G. Madec)  F90: Free form
+   !!             9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
+   !!             9.0  !  06-07  (S. Masson)  use IOM for restart
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   rst_opn    : open the ocean restart file
+   !!   rst_write  : write the ocean restart file
+   !!   rst_read   : read the ocean restart file
+   !!----------------------------------------------------------------------
    USE dom_oce         ! ocean space and time domain
    USE oce             ! ocean dynamics and tracers 
    USE phycst          ! physical constants
-   USE in_out_manager  ! I/O manager
    USE daymod          ! calendar
-   USE sol_oce         ! ocean elliptic solver
-   USE zdf_oce         ! ???
-   USE zdftke          ! turbulent kinetic energy scheme
    USE ice_oce         ! ice variables
    USE blk_oce         ! bulk variables
-   USE flx_oce         ! sea-ice/ocean forcings variables
-   USE dynspg_oce      ! free surface time splitting scheme variables
    USE cpl_oce, ONLY : lk_cpl              !
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O module
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC rst_write  ! routine called by step.F90
-   PUBLIC rst_read   ! routine called by inidtr.F90
-
-   !! * Module variables
-   CHARACTER (len=48) ::   &
-      crestart = 'initial.nc'   ! restart file name
-   !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   PUBLIC   rst_opn    ! routine called by step module
+   PUBLIC   rst_write  ! routine called by step module
+   PUBLIC   rst_read   ! routine called by opa  module
+
+   LOGICAL, PUBLIC ::   lrst_oce         !: logical to control the oce restart write 
+   INTEGER, PUBLIC ::   nitrst           !: time step at which restart file should be written
+   INTEGER, PUBLIC ::   numror, numrow   !: logical unit for cean restart (read and write)
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2006) 
    !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence (modipsl/doc/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   ! ice output restart file name
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 ) THEN   ! default initialization, to do: should be read in the namelist...
+         nitrst = nitend        ! to do: should be read in the namelist in a cleaver way...
+         lrst_oce = .FALSE.
+      ENDIF
+      
+      IF    ( kt == nitrst-1 .AND. lrst_oce         ) THEN
+         CALL ctl_stop( 'rst_opn: we cannot create an ocean restart at every time step' )
+         numrow = 0
+      ELSEIF( kt == nitrst-1 .OR.  nitend == nit000 ) THEN   ! beware if model runs only one time step
+         ! beware of the format used to write kt (default is i8.8, that should be large enough)
+         IF( nitrst > 1.0e9 ) THEN   
+            WRITE(clkt,*) nitrst
+         ELSE
+            WRITE(clkt,'(i8.8)') nitrst
+         ENDIF
+         ! create the file
+         IF(lwp) WRITE(numout,*)
+         clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart"
+         IF(lwp) WRITE(numout,*) '             open ocean restart.output NetCDF file: '//clname
+         CALL iom_open( clname, numrow, ldwrt = .TRUE. )
+         lrst_oce = .TRUE.
+      ENDIF
+      !
+   END SUBROUTINE rst_opn
+
 
 #if  ( defined key_mpp_mpi   ||   defined key_mpp_shmem ) && defined key_dimgout
@@ -66,153 +111,56 @@
       !! ** Purpose :   Write restart fields in NetCDF format
       !!
-      !! ** Method  :   Write in numwrs file each nstock time step in NetCDF
+      !! ** Method  :   Write in numrow when kt == nitrst in NetCDF
       !!      file, save fields which are necessary for restart
-      !!
-      !! History :
-      !!        !  99-11  (M. Imbard)  Original code
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form
-      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE ioipsl
-
-      !! * Arguments 
-      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
-
-      !! * Local declarations
-      LOGICAL ::   llbon
-      CHARACTER (len=50) ::   clname, cln
-      INTEGER ::   ic, jc, itime
-      INTEGER ::   inumwrs
-      REAL(wp) ::   zdate0
-      REAL(wp), DIMENSION( 1) ::   zfice, zfblk   ! used only in case of ice & bulk
-      REAL(wp), DIMENSION(10) ::   zinfo(10)
-      REAL(wp), DIMENSION(jpi,jpj) :: ztab 
-#if defined key_agrif
-       Integer :: knum
-#endif
-      !!----------------------------------------------------------------------
-
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'rst_wri : write restart.output NetCDF file'
-         IF(lwp) WRITE(numout,*) '~~~~~~~'
-         zfice(1) = 1.e0   ;   zfblk(1) = 1.e0
-      ENDIF
-
-
-      IF( MOD( kt, nstock ) == 0 .OR. kt == nitend ) THEN
-         
-         ! 0. Initializations
-         ! ------------------
-
-         IF(lwp) WRITE(numout,*) ' '
-         IF(lwp) WRITE(numout,*) 'rst_write : write the restart file in NetCDF format ',   &
-                                              'at it= ',kt,' date= ',ndastp
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~'
-
-         ! Job informations 
-         zinfo(:) = 0.e0 
-         zinfo(1) = FLOAT( no        )   ! job number
-         zinfo(2) = FLOAT( kt        )   ! time-step
-         zinfo(3) = FLOAT( 2 - nsolv )   ! pcg solver
-         zinfo(4) = FLOAT( nsolv - 1 )   ! sor solver
-         IF( lk_zdftke ) THEN
-            zinfo(5) = 1.e0              ! TKE 
-         ELSE
-            zinfo(5) = 0.e0              ! no TKE 
-         ENDIF
-         zinfo(6) = FLOAT( ndastp )      ! date
-         zinfo(7) = adatrj               ! ???
-
-         ! delete the restart file if it exists 
-         INQUIRE( FILE=crestart, EXIST=llbon )
-         IF(llbon) THEN
-#if defined key_agrif
-       knum =Agrif_Get_Unit()
-            OPEN( UNIT=knum, FILE=crestart, STATUS='old' )
-            CLOSE( knum, STATUS='delete' )
-#else            
-            OPEN( UNIT=inumwrs, FILE=crestart, STATUS='old' )
-            CLOSE( inumwrs, STATUS='delete' )
-#endif
-         ENDIF
-
-         ! Name of the new restart file
-         ic     = 1
-         DO jc = 1, 16
-            IF( cexper(jc:jc) /= ' ' )   ic = jc
-         END DO
-         WRITE(cln,'("_",i4.4,i2.2,i2.2,"_restart")') nyear, nmonth, nday
-         clname = cexper(1:ic)//cln
-         ic = 1
-         DO jc = 1, 48
-            IF( clname(jc:jc) /= ' ' ) ic = jc
-         END DO
-         crestart = clname(1:ic)//".nc"
-         itime = 0
-         CALL ymds2ju( nyear, nmonth, nday, 0.e0, zdate0 )
-         CALL restini( 'NONE', jpi, jpj, glamt, gphit, jpk, gdept_0, clname,   &
-                        itime, zdate0, rdt*nstock ,inumwrs, domain_id=nidom )
-
-         CALL restput( inumwrs, 'info'   , 1  , 1  , 10 , 0, zinfo   )   ! restart informations
-         
-         CALL restput( inumwrs, 'ub'     , jpi, jpj, jpk, 0, ub      )   ! prognostic variables
-         CALL restput( inumwrs, 'vb'     , jpi, jpj, jpk, 0, vb      )
-         CALL restput( inumwrs, 'tb'     , jpi, jpj, jpk, 0, tb      )
-         CALL restput( inumwrs, 'sb'     , jpi, jpj, jpk, 0, sb      )
-         CALL restput( inumwrs, 'rotb'   , jpi, jpj, jpk, 0, rotb    )
-         CALL restput( inumwrs, 'hdivb'  , jpi, jpj, jpk, 0, hdivb   )
-         CALL restput( inumwrs, 'un'     , jpi, jpj, jpk, 0, un      )
-         CALL restput( inumwrs, 'vn'     , jpi, jpj, jpk, 0, vn      )
-         CALL restput( inumwrs, 'tn'     , jpi, jpj, jpk, 0, tn      )
-         CALL restput( inumwrs, 'sn'     , jpi, jpj, jpk, 0, sn      )
-         CALL restput( inumwrs, 'rotn'   , jpi, jpj, jpk, 0, rotn    )
-         CALL restput( inumwrs, 'hdivn'  , jpi, jpj, jpk, 0, hdivn   )
-
-         ztab(:,:) = gcx(1:jpi,1:jpj)
-         CALL restput( inumwrs, 'gcx'    , jpi, jpj, 1  , 0, ztab    )   ! Read elliptic solver arrays
-         ztab(:,:) = gcxb(1:jpi,1:jpj)
-         CALL restput( inumwrs, 'gcxb'   , jpi, jpj, 1  , 0, ztab    )
-# if defined key_dynspg_rl
-         CALL restput( inumwrs, 'bsfb'   , jpi, jpj, 1  , 0, bsfb    )   ! Rigid-lid formulation (bsf)
-         CALL restput( inumwrs, 'bsfn'   , jpi, jpj, 1  , 0, bsfn    )
-         CALL restput( inumwrs, 'bsfd'   , jpi, jpj, 1  , 0, bsfd    )
-# else
-         CALL restput( inumwrs, 'sshb'   , jpi, jpj, 1  , 0, sshb    )   ! free surface formulation (ssh)
-         CALL restput( inumwrs, 'sshn'   , jpi, jpj, 1  , 0, sshn    )
-#  if defined key_dynspg_ts
-         CALL restput( inumwrs, 'sshb_b' , jpi, jpj, 1  , 0, sshb_b  )   ! free surface formulation (ssh)
-         CALL restput( inumwrs, 'sshn_b' , jpi, jpj, 1  , 0, sshn_b  )   ! issued from barotropic loop
-         CALL restput( inumwrs, 'un_b'   , jpi, jpj, 1  , 0, un_b    )   ! horizontal transports
-         CALL restput( inumwrs, 'vn_b'   , jpi, jpj, 1  , 0, vn_b    )   ! issued from barotropic loop
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step
+      !!----------------------------------------------------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'rst_write : write ocean NetCDF restart file  kt =', kt,' date= ', ndastp
+         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]
+
+      ! prognostic variables
+      CALL iom_rstput( kt, nitrst, numrow, 'ub'     , ub      )   
+      CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vb      )
+      CALL iom_rstput( kt, nitrst, numrow, 'tb'     , tb      )
+      CALL iom_rstput( kt, nitrst, numrow, 'sb'     , sb      )
+      CALL iom_rstput( kt, nitrst, numrow, 'rotb'   , rotb    )
+      CALL iom_rstput( kt, nitrst, numrow, 'hdivb'  , hdivb   )
+      CALL iom_rstput( kt, nitrst, numrow, 'un'     , un      )
+      CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vn      )
+      CALL iom_rstput( kt, nitrst, numrow, 'tn'     , tn      )
+      CALL iom_rstput( kt, nitrst, numrow, 'sn'     , sn      )
+      CALL iom_rstput( kt, nitrst, numrow, 'rotn'   , rotn    )
+      CALL iom_rstput( kt, nitrst, numrow, 'hdivn'  , hdivn   )
+
+# if defined key_ice_lim        
+      CALL iom_rstput( kt, nitrst, numrow, 'nfice'  , REAL( nfice, wp) )   !  ice computation frequency
+      CALL iom_rstput( kt, nitrst, numrow, 'sst_io' , sst_io  )
+      CALL iom_rstput( kt, nitrst, numrow, 'sss_io' , sss_io  )
+      CALL iom_rstput( kt, nitrst, numrow, 'u_io'   , u_io    )
+      CALL iom_rstput( kt, nitrst, numrow, 'v_io'   , v_io    )
+#  if defined key_coupled
+      CALL iom_rstput( kt, nitrst, numrow, 'alb_ice', alb_ice )
 #  endif
 # endif
-# if defined key_zdftke   ||   defined key_esopa
-         IF( lk_zdftke ) THEN
-            CALL restput( inumwrs, 'en'     , jpi, jpj, jpk, 0, en      )   ! TKE arrays
-         ENDIF
-# endif
-# if defined key_ice_lim
-         zfice(1) = FLOAT( nfice )                                      ! Louvain La Neuve Sea Ice Model
-         CALL restput( inumwrs, 'nfice'  ,   1,   1, 1  , 0, zfice   )
-         CALL restput( inumwrs, 'sst_io' , jpi, jpj, 1  , 0, sst_io  )
-         CALL restput( inumwrs, 'sss_io' , jpi, jpj, 1  , 0, sss_io  )
-         CALL restput( inumwrs, 'u_io'   , jpi, jpj, 1  , 0, u_io    )
-         CALL restput( inumwrs, 'v_io'   , jpi, jpj, 1  , 0, v_io    )
-# if defined key_coupled
-         CALL restput( inumwrs, 'alb_ice', jpi, jpj, 1  , 0, alb_ice )
-# endif
-# endif
 # if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
-         zfblk(1) = FLOAT( nfbulk )                                 ! Bulk
-         CALL restput( inumwrs, 'nfbulk' ,   1,   1, 1  , 0, zfblk   )
-         CALL restput( inumwrs, 'gsst'   , jpi, jpj, 1  , 0, gsst    )
-# endif
-
-         CALL restclo( inumwrs )                                         ! close the restart file
-         
-      ENDIF
-
+      CALL iom_rstput( kt, nitrst, numrow, 'nfbulk' , REAL( nfbulk, wp) )   !  bulk computation frequency
+      CALL iom_rstput( kt, nitrst, numrow, 'gsst'   , gsst    )
+# endif
+
+      IF( kt == nitrst ) THEN
+         CALL iom_close( numrow )     ! close the restart file (only at last time step)
+         lrst_oce = .FALSE.
+      ENDIF
+      !
    END SUBROUTINE rst_write
 
@@ -246,154 +194,126 @@
       !!       nrstdt = 2  the duration of the experiment in days (adatrj)
       !!                    has been stored in the restart file.
-      !!
-      !! History :
-      !!        !  99-05  (M. Imbard)  Original code
-      !!   8.5  !  02-09  (G. Madec)  F90: Free form
-      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE iom
-
-      !! * Local declarations
-      INTEGER  ::   &
-         inum                 ! temporary logical unit
-      REAL(wp), DIMENSION(1, 1, 10)  ::   zinfo
-      REAL(wp), DIMENSION(1, 1, 1)   ::   zzz 
-      INTEGER  ::   ios
-#   if defined key_ice_lim
+      !!----------------------------------------------------------------------
+      REAL(wp) ::   zcoef, zkt, zndastp, znfice, znfbulk
+# if defined key_ice_lim
       INTEGER  ::   ji, jj
-#   endif
-      !!----------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'rst_read : read the NetCDF restart file'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~'
-
-      IF(lwp) WRITE(numout,*) ' Info on the present job : '
-      IF(lwp) WRITE(numout,*) '   job number          : ', no
-      IF(lwp) WRITE(numout,*) '   time-step           : ', nit000
-      IF(lwp) WRITE(numout,*) '   solver type         : ', nsolv
-      IF( lk_zdftke ) THEN
-         IF(lwp) WRITE(numout,*) '   tke option          : 1 '
-      ELSE
-         IF(lwp) WRITE(numout,*) '   tke option          : 0 '
-      ENDIF
-      IF(lwp) WRITE(numout,*) '   date ndastp         : ', ndastp
-      IF(lwp) WRITE(numout,*)
-
-      ! Time domain : restart
-      ! -------------------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' *** restart option'
-      SELECT CASE ( nrstdt )
-      CASE ( 0 ) 
-         IF(lwp) WRITE(numout,*) ' nrstdt = 0 no control of nit000'
-      CASE ( 1 ) 
-         IF(lwp) WRITE(numout,*) ' nrstdt = 1 we control the date of nit000'
-      CASE ( 2 )
-         IF(lwp) WRITE(numout,*) ' nrstdt = 2 the date adatrj is read in restart file'
-      CASE DEFAULT
-         IF(lwp) WRITE(numout,*) '  ===>>>> nrstdt not equal 0, 1 or 2 : no control of the date'
-         IF(lwp) WRITE(numout,*) ' =======                   ========='
-      END SELECT
-
-      CALL iom_open ( 'restart', inum )
-      
-      CALL iom_get ( inum, jpdom_unknown, 'info', zinfo )
-      
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' Info on the restart file read : '
-      IF(lwp) WRITE(numout,*) '   job number          : ', NINT( zinfo(1, 1, 1) )
-      IF(lwp) WRITE(numout,*) '   time-step           : ', NINT( zinfo(1, 1, 2) )
-      IF(lwp) WRITE(numout,*) '   solver type         : ', NINT( zinfo(1, 1, 4) ) + 1
-      IF(lwp) WRITE(numout,*) '   tke option          : ', NINT( zinfo(1, 1, 5) )
-      IF(lwp) WRITE(numout,*) '   date ndastp         : ', NINT( zinfo(1, 1, 6) )
-      IF(lwp) WRITE(numout,*)
-
+# endif
+      !!----------------------------------------------------------------------
+
+      IF(lwp) THEN                                             ! Contol prints
+         WRITE(numout,*)
+         WRITE(numout,*) 'rst_read : read oce NetCDF restart file'
+         WRITE(numout,*) '~~~~~~~~'
+         
+         WRITE(numout,*) ' *** Info on the present job : '
+         WRITE(numout,*) '   time-step           : ', nit000
+!!$         WRITE(numout,*) '   solver type         : ', nsolv
+!!$         IF( lk_zdftke ) THEN
+!!$            WRITE(numout,*) '   tke option          : 1 '
+!!$         ELSE
+!!$            WRITE(numout,*) '   tke option          : 0 '
+!!$         ENDIF
+         WRITE(numout,*) '   date ndastp         : ', ndastp
+         WRITE(numout,*)
+         WRITE(numout,*) ' *** restart option'
+         SELECT CASE ( nrstdt )
+         CASE ( 0 ) 
+            WRITE(numout,*) ' nrstdt = 0 no control of nit000'
+         CASE ( 1 ) 
+            WRITE(numout,*) ' nrstdt = 1 we control the date of nit000'
+         CASE ( 2 )
+            WRITE(numout,*) ' nrstdt = 2 the date adatrj is read in restart file'
+         CASE DEFAULT
+            WRITE(numout,*) '  ===>>>> nrstdt not equal 0, 1 or 2 : no control of the date'
+            WRITE(numout,*) '  =======                  ========='
+         END SELECT
+         WRITE(numout,*)
+      ENDIF
+
+      CALL iom_open( 'restart', numror )                       ! Open
+
+      ! Calendar informations
+      CALL iom_get( numror, 'kt'    , zkt     )   ! time-step 
+      CALL iom_get( numror, 'ndastp', zndastp )   ! date
+      ! Additional contol prints
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' *** Info on the restart file read : '
+         WRITE(numout,*) '   time-step           : ', NINT( zkt )
+!!$         WRITE(numout,*) '   solver type         : ', +++
+!!$         WRITE(numout,*) '   tke option          : ', +++
+         WRITE(numout,*) '   date ndastp         : ', NINT( zndastp )
+         WRITE(numout,*)
+      ENDIF
       ! Control of date
-      IF( nit000 - NINT( zinfo(1, 1, 2) )  /= 1 .AND. nrstdt /= 0 ) &
+      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)' )
-
       ! re-initialisation of  adatrj0
-      adatrj0 =  ( FLOAT( nit000-1 ) * rdttra(1) ) / rday
-
+      adatrj0 = ( REAL( nit000-1, wp ) * rdttra(1) ) / rday
       IF ( nrstdt == 2 ) THEN
 !                             by default ndatsp has been set to ndate0 in dom_nam
 !                             ndate0 has been read in the namelist (standard OPA 8)
 !                             here when nrstdt=2 we keep the  final date of previous run
-        ndastp = NINT( zinfo(1, 1, 6) )
-        adatrj0 =  zinfo(1, 1, 7)
-      ENDIF
-
-      CALL iom_get( inum, jpdom_local, 'ub'   , ub    )   ! Read prognostic variables
-      CALL iom_get( inum, jpdom_local, 'vb'   , vb    )
-      CALL iom_get( inum, jpdom_local, 'tb'   , tb    )
-      CALL iom_get( inum, jpdom_local, 'sb'   , sb    )
-      CALL iom_get( inum, jpdom_local, 'rotb' , rotb  )
-      CALL iom_get( inum, jpdom_local, 'hdivb', hdivb )
-      CALL iom_get( inum, jpdom_local, 'un'   , un    )
-      CALL iom_get( inum, jpdom_local, 'vn'   , vn    )
-      CALL iom_get( inum, jpdom_local, 'tn'   , tn    )
-      CALL iom_get( inum, jpdom_local, 'sn'   , sn    )
-      CALL iom_get( inum, jpdom_local, 'rotn' , rotn  )
-      CALL iom_get( inum, jpdom_local, 'hdivn', hdivn )
-! Caution : extrahallow 
-! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
-      CALL iom_get( inum, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) )
-      CALL iom_get( inum, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) )     ! Read elliptic solver arrays
-# if defined key_dynspg_rl
-      CALL iom_get( inum, jpdom_local, 'bsfb', bsfb )     ! Rigid-lid formulation (bsf)
-      CALL iom_get( inum, jpdom_local, 'bsfn', bsfn )
-      CALL iom_get( inum, jpdom_local, 'bsfd', bsfd )
-# else
-      CALL iom_get( inum, jpdom_local, 'sshb', sshb )     ! free surface formulation (ssh)
-      CALL iom_get( inum, jpdom_local, 'sshn', sshn )
-#  if defined key_dynspg_ts
-      CALL iom_get( inum, jpdom_local, 'sshb_b', sshb_b ) ! free surface formulation (ssh)
-      CALL iom_get( inum, jpdom_local, 'sshn_b', sshn_b ) ! issued from barotropic loop
-      CALL iom_get( inum, jpdom_local, 'un_b'  , un_b )   ! horizontal transports
-      CALL iom_get( inum, jpdom_local, 'vn_b'  , vn_b )   ! issued from barotropic loop
-#  endif
-# endif
-# if defined key_zdftke   ||   defined key_esopa
-      IF( lk_zdftke ) THEN
-         IF( NINT( zinfo(1, 1, 5) ) == 1 ) THEN                                ! Read tke arrays
-            CALL iom_get( inum, jpdom_local, 'en', en )
-            ln_rstke = .FALSE.
-         ELSE
-            IF(lwp) WRITE(numout,*) ' ===>>>> : the previous restart file did not used  tke scheme'
-            IF(lwp) WRITE(numout,*) ' =======                ======='
-            nrstdt = 2
-            ln_rstke = .TRUE.
-         ENDIF
-      ENDIF
-# endif
+         ndastp = NINT( zndastp )
+        CALL iom_get( numror, 'adatrj', adatrj )   ! number of elapsed days since the begining of last run
+      ENDIF
+
+      !                                                       ! Read prognostic variables
+      CALL iom_get( numror, jpdom_local, 'ub'   , ub    )        ! before i-component velocity
+      CALL iom_get( numror, jpdom_local, 'vb'   , vb    )        ! before j-component velocity
+      CALL iom_get( numror, jpdom_local, 'tb'   , tb    )        ! before temperature
+      CALL iom_get( numror, jpdom_local, 'sb'   , sb    )        ! before salinity
+      CALL iom_get( numror, jpdom_local, 'rotb' , rotb  )        ! before curl
+      CALL iom_get( numror, jpdom_local, 'hdivb', hdivb )        ! before horizontal divergence
+      CALL iom_get( numror, jpdom_local, 'un'   , un    )        ! now    i-component velocity
+      CALL iom_get( numror, jpdom_local, 'vn'   , vn    )        ! now    j-component velocity
+      CALL iom_get( numror, jpdom_local, 'tn'   , tn    )        ! now    temperature
+      CALL iom_get( numror, jpdom_local, 'sn'   , sn    )        ! now    salinity
+      CALL iom_get( numror, jpdom_local, 'rotn' , rotn  )        ! now    curl
+      CALL iom_get( numror, jpdom_local, 'hdivn', hdivn )        ! now    horizontal divergence
+
+
+      IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0)
+         tb   (:,:,:) = tn   (:,:,:)                             ! all before fields set to now field values
+         sb   (:,:,:) = sn   (:,:,:)
+         ub   (:,:,:) = un   (:,:,:)
+         vb   (:,:,:) = vn   (:,:,:)
+         rotb (:,:,:) = rotn (:,:,:)
+         hdivb(:,:,:) = hdivn(:,:,:)
+      ENDIF
+
+      !!sm: TO BE MOVED IN NEW SURFACE MODULE...
+
 # if defined key_ice_lim
       ! Louvain La Neuve Sea Ice Model
-      ios = iom_varid( inum, 'nfice' )
-      IF( ios > 0 ) then 
-         CALL iom_get( inum, jpdom_unknown, 'nfice' , zzz )
-         zinfo(1, 1, 8) = zzz(1, 1, 1)
-         CALL iom_get( inum, jpdom_local, 'sst_io', sst_io )
-         CALL iom_get( inum, jpdom_local, 'sss_io', sss_io )
-         CALL iom_get( inum, jpdom_local, 'u_io'  , u_io )
-         CALL iom_get( inum, jpdom_local, 'v_io'  , v_io )
+      IF( iom_varid( numror, 'nfice' ) > 0 ) then 
+         CALL iom_get( numror             , 'nfice'  , znfice  )   ! ice computation frequency
+         CALL iom_get( numror, jpdom_local, 'sst_io' , sst_io  )
+         CALL iom_get( numror, jpdom_local, 'sss_io' , sss_io  )
+         CALL iom_get( numror, jpdom_local, 'u_io'   , u_io    )
+         CALL iom_get( numror, jpdom_local, 'v_io'   , v_io    )
 #if defined key_coupled
-         CALL iom_get( inum, jpdom_local, 'alb_ice', alb_ice )
+         CALL iom_get( numror, jpdom_local, 'alb_ice', alb_ice )
 #endif
-      ENDIF
-      IF( zinfo(1, 1, 8) /= FLOAT(nfice) .OR. ios == 0 ) THEN
+         IF( znfice /= REAL( nfice, wp ) ) THEN      ! if nfice changed between 2 runs
+            zcoef = REAL( nfice-1, wp ) / znfice
+            sst_io(:,:) = zcoef * sst_io(:,:)
+            sss_io(:,:) = zcoef * sss_io(:,:)
+            u_io  (:,:) = zcoef * u_io  (:,:)
+            v_io  (:,:) = zcoef * v_io  (:,:)
+         ENDIF
+      ELSE
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'rst_read :  LLN sea Ice Model => Ice initialization'
          IF(lwp) WRITE(numout,*)
-         sst_io(:,:) = ( nfice-1 )*( tn(:,:,1) + rt0 )          !!bug a explanation is needed here!
-         sss_io(:,:) = ( nfice-1 )*  sn(:,:,1)
+         zcoef = REAL( nfice-1, wp )
+         sst_io(:,:) = zcoef *( tn(:,:,1) + rt0 )          !!bug a explanation is needed here!
+         sss_io(:,:) = zcoef *  sn(:,:,1)
+         zcoef = 0.5 * REAL( nfice-1, wp )
          DO jj = 2, jpj
-            DO ji = 2, jpi
-               u_io(ji,jj) = ( nfice-1 ) * 0.5 * ( un(ji-1,jj  ,1) + un(ji-1,jj-1,1) )
-               v_io(ji,jj) = ( nfice-1 ) * 0.5 * ( vn(ji  ,jj-1,1) + vn(ji-1,jj-1,1) )
+            DO ji = fs_2, jpi   ! vector opt.
+               u_io(ji,jj) = zcoef * ( un(ji-1,jj  ,1) + un(ji-1,jj-1,1) )
+               v_io(ji,jj) = zcoef * ( vn(ji  ,jj-1,1) + vn(ji-1,jj-1,1) )
             END DO
          END DO
@@ -405,41 +325,25 @@
 # if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
       ! Louvain La Neuve Sea Ice Model
-      ios = iom_varid( inum, 'nfbulk' )
-      IF( ios > 0 ) then 
-         CALL iom_get( inum, jpdom_unknown, 'nfbulk' , zzz )
-         CALL iom_get( inum, jpdom_local, 'gsst' , gsst )
-         zinfo(1, 1, 9) = zzz(1, 1, 1)
-      ENDIF
-      IF( zinfo(1, 1, 9) /= FLOAT(nfbulk) .OR. ios == 0 ) THEN
+      IF( iom_varid( numror, 'nfbulk' ) > 0 ) THEN 
+         CALL iom_get( numror             , 'nfbulk', znfbulk )   ! bulk computation frequency
+         CALL iom_get( numror, jpdom_local, 'gsst'  , gsst    )
+         IF( znfbulk /= REAL(nfbulk, wp) ) THEN      ! if you change nfbulk between 2 runs
+            zcoef = REAL( nfbulk-1, wp ) / znfbulk
+            gsst(:,:) = zcoef * gsst(:,:)
+         ENDIF
+      ELSE
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'rst_read :  LLN sea Ice Model => Ice initialization'
          IF(lwp) WRITE(numout,*)
-         gsst(:,:) = 0.
-         gsst(:,:) = gsst(:,:) + ( nfbulk-1 )*( tn(:,:,1) + rt0 )
+         gsst(:,:) = REAL( nfbulk - 1, wp )*( tn(:,:,1) + rt0 )
       ENDIF
 # endif
       
-      CALL iom_close( inum )
-
-  ! In case of restart with neuler = 0 then put all before fields = to now fields
-    IF ( neuler == 0 ) THEN
-      tb(:,:,:)=tn(:,:,:)
-      sb(:,:,:)=sn(:,:,:)
-      ub(:,:,:)=un(:,:,:)
-      vb(:,:,:)=vn(:,:,:)
-      rotb(:,:,:)=rotn(:,:,:)
-      hdivb(:,:,:)=hdivn(:,:,:)
-#if defined key_dynspg_rl
-    ! rigid lid
-      bsfb(:,:)=bsfn(:,:)
-#else
-    ! free surface formulation (eta)
-      sshb(:,:)=sshn(:,:)
+      !!sm: end of TO BE MOVED IN NEW SURFACE MODULE...
+      !
+   END SUBROUTINE rst_read
+
 #endif
-    ENDIF
-
-   END SUBROUTINE rst_read
-
-#endif
+
    !!=====================================================================
 END MODULE restart
Index: trunk/NEMO/OPA_SRC/step.F90
===================================================================
--- trunk/NEMO/OPA_SRC/step.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/step.F90	(revision 508)
@@ -4,23 +4,23 @@
    !! Time-stepping    : manager of the ocean, tracer and ice time stepping
    !!======================================================================
-   !! History :
-   !!        !  91-03  ()  Original code
-   !!        !  91-11  (G. Madec)
-   !!        !  92-06  (M. Imbard)  add a first output record
-   !!        !  96-04  (G. Madec)  introduction of dynspg
-   !!        !  96-04  (M.A. Foujols)  introduction of passive tracer
-   !!   8.0  !  97-06  (G. Madec)  new architecture of call
-   !!   8.2  !  97-06  (G. Madec, M. Imbard, G. Roullet)  free surface
-   !!   8.2  !  99-02  (G. Madec, N. Grima)  hpg implicit
-   !!   8.2  !  00-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
-   !!   9.0  !  02-06  (G. Madec)  free form, suppress macro-tasking
-   !!    "   !  04-08  (C. Talandier) New trends organization
-   !!    "   !  05-01  (C. Ethe) Add the KPP closure scheme
-   !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
-   !!    "   !  05-11  (G. Madec)  Reorganisation of tra and dyn calls
+   !! History :        !  91-03  ()  Original code
+   !!                  !  91-11  (G. Madec)
+   !!                  !  92-06  (M. Imbard)  add a first output record
+   !!                  !  96-04  (G. Madec)  introduction of dynspg
+   !!                  !  96-04  (M.A. Foujols)  introduction of passive tracer
+   !!             8.0  !  97-06  (G. Madec)  new architecture of call
+   !!             8.2  !  97-06  (G. Madec, M. Imbard, G. Roullet)  free surface
+   !!             8.2  !  99-02  (G. Madec, N. Grima)  hpg implicit
+   !!             8.2  !  00-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
+   !!             9.0  !  02-06  (G. Madec)  free form, suppress macro-tasking
+   !!             " "  !  04-08  (C. Talandier) New trends organization
+   !!             " "  !  05-01  (C. Ethe) Add the KPP closure scheme
+   !!             " "  !  05-11  (V. Garnier) Surface pressure gradient organization
+   !!             " "  !  05-11  (G. Madec)  Reorganisation of tra and dyn calls
+   !!             " "  !  06-01  (L. Debreu, C. Mazauric)  Agrif implementation
+   !!             " "  !  06-07  (S. Masson)  restart using iom
    !!----------------------------------------------------------------------
    !!   stp            : OPA system time-stepping
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain variables 
@@ -30,4 +30,5 @@
    USE cpl_oce         ! coupled ocean-atmosphere variables
    USE in_out_manager  ! I/O manager
+   USE iom
    USE lbclnk
 
@@ -126,5 +127,4 @@
    PRIVATE
 
-   !! * Routine accessibility
    PUBLIC stp            ! called by opa.F90
 
@@ -135,13 +135,13 @@
    !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    !!----------------------------------------------------------------------
 
 CONTAINS
 
-#if !defined key_agrif
+#if defined key_agrif
+   SUBROUTINE stp( )
+#else
    SUBROUTINE stp( kstp )
-#else
-   SUBROUTINE stp( )
 #endif
       !!----------------------------------------------------------------------
@@ -160,11 +160,10 @@
       !!              -7- Compute the diagnostics variables (rd,N2, div,cur,w)
       !!              -8- Outputs and diagnostics
-      !!
       !!----------------------------------------------------------------------
       !! * Arguments
-#if !defined key_agrif   
+#if defined key_agrif   
+      INTEGER               :: kstp   ! ocean time-step index
+#else
       INTEGER, INTENT( in ) :: kstp   ! ocean time-step index
-#else
-      INTEGER               :: kstp   ! ocean time-step index
 #endif      
 
@@ -176,10 +175,13 @@
       kstp = nit000 + Agrif_Nb_Step()
 !      IF ( Agrif_Root() .and. lwp) Write(*,*) '---'
-!      IF (lwp) Write(*,*) 'Grid N°',Agrif_Fixed(),' time step ',kstp
+!      IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp
 #endif   
       indic = 1                    ! reset to no error condition
+
       adatrj = adatrj + rdt/86400._wp
 
       CALL day( kstp )             ! Calendar
+
+      CALL rst_opn( kstp )         ! Open the restart file
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -326,5 +328,5 @@
                              CALL tra_ldf    ( kstp )       ! lateral mixing
 #if defined key_agrif
-      IF (.NOT. Agrif_Root())  CALL Agrif_Sponge_tra( kstp )          ! tracers sponge
+      IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra( kstp )  ! tracers sponge
 #endif
                              CALL tra_zdf    ( kstp )       ! vertical mixing
@@ -363,5 +365,5 @@
                                CALL dyn_ldf( kstp )           ! lateral mixing
 #if defined key_agrif
-      IF (.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn( kstp )         ! momemtum sponge
+      IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn( kstp )  ! momemtum sponge
 #endif
                                CALL dyn_hpg( kstp )           ! horizontal gradient of Hydrostatic pressure
@@ -400,15 +402,25 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Control, diagnostics and outputs
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
-      !-----------------------------------------------------------------------
-
-      !                                            ! Time loop: control and print
-                       CALL stp_ctl( kstp, indic )
-                       IF ( indic < 0 ) CALL ctl_stop( 'step: indic < 0' ) 
-
-      IF ( nstop == 0 ) THEN
-         !                                         ! Diagnostics:
+      ! Control, and restarts
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
+      !-----------------------------------------------------------------------
+      !                                                     ! Time loop: control and print
+                                 CALL stp_ctl( kstp, indic )
+      IF( indic < 0          )   CALL ctl_stop( 'step: indic < 0' )
+
+      IF( kstp == nit000     )   CALL iom_close( numror )             ! close input  ocean restart file
+      IF( lrst_oce           )   CALL rst_write  ( kstp )             ! write output ocean restart file
+      IF( lk_obc             )   CALL obc_rst_wri( kstp )             ! write open boundary restart file
+      IF( lk_trdmld          )   CALL trd_mld_rst_write( kstp )       ! ocean model: restart file output for trends diagnostics
+
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! diagnostics and outputs
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
+      !-----------------------------------------------------------------------
+
+      IF ( nstop == 0 ) THEN                                ! Diagnostics
          IF( lk_floats  )   CALL flo_stp( kstp )                 ! drifting Floats
          IF( lk_trddyn  )   CALL trd_dwr( kstp )                 ! trends: dynamics 
@@ -423,10 +435,6 @@
          IF( ln_diaptr  )   CALL dia_ptr( kstp )                 ! Poleward TRansports diagnostics
 
-         !                                         ! save and outputs
-                            CALL rst_write  ( kstp )             ! ocean model: restart file output
-         IF( lk_obc     )   CALL obc_rst_wri( kstp )             ! ocean model: open boundary restart file output
-         IF( lk_trdmld  )   CALL trd_mld_rst_write( kstp )       ! ocean model: restart file output for trends diagnostics
+         !                                                 ! Outputs
                             CALL dia_wri    ( kstp, indic )      ! ocean model: outputs
-
       ENDIF
 
@@ -436,5 +444,5 @@
 
       IF( lk_cpl    )   CALL cpl_stp( kstp )                 ! coupled mode : field exchanges
-
+      !
    END SUBROUTINE stp
 
