Index: /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90
===================================================================
--- /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90	(revision 8816)
+++ /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90	(revision 8817)
@@ -12,5 +12,8 @@
    !!   'key_lim3'                                       ESIM sea-ice model
    !!----------------------------------------------------------------------
-   !!   ice_dyn_adv_pra   : advection of sea ice using Prather scheme
+   !!   ice_dyn_adv_pra : advection of sea ice using Prather scheme
+   !!   adv_x, adv_y    : Prather scheme applied in i- and j-direction, resp.
+   !!   adv_pra_init    : initialisation of the Prather scheme
+   !!   adv_pra_rst     : read/write Prather field in ice restart file, or initialized to zero
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean domain
@@ -32,5 +35,4 @@
 
    ! Moments for advection
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   sxopw, syopw, sxxopw, syyopw, sxyopw   ! open water in sea ice
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxice, syice, sxxice, syyice, sxyice   ! ice thickness 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsn , sysn , sxxsn , syysn , sxysn    ! snow thickness
@@ -39,8 +41,9 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   ! ice salinity
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   ! ice age
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   sxopw, syopw, sxxopw, syyopw, sxyopw   ! open water in sea ice
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     ! ice layers heat content
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     ! ice layers heat content
-   !
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -262,4 +265,5 @@
    END SUBROUTINE ice_dyn_adv_pra
    
+   
    SUBROUTINE adv_x( pdf, put , pcrh, psm , ps0 ,   &
       &              psx, psxx, psy , psyy, psxy )
@@ -608,4 +612,5 @@
    END SUBROUTINE adv_y
 
+
    SUBROUTINE adv_pra_init
       !!-------------------------------------------------------------------
@@ -616,4 +621,6 @@
       INTEGER ::   ierr
       !!-------------------------------------------------------------------
+      !
+      !                             !* allocate prather fields
       ALLOCATE( sxopw(jpi,jpj)     , syopw(jpi,jpj)     , sxxopw(jpi,jpj)     , syyopw(jpi,jpj)     , sxyopw(jpi,jpj)     ,   &
          &      sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
@@ -632,7 +639,8 @@
       IF( ierr /= 0 )   CALL ctl_stop('STOP', 'adv_pra_init : unable to allocate ice arrays for Prather advection scheme')
       !
-      CALL adv_pra_rst( 'READ' )  !* read or initialize all required files
+      CALL adv_pra_rst( 'READ' )    !* read or initialize all required files
       !
    END SUBROUTINE adv_pra_init
+
 
    SUBROUTINE adv_pra_rst( cdrw, kt )
@@ -652,376 +660,175 @@
       CHARACTER(len=25) ::   znam
       CHARACTER(len=2)  ::   zchar, zchar1
-      REAL(wp), DIMENSION(jpi,jpj) :: z2d
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z3d   ! 3D workspace
       !!----------------------------------------------------------------------
       !
-      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
-         !                                   ! ---------------
-         IF( ln_rstart ) THEN                   !* Read the restart file
-            !
-            id1 = iom_varid( numrir, 'sxopw' , ldstop = .FALSE. )
-            !
-            IF( id1 > 0 ) THEN      ! fields exist
-               DO jl = 1, jpl 
-                  WRITE(zchar,'(I2.2)') jl
-                  znam = 'sxice'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxice(:,:,jl) = z2d(:,:)
-                  znam = 'syice'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  syice(:,:,jl) = z2d(:,:)
-                  znam = 'sxxice'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxxice(:,:,jl) = z2d(:,:)
-                  znam = 'syyice'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  syyice(:,:,jl) = z2d(:,:)
-                  znam = 'sxyice'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxyice(:,:,jl) = z2d(:,:)
-                  znam = 'sxsn'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxsn(:,:,jl) = z2d(:,:)
-                  znam = 'sysn'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sysn(:,:,jl) = z2d(:,:)
-                  znam = 'sxxsn'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxxsn(:,:,jl) = z2d(:,:)
-                  znam = 'syysn'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  syysn(:,:,jl) = z2d(:,:)
-                  znam = 'sxysn'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxysn(:,:,jl) = z2d(:,:)
-                  znam = 'sxa'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxa(:,:,jl) = z2d(:,:)
-                  znam = 'sya'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sya(:,:,jl) = z2d(:,:)
-                  znam = 'sxxa'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxxa(:,:,jl) = z2d(:,:)
-                  znam = 'syya'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  syya(:,:,jl) = z2d(:,:)
-                  znam = 'sxya'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxya(:,:,jl) = z2d(:,:)
-                  znam = 'sxc0'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxc0(:,:,jl) = z2d(:,:)
-                  znam = 'syc0'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  syc0(:,:,jl) = z2d(:,:)
-                  znam = 'sxxc0'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxxc0(:,:,jl) = z2d(:,:)
-                  znam = 'syyc0'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  syyc0(:,:,jl) = z2d(:,:)
-                  znam = 'sxyc0'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxyc0(:,:,jl) = z2d(:,:)
-                  znam = 'sxsal'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxsal(:,:,jl) = z2d(:,:)
-                  znam = 'sysal'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sysal(:,:,jl) = z2d(:,:)
-                  znam = 'sxxsal'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxxsal(:,:,jl) = z2d(:,:)
-                  znam = 'syysal'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  syysal(:,:,jl) = z2d(:,:)
-                  znam = 'sxysal'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxysal(:,:,jl) = z2d(:,:)
-                  znam = 'sxage'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxage(:,:,jl) = z2d(:,:)
-                  znam = 'syage'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  syage(:,:,jl) = z2d(:,:)
-                  znam = 'sxxage'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxxage(:,:,jl) = z2d(:,:)
-                  znam = 'syyage'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  syyage(:,:,jl) = z2d(:,:)
-                  znam = 'sxyage'//'_htc'//zchar
-                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                  sxyage(:,:,jl)= z2d(:,:)
-               END DO
-               IF ( ln_pnd_H12 ) THEN
-                  DO jl = 1, jpl 
-                     WRITE(zchar,'(I2.2)') jl
-                     znam = 'sxap'//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     sxap(:,:,jl) = z2d(:,:)
-                     znam = 'syap'//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     syap(:,:,jl) = z2d(:,:)
-                     znam = 'sxxap'//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     sxxap(:,:,jl) = z2d(:,:)
-                     znam = 'syyap'//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     syyap(:,:,jl) = z2d(:,:)
-                     znam = 'sxyap'//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     sxyap(:,:,jl) = z2d(:,:)
-
-                     znam = 'sxvp'//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     sxvp(:,:,jl) = z2d(:,:)
-                     znam = 'syvp'//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     syvp(:,:,jl) = z2d(:,:)
-                     znam = 'sxxvp'//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     sxxvp(:,:,jl) = z2d(:,:)
-                     znam = 'syyvp'//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     syyvp(:,:,jl) = z2d(:,:)
-                     znam = 'sxyvp'//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     sxyvp(:,:,jl) = z2d(:,:)
-                  END DO
-               ENDIF
-
-               CALL iom_get( numrir, jpdom_autoglo, 'sxopw ' ,  sxopw  )
-               CALL iom_get( numrir, jpdom_autoglo, 'syopw ' ,  syopw  )
-               CALL iom_get( numrir, jpdom_autoglo, 'sxxopw' ,  sxxopw )
-               CALL iom_get( numrir, jpdom_autoglo, 'syyopw' ,  syyopw )
-               CALL iom_get( numrir, jpdom_autoglo, 'sxyopw' ,  sxyopw )
-
-               DO jl = 1, jpl 
-                  WRITE(zchar,'(I2.2)') jl
-                  DO jk = 1, nlay_i 
-                     WRITE(zchar1,'(I2.2)') jk
-                     znam = 'sxe'//'_il'//zchar1//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     sxe(:,:,jk,jl) = z2d(:,:)
-                     znam = 'sye'//'_il'//zchar1//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     sye(:,:,jk,jl) = z2d(:,:)
-                     znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     sxxe(:,:,jk,jl) = z2d(:,:)
-                     znam = 'syye'//'_il'//zchar1//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     syye(:,:,jk,jl) = z2d(:,:)
-                     znam = 'sxye'//'_il'//zchar1//'_htc'//zchar
-                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-                     sxye(:,:,jk,jl) = z2d(:,:)
-                  END DO
-               END DO
-               !
-            ELSE                                     ! start rheology from rest
-               IF(lwp) WRITE(numout,*) '   ==>>   previous run without Prather, set moments to 0'
-               sxopw (:,:) = 0._wp   ;   sxice (:,:,:) = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:) = 0._wp
-               syopw (:,:) = 0._wp   ;   syice (:,:,:) = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:) = 0._wp
-               sxxopw(:,:) = 0._wp   ;   sxxice(:,:,:) = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:) = 0._wp
-               syyopw(:,:) = 0._wp   ;   syyice(:,:,:) = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:) = 0._wp
-               sxyopw(:,:) = 0._wp   ;   sxyice(:,:,:) = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:) = 0._wp
-               !
-               sxc0  (:,:,:) = 0._wp   ;   sxe  (:,:,:,:) = 0._wp   ;   sxsal  (:,:,:) = 0._wp   ;   sxage  (:,:,:) = 0._wp
-               syc0  (:,:,:) = 0._wp   ;   sye  (:,:,:,:) = 0._wp   ;   sysal  (:,:,:) = 0._wp   ;   syage  (:,:,:) = 0._wp
-               sxxc0 (:,:,:) = 0._wp   ;   sxxe (:,:,:,:) = 0._wp   ;   sxxsal (:,:,:) = 0._wp   ;   sxxage (:,:,:) = 0._wp
-               syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp
-               sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp
-               IF ( ln_pnd_H12 ) THEN
-                  sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
-                  syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
-                  sxxap (:,:,:) = 0._wp    ; sxxvp (:,:,:) = 0._wp 
-                  syyap (:,:,:) = 0._wp    ; syyvp (:,:,:) = 0._wp 
-                  sxyap (:,:,:) = 0._wp    ; sxyvp (:,:,:) = 0._wp
-               ENDIF
+      !                                      !==========================!
+      IF( TRIM(cdrw) == 'READ' ) THEN        !==  Read or initialize  ==!
+         !                                   !==========================!
+         !
+         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxopw' , ldstop = .FALSE. )    ! file exist: id1>0
+         ELSE                   ;   id1 = 0                                                  ! no restart: id1=0
+         ENDIF
+         !
+         IF( id1 > 0 ) THEN                     !**  Read the restart file  **!
+            !
+            !                                                        ! ice thickness
+            CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice  )
+            CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice  )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice )
+            CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice )
+            !                                                        ! snow thickness
+            CALL iom_get( numrir, jpdom_autoglo, 'sxsn'  , sxsn   )
+            CALL iom_get( numrir, jpdom_autoglo, 'sysn'  , sysn   )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn  )
+            CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  )
+            !                                                        ! lead fraction
+            CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    )
+            CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxxa'  , sxxa   )
+            CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   )
+            !                                                        ! snow thermal content
+            CALL iom_get( numrir, jpdom_autoglo, 'sxc0'  , sxc0   )
+            CALL iom_get( numrir, jpdom_autoglo, 'syc0'  , syc0   )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxxc0' , sxxc0  )
+            CALL iom_get( numrir, jpdom_autoglo, 'syyc0' , syyc0  )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxyc0' , sxyc0  )
+            !                                                        ! ice salinity
+            CALL iom_get( numrir, jpdom_autoglo, 'sxsal' , sxsal  )
+            CALL iom_get( numrir, jpdom_autoglo, 'sysal' , sysal  )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxxsal', sxxsal )
+            CALL iom_get( numrir, jpdom_autoglo, 'syysal', syysal )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxysal', sxysal )
+            !                                                        ! ice age
+            CALL iom_get( numrir, jpdom_autoglo, 'sxage' , sxage  )
+            CALL iom_get( numrir, jpdom_autoglo, 'syage' , syage  )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxxage', sxxage )
+            CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage )
+            !                                                        ! open water in sea ice
+            CALL iom_get( numrir, jpdom_autoglo, 'sxopw ', sxopw  )
+            CALL iom_get( numrir, jpdom_autoglo, 'syopw ', syopw  )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxxopw', sxxopw )
+            CALL iom_get( numrir, jpdom_autoglo, 'syyopw', syyopw )
+            CALL iom_get( numrir, jpdom_autoglo, 'sxyopw', sxyopw )
+            !                                                        ! ice layers heat content
+            DO jk = 1, nlay_i 
+               WRITE(zchar1,'(I2.2)') jk
+               znam = 'sxe'//'_il'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
+               znam = 'sye'//'_il'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:)
+               znam = 'sxxe'//'_il'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
+               znam = 'syye'//'_il'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
+               znam = 'sxye'//'_il'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:)
+            END DO
+            !
+            IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction
+               CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  )
+               CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap  )
+               CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap )
+               CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap )
+               CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap )
+               !                                                     ! melt pond volume
+               CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp  )
+               CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp  )
+               CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp )
+               CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp )
+               CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp )
             ENDIF
-         ELSE                                   !* Start from rest
-            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set moments to 0'
-            sxopw (:,:) = 0._wp   ;   sxice (:,:,:) = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:) = 0._wp
-            syopw (:,:) = 0._wp   ;   syice (:,:,:) = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:) = 0._wp
-            sxxopw(:,:) = 0._wp   ;   sxxice(:,:,:) = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:) = 0._wp
-            syyopw(:,:) = 0._wp   ;   syyice(:,:,:) = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:) = 0._wp
-            sxyopw(:,:) = 0._wp   ;   sxyice(:,:,:) = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:) = 0._wp
-            !
-            sxc0  (:,:,:) = 0._wp   ;   sxe  (:,:,:,:) = 0._wp   ;   sxsal  (:,:,:) = 0._wp   ;   sxage  (:,:,:) = 0._wp
-            syc0  (:,:,:) = 0._wp   ;   sye  (:,:,:,:) = 0._wp   ;   sysal  (:,:,:) = 0._wp   ;   syage  (:,:,:) = 0._wp
-            sxxc0 (:,:,:) = 0._wp   ;   sxxe (:,:,:,:) = 0._wp   ;   sxxsal (:,:,:) = 0._wp   ;   sxxage (:,:,:) = 0._wp
-            syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp
-            sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp
-            IF ( ln_pnd_H12 ) THEN
-               sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
-               syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
-               sxxap (:,:,:) = 0._wp    ; sxxvp (:,:,:) = 0._wp 
-               syyap (:,:,:) = 0._wp    ; syyvp (:,:,:) = 0._wp 
-               sxyap (:,:,:) = 0._wp    ; sxyvp (:,:,:) = 0._wp
+            !
+         ELSE                                   !**  start rheology from rest  **!
+            !
+            IF(lwp) WRITE(numout,*) '   ==>>   start from rest OR previous run without Prather, set moments to 0'
+            !
+            sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness
+            sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness
+            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! lead fraction
+            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow thermal content
+            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
+            sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
+            sxopw = 0._wp   ;   syopw = 0._wp   ;   sxxopw = 0._wp   ;   syyopw = 0._wp   ;   sxyopw = 0._wp      ! open water in sea ice
+            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
+            IF( ln_pnd_H12 ) THEN
+               sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction
+               sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume
             ENDIF
          ENDIF
          !
-      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
-         !                                   ! -------------------
-         IF(lwp) WRITE(numout,*) '---- adv-rst ----'
+         !                                   !=====================================!
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   !==  write in the ice restart file  ==!
+         !                                   !=====================================!
+         IF(lwp) WRITE(numout,*) '----  ice-adv-rst  ----'
          iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
          !
-         DO jl = 1, jpl 
-            WRITE(zchar,'(I2.2)') jl
-            znam = 'sxice'//'_htc'//zchar
-            z2d(:,:) = sxice(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syice'//'_htc'//zchar
-            z2d(:,:) = syice(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxxice'//'_htc'//zchar
-            z2d(:,:) = sxxice(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syyice'//'_htc'//zchar
-            z2d(:,:) = syyice(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxyice'//'_htc'//zchar
-            z2d(:,:) = sxyice(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxsn'//'_htc'//zchar
-            z2d(:,:) = sxsn(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sysn'//'_htc'//zchar
-            z2d(:,:) = sysn(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxxsn'//'_htc'//zchar
-            z2d(:,:) = sxxsn(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syysn'//'_htc'//zchar
-            z2d(:,:) = syysn(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxysn'//'_htc'//zchar
-            z2d(:,:) = sxysn(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxa'//'_htc'//zchar
-            z2d(:,:) = sxa(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sya'//'_htc'//zchar
-            z2d(:,:) = sya(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxxa'//'_htc'//zchar
-            z2d(:,:) = sxxa(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syya'//'_htc'//zchar
-            z2d(:,:) = syya(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxya'//'_htc'//zchar
-            z2d(:,:) = sxya(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxc0'//'_htc'//zchar
-            z2d(:,:) = sxc0(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syc0'//'_htc'//zchar
-            z2d(:,:) = syc0(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxxc0'//'_htc'//zchar
-            z2d(:,:) = sxxc0(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syyc0'//'_htc'//zchar
-            z2d(:,:) = syyc0(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxyc0'//'_htc'//zchar
-            z2d(:,:) = sxyc0(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxsal'//'_htc'//zchar
-            z2d(:,:) = sxsal(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sysal'//'_htc'//zchar
-            z2d(:,:) = sysal(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxxsal'//'_htc'//zchar
-            z2d(:,:) = sxxsal(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syysal'//'_htc'//zchar
-            z2d(:,:) = syysal(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxysal'//'_htc'//zchar
-            z2d(:,:) = sxysal(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxage'//'_htc'//zchar
-            z2d(:,:) = sxage(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syage'//'_htc'//zchar
-            z2d(:,:) = syage(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxxage'//'_htc'//zchar
-            z2d(:,:) = sxxage(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syyage'//'_htc'//zchar
-            z2d(:,:) = syyage(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxyage'//'_htc'//zchar
-            z2d(:,:) = sxyage(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-         END DO
-
-         CALL iom_rstput( iter, nitrst, numriw, 'sxopw ' ,  sxopw  )
-         CALL iom_rstput( iter, nitrst, numriw, 'syopw ' ,  syopw  )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxxopw' ,  sxxopw )
-         CALL iom_rstput( iter, nitrst, numriw, 'syyopw' ,  syyopw )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxyopw' ,  sxyopw )
-         
-         DO jl = 1, jpl 
-            WRITE(zchar,'(I2.2)') jl
-            DO jk = 1, nlay_i 
-               WRITE(zchar1,'(I2.2)') jk
-               znam = 'sxe'//'_il'//zchar1//'_htc'//zchar
-               z2d(:,:) = sxe(:,:,jk,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sye'//'_il'//zchar1//'_htc'//zchar
-               z2d(:,:) = sye(:,:,jk,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar
-               z2d(:,:) = sxxe(:,:,jk,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'syye'//'_il'//zchar1//'_htc'//zchar
-               z2d(:,:) = syye(:,:,jk,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sxye'//'_il'//zchar1//'_htc'//zchar
-               z2d(:,:) = sxye(:,:,jk,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            END DO
-         END DO
-         IF ( ln_pnd_H12 ) THEN
-            DO jl = 1, jpl 
-               WRITE(zchar,'(I2.2)') jl
-               znam = 'sxap'//'_htc'//zchar
-               z2d(:,:) = sxap(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'syap'//'_htc'//zchar
-               z2d(:,:) = syap(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sxxap'//'_htc'//zchar
-               z2d(:,:) = sxxap(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'syyap'//'_htc'//zchar
-               z2d(:,:) = syyap(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sxyap'//'_htc'//zchar
-               z2d(:,:) = sxyap(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-   
-               znam = 'sxvp'//'_htc'//zchar
-               z2d(:,:) = sxvp(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'syvp'//'_htc'//zchar
-               z2d(:,:) = syvp(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sxxvp'//'_htc'//zchar
-               z2d(:,:) = sxxvp(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'syyvp'//'_htc'//zchar
-               z2d(:,:) = syyvp(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sxyvp'//'_htc'//zchar
-               z2d(:,:) = sxyvp(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            END DO
+         !
+         ! In case Prather scheme is used for advection, write second order moments
+         ! ------------------------------------------------------------------------
+         !
+         !                                                           ! ice thickness
+         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice  )
+         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice  )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
+         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
+         !                                                           ! snow thickness
+         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn   )
+         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn   )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn  )
+         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  )
+         !                                                           ! lead fraction
+         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    )
+         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa   )
+         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   )
+         !                                                           ! snow thermal content
+         CALL iom_rstput( iter, nitrst, numriw, 'sxc0'  , sxc0   )
+         CALL iom_rstput( iter, nitrst, numriw, 'syc0'  , syc0   )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxc0' , sxxc0  )
+         CALL iom_rstput( iter, nitrst, numriw, 'syyc0' , syyc0  )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxyc0' , sxyc0  )
+         !                                                           ! ice salinity
+         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  )
+         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal  )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
+         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
+         !                                                           ! ice age
+         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage  )
+         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage  )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
+         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
+         !                                                           ! open water in sea ice
+         CALL iom_rstput( iter, nitrst, numriw, 'sxopw ', sxopw  )
+         CALL iom_rstput( iter, nitrst, numriw, 'syopw ', syopw  )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxopw', sxxopw )
+         CALL iom_rstput( iter, nitrst, numriw, 'syyopw', syyopw )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxyopw', sxyopw )
+         !                                                           ! ice layers heat content
+         DO jk = 1, nlay_i 
+            WRITE(zchar1,'(I2.2)') jk
+            znam = 'sxe'//'_il'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
+            znam = 'sye'//'_il'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
+            znam = 'sxxe'//'_il'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
+            znam = 'syye'//'_il'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
+            znam = 'sxye'//'_il'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
+         END DO
+         !
+         IF( ln_pnd_H12 ) THEN                                       ! melt pond fraction
+            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
+            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
+            CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
+            CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
+            CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
+            !                                                        ! melt pond volume
+            CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp  )
+            CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp  )
+            CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
+            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
+            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
          ENDIF
          !
Index: /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90
===================================================================
--- /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 8816)
+++ /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 8817)
@@ -6,5 +6,6 @@
    !! History:   3.0  ! 2005-04 (M. Vancoppenolle) Original code
    !!             -   ! 2008-03 (C. Ethe) restart files in using IOM interface
-   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
+   !!            4.0  ! 2017-11 (M. Andrejczuk) Extend IOM interface to write 3D ice fields
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -82,5 +83,5 @@
             ENDIF
             !
-            CALL iom_open( TRIM(clpath)//TRIM(clname), numriw, ldwrt = .TRUE., kiolib = jprstlib )
+            CALL iom_open( TRIM(clpath)//TRIM(clname), numriw, ldwrt = .TRUE., kiolib = jprstlib, kdlev = jpl )
             lrst_ice = .TRUE.
          ENDIF
@@ -88,4 +89,5 @@
       !
       IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - Beginning the time step - ' )   ! control print
+      !
    END SUBROUTINE ice_rst_opn
 
@@ -103,5 +105,5 @@
       CHARACTER(len=25) ::   znam
       CHARACTER(len=2)  ::   zchar, zchar1
-      REAL(wp), DIMENSION(jpi,jpj) :: z2d
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z3d   ! 3D workspace
       !!----------------------------------------------------------------------
 
@@ -120,54 +122,30 @@
       CALL iom_rstput( iter, nitrst, numriw, 'kt_ice' , REAL( iter   , wp ) )      ! date
 
-!!gm   It is possible and easy to define a 3D domain size (jpi,jpj,jpl) or use a SIZE( tab, DIM=3) in iom_rtput )
-!!gm         ===>>> just a simple   iom_rstput( iter, nitrst, numriw, 'v_i', v_i )  etc...
-!!gm   "just" ask Sebatien
-
-      ! Prognostic variables 
-      DO jl = 1, jpl 
-         WRITE(zchar,'(I2.2)') jl
-         znam = 'v_i'//'_htc'//zchar
-         z2d(:,:) = v_i(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_i
-         znam = 'v_s'//'_htc'//zchar
-         z2d(:,:) = v_s(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_s
-         znam = 'sv_i'//'_htc'//zchar
-         z2d(:,:) = sv_i(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! sv_i
-         znam = 'oa_i'//'_htc'//zchar
-         z2d(:,:) = oa_i(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! oa_i
-         znam = 'a_i'//'_htc'//zchar
-         z2d(:,:) = a_i(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_i
-         znam = 't_su'//'_htc'//zchar
-         z2d(:,:) = t_su(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! t_su
+      ! Prognostic variables
+      CALL iom_rstput( iter, nitrst, numriw, 'v_i' , v_i  )
+      CALL iom_rstput( iter, nitrst, numriw, 'v_s' , v_s  )
+      CALL iom_rstput( iter, nitrst, numriw, 'sv_i', sv_i )
+      CALL iom_rstput( iter, nitrst, numriw, 'oa_i', oa_i )
+      CALL iom_rstput( iter, nitrst, numriw, 'a_i' , a_i  )
+      CALL iom_rstput( iter, nitrst, numriw, 't_su', t_su )
+      !
+      ! Melt ponds
+      CALL iom_rstput( iter, nitrst, numriw, 'a_ip', a_ip )
+      CALL iom_rstput( iter, nitrst, numriw, 'v_ip', v_ip )
+      !
+!!gm dangerous !!!!!  ===>>>> better reading writing all snow layers !
+      ! Snow enthalpy (1st snow layer only)
+      z3d = e_s(:,:,1,:)
+      CALL iom_rstput( iter, nitrst, numriw, 'tempt_sl1' , z3d )
+      !
+      ! Ice enthalpy (all ice layers)
+      DO jk = 1, nlay_i 
+         WRITE(zchar1,'(I2.2)') jk
+         znam = 'tempt'//'_il'//zchar1
+         z3d(:,:,:) = e_i(:,:,jk,:)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
       END DO
-
-      DO jl = 1, jpl 
-         WRITE(zchar,'(I2.2)') jl
-         znam = 'a_ip'//'_htc'//zchar
-         z2d(:,:) = a_ip(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_ip
-         znam = 'v_ip'//'_htc'//zchar
-         z2d(:,:) = v_ip(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_ip
-      END DO
-
-      DO jl = 1, jpl 
-         WRITE(zchar,'(I2.2)') jl
-         znam = 'tempt_sl1'//'_htc'//zchar
-         z2d(:,:) = e_s(:,:,1,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! e_s
-         DO jk = 1, nlay_i 
-            WRITE(zchar1,'(I2.2)') jk
-            znam = 'tempt'//'_il'//zchar1//'_htc'//zchar
-            z2d(:,:) = e_i(:,:,jk,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! e_i
-         END DO
-      END DO
-
+      !
+      ! ice velocity
       CALL iom_rstput( iter, nitrst, numriw, 'u_ice', u_ice ) ! u_ice
       CALL iom_rstput( iter, nitrst, numriw, 'v_ice', v_ice ) ! v_ice
@@ -189,12 +167,12 @@
       !! ** purpose  :   read restart file
       !!----------------------------------------------------------------------
-      INTEGER  ::   jk, jl
-      INTEGER  ::   id1            ! local integer
-      REAL(wp) ::   zfice, ziter
-      REAL(wp), DIMENSION(jpi,jpj) ::   z2d
+      INTEGER           ::   jk, jl
+      LOGICAL           ::   llok
+      INTEGER           ::   id1            ! local integer
+      INTEGER           ::   jlibalt = jprstlib
       CHARACTER(len=25) ::   znam
       CHARACTER(len=2)  ::   zchar, zchar1
-      INTEGER           ::   jlibalt = jprstlib
-      LOGICAL           ::   llok
+      REAL(wp)          ::   zfice, ziter
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z3d   ! 3D workspace
       !!----------------------------------------------------------------------
 
@@ -205,5 +183,5 @@
       ENDIF
 
-      CALL iom_open ( TRIM(cn_icerst_indir)//'/'//cn_icerst_in, numrir, kiolib = jprstlib )
+      CALL iom_open ( TRIM(cn_icerst_indir)//'/'//cn_icerst_in, numrir, kiolib = jprstlib, kdlev = jpl )
 
       CALL iom_get( numrir, 'nn_fsbc', zfice )
@@ -223,56 +201,36 @@
 
       ! Prognostic variables 
-      DO jl = 1, jpl 
-         WRITE(zchar,'(I2.2)') jl
-         znam = 'v_i'//'_htc'//zchar
-         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-         v_i(:,:,jl) = z2d(:,:)
-         znam = 'v_s'//'_htc'//zchar
-         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-         v_s(:,:,jl) = z2d(:,:) 
-         znam = 'sv_i'//'_htc'//zchar
-         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-         sv_i(:,:,jl) = z2d(:,:)
-         znam = 'oa_i'//'_htc'//zchar
-         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-         oa_i(:,:,jl) = z2d(:,:)
-         znam = 'a_i'//'_htc'//zchar
-         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-         a_i(:,:,jl) = z2d(:,:)
-         znam = 't_su'//'_htc'//zchar
-         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-         t_su(:,:,jl) = z2d(:,:)
-      END DO
-
+      CALL iom_get( numrir, jpdom_autoglo, 'v_i' , v_i  )
+      CALL iom_get( numrir, jpdom_autoglo, 'v_s' , v_s  )
+      CALL iom_get( numrir, jpdom_autoglo, 'sv_i', sv_i )
+      CALL iom_get( numrir, jpdom_autoglo, 'oa_i', oa_i )
+      CALL iom_get( numrir, jpdom_autoglo, 'a_i' , a_i  )
+      CALL iom_get( numrir, jpdom_autoglo, 't_su', t_su )
+      !
+      ! Melt ponds
       id1 = iom_varid( numrir, 'a_ip_htc01' , ldstop = .FALSE. )
       IF( id1 > 0 ) THEN                       ! fields exist (melt ponds)
-         DO jl = 1, jpl 
-            WRITE(zchar,'(I2.2)') jl
-            znam = 'a_ip'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            a_ip(:,:,jl) = z2d(:,:)
-            znam = 'v_ip'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            v_ip(:,:,jl) = z2d(:,:)
-         END DO
+         CALL iom_get( numrir, jpdom_autoglo, 'a_ip' , a_ip )
+         CALL iom_get( numrir, jpdom_autoglo, 'v_ip' , v_ip )
       ELSE                                     ! start from rest
-         IF(lwp) WRITE(numout,*) '   ==>>   previous run without melt ponds output then set it'
+         IF(lwp) WRITE(numout,*) '   ==>>   previous run without melt ponds output then set it to zero'
          a_ip(:,:,:) = 0._wp
          v_ip(:,:,:) = 0._wp
       ENDIF
-
-      DO jl = 1, jpl 
-         WRITE(zchar,'(I2.2)') jl
-         znam = 'tempt_sl1'//'_htc'//zchar
-         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-         e_s(:,:,1,jl) = z2d(:,:)
-         DO jk = 1, nlay_i 
-            WRITE(zchar1,'(I2.2)') jk
-            znam = 'tempt'//'_il'//zchar1//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            e_i(:,:,jk,jl) = z2d(:,:)
-         END DO
+      !
+!!gm dangerous !!!!!  ===>>>> better reading writing all snow layers !
+      ! Snow enthalpy (1st snow layer only)
+      CALL iom_get( numrir, jpdom_autoglo, 'tempt_sl1' , z3d )
+      e_s(:,:,1,:) = z3d
+      !
+      ! Ice enthalpy (all ice layers)
+      DO jk = 1, nlay_i 
+         WRITE(zchar1,'(I2.2)') jk
+         znam = 'tempt'//'_il'//zchar1
+         CALL iom_get( numrir, jpdom_autoglo, znam , z3d )
+         e_i(:,:,jk,:) = z3d(:,:,:)
       END DO
-
+      !
+      ! ice velocity
       CALL iom_get( numrir, jpdom_autoglo, 'u_ice', u_ice )
       CALL iom_get( numrir, jpdom_autoglo, 'v_ice', v_ice )
Index: /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 8816)
+++ /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 8817)
@@ -1,7 +1,7 @@
 MODULE iom
-   !!=====================================================================
+   !!======================================================================
    !!                    ***  MODULE  iom ***
    !! Input/Output manager :  Library to read input files
-   !!====================================================================
+   !!======================================================================
    !! History :  2.0  ! 2005-12  (J. Belier) Original code
    !!            2.0  ! 2006-02  (S. Masson) Adaptation to NEMO
@@ -10,7 +10,8 @@
    !!            3.6  ! 2014-15  DIMG format removed
    !!            3.6  ! 2015-15  (J. Harle) Added procedure to read REAL attributes
-   !!--------------------------------------------------------------------
-
-   !!--------------------------------------------------------------------
+   !!            4.0  ! 2017-11  (M. Andrejczuk) Extend IOM interface to write any 3D fields
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
    !!   iom_open       : open a file read only
    !!   iom_close      : close a file or all files opened by iom
@@ -19,5 +20,5 @@
    !!   iom_varid      : get the id of a variable in a file
    !!   iom_rstput     : write a field in a restart file (interfaced to several routines)
-   !!--------------------------------------------------------------------
+   !!----------------------------------------------------------------------
    USE dom_oce         ! ocean space and time domain
    USE c1d             ! 1D vertical configuration
@@ -29,9 +30,9 @@
    USE lib_mpp           ! MPP library
 #if defined key_iomput
-   USE sbc_oce, ONLY :   nn_fsbc         ! ocean space and time domain
-   USE trc_oce, ONLY :   nn_dttrc        !  !: frequency of step on passive tracers
-   USE icb_oce, ONLY :   nclasses, class_num       !  !: iceberg classes
+   USE sbc_oce  , ONLY :   nn_fsbc         ! ocean space and time domain
+   USE trc_oce  , ONLY :   nn_dttrc        !  !: frequency of step on passive tracers
+   USE icb_oce  , ONLY :   nclasses, class_num       !  !: iceberg classes
 #if defined key_lim3
-   USE ice    , ONLY :   jpl
+   USE ice      , ONLY :   jpl
 #endif
    USE domngb          ! ocean space and time domain
@@ -80,9 +81,8 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -95,6 +95,7 @@
       !!----------------------------------------------------------------------
       CHARACTER(len=*), INTENT(in)  :: cdname
+      !
 #if defined key_iomput
-
+      !
       TYPE(xios_duration) :: dtime    = xios_duration(0, 0, 0, 0, 0, 0)
       TYPE(xios_date)     :: start_date
@@ -104,7 +105,7 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zt_bnds, zw_bnds
       !!----------------------------------------------------------------------
-
+      !
       ALLOCATE( zt_bnds(2,jpk), zw_bnds(2,jpk) )
-
+      !
       clname = cdname
       IF( TRIM(Agrif_CFixed()) /= '0' )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(cdname)
@@ -125,5 +126,5 @@
       ! horizontal grid definition
       CALL set_scalar
-
+      !
       IF( TRIM(cdname) == TRIM(cxios_context) ) THEN  
          CALL set_grid( "T", glamt, gphit ) 
@@ -144,5 +145,5 @@
          ENDIF
       ENDIF
-
+      !
       IF( TRIM(cdname) == TRIM(cxios_context)//"_crs" ) THEN  
          CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
@@ -167,5 +168,5 @@
          ENDIF
       ENDIF
-
+      !
       ! vertical grid definition
       CALL iom_set_axis_attr( "deptht", gdept_1d )
@@ -173,5 +174,5 @@
       CALL iom_set_axis_attr( "depthv", gdept_1d )
       CALL iom_set_axis_attr( "depthw", gdepw_1d )
-
+      !
       ! Add vertical grid bounds
       jkmin = MIN(2,jpk)  ! in case jpk=1 (i.e. sas2D)
@@ -186,14 +187,13 @@
       CALL iom_set_axis_attr( "depthv", bounds=zt_bnds )
       CALL iom_set_axis_attr( "depthw", bounds=zw_bnds )
-
-
+      !
 # if defined key_floats
       CALL iom_set_axis_attr( "nfloat", (/ (REAL(ji,wp), ji=1,nfloat) /) )
 # endif
-#if defined key_lim3
+# if defined key_lim3
       CALL iom_set_axis_attr( "ncatice", (/ (REAL(ji,wp), ji=1,jpl) /) )
       ! SIMIP diagnostics (4 main arctic straits)
       CALL iom_set_axis_attr( "nstrait", (/ (REAL(ji,wp), ji=1,4) /) )
-#endif
+# endif
       CALL iom_set_axis_attr( "icbcla", class_num )
       CALL iom_set_axis_attr( "iax_20C", (/ REAL(20,wp) /) )
@@ -202,5 +202,5 @@
       ! automatic definitions of some of the xml attributs
       CALL set_xmlatt
-
+      !
       ! end file definition
       dtime%second = rdt
@@ -209,9 +209,9 @@
       
       CALL xios_update_calendar(0)
-
+      !
       DEALLOCATE( zt_bnds, zw_bnds )
-
+      !
 #endif
-      
+      !
    END SUBROUTINE iom_init
 
@@ -239,5 +239,5 @@
 
 
-   SUBROUTINE iom_open( cdname, kiomid, ldwrt, kdom, kiolib, ldstop, ldiof )
+   SUBROUTINE iom_open( cdname, kiomid, ldwrt, kdom, kiolib, ldstop, ldiof, kdlev )
       !!---------------------------------------------------------------------
       !!                   ***  SUBROUTINE  iom_open  ***
@@ -252,4 +252,5 @@
       LOGICAL         , INTENT(in   ), OPTIONAL ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.)
       LOGICAL         , INTENT(in   ), OPTIONAL ::   ldiof    ! Interp On the Fly, needed for AGRIF (default = .FALSE.)
+      INTEGER         , INTENT(in   ), OPTIONAL ::   kdlev    ! number of vertical levels
 
       CHARACTER(LEN=256)    ::   clname    ! the name of the file based on cdname [[+clcpu]+clcpu]
@@ -405,5 +406,5 @@
       IF( istop == nstop ) THEN   ! no error within this routine
          SELECT CASE (iolib)
-         CASE (jpnf90   )   ;   CALL iom_nf90_open(    clname, kiomid, llwrt, llok, idompar )
+         CASE (jpnf90   )   ;   CALL iom_nf90_open(    clname, kiomid, llwrt, llok, idompar, kdlev = kdlev )
          CASE DEFAULT
             CALL ctl_stop( TRIM(clinfo)//' accepted IO library is only jpnf90 (jpioipsl option has been removed) ' )
@@ -672,6 +673,9 @@
       CHARACTER(LEN=1)               ::   clrankpv, cldmspc      ! 
       LOGICAL                        ::   ll_depth_spec ! T => if kstart, kcount present then *only* use values for 3rd spatial dimension.
+      INTEGER                        ::   inlev       ! number of levels for 3D data
       !---------------------------------------------------------------------
       !
+      inlev = -1
+      IF( PRESENT(pv_r3d) )   inlev = SIZE(pv_r3d, 3)
       clname = iom_file(kiomid)%name   !   esier to read
       clinfo = '          iom_get_123d, file: '//trim(clname)//', var: '//trim(cdvar)
@@ -774,7 +778,10 @@
          istart(idmspc+1) = itime
 
-         IF( PRESENT(kstart) .AND. .NOT. ll_depth_spec ) THEN ; istart(1:idmspc) = kstart(1:idmspc) ; icnt(1:idmspc) = kcount(1:idmspc)
+         IF( PRESENT(kstart) .AND. .NOT. ll_depth_spec ) THEN 
+            istart(1:idmspc) = kstart(1:idmspc) 
+            icnt(1:idmspc) = kcount(1:idmspc)
          ELSE
-            IF(           idom == jpdom_unknown ) THEN                                                ; icnt(1:idmspc) = idimsz(1:idmspc)
+            IF(idom == jpdom_unknown ) THEN
+               icnt(1:idmspc) = idimsz(1:idmspc)
             ELSE 
                IF( .NOT. PRESENT(pv_r1d) ) THEN   !   not a 1D array
@@ -799,7 +806,7 @@
                   ENDIF
                   IF( PRESENT(pv_r3d) ) THEN
-                     IF( idom == jpdom_data ) THEN                                  ; icnt(3) = jpkglo
+                     IF( idom == jpdom_data ) THEN                                  ; icnt(3) = inlev
                      ELSE IF( ll_depth_spec .AND. PRESENT(kstart) ) THEN            ; istart(3) = kstart(3); icnt(3) = kcount(3)
-                     ELSE                                                           ; icnt(3) = jpk
+                     ELSE                                                           ; icnt(3) = inlev
                      ENDIF
                   ENDIF
@@ -884,5 +891,5 @@
             ELSEIF( PRESENT(pv_r3d) .AND. idom /= jpdom_unknown ) THEN
                ! this if could be simplified with the new lbc_lnk that works with any size of the 3rd dimension
-               IF( icnt(3) == jpk ) THEN
+               IF( icnt(3) == inlev ) THEN
                   CALL lbc_lnk( pv_r3d,'Z',-999.,'no0' )
                ELSE   ! put some arbitrary value (a call to lbc_lnk will be done later...)
@@ -1133,4 +1140,5 @@
    END SUBROUTINE iom_rp0d
 
+
    SUBROUTINE iom_rp1d( kt, kwrite, kiomid, cdvar, pvar, ktype )
       INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
@@ -1153,4 +1161,5 @@
    END SUBROUTINE iom_rp1d
 
+
    SUBROUTINE iom_rp2d( kt, kwrite, kiomid, cdvar, pvar, ktype )
       INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
@@ -1173,4 +1182,5 @@
    END SUBROUTINE iom_rp2d
 
+
    SUBROUTINE iom_rp3d( kt, kwrite, kiomid, cdvar, pvar, ktype )
       INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
@@ -1234,26 +1244,31 @@
       REAL(wp),       DIMENSION(:,:,:), INTENT(in) ::   pfield3d
 #if defined key_iomput
-      CALL xios_send_field(cdname, pfield3d)
+      CALL xios_send_field( cdname, pfield3d )
 #else
       IF( .FALSE. )   WRITE(numout,*) cdname, pfield3d   ! useless test to avoid compilation warnings
 #endif
    END SUBROUTINE iom_p3d
+
+#if defined key_iomput
+
    !!----------------------------------------------------------------------
-
-#if defined key_iomput
+   !!   'key_iomput'                                         IOM  interface
+   !!----------------------------------------------------------------------
 
    SUBROUTINE iom_set_domain_attr( cdid, ni_glo, nj_glo, ibegin, jbegin, ni, nj, zoom_ibegin, zoom_jbegin, zoom_ni, zoom_nj,   &
       &                                    data_dim, data_ibegin, data_ni, data_jbegin, data_nj, lonvalue, latvalue, mask,     &
       &                                    nvertex, bounds_lon, bounds_lat, area )
-      CHARACTER(LEN=*)                   , INTENT(in) ::   cdid
-      INTEGER                  , OPTIONAL, INTENT(in) ::   ni_glo, nj_glo, ibegin, jbegin, ni, nj
-      INTEGER                  , OPTIONAL, INTENT(in) ::   data_dim, data_ibegin, data_ni, data_jbegin, data_nj
-      INTEGER                  , OPTIONAL, INTENT(in) ::   zoom_ibegin, zoom_jbegin, zoom_ni, zoom_nj, nvertex
-      REAL(wp), DIMENSION(:)   , OPTIONAL, INTENT(in) ::   lonvalue, latvalue
-      REAL(wp), DIMENSION(:,:) , OPTIONAL, INTENT(in) ::   bounds_lon, bounds_lat, area
-      LOGICAL, DIMENSION(:)   , OPTIONAL, INTENT(in) ::   mask
-
-
-      IF ( xios_is_valid_domain     (cdid) ) THEN
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      CHARACTER(LEN=*)                  , INTENT(in) ::   cdid
+      INTEGER                 , OPTIONAL, INTENT(in) ::   ni_glo, nj_glo, ibegin, jbegin, ni, nj
+      INTEGER                 , OPTIONAL, INTENT(in) ::   data_dim, data_ibegin, data_ni, data_jbegin, data_nj
+      INTEGER                 , OPTIONAL, INTENT(in) ::   zoom_ibegin, zoom_jbegin, zoom_ni, zoom_nj, nvertex
+      REAL(wp), DIMENSION(:)  , OPTIONAL, INTENT(in) ::   lonvalue, latvalue
+      REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT(in) ::   bounds_lon, bounds_lat, area
+      LOGICAL , DIMENSION(:)  , OPTIONAL, INTENT(in) ::   mask
+      !!----------------------------------------------------------------------
+      !
+      IF( xios_is_valid_domain     (cdid) ) THEN
          CALL xios_set_domain_attr     ( cdid, ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj,   &
             &    data_dim=data_dim, data_ibegin=data_ibegin, data_ni=data_ni, data_jbegin=data_jbegin, data_nj=data_nj ,   &
@@ -1261,5 +1276,5 @@
             &    bounds_lat_1D=bounds_lat, area=area, type='curvilinear')
      ENDIF
-      IF ( xios_is_valid_domaingroup(cdid) ) THEN
+      IF( xios_is_valid_domaingroup(cdid) ) THEN
          CALL xios_set_domaingroup_attr( cdid, ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj,   &
             &    data_dim=data_dim, data_ibegin=data_ibegin, data_ni=data_ni, data_jbegin=data_jbegin, data_nj=data_nj ,   &
@@ -1267,32 +1282,33 @@
             &    bounds_lat_1D=bounds_lat, area=area, type='curvilinear' )
       ENDIF
-
+      !
       CALL xios_solve_inheritance()
-
+      !
    END SUBROUTINE iom_set_domain_attr
 
 
-   SUBROUTINE iom_set_zoom_domain_attr( cdid, ibegin, jbegin, ni, nj)
-      CHARACTER(LEN=*)                   , INTENT(in) ::   cdid
-      INTEGER                  , OPTIONAL, INTENT(in) ::   ibegin, jbegin, ni, nj
-
-      IF ( xios_is_valid_zoom_domain     (cdid) ) THEN
-          CALL xios_set_zoom_domain_attr     ( cdid, ibegin=ibegin, jbegin=jbegin, ni=ni,    &
-            &   nj=nj)
-     ENDIF
+   SUBROUTINE iom_set_zoom_domain_attr( cdid, ibegin, jbegin, ni, nj )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      CHARACTER(LEN=*)          , INTENT(in) ::   cdid
+      INTEGER         , OPTIONAL, INTENT(in) ::   ibegin, jbegin, ni, nj
+      !!----------------------------------------------------------------------
+      IF( xios_is_valid_zoom_domain(cdid) )   CALL xios_set_zoom_domain_attr( cdid, ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj )
    END SUBROUTINE iom_set_zoom_domain_attr
 
 
    SUBROUTINE iom_set_axis_attr( cdid, paxis, bounds )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
       CHARACTER(LEN=*)      , INTENT(in) ::   cdid
       REAL(wp), DIMENSION(:)  , OPTIONAL, INTENT(in) ::   paxis
       REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT(in) ::   bounds
-
-      IF ( PRESENT(paxis) ) THEN
-         IF ( xios_is_valid_axis     (cdid) )   CALL xios_set_axis_attr     ( cdid, n_glo=SIZE(paxis), value=paxis )
-         IF ( xios_is_valid_axisgroup(cdid) )   CALL xios_set_axisgroup_attr( cdid, n_glo=SIZE(paxis), value=paxis )
-      ENDIF
-      IF ( xios_is_valid_axis     (cdid) )   CALL xios_set_axis_attr     ( cdid, bounds=bounds )
-      IF ( xios_is_valid_axisgroup(cdid) )   CALL xios_set_axisgroup_attr( cdid, bounds=bounds )
+      !!----------------------------------------------------------------------
+      IF( PRESENT(paxis) ) THEN
+         IF( xios_is_valid_axis     (cdid) )   CALL xios_set_axis_attr     ( cdid, n_glo=SIZE(paxis), value=paxis )
+         IF( xios_is_valid_axisgroup(cdid) )   CALL xios_set_axisgroup_attr( cdid, n_glo=SIZE(paxis), value=paxis )
+      ENDIF
+      IF( xios_is_valid_axis     (cdid) )   CALL xios_set_axis_attr     ( cdid, bounds=bounds )
+      IF( xios_is_valid_axisgroup(cdid) )   CALL xios_set_axisgroup_attr( cdid, bounds=bounds )
       CALL xios_solve_inheritance()
    END SUBROUTINE iom_set_axis_attr
@@ -1300,11 +1316,12 @@
 
    SUBROUTINE iom_set_field_attr( cdid, freq_op, freq_offset )
-      CHARACTER(LEN=*)          , INTENT(in) ::   cdid
-      TYPE(xios_duration),OPTIONAL , INTENT(in) ::   freq_op
-      TYPE(xios_duration),OPTIONAL , INTENT(in) ::   freq_offset
-      IF ( xios_is_valid_field     (cdid) )   CALL xios_set_field_attr       &
-    &     ( cdid, freq_op=freq_op, freq_offset=freq_offset )
-      IF ( xios_is_valid_fieldgroup(cdid) )   CALL xios_set_fieldgroup_attr  &
-    &                    ( cdid, freq_op=freq_op, freq_offset=freq_offset )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      CHARACTER(LEN=*)             , INTENT(in) ::   cdid
+      TYPE(xios_duration), OPTIONAL, INTENT(in) ::   freq_op
+      TYPE(xios_duration), OPTIONAL, INTENT(in) ::   freq_offset
+      !!----------------------------------------------------------------------
+      IF( xios_is_valid_field     (cdid) )   CALL xios_set_field_attr     ( cdid, freq_op=freq_op, freq_offset=freq_offset )
+      IF( xios_is_valid_fieldgroup(cdid) )   CALL xios_set_fieldgroup_attr( cdid, freq_op=freq_op, freq_offset=freq_offset )
       CALL xios_solve_inheritance()
    END SUBROUTINE iom_set_field_attr
@@ -1312,8 +1329,11 @@
 
    SUBROUTINE iom_set_file_attr( cdid, name, name_suffix )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
       CHARACTER(LEN=*)          , INTENT(in) ::   cdid
       CHARACTER(LEN=*),OPTIONAL , INTENT(in) ::   name, name_suffix
-      IF ( xios_is_valid_file     (cdid) )   CALL xios_set_file_attr     ( cdid, name=name, name_suffix=name_suffix )
-      IF ( xios_is_valid_filegroup(cdid) )   CALL xios_set_filegroup_attr( cdid, name=name, name_suffix=name_suffix )
+      !!----------------------------------------------------------------------
+      IF( xios_is_valid_file     (cdid) )   CALL xios_set_file_attr     ( cdid, name=name, name_suffix=name_suffix )
+      IF( xios_is_valid_filegroup(cdid) )   CALL xios_set_filegroup_attr( cdid, name=name, name_suffix=name_suffix )
       CALL xios_solve_inheritance()
    END SUBROUTINE iom_set_file_attr
@@ -1321,4 +1341,6 @@
 
    SUBROUTINE iom_get_file_attr( cdid, name, name_suffix, output_freq )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
       CHARACTER(LEN=*)          , INTENT(in ) ::   cdid
       CHARACTER(LEN=*),OPTIONAL , INTENT(out) ::   name, name_suffix
@@ -1329,5 +1351,5 @@
       IF( PRESENT( name_suffix ) )   name_suffix = ''
       IF( PRESENT( output_freq ) )   output_freq = xios_duration(0,0,0,0,0,0)
-      IF ( xios_is_valid_file     (cdid) ) THEN
+      IF( xios_is_valid_file     (cdid) ) THEN
          CALL xios_solve_inheritance()
          CALL xios_is_defined_file_attr     ( cdid, name = llexist1, name_suffix = llexist2, output_freq = llexist3)
@@ -1336,5 +1358,5 @@
          IF(llexist3)   CALL xios_get_file_attr     ( cdid, output_freq = output_freq )
       ENDIF
-      IF ( xios_is_valid_filegroup(cdid) ) THEN
+      IF( xios_is_valid_filegroup(cdid) ) THEN
          CALL xios_solve_inheritance()
          CALL xios_is_defined_filegroup_attr( cdid, name = llexist1, name_suffix = llexist2, output_freq = llexist3)
@@ -1347,24 +1369,30 @@
 
    SUBROUTINE iom_set_grid_attr( cdid, mask )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
       CHARACTER(LEN=*)                   , INTENT(in) ::   cdid
       LOGICAL, DIMENSION(:,:,:), OPTIONAL, INTENT(in) ::   mask
-      IF ( xios_is_valid_grid     (cdid) )   CALL xios_set_grid_attr     ( cdid, mask_3D=mask )
-      IF ( xios_is_valid_gridgroup(cdid) )   CALL xios_set_gridgroup_attr( cdid, mask_3D=mask )
+      !!----------------------------------------------------------------------
+      IF( xios_is_valid_grid     (cdid) )   CALL xios_set_grid_attr     ( cdid, mask_3D=mask )
+      IF( xios_is_valid_gridgroup(cdid) )   CALL xios_set_gridgroup_attr( cdid, mask_3D=mask )
       CALL xios_solve_inheritance()
    END SUBROUTINE iom_set_grid_attr
 
    SUBROUTINE iom_setkt( kt, cdname )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
       INTEGER         , INTENT(in) ::   kt 
       CHARACTER(LEN=*), INTENT(in) ::   cdname
-      !     
+      !!----------------------------------------------------------------------
       CALL iom_swap( cdname )   ! swap to cdname context
       CALL xios_update_calendar(kt)
-      IF( cdname /= TRIM(cxios_context) ) CALL iom_swap( TRIM(cxios_context) )   ! return back to nemo context
-      !
+      IF( cdname /= TRIM(cxios_context) )   CALL iom_swap( TRIM(cxios_context) )   ! return back to nemo context
    END SUBROUTINE iom_setkt
 
    SUBROUTINE iom_context_finalize( cdname )
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
       CHARACTER(LEN=*), INTENT(in) :: cdname
-      !
+      !!----------------------------------------------------------------------
       IF( xios_is_valid_context(cdname) ) THEN
          CALL iom_swap( cdname )   ! swap to cdname context
@@ -1372,5 +1400,4 @@
          IF( cdname /= TRIM(cxios_context) ) CALL iom_swap( TRIM(cxios_context) )   ! return back to nemo context
       ENDIF
-      !
    END SUBROUTINE iom_context_finalize
 
@@ -1381,5 +1408,4 @@
       !!
       !! ** Purpose :   define horizontal grids
-      !!
       !!----------------------------------------------------------------------
       CHARACTER(LEN=1)            , INTENT(in) ::   cdgrd
@@ -1387,14 +1413,15 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   plat
       !
+      INTEGER  :: ni,nj
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmask
-      INTEGER  :: ni,nj
-      
+      !!----------------------------------------------------------------------
+      !
       ni=nlei-nldi+1 ; nj=nlej-nldj+1
-
+      !
       CALL iom_set_domain_attr("grid_"//cdgrd, ni_glo=jpiglo, nj_glo=jpjglo, ibegin=nimpp+nldi-2, jbegin=njmpp+nldj-2, ni=ni, nj=nj)
       CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = 1-nldi, data_ni = jpi, data_jbegin = 1-nldj, data_nj = jpj)
       CALL iom_set_domain_attr("grid_"//cdgrd, lonvalue = RESHAPE(plon(nldi:nlei, nldj:nlej),(/ ni*nj /)),   &
          &                                     latvalue = RESHAPE(plat(nldi:nlei, nldj:nlej),(/ ni*nj /)))  
-
+      !
       IF ( ln_mskland ) THEN
          ! mask land points, keep values on coast line -> specific mask for U, V and W points
@@ -1409,5 +1436,5 @@
          CALL iom_set_grid_attr  ( "grid_"//cdgrd//"_3D", mask = RESHAPE(zmask(nldi:nlei,nldj:nlej,:),(/ni,nj,jpk/)) /= 0. )
       ENDIF
-      
+      !
    END SUBROUTINE set_grid
 
@@ -1420,34 +1447,30 @@
       !!
       !!----------------------------------------------------------------------
-      CHARACTER(LEN=1) , INTENT(in) :: cdgrd
-      !
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)           :: plon_cnr, plat_cnr  ! Lat/lon coordinates of a contiguous vertex of cell (i,j)
-      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: plon_pnt, plat_pnt  ! Lat/lon coordinates of the point of cell (i,j)
-      !
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)   :: z_bnds      ! Lat/lon coordinates of the vertices of cell (i,j)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)       :: z_fld       ! Working array to determine where to rotate cells
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)       :: z_rot       ! Lat/lon working array for rotation of cells
-      !
-      INTEGER :: icnr, jcnr                                      ! Offset such that the vertex coordinate (i+icnr,j+jcnr)
-      !                                                          ! represents the bottom-left corner of cell (i,j)
+      CHARACTER(LEN=1)                      , INTENT(in) :: cdgrd
+      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) :: plon_cnr, plat_cnr  ! Lat/lon coord. of a contiguous vertex of cell (i,j)
+      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: plon_pnt, plat_pnt  ! Lat/lon coord. of the point of cell (i,j)
+      !
       INTEGER :: ji, jj, jn, ni, nj
-
+      INTEGER :: icnr, jcnr                                    ! Offset such that the vertex coordinate (i+icnr,j+jcnr)
+      !                                                        ! represents the bottom-left corner of cell (i,j)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z_bnds      ! Lat/lon coordinates of the vertices of cell (i,j)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z_fld       ! Working array to determine where to rotate cells
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z_rot       ! Lat/lon working array for rotation of cells
+      !!----------------------------------------------------------------------
+      !
       ALLOCATE( z_bnds(4,jpi,jpj,2), z_fld(jpi,jpj), z_rot(4,2)  )
-
+      !
       ! Offset of coordinate representing bottom-left corner
       SELECT CASE ( TRIM(cdgrd) )
-         CASE ('T', 'W')
-            icnr = -1 ; jcnr = -1
-         CASE ('U')
-            icnr =  0 ; jcnr = -1
-         CASE ('V')
-            icnr = -1 ; jcnr =  0
+      CASE ('T', 'W')   ;   icnr = -1   ;   jcnr = -1
+      CASE ('U')        ;   icnr =  0   ;   jcnr = -1
+      CASE ('V')        ;   icnr = -1   ;   jcnr =  0
       END SELECT
-
+      !
       ni = nlei-nldi+1 ; nj = nlej-nldj+1  ! Dimensions of subdomain interior
-
+      !
       z_fld(:,:) = 1._wp
       CALL lbc_lnk( z_fld, cdgrd, -1. )    ! Working array for location of northfold
-
+      !
       ! Cell vertices that can be defined
       DO jj = 2, jpjm1
@@ -1463,5 +1486,5 @@
          END DO
       END DO
-
+      !
       ! Cell vertices on boundries
       DO jn = 1, 4
@@ -1469,5 +1492,5 @@
          CALL lbc_lnk( z_bnds(jn,:,:,2), cdgrd, 1., pval=999._wp )
       END DO
-
+      !
       ! Zero-size cells at closed boundaries if cell points provided,
       ! otherwise they are closed cells with unrealistic bounds
@@ -1494,7 +1517,6 @@
          ENDIF
       ENDIF
-
-      ! Rotate cells at the north fold
-      IF( (nbondj == 1 .OR. nbondj == 2) .AND. jperio >= 3 ) THEN
+      !
+      IF( (nbondj == 1 .OR. nbondj == 2) .AND. jperio >= 3 ) THEN    ! Rotate cells at the north fold
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -1506,7 +1528,5 @@
             END DO
          END DO
-
-      ! Invert cells at the symmetric equator
-      ELSE IF( nbondj == 2 .AND. jperio == 2 ) THEN
+      ELSE IF( nbondj == 2 .AND. jperio == 2 ) THEN                  ! Invert cells at the symmetric equator
          DO ji = 1, jpi
             z_rot(1:2,:) = z_bnds(3:4,ji,1,:)
@@ -1515,10 +1535,10 @@
          END DO
       ENDIF
-
+      !
       CALL iom_set_domain_attr("grid_"//cdgrd, bounds_lat = RESHAPE(z_bnds(:,nldi:nlei,nldj:nlej,1),(/ 4,ni*nj /)),           &
-                                               bounds_lon = RESHAPE(z_bnds(:,nldi:nlei,nldj:nlej,2),(/ 4,ni*nj /)), nvertex=4 )
-
+          &                                    bounds_lon = RESHAPE(z_bnds(:,nldi:nlei,nldj:nlej,2),(/ 4,ni*nj /)), nvertex=4 )
+      !
       DEALLOCATE( z_bnds, z_fld, z_rot ) 
-
+      !
    END SUBROUTINE set_grid_bounds
 
@@ -1535,10 +1555,11 @@
       REAL(wp), DIMENSION(:), ALLOCATABLE  ::   zlon
       INTEGER  :: ni,nj, ix, iy
-
-      
+      !!----------------------------------------------------------------------
+      !
       ni=nlei-nldi+1 ; nj=nlej-nldj+1            ! define zonal mean domain (jpj*jpk)
       ALLOCATE( zlon(ni*nj) )       ;       zlon(:) = 0.
-
-      CALL dom_ngb( 180., 90., ix, iy, 'T' ) !  i-line that passes near the North Pole : Reference latitude (used in plots)
+      !
+      CALL dom_ngb( -168.53, 65.03, ix, iy, 'T' ) !  i-line that passes through Bering Strait: Reference latitude (used in plots)
+!      CALL dom_ngb( 180., 90., ix, iy, 'T' ) !  i-line that passes near the North Pole : Reference latitude (used in plots)
       CALL iom_set_domain_attr("gznl", ni_glo=jpiglo, nj_glo=jpjglo, ibegin=nimpp+nldi-2, jbegin=njmpp+nldj-2, ni=ni, nj=nj)
       CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = 1-nldi, data_ni = jpi, data_jbegin = 1-nldj, data_nj = jpj)
@@ -1551,4 +1572,5 @@
    END SUBROUTINE set_grid_znl
 
+
    SUBROUTINE set_scalar
       !!----------------------------------------------------------------------
@@ -1560,11 +1582,11 @@
       REAL(wp), DIMENSION(1)   ::   zz = 1.
       !!----------------------------------------------------------------------
-
+      !
       CALL iom_set_domain_attr('scalarpoint', ni_glo=jpnij, nj_glo=1, ibegin=narea-1, jbegin=0, ni=1, nj=1)
       CALL iom_set_domain_attr('scalarpoint', data_dim=2, data_ibegin = 1, data_ni = 1, data_jbegin = 1, data_nj = 1)
-      
-      zz=REAL(narea,wp)
+      !
+      zz = REAL( narea, wp )
       CALL iom_set_domain_attr('scalarpoint', lonvalue=zz, latvalue=zz)
-
+      !
    END SUBROUTINE set_scalar
 
@@ -1637,10 +1659,9 @@
       zlatpira = (/ -19.0, -14.0,  -8.0, 0.0, 4.0, 8.0, 12.0, 15.0, 20.0 /)
       CALL set_mooring( zlonpira, zlatpira )
-
-      
+      !
    END SUBROUTINE set_xmlatt
 
 
-   SUBROUTINE set_mooring( plon, plat)
+   SUBROUTINE set_mooring( plon, plat )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE set_mooring  ***
@@ -1649,5 +1670,5 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:), INTENT(in) ::  plon, plat           ! longitudes/latitudes oft the mooring
+      REAL(wp), DIMENSION(:), INTENT(in) ::   plon, plat   ! longitudes/latitudes oft the mooring
       !
 !!$      CHARACTER(len=1),DIMENSION(4) ::   clgrd = (/ 'T', 'U', 'V', 'W' /)   ! suffix name
@@ -1798,5 +1819,5 @@
                idx = INDEX(clname,'@enddate@') + INDEX(clname,'@ENDDATE@')
             END DO
-
+            !
             idx = INDEX(clname,'@enddatefull@') + INDEX(clname,'@ENDDATEFULL@')
             DO WHILE ( idx /= 0 ) 
@@ -1805,13 +1826,13 @@
                idx = INDEX(clname,'@enddatefull@') + INDEX(clname,'@ENDDATEFULL@')
             END DO
-
+            !
             IF( jn == 1 .AND. TRIM(Agrif_CFixed()) /= '0' )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
             IF( jn == 1 )   CALL iom_set_file_attr( cdid, name        = clname )
             IF( jn == 2 )   CALL iom_set_file_attr( cdid, name_suffix = clname )
-
-         ENDIF
-
+            !
+         ENDIF
+         !
       END DO
-
+      !
    END SUBROUTINE iom_update_file_name
 
@@ -1822,5 +1843,4 @@
       !!
       !! ** Purpose :   send back the date corresponding to the given julian day
-      !!
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   )           ::   pjday         ! julian day
@@ -1833,27 +1853,28 @@
       REAL(wp)          ::   zsec
       LOGICAL           ::   ll24, llfull
+      !!----------------------------------------------------------------------
       !
       IF( PRESENT(ld24) ) THEN   ;   ll24 = ld24
       ELSE                       ;   ll24 = .FALSE.
       ENDIF
-
+      !
       IF( PRESENT(ldfull) ) THEN   ;   llfull = ldfull
       ELSE                         ;   llfull = .FALSE.
       ENDIF
-
+      !
       CALL ju2ymds( pjday, iyear, imonth, iday, zsec )
       isec = NINT(zsec)
-
+      !
       IF ( ll24 .AND. isec == 0 ) THEN   ! 00:00 of the next day -> move to 24:00 of the current day
          CALL ju2ymds( pjday - 1., iyear, imonth, iday, zsec )
          isec = 86400
       ENDIF
-
+      !
       IF( iyear < 10000 ) THEN   ;   clfmt = "i4.4,2i2.2"                ! format used to write the date 
       ELSE                       ;   WRITE(clfmt, "('i',i1,',2i2.2')") INT(LOG10(REAL(iyear,wp))) + 1
       ENDIF
-      
+      !
 !$AGRIF_DO_NOT_TREAT      
-! Should be fixed in the conv
+      ! needed in the conv
       IF( llfull ) THEN 
          clfmt = TRIM(clfmt)//",'_',i2.2,':',i2.2,':',i2.2"
@@ -1867,9 +1888,8 @@
       ENDIF
 !$AGRIF_END_DO_NOT_TREAT      
-
+      !
    END FUNCTION iom_sdate
 
 #else
-
 
    SUBROUTINE iom_setkt( kt, cdname )
Index: /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90
===================================================================
--- /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90	(revision 8816)
+++ /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90	(revision 8817)
@@ -1,16 +1,11 @@
 MODULE iom_def
-   !!=====================================================================
+   !!======================================================================
    !!                    ***  MODULE  iom_def ***
    !! IOM variables definitions
-   !!====================================================================
-   !! History :  9.0  ! 06 09  (S. Masson) Original code
-   !!             "   ! 07 07  (D. Storkey) Add uldname
-   !!--------------------------------------------------------------------
-   !!---------------------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!---------------------------------------------------------------------------------
-
+   !!======================================================================
+   !! History :  9.0  ! 2006 09  (S. Masson) Original code
+   !!             -   ! 2007 07  (D. Storkey) Add uldname
+   !!            4.0  ! 2017-11 (M. Andrejczuk) Extend IOM interface to write any 3D fields
+   !!----------------------------------------------------------------------
    USE par_kind
 
@@ -64,8 +59,13 @@
       REAL(kind=wp), DIMENSION(jpmax_vars)      ::   scf      !: scale_factor of the variables
       REAL(kind=wp), DIMENSION(jpmax_vars)      ::   ofs      !: add_offset of the variables
+      INTEGER                                   ::   nlev     ! number of vertical levels
    END TYPE file_descriptor
    TYPE(file_descriptor), DIMENSION(jpmax_files), PUBLIC ::   iom_file !: array containing the info for all opened files
 !$AGRIF_END_DO_NOT_TREAT
 
-   !!=====================================================================
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!======================================================================
 END MODULE iom_def
Index: /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90
===================================================================
--- /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90	(revision 8816)
+++ /branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90	(revision 8817)
@@ -1,15 +1,15 @@
 MODULE iom_nf90
-   !!=====================================================================
+   !!======================================================================
    !!                    ***  MODULE  iom_nf90 ***
    !! Input/Output manager :  Library to read input files with NF90 (only fliocom module)
-   !!====================================================================
+   !!======================================================================
    !! History :  9.0  ! 05 12  (J. Belier) Original code
    !!            9.0  ! 06 02  (S. Masson) Adaptation to NEMO
    !!             "   ! 07 07  (D. Storkey) Changes to iom_nf90_gettime
    !!            3.6  ! 2015-15  (J. Harle) Added procedure to read REAL attributes
-   !!--------------------------------------------------------------------
-   !!gm  caution add !DIR nec: improved performance to be checked as well as no result changes
-
-   !!--------------------------------------------------------------------
+   !!            4.0  ! 2017-11 (M. Andrejczuk) Extend IOM interface to write any 3D fields
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
    !!   iom_open       : open a file read only
    !!   iom_close      : close a file or all files opened by iom
@@ -18,5 +18,5 @@
    !!   iom_varid      : get the id of a variable in a file
    !!   iom_rstput     : write a field in a restart file (interfaced to several routines)
-   !!--------------------------------------------------------------------
+   !!----------------------------------------------------------------------
    USE dom_oce         ! ocean space and time domain
    USE lbclnk          ! lateal boundary condition / mpp exchanges
@@ -29,5 +29,5 @@
    PRIVATE
 
-   PUBLIC iom_nf90_open, iom_nf90_close, iom_nf90_varid, iom_nf90_get, iom_nf90_gettime, iom_nf90_rstput
+   PUBLIC iom_nf90_open  , iom_nf90_close, iom_nf90_varid, iom_nf90_get, iom_nf90_gettime, iom_nf90_rstput
    PUBLIC iom_nf90_getatt, iom_nf90_putatt
 
@@ -46,12 +46,11 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
-   SUBROUTINE iom_nf90_open( cdname, kiomid, ldwrt, ldok, kdompar )
+   SUBROUTINE iom_nf90_open( cdname, kiomid, ldwrt, ldok, kdompar, kdlev )
       !!---------------------------------------------------------------------
       !!                   ***  SUBROUTINE  iom_open  ***
@@ -64,4 +63,5 @@
       LOGICAL                , INTENT(in   )           ::   ldok        ! check the existence 
       INTEGER, DIMENSION(2,5), INTENT(in   ), OPTIONAL ::   kdompar     ! domain parameters: 
+      INTEGER                , INTENT(in   ), OPTIONAL ::   kdlev       ! size of the third dimension
 
       CHARACTER(LEN=256) ::   clinfo           ! info character
@@ -76,8 +76,15 @@
       INTEGER            ::   ihdf5            ! local variable for retrieval of value for NF90_HDF5
       LOGICAL            ::   llclobber        ! local definition of ln_clobber
-      !---------------------------------------------------------------------
-
+      INTEGER            ::   ilevels           ! vertical levels
+      !---------------------------------------------------------------------
+      !
       clinfo = '                    iom_nf90_open ~~~  '
-      istop = nstop   ! store the actual value of nstop
+      istop = nstop     ! store the actual value of nstop
+      !
+      !                 !number of vertical levels
+      IF( PRESENT(kdlev) ) THEN   ;   ilevels = kdlev    ! use input value (useful for sea-ice)
+      ELSE                        ;   ilevels = jpk      ! by default jpk
+      ENDIF
+      !
       IF( nn_chunksz > 0 ) THEN   ;   ichunk = nn_chunksz
       ELSE                        ;   ichunk = NF90_SIZEHINT_DEFAULT
@@ -85,6 +92,6 @@
       !
       llclobber = ldwrt .AND. ln_clobber
-      IF( ldok .AND. .NOT. llclobber ) THEN      ! Open existing file...
-         !                 ! =============
+      IF( ldok .AND. .NOT. llclobber ) THEN      !==  Open existing file ==!
+         !                                       !=========================!
          IF( ldwrt ) THEN  ! ... in write mode
             IF(lwp) WRITE(numout,*) TRIM(clinfo)//' open existing file: '//TRIM(cdname)//' in WRITE mode'
@@ -99,8 +106,8 @@
             CALL iom_nf90_check(NF90_OPEN( TRIM(cdname), NF90_NOWRITE, if90id, chunksize = ichunk ), clinfo)
          ENDIF
-      ELSE                                       ! the file does not exist (or we overwrite it)
-         !                 ! =============
+      ELSE                                       !== the file doesn't exist ==!   (or we overwrite it)
+         !                                       !============================!
          iln = INDEX( cdname, '.nc' )
-         IF( ldwrt ) THEN  ! the file should be open in write mode so we create it...
+         IF( ldwrt ) THEN              !* the file should be open in write mode so we create it...
             IF( jpnij > 1 ) THEN
                WRITE(cltmp,'(a,a,i4.4,a)') cdname(1:iln-1), '_', narea-1, '.nc'
@@ -126,5 +133,5 @@
             CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'x', kdompar(1,1)  , idmy ), clinfo)
             CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'y', kdompar(2,1)  , idmy ), clinfo)
-            CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'z', jpk           , idmy ), clinfo)
+            CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'z', ilevels       , idmy ), clinfo)
             CALL iom_nf90_check(NF90_DEF_DIM( if90id, 't', NF90_UNLIMITED, idmy ), clinfo)
             ! global attributes
@@ -139,8 +146,9 @@
             CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_halo_size_end'  , kdompar(:,5)       ), clinfo)
             CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_type'           , 'BOX'              ), clinfo)
-         ELSE              ! the file should be open for read mode so it must exist...
+         ELSE                          !* the file should be open for read mode so it must exist...
             CALL ctl_stop( TRIM(clinfo), ' should be impossible case...' )
          ENDIF
       ENDIF
+      !
       ! start to fill file informations
       ! =============
@@ -156,9 +164,10 @@
          iom_file(kiomid)%nvars  = 0
          iom_file(kiomid)%irec   = -1   ! useless for NetCDF files, used to know if the file is in define mode 
+         iom_file(kiomid)%nlev   = ilevels
          CALL iom_nf90_check(NF90_Inquire(if90id, unlimitedDimId = iom_file(kiomid)%iduld), clinfo)
-         IF ( iom_file(kiomid)%iduld .GE. 0 ) THEN
-           CALL iom_nf90_check(NF90_Inquire_Dimension(if90id, iom_file(kiomid)%iduld,     & 
-        &                                               name = iom_file(kiomid)%uldname,  &
-        &                                               len  = iom_file(kiomid)%lenuld ), clinfo )
+         IF( iom_file(kiomid)%iduld .GE. 0 ) THEN
+            CALL iom_nf90_check(NF90_Inquire_Dimension(if90id, iom_file(kiomid)%iduld,    & 
+               &                                       name = iom_file(kiomid)%uldname,   &
+               &                                       len  = iom_file(kiomid)%lenuld ), clinfo )
          ENDIF
          IF(lwp) WRITE(numout,*) '                   ---> '//TRIM(cdname)//' OK'
@@ -179,8 +188,6 @@
       CHARACTER(LEN=100)  ::   clinfo   ! info character
       !---------------------------------------------------------------------
-      !
       clinfo = '      iom_nf90_close    , file: '//TRIM(iom_file(kiomid)%name)
       CALL iom_nf90_check(NF90_CLOSE(iom_file(kiomid)%nfid), clinfo)
-      !    
    END SUBROUTINE iom_nf90_close
 
@@ -275,5 +282,4 @@
       clinfo = 'iom_nf90_g0d , file: '//TRIM(iom_file(kiomid)%name)//', var: '//TRIM(iom_file(kiomid)%cn_var(kvid))
       CALL iom_nf90_check(NF90_GET_VAR(iom_file(kiomid)%nfid, iom_file(kiomid)%nvid(kvid), pvar, start = kstart), clinfo )
-      ! 
    END SUBROUTINE iom_nf90_g0d
 
@@ -357,5 +363,5 @@
          ivarid = NF90_GLOBAL
       ENDIF
-!
+      !
       IF( llok) THEN
          clinfo = 'iom_nf90_getatt, file: '//TRIM(iom_file(kiomid)%name)//', giatt: '//TRIM(cdatt)
@@ -368,5 +374,6 @@
    END SUBROUTINE iom_nf90_giatt
 
-   SUBROUTINE iom_nf90_gratt( kiomid, cdatt, pv_r0d, cdvar)
+
+   SUBROUTINE iom_nf90_gratt( kiomid, cdatt, pv_r0d, cdvar )
       !!-----------------------------------------------------------------------
       !!                  ***  ROUTINE  iom_nf90_gratt  ***
@@ -376,14 +383,13 @@
       !!               attribute if optional variable name is supplied (cdvar))
       !!-----------------------------------------------------------------------
-      INTEGER         , INTENT(in   ) ::   kiomid   ! Identifier of the file
-      CHARACTER(len=*), INTENT(in   ) ::   cdatt    ! attribute name
-      REAL(wp)        , INTENT(  out) ::   pv_r0d   ! read field
-      CHARACTER(len=*), INTENT(in   ), OPTIONAL     &
-                      &               ::   cdvar    ! name of the variable
-      !
-      INTEGER                         ::   if90id   ! temporary integer
-      INTEGER                         ::   ivarid   ! NetCDF variable Id
-      LOGICAL                         ::   llok     ! temporary logical
-      CHARACTER(LEN=100)              ::   clinfo   ! info character
+      INTEGER                   , INTENT(in   ) ::   kiomid   ! Identifier of the file
+      CHARACTER(len=*)          , INTENT(in   ) ::   cdatt    ! attribute name
+      REAL(wp)                  , INTENT(  out) ::   pv_r0d   ! read field
+      CHARACTER(len=*), OPTIONAL, INTENT(in   ) ::   cdvar    ! name of the variable
+      !
+      INTEGER            ::   if90id   ! temporary integer
+      INTEGER            ::   ivarid   ! NetCDF variable Id
+      LOGICAL            ::   llok     ! temporary logical
+      CHARACTER(LEN=100) ::   clinfo   ! info character
       !---------------------------------------------------------------------
       !
@@ -402,5 +408,5 @@
          ivarid = NF90_GLOBAL
       ENDIF
-!
+      !
       IF( llok) THEN
          clinfo = 'iom_nf90_getatt, file: '//TRIM(iom_file(kiomid)%name)//', gratt: '//TRIM(cdatt)
@@ -413,5 +419,6 @@
    END SUBROUTINE iom_nf90_gratt
 
-   SUBROUTINE iom_nf90_gcatt( kiomid, cdatt, pv_c0d, cdvar)
+
+   SUBROUTINE iom_nf90_gcatt( kiomid, cdatt, pv_c0d, cdvar )
       !!-----------------------------------------------------------------------
       !!                  ***  ROUTINE  iom_nf90_gcatt  ***
@@ -421,14 +428,13 @@
       !!               attribute if optional variable name is supplied (cdvar))
       !!-----------------------------------------------------------------------
-      INTEGER         , INTENT(in   ) ::   kiomid   ! Identifier of the file
-      CHARACTER(len=*), INTENT(in   ) ::   cdatt    ! attribute name
-      CHARACTER(len=*), INTENT(  out) ::   pv_c0d   ! read field
-      CHARACTER(len=*), INTENT(in   ), OPTIONAL     &
-                      &               ::   cdvar    ! name of the variable
-      !
-      INTEGER                         ::   if90id   ! temporary integer
-      INTEGER                         ::   ivarid   ! NetCDF variable Id
-      LOGICAL                         ::   llok     ! temporary logical
-      CHARACTER(LEN=100)              ::   clinfo   ! info character
+      INTEGER                   , INTENT(in   ) ::   kiomid   ! Identifier of the file
+      CHARACTER(len=*)          , INTENT(in   ) ::   cdatt    ! attribute name
+      CHARACTER(len=*)          , INTENT(  out) ::   pv_c0d   ! read field
+      CHARACTER(len=*), OPTIONAL, INTENT(in   ) ::   cdvar    ! name of the variable
+      !
+      INTEGER            ::   if90id   ! temporary integer
+      INTEGER            ::   ivarid   ! NetCDF variable Id
+      LOGICAL            ::   llok     ! temporary logical
+      CHARACTER(LEN=100) ::   clinfo   ! info character
       !---------------------------------------------------------------------
       !
@@ -458,4 +464,5 @@
    END SUBROUTINE iom_nf90_gcatt
 
+
    !!----------------------------------------------------------------------
    !!                   INTERFACE iom_nf90_putatt
@@ -495,5 +502,5 @@
          ivarid = NF90_GLOBAL
       ENDIF
-!
+      !
       IF( llok) THEN
          clinfo = 'iom_nf90_putatt, file: '//TRIM(iom_file(kiomid)%name)//', piatt: '//TRIM(cdatt)
@@ -517,5 +524,6 @@
    END SUBROUTINE iom_nf90_piatt
 
-   SUBROUTINE iom_nf90_pratt( kiomid, cdatt, pv_r0d, cdvar)
+
+   SUBROUTINE iom_nf90_pratt( kiomid, cdatt, pv_r0d, cdvar )
       !!-----------------------------------------------------------------------
       !!                  ***  ROUTINE  iom_nf90_pratt  ***
@@ -525,15 +533,14 @@
       !!               attribute if optional variable name is supplied (cdvar))
       !!-----------------------------------------------------------------------
-      INTEGER         , INTENT(in   ) ::   kiomid   ! Identifier of the file
-      CHARACTER(len=*), INTENT(in   ) ::   cdatt    ! attribute name
-      REAL(wp)        , INTENT(in   ) ::   pv_r0d   ! write field
-      CHARACTER(len=*), INTENT(in   ), OPTIONAL     &
-                      &               ::   cdvar    ! name of the variable
-      !
-      INTEGER                         ::   if90id   ! temporary integer
-      INTEGER                         ::   ivarid   ! NetCDF variable Id
-      LOGICAL                         ::   llok     ! temporary logical
-      LOGICAL                         ::   lenddef  ! temporary logical
-      CHARACTER(LEN=100)              ::   clinfo   ! info character
+      INTEGER                   , INTENT(in   ) ::   kiomid   ! Identifier of the file
+      CHARACTER(len=*)          , INTENT(in   ) ::   cdatt    ! attribute name
+      REAL(wp)                  , INTENT(in   ) ::   pv_r0d   ! write field
+      CHARACTER(len=*), OPTIONAL, INTENT(in   ) ::   cdvar    ! name of the variable
+      !
+      INTEGER            ::   if90id   ! temporary integer
+      INTEGER            ::   ivarid   ! NetCDF variable Id
+      LOGICAL            ::   llok     ! temporary logical
+      LOGICAL            ::   lenddef  ! temporary logical
+      CHARACTER(LEN=100) ::   clinfo   ! info character
       !---------------------------------------------------------------------
       !
@@ -550,5 +557,5 @@
          ivarid = NF90_GLOBAL
       ENDIF
-!
+      !
       IF( llok) THEN
          clinfo = 'iom_nf90_putatt, file: '//TRIM(iom_file(kiomid)%name)//', pratt: '//TRIM(cdatt)
@@ -572,5 +579,6 @@
    END SUBROUTINE iom_nf90_pratt
 
-   SUBROUTINE iom_nf90_pcatt( kiomid, cdatt, pv_c0d, cdvar)
+
+   SUBROUTINE iom_nf90_pcatt( kiomid, cdatt, pv_c0d, cdvar )
       !!-----------------------------------------------------------------------
       !!                  ***  ROUTINE  iom_nf90_pcatt  ***
@@ -580,15 +588,14 @@
       !!               attribute if optional variable name is supplied (cdvar))
       !!-----------------------------------------------------------------------
-      INTEGER         , INTENT(in   ) ::   kiomid   ! Identifier of the file
-      CHARACTER(len=*), INTENT(in   ) ::   cdatt    ! attribute name
-      CHARACTER(len=*), INTENT(in   ) ::   pv_c0d   ! write field
-      CHARACTER(len=*), INTENT(in   ), OPTIONAL     &
-                      &               ::   cdvar    ! name of the variable
-      !
-      INTEGER                         ::   if90id   ! temporary integer
-      INTEGER                         ::   ivarid   ! NetCDF variable Id
-      LOGICAL                         ::   llok     ! temporary logical
-      LOGICAL                         ::   lenddef  ! temporary logical
-      CHARACTER(LEN=100)              ::   clinfo   ! info character
+      INTEGER                   , INTENT(in   ) ::   kiomid   ! Identifier of the file
+      CHARACTER(len=*)          , INTENT(in   ) ::   cdatt    ! attribute name
+      CHARACTER(len=*)          , INTENT(in   ) ::   pv_c0d   ! write field
+      CHARACTER(len=*), OPTIONAL, INTENT(in   ) ::   cdvar    ! name of the variable
+      !
+      INTEGER            ::   if90id   ! temporary integer
+      INTEGER            ::   ivarid   ! NetCDF variable Id
+      LOGICAL            ::   llok     ! temporary logical
+      LOGICAL            ::   lenddef  ! temporary logical
+      CHARACTER(LEN=100) ::   clinfo   ! info character
       !---------------------------------------------------------------------
       !
@@ -605,5 +612,5 @@
          ivarid = NF90_GLOBAL
       ENDIF
-!
+      !
       IF( llok) THEN
          clinfo = 'iom_nf90_putatt, file: '//TRIM(iom_file(kiomid)%name)//', pcatt: '//TRIM(cdatt)
@@ -658,5 +665,5 @@
 
    SUBROUTINE iom_nf90_rp0123d( kt, kwrite, kiomid, cdvar , kvid  , ktype,   &
-         &                               pv_r0d, pv_r1d, pv_r2d, pv_r3d )
+         &                                  pv_r0d, pv_r1d, pv_r2d, pv_r3d )
       !!--------------------------------------------------------------------
       !!                   ***  SUBROUTINE  iom_nf90_rstput  ***
@@ -687,10 +694,10 @@
       INTEGER               :: itype                ! variable type
       INTEGER, DIMENSION(4) :: ichunksz             ! NetCDF4 chunk sizes. Will be computed using
-                                                    ! nn_nchunks_[i,j,k,t] namelist parameters
-      INTEGER               :: ichunkalg, ishuffle,&
-                               ideflate, ideflate_level
-                                                    ! NetCDF4 internally fixed parameters
+      !                                             ! nn_nchunks_[i,j,k,t] namelist parameters
+      INTEGER               :: ichunkalg, ishuffle, ideflate, ideflate_level
+      !                                             ! NetCDF4 internally fixed parameters
       LOGICAL               :: lchunk               ! logical switch to activate chunking and compression
-                                                    ! when appropriate (currently chunking is applied to 4d fields only)
+      !                                             ! when appropriate (currently chunking is applied to 4d fields only)
+      INTEGER               :: idlv                 ! local variable
       !---------------------------------------------------------------------
       !
@@ -706,5 +713,9 @@
          ENDIF
          ! define the dimension variables if it is not already done
-         cltmp = (/ 'nav_lon     ', 'nav_lat     ', 'nav_lev     ', 'time_counter' /)
+         IF(iom_file(kiomid)%nlev == jpk ) THEN
+          cltmp = (/ 'nav_lon     ', 'nav_lat     ', 'nav_lev     ', 'time_counter' /)
+         ELSE
+          cltmp = (/ 'nav_lon     ', 'nav_lat     ', 'numcat      ', 'time_counter' /)
+         ENDIF
          CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(1)), NF90_FLOAT , (/ 1, 2 /), iom_file(kiomid)%nvid(1) ), clinfo)
          CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(2)), NF90_FLOAT , (/ 1, 2 /), iom_file(kiomid)%nvid(2) ), clinfo)
@@ -755,11 +766,11 @@
          IF( PRESENT(pv_r0d) ) THEN
             CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype,                    &
-                 &                            iom_file(kiomid)%nvid(idvar) ), clinfo)
+               &                              iom_file(kiomid)%nvid(idvar) ), clinfo )
          ELSE
             CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, idimid(1:idims),   &
-                 &                            iom_file(kiomid)%nvid(idvar) ), clinfo)
+               &                              iom_file(kiomid)%nvid(idvar) ), clinfo )
          ENDIF
          lchunk = .false.
-         IF( snc4set%luse .AND. idims.eq.4 ) lchunk = .true.
+         IF( snc4set%luse .AND. idims == 4 )   lchunk = .true.
          ! update informations structure related the new variable we want to add...
          iom_file(kiomid)%nvars         = idvar
@@ -782,5 +793,5 @@
             ichunksz(3) = MIN( ichunksz(3),MAX( (ichunksz(3)-1)/snc4set%nk + 1 , 1 ) ) ! Suggested default nc4set%nk=6
             ichunksz(4) = 1                                                            ! Do not allow chunks to span the
-                                                                                       ! unlimited dimension
+            !                                                                          ! unlimited dimension
             CALL iom_nf90_check(SET_NF90_DEF_VAR_CHUNKING(if90id, idvar, ichunkalg, ichunksz), clinfo)
             CALL iom_nf90_check(SET_NF90_DEF_VAR_DEFLATE(if90id, idvar, ishuffle, ideflate, ideflate_level), clinfo)
@@ -791,5 +802,5 @@
          idvar = kvid
       ENDIF
-
+      !
       ! time step kwrite : write the variable
       IF( kt == kwrite ) THEN
@@ -815,18 +826,24 @@
             ! trick: is defined to 0 => dimension variable are defined but not yet written
             IF( iom_file(kiomid)%dimsz(1, 1) == 0 ) THEN
-               CALL iom_nf90_check(NF90_INQ_VARID( if90id, 'nav_lon'     , idmy ), clinfo)
-               CALL iom_nf90_check(NF90_PUT_VAR( if90id, idmy, glamt(ix1:ix2, iy1:iy2) ), clinfo)
-               CALL iom_nf90_check(NF90_INQ_VARID( if90id, 'nav_lat'     , idmy ), clinfo)
-               CALL iom_nf90_check(NF90_PUT_VAR( if90id, idmy, gphit(ix1:ix2, iy1:iy2) ), clinfo)
-               CALL iom_nf90_check(NF90_INQ_VARID( if90id, 'nav_lev'     , idmy ), clinfo)
-               CALL iom_nf90_check(NF90_PUT_VAR( if90id, idmy, gdept_1d                ), clinfo)
+               CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lon'     , idmy )         , clinfo )
+               CALL iom_nf90_check( NF90_PUT_VAR  ( if90id, idmy, glamt(ix1:ix2, iy1:iy2) ), clinfo )
+               CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lat'     , idmy )         , clinfo )
+               CALL iom_nf90_check( NF90_PUT_VAR  ( if90id, idmy, gphit(ix1:ix2, iy1:iy2) ), clinfo )
+               IF(iom_file(kiomid)%nlev == jpk ) THEN 
+                  !NEMO
+                  CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lev'     , idmy ), clinfo )
+                  CALL iom_nf90_check( NF90_PUT_VAR  ( if90id, idmy, gdept_1d       ), clinfo )
+               ELSE
+                  CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'numcat'     , idmy ), clinfo)
+                  CALL iom_nf90_check( NF90_PUT_VAR  ( if90id, idmy, (/ (idlv, idlv = 1,iom_file(kiomid)%nlev) /)), clinfo )
+               ENDIF
                ! +++ WRONG VALUE: to be improved but not really useful...
-               CALL iom_nf90_check(NF90_INQ_VARID( if90id, 'time_counter', idmy ), clinfo)
-               CALL iom_nf90_check(NF90_PUT_VAR( if90id, idmy, kt                      ), clinfo)   
+               CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'time_counter', idmy ), clinfo )
+               CALL iom_nf90_check( NF90_PUT_VAR( if90id, idmy, kt                      ), clinfo )   
                ! update the values of the variables dimensions size
-               CALL iom_nf90_check(NF90_INQUIRE_DIMENSION( if90id, 1, len = iom_file(kiomid)%dimsz(1,1) ), clinfo)
-               CALL iom_nf90_check(NF90_INQUIRE_DIMENSION( if90id, 2, len = iom_file(kiomid)%dimsz(2,1) ), clinfo)
+               CALL iom_nf90_check( NF90_INQUIRE_DIMENSION( if90id, 1, len = iom_file(kiomid)%dimsz(1,1) ), clinfo )
+               CALL iom_nf90_check( NF90_INQUIRE_DIMENSION( if90id, 2, len = iom_file(kiomid)%dimsz(2,1) ), clinfo )
                iom_file(kiomid)%dimsz(1:2, 2) = iom_file(kiomid)%dimsz(1:2, 1)
-               CALL iom_nf90_check(NF90_INQUIRE_DIMENSION( if90id, 3, len = iom_file(kiomid)%dimsz(1,3) ), clinfo)
+               CALL iom_nf90_check( NF90_INQUIRE_DIMENSION( if90id, 3, len = iom_file(kiomid)%dimsz(1,3) ), clinfo )
                iom_file(kiomid)%dimsz(1  , 4) = 1   ! unlimited dimension
                IF(lwp) WRITE(numout,*) TRIM(clinfo)//' write dimension variables done'
@@ -837,11 +854,11 @@
          ! =============
          IF(     PRESENT(pv_r0d) ) THEN
-            CALL iom_nf90_check(NF90_PUT_VAR( if90id, idvar, pv_r0d                      ), clinfo)
+            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r0d                    ), clinfo )
          ELSEIF( PRESENT(pv_r1d) ) THEN
-            CALL iom_nf90_check(NF90_PUT_VAR( if90id, idvar, pv_r1d(                  :) ), clinfo)
+            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r1d(:)                 ), clinfo )
          ELSEIF( PRESENT(pv_r2d) ) THEN
-            CALL iom_nf90_check(NF90_PUT_VAR( if90id, idvar, pv_r2d(ix1:ix2, iy1:iy2   ) ), clinfo)
+            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r2d(ix1:ix2,iy1:iy2)   ), clinfo )
          ELSEIF( PRESENT(pv_r3d) ) THEN
-            CALL iom_nf90_check(NF90_PUT_VAR( if90id, idvar, pv_r3d(ix1:ix2, iy1:iy2, :) ), clinfo)
+            CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r3d(ix1:ix2,iy1:iy2,:) ), clinfo )
          ENDIF
          ! add 1 to the size of the temporal dimension (not really useful...)
