Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/IOM/iom.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/IOM/iom.F90	(revision 13998)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/IOM/iom.F90	(revision 14018)
@@ -46,4 +46,6 @@
    USE lib_fortran 
    USE diu_bulk, ONLY : ln_diurnal_only, ln_diurnal
+   USE iom_nf90
+   USE netcdf
 
    IMPLICIT NONE
@@ -58,4 +60,5 @@
    PUBLIC iom_chkatt, iom_getatt, iom_putatt, iom_getszuld, iom_rstput, iom_delay_rst, iom_put
    PUBLIC iom_use, iom_context_finalize, iom_update_file_name, iom_miss_val
+   PUBLIC iom_xios_setid
 
    PRIVATE iom_rp0d_sp, iom_rp1d_sp, iom_rp2d_sp, iom_rp3d_sp
@@ -69,7 +72,8 @@
    PRIVATE iom_set_domain_attr, iom_set_axis_attr, iom_set_field_attr, iom_set_file_attr, iom_get_file_attr, iom_set_grid_attr
    PRIVATE set_grid, set_grid_bounds, set_scalar, set_xmlatt, set_mooring, iom_sdate
-   PRIVATE iom_set_rst_context, iom_set_rstw_active, iom_set_rstr_active
+   PRIVATE iom_set_rst_context, iom_set_vars_active
 # endif
-   PUBLIC iom_set_rstw_var_active, iom_set_rstw_core, iom_set_rst_vars
+   PRIVATE set_xios_context
+   PRIVATE iom_set_rstw_active
 
    INTERFACE iom_get
@@ -101,5 +105,5 @@
 CONTAINS
 
-   SUBROUTINE iom_init( cdname, fname, ld_closedef ) 
+   SUBROUTINE iom_init( cdname, kdid, ld_closedef ) 
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE   ***
@@ -109,5 +113,5 @@
       !!----------------------------------------------------------------------
       CHARACTER(len=*),           INTENT(in)  :: cdname
-      CHARACTER(len=*), OPTIONAL, INTENT(in)  :: fname
+      INTEGER         , OPTIONAL, INTENT(in)  :: kdid          
       LOGICAL         , OPTIONAL, INTENT(in)  :: ld_closedef
 #if defined key_iomput
@@ -118,12 +122,15 @@
       INTEGER             :: irefyear, irefmonth, irefday
       INTEGER           :: ji
-      LOGICAL :: llrst_context              ! is context related to restart
+      LOGICAL           :: llrst_context              ! is context related to restart
+      LOGICAL           :: llrstr, llrstw 
+      INTEGER           :: inum
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zt_bnds, zw_bnds
       REAL(wp), DIMENSION(2,jpkam1)         :: za_bnds   ! ABL vertical boundaries
-      LOGICAL ::   ll_closedef = .TRUE.
+      LOGICAL ::   ll_closedef
       LOGICAL ::   ll_exist
       !!----------------------------------------------------------------------
       !
+      ll_closedef = .TRUE.
       IF ( PRESENT(ld_closedef) ) ll_closedef = ld_closedef
       !
@@ -134,5 +141,14 @@
       CALL xios_context_initialize(TRIM(clname), mpi_comm_oce)
       CALL iom_swap( cdname )
-      llrst_context =  (TRIM(cdname) == TRIM(crxios_context) .OR. TRIM(cdname) == TRIM(cwxios_context))
+
+      llrstr = (cdname == cr_ocerst_cxt) .OR. (cdname == cr_icerst_cxt)
+      llrstr = llrstr .OR. (cdname == cr_toprst_cxt)
+      llrstr = llrstr .OR. (cdname == cr_sedrst_cxt)
+
+      llrstw = (cdname == cw_ocerst_cxt) .OR. (cdname == cw_icerst_cxt)
+      llrstw = llrstw .OR. (cdname == cw_toprst_cxt)
+      llrstw = llrstw .OR. (cdname == cw_sedrst_cxt)
+
+      llrst_context = llrstr .OR. llrstw
 
       ! Calendar type is now defined in xml file 
@@ -153,5 +169,5 @@
       IF(.NOT.llrst_context) CALL set_scalar
       !
-      IF( TRIM(cdname) == TRIM(cxios_context) ) THEN  
+      IF( cdname == cxios_context ) THEN  
          CALL set_grid( "T", glamt, gphit, .FALSE., .FALSE. ) 
          CALL set_grid( "U", glamu, gphiu, .FALSE., .FALSE. )
@@ -200,74 +216,72 @@
       ! vertical grid definition
       IF(.NOT.llrst_context) THEN
-          CALL iom_set_axis_attr(  "deptht", paxis = gdept_1d )
-          CALL iom_set_axis_attr(  "depthu", paxis = gdept_1d )
-          CALL iom_set_axis_attr(  "depthv", paxis = gdept_1d )
-          CALL iom_set_axis_attr(  "depthw", paxis = gdepw_1d )
+         CALL iom_set_axis_attr(  "deptht", paxis = gdept_1d )
+         CALL iom_set_axis_attr(  "depthu", paxis = gdept_1d )
+         CALL iom_set_axis_attr(  "depthv", paxis = gdept_1d )
+         CALL iom_set_axis_attr(  "depthw", paxis = gdepw_1d )
           CALL iom_set_axis_attr(  "depthf", paxis = gdept_1d )
 
           ! ABL
-          IF( .NOT. ALLOCATED(ght_abl) ) THEN   ! force definition for xml files (xios) 
-             ALLOCATE( ght_abl(jpka), ghw_abl(jpka), e3t_abl(jpka), e3w_abl(jpka) )   ! default allocation needed by iom
-             ght_abl(:) = -1._wp   ;   ghw_abl(:) = -1._wp
-             e3t_abl(:) = -1._wp   ;   e3w_abl(:) = -1._wp
-          ENDIF
-          CALL iom_set_axis_attr( "ght_abl", ght_abl(2:jpka) )
-          CALL iom_set_axis_attr( "ghw_abl", ghw_abl(2:jpka) )
+         IF( .NOT. ALLOCATED(ght_abl) ) THEN   ! force definition for xml files (xios) 
+            ALLOCATE( ght_abl(jpka), ghw_abl(jpka), e3t_abl(jpka), e3w_abl(jpka) )   ! default allocation needed by iom
+            ght_abl(:) = -1._wp   ;   ghw_abl(:) = -1._wp
+            e3t_abl(:) = -1._wp   ;   e3w_abl(:) = -1._wp
+         ENDIF
+         CALL iom_set_axis_attr( "ght_abl", ght_abl(2:jpka) )
+         CALL iom_set_axis_attr( "ghw_abl", ghw_abl(2:jpka) )
           
-          ! Add vertical grid bounds
-          zt_bnds(2,:      ) = gdept_1d(:)
-          zt_bnds(1,2:jpk  ) = gdept_1d(1:jpkm1)
-          zt_bnds(1,1      ) = gdept_1d(1) - e3w_1d(1)
-          zw_bnds(1,:      ) = gdepw_1d(:)
-          zw_bnds(2,1:jpkm1) = gdepw_1d(2:jpk)
-          zw_bnds(2,jpk:   ) = gdepw_1d(jpk) + e3t_1d(jpk)
-          CALL iom_set_axis_attr(  "deptht", bounds=zw_bnds )
-          CALL iom_set_axis_attr(  "depthu", bounds=zw_bnds )
-          CALL iom_set_axis_attr(  "depthv", bounds=zw_bnds )
-          CALL iom_set_axis_attr(  "depthw", bounds=zt_bnds )
+         ! Add vertical grid bounds
+         zt_bnds(2,:      ) = gdept_1d(:)
+         zt_bnds(1,2:jpk  ) = gdept_1d(1:jpkm1)
+         zt_bnds(1,1      ) = gdept_1d(1) - e3w_1d(1)
+         zw_bnds(1,:      ) = gdepw_1d(:)
+         zw_bnds(2,1:jpkm1) = gdepw_1d(2:jpk)
+         zw_bnds(2,jpk:   ) = gdepw_1d(jpk) + e3t_1d(jpk)
+         CALL iom_set_axis_attr(  "deptht", bounds=zw_bnds )
+         CALL iom_set_axis_attr(  "depthu", bounds=zw_bnds )
+         CALL iom_set_axis_attr(  "depthv", bounds=zw_bnds )
+         CALL iom_set_axis_attr(  "depthw", bounds=zt_bnds )
           CALL iom_set_axis_attr(  "depthf", bounds=zw_bnds )
 
-          ! ABL
-          za_bnds(1,:) = ghw_abl(1:jpkam1)
-          za_bnds(2,:) = ghw_abl(2:jpka  )
-          CALL iom_set_axis_attr( "ght_abl", bounds=za_bnds )
-          za_bnds(1,:) = ght_abl(2:jpka  )
-          za_bnds(2,:) = ght_abl(2:jpka  ) + e3w_abl(2:jpka)
-          CALL iom_set_axis_attr( "ghw_abl", bounds=za_bnds )
-
-          CALL iom_set_axis_attr(  "nfloat", (/ (REAL(ji,wp), ji=1,jpnfl) /) )
+         ! ABL
+         za_bnds(1,:) = ghw_abl(1:jpkam1)
+         za_bnds(2,:) = ghw_abl(2:jpka  )
+         CALL iom_set_axis_attr( "ght_abl", bounds=za_bnds )
+         za_bnds(1,:) = ght_abl(2:jpka  )
+         za_bnds(2,:) = ght_abl(2:jpka  ) + e3w_abl(2:jpka)
+         CALL iom_set_axis_attr( "ghw_abl", bounds=za_bnds )
+
+         CALL iom_set_axis_attr(  "nfloat", (/ (REAL(ji,wp), ji=1,jpnfl) /) )
 # if defined key_si3
-          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) /) )
+         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
 #if defined key_top
-          IF( ALLOCATED(profsed) ) CALL iom_set_axis_attr( "profsed", paxis = profsed )
-#endif
-          CALL iom_set_axis_attr( "icbcla", class_num )
-          CALL iom_set_axis_attr( "iax_20C", (/ REAL(20,wp) /) )   ! strange syntaxe and idea...
-          CALL iom_set_axis_attr( "iax_26C", (/ REAL(26,wp) /) )   ! strange syntaxe and idea...
-          CALL iom_set_axis_attr( "iax_28C", (/ REAL(28,wp) /) )   ! strange syntaxe and idea...
-          ! for diaprt, we need to define an axis which size can be 1 (default) or 5 (if the file subbasins.nc exists)
-          INQUIRE( FILE = 'subbasins.nc', EXIST = ll_exist )
-          nbasin = 1 + 4 * COUNT( (/ll_exist/) )
-          CALL iom_set_axis_attr( "basin"  , (/ (REAL(ji,wp), ji=1,nbasin) /) )
+         IF( ALLOCATED(profsed) ) CALL iom_set_axis_attr( "profsed", paxis = profsed )
+#endif
+         CALL iom_set_axis_attr( "icbcla", class_num )
+         CALL iom_set_axis_attr( "iax_20C", (/ REAL(20,wp) /) )   ! strange syntaxe and idea...
+         CALL iom_set_axis_attr( "iax_26C", (/ REAL(26,wp) /) )   ! strange syntaxe and idea...
+         CALL iom_set_axis_attr( "iax_28C", (/ REAL(28,wp) /) )   ! strange syntaxe and idea...
+         ! for diaprt, we need to define an axis which size can be 1 (default) or 5 (if the file subbasins.nc exists)
+         INQUIRE( FILE = 'subbasins.nc', EXIST = ll_exist )
+         nbasin = 1 + 4 * COUNT( (/ll_exist/) )
+         CALL iom_set_axis_attr( "basin"  , (/ (REAL(ji,wp), ji=1,nbasin) /) )
       ENDIF
       !
       ! automatic definitions of some of the xml attributs
-      IF( TRIM(cdname) == TRIM(crxios_context) ) THEN
-!set names of the fields in restart file IF using XIOS to read data
-          CALL iom_set_rst_context(.TRUE.)
-          CALL iom_set_rst_vars(rst_rfields)
-!set which fields are to be read from restart file
-          CALL iom_set_rstr_active()
-      ELSE IF( TRIM(cdname) == TRIM(cwxios_context) ) THEN
-!set names of the fields in restart file IF using XIOS to write data
-          CALL iom_set_rst_context(.FALSE.)
-          CALL iom_set_rst_vars(rst_wfields)
-!set which fields are to be written to a restart file
-          CALL iom_set_rstw_active(fname)
+      IF(llrstr) THEN
+         IF(PRESENT(kdid)) THEN
+            CALL iom_set_rst_context(.TRUE.)
+!set which fields will be read from restart file
+            CALL iom_set_vars_active(kdid)
+         ELSE
+            CALL ctl_stop( 'iom_init:', 'restart read with XIOS: missing pointer to NETCDF file' )
+         ENDIF
+      ELSE IF(llrstw) THEN
+         CALL iom_set_rstw_file(iom_file(kdid)%name)
       ELSE
-          CALL set_xmlatt
+         CALL set_xmlatt
       ENDIF
       !
@@ -285,5 +299,5 @@
    END SUBROUTINE iom_init
 
-   SUBROUTINE iom_init_closedef
+   SUBROUTINE iom_init_closedef(cdname)
       !!----------------------------------------------------------------------
       !!            ***  SUBROUTINE iom_init_closedef  ***
@@ -293,8 +307,24 @@
       !!
       !!----------------------------------------------------------------------
-
+      CHARACTER(len=*), OPTIONAL, INTENT(IN) :: cdname
 #if defined key_iomput
-      CALL xios_close_context_definition()
-      CALL xios_update_calendar( 0 )
+      LOGICAL :: llrstw
+
+      llrstw = .FALSE.
+      IF(PRESENT(cdname)) THEN
+         llrstw = (cdname == cw_ocerst_cxt)
+         llrstw = llrstw .OR. (cdname == cw_icerst_cxt)
+         llrstw = llrstw .OR. (cdname == cw_toprst_cxt)
+         llrstw = llrstw .OR. (cdname == cw_sedrst_cxt)
+      ENDIF
+
+      IF( llrstw ) THEN
+!set names of the fields in restart file IF using XIOS to write data
+         CALL iom_set_rst_context(.FALSE.)
+         CALL xios_close_context_definition()
+      ELSE
+         CALL xios_close_context_definition()
+         CALL xios_update_calendar( 0 )
+      ENDIF
 #else
       IF( .FALSE. )   WRITE(numout,*) 'iom_init_closedef: should not see this'   ! useless statement to avoid compilation warnings
@@ -303,252 +333,143 @@
    END SUBROUTINE iom_init_closedef
 
-   SUBROUTINE iom_set_rstw_var_active(field)
+   SUBROUTINE iom_set_vars_active(idnum)
       !!---------------------------------------------------------------------
-      !!                   ***  SUBROUTINE  iom_set_rstw_var_active  ***
-      !!
-      !! ** Purpose :  enable variable in restart file when writing with XIOS 
+      !!                   ***  SUBROUTINE  iom_set_vars_active  ***
+      !!
+      !! ** Purpose :  define filename in XIOS context for reading file,
+      !!               enable variables present in a file for reading with XIOS 
+      !!               id of the file is assumed to be rrestart.
       !!---------------------------------------------------------------------
-   CHARACTER(len = *), INTENT(IN) :: field
-   INTEGER :: i
-   LOGICAL :: llis_set
-   CHARACTER(LEN=256) :: clinfo    ! info character
-
+      INTEGER, INTENT(IN) :: idnum 
+      
 #if defined key_iomput
-   llis_set = .FALSE.
-
-   DO i = 1, max_rst_fields
-       IF(TRIM(rst_wfields(i)%vname) == field) THEN 
-          rst_wfields(i)%active = .TRUE.
-          llis_set = .TRUE.
-          EXIT
-       ENDIF
-   ENDDO
-!Warn if variable is not in defined in rst_wfields
-   IF(.NOT.llis_set) THEN
-      WRITE(ctmp1,*) 'iom_set_rstw_var_active: variable ', field ,' is available for writing but not defined' 
-      CALL ctl_stop( 'iom_set_rstw_var_active:', ctmp1 )
-   ENDIF
-#else
-        clinfo = 'iom_set_rstw_var_active: key_iomput is needed to use XIOS restart read/write functionality'
-        CALL ctl_stop('STOP', TRIM(clinfo))
-#endif
-
-   END SUBROUTINE iom_set_rstw_var_active
-
-   SUBROUTINE iom_set_rstr_active()
+      INTEGER                                    :: ndims, nvars, natts, unlimitedDimId, dimlen, xtype,mdims
+      TYPE(xios_field)                           :: field_hdl
+      TYPE(xios_file)                            :: file_hdl
+      TYPE(xios_filegroup)                       :: filegroup_hdl
+      INTEGER                                    :: dimids(4), jv,i, idim
+      CHARACTER(LEN=256)                         :: clinfo               ! info character
+      INTEGER, ALLOCATABLE                       :: indimlens(:)
+      CHARACTER(LEN=nf90_max_name), ALLOCATABLE  :: indimnames(:)
+      CHARACTER(LEN=nf90_max_name)               :: dimname, varname
+      INTEGER                                    :: iln
+      CHARACTER(LEN=lc)                          :: fname
+      LOGICAL                                    :: lmeta
+!metadata in restart file for restart read with XIOS
+      INTEGER, PARAMETER                         :: NMETA = 10
+      CHARACTER(LEN=lc)                          :: meta(NMETA)
+
+
+      meta(1) = "nav_lat"
+      meta(2) = "nav_lon"
+      meta(3) = "nav_lev"
+      meta(4) = "time_instant"
+      meta(5) = "time_instant_bounds"
+      meta(6) = "time_counter"
+      meta(7) = "time_counter_bounds"
+      meta(8) = "x"
+      meta(9) = "y"
+      meta(10) = "numcat"
+
+      clinfo = '          iom_set_vars_active, file: '//TRIM(iom_file(idnum)%name)
+
+      iln = INDEX( iom_file(idnum)%name, '.nc' )
+!XIOS doee not need .nc
+      IF(iln > 0) THEN
+        fname =  iom_file(idnum)%name(1:iln-1)
+      ELSE
+        fname =  iom_file(idnum)%name
+      ENDIF
+
+!set name of the restart file and enable available fields
+      CALL xios_get_handle("file_definition", filegroup_hdl )
+      CALL xios_add_child(filegroup_hdl, file_hdl, 'rrestart')
+      CALL xios_set_file_attr( "rrestart", name=fname, type="one_file",      &
+           par_access="collective", enabled=.TRUE., mode="read",              &
+                                                    output_freq=xios_timestep )
+
+      CALL iom_nf90_check( nf90_inquire(iom_file(idnum)%nfid, ndims, nvars, natts ), clinfo )
+      ALLOCATE(indimlens(ndims), indimnames(ndims))
+      CALL iom_nf90_check( nf90_inquire(iom_file(idnum)%nfid, unlimitedDimId = unlimitedDimId ), clinfo )
+
+      DO idim = 1, ndims
+         CALL iom_nf90_check( nf90_inquire_dimension(iom_file(idnum)%nfid, idim, dimname, dimlen ), clinfo )
+         indimlens(idim) = dimlen
+         indimnames(idim) = dimname
+      ENDDO
+
+      DO jv =1, nvars
+         lmeta = .FALSE.
+         CALL iom_nf90_check( nf90_inquire_variable(iom_file(idnum)%nfid, jv, varname, xtype, ndims, dimids, natts ), clinfo )
+         DO i = 1, NMETA
+           IF(varname == meta(i)) THEN
+             lmeta = .TRUE.
+           ENDIF
+         ENDDO
+         IF(.NOT.lmeta) THEN
+            CALL xios_add_child(file_hdl, field_hdl, varname)
+            mdims = ndims
+
+            IF(ANY(dimids(1:ndims) == unlimitedDimId)) THEN
+               mdims = mdims - 1
+            ENDIF
+
+            IF(mdims == 3) THEN
+               CALL xios_set_attr (field_hdl, enabled = .TRUE., name = varname,   &
+                                   domain_ref="grid_N",                           &
+                                   axis_ref=iom_axis(indimlens(dimids(mdims))),   &
+                                   prec = 8, operation = "instant"                )
+            ELSEIF(mdims == 2) THEN
+               CALL xios_set_attr (field_hdl, enabled = .TRUE., name = varname,  &
+                                   domain_ref="grid_N", prec = 8,                &
+                                   operation = "instant"                         ) 
+            ELSEIF(mdims == 1) THEN
+               CALL xios_set_attr (field_hdl, enabled = .TRUE., name = varname, &
+                                   axis_ref=iom_axis(indimlens(dimids(mdims))), &
+                                   prec = 8, operation = "instant"              )
+            ELSEIF(mdims == 0) THEN
+               CALL xios_set_attr (field_hdl, enabled = .TRUE., name = varname, &
+                                   scalar_ref = "grid_scalar", prec = 8,        &
+                                   operation = "instant"                        )
+            ELSE
+               WRITE(ctmp1,*) 'iom_set_vars_active: variable ', TRIM(varname) ,' incorrect number of dimensions' 
+               CALL ctl_stop( 'iom_set_vars_active:', ctmp1 )
+            ENDIF
+         ENDIF
+      ENDDO
+      DEALLOCATE(indimlens, indimnames)
+#endif
+   END SUBROUTINE iom_set_vars_active
+
+   SUBROUTINE iom_set_rstw_file(cdrst_file)
       !!---------------------------------------------------------------------
-      !!                   ***  SUBROUTINE  iom_set_rstr_active  ***
-      !!
-      !! ** Purpose :  define file name in XIOS context for reading restart file,
-      !!               enable variables present in restart file for reading with XIOS 
+      !!                   ***  SUBROUTINE iom_set_rstw_file   ***
+      !!
+      !! ** Purpose :  define file name in XIOS context for writing restart
       !!---------------------------------------------------------------------
-
-!sets enabled = .TRUE. for each field in restart file
-   CHARACTER(len=256) :: rst_file
-
+      CHARACTER(len=*) :: cdrst_file
 #if defined key_iomput
-   TYPE(xios_field) :: field_hdl
-   TYPE(xios_file) :: file_hdl
-   TYPE(xios_filegroup) :: filegroup_hdl
-   INTEGER :: i
-   CHARACTER(lc)  ::   clpath
-
-        clpath = TRIM(cn_ocerst_indir)
-        IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
-        IF( TRIM(Agrif_CFixed()) == '0' ) THEN
-           rst_file = TRIM(clpath)//TRIM(cn_ocerst_in)
-        ELSE
-           rst_file = TRIM(clpath)//TRIM(Agrif_CFixed())//'_'//TRIM(cn_ocerst_in)
-        ENDIF
+      TYPE(xios_file) :: file_hdl
+      TYPE(xios_filegroup) :: filegroup_hdl
+
 !set name of the restart file and enable available fields
-        if(lwp) WRITE(numout,*) 'Setting restart filename (for XIOS) to: ',rst_file
-        CALL xios_get_handle("file_definition", filegroup_hdl )
-        CALL xios_add_child(filegroup_hdl, file_hdl, 'rrestart')
-        CALL xios_set_file_attr( "rrestart", name=trim(rst_file), type="one_file", &
-             par_access="collective", enabled=.TRUE., mode="read",                 &
-             output_freq=xios_timestep)
-!define variables for restart context
-        DO i = 1, max_rst_fields
-         IF( TRIM(rst_rfields(i)%vname) /= "NO_NAME") THEN
-           IF( iom_varid( numror, TRIM(rst_rfields(i)%vname), ldstop = .FALSE. ) > 0 ) THEN
-                CALL xios_add_child(file_hdl, field_hdl, TRIM(rst_rfields(i)%vname))
-                SELECT CASE (TRIM(rst_rfields(i)%grid))
-                 CASE ("grid_N_3D")
-                    CALL xios_set_attr (field_hdl, enabled = .TRUE., name = TRIM(rst_rfields(i)%vname), &
-                        domain_ref="grid_N", axis_ref="nav_lev", operation = "instant")
-                 CASE ("grid_N")
-                    CALL xios_set_attr (field_hdl, enabled = .TRUE., name = TRIM(rst_rfields(i)%vname), &
-                        domain_ref="grid_N", operation = "instant") 
-                CASE ("grid_vector")
-                    CALL xios_set_attr (field_hdl, enabled = .TRUE., name = TRIM(rst_rfields(i)%vname), &
-                         axis_ref="nav_lev", operation = "instant")
-                 CASE ("grid_scalar")
-                    CALL xios_set_attr (field_hdl, enabled = .TRUE., name = TRIM(rst_rfields(i)%vname), &
-                        scalar_ref = "grid_scalar", operation = "instant")
-                END SELECT
-                IF(lwp) WRITE(numout,*) 'XIOS read: ', TRIM(rst_rfields(i)%vname), ' enabled in ', TRIM(rst_file)
-           ENDIF
-         ENDIF
-        END DO
-#endif
-   END SUBROUTINE iom_set_rstr_active
-
-   SUBROUTINE iom_set_rstw_core(cdmdl)
-      !!---------------------------------------------------------------------
-      !!                   ***  SUBROUTINE  iom_set_rstw_core  ***
-      !!
-      !! ** Purpose :  set variables which are always in restart file 
-      !!---------------------------------------------------------------------
-   CHARACTER (len=*), INTENT (IN) :: cdmdl ! model OPA or SAS
-   CHARACTER(LEN=256)             :: clinfo    ! info character
-#if defined key_iomput
-   IF(cdmdl == "OPA") THEN
-!from restart.F90
-   CALL iom_set_rstw_var_active("rn_Dt")
-   IF ( .NOT. ln_diurnal_only ) THEN
-        CALL iom_set_rstw_var_active('ub'  )
-        CALL iom_set_rstw_var_active('vb'  )
-        CALL iom_set_rstw_var_active('tb'  )
-        CALL iom_set_rstw_var_active('sb'  )
-        CALL iom_set_rstw_var_active('sshb')
-        !
-        CALL iom_set_rstw_var_active('un'  )
-        CALL iom_set_rstw_var_active('vn'  )
-        CALL iom_set_rstw_var_active('tn'  )
-        CALL iom_set_rstw_var_active('sn'  )
-        CALL iom_set_rstw_var_active('sshn')
-        CALL iom_set_rstw_var_active('rhop')
-      ENDIF
-      IF(ln_diurnal) CALL iom_set_rstw_var_active('Dsst')
-!from trasbc.F90
-         CALL iom_set_rstw_var_active('sbc_hc_b')
-         CALL iom_set_rstw_var_active('sbc_sc_b')
-   ENDIF
-#else
-        clinfo = 'iom_set_rstw_core: key_iomput is needed to use XIOS restart read/write functionality'
-        CALL ctl_stop('STOP', TRIM(clinfo))
-#endif
-   END SUBROUTINE iom_set_rstw_core
-
-   SUBROUTINE iom_set_rst_vars(fields)
-      !!---------------------------------------------------------------------
-      !!                   ***  SUBROUTINE iom_set_rst_vars   ***
-      !!
-      !! ** Purpose :  Fill array fields with the information about all 
-      !!               possible variables and corresponding grids definition 
-      !!               for reading/writing restart with XIOS
-      !!---------------------------------------------------------------------
-   TYPE(RST_FIELD), INTENT(INOUT) :: fields(max_rst_fields)
-   INTEGER :: i
-
-        i = 0
-        i = i + 1; fields(i)%vname="rn_Dt";            fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="un";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="ub";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="vn";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="vb";             fields(i)%grid="grid_N_3D"  
-        i = i + 1; fields(i)%vname="tn";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="tb";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="sn";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="sb";             fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="sshn";           fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sshb";           fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="rhop";           fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="kt";             fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="ndastp";         fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="adatrj";         fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="utau_b";         fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="vtau_b";         fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="qns_b";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="emp_b";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sfx_b";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="en" ;            fields(i)%grid="grid_N_3D" 
-        i = i + 1; fields(i)%vname="avt_k";            fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="avm_k";            fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="dissl";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="sbc_hc_b";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sbc_sc_b";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="qsr_hc_b";       fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="fraqsr_1lev";    fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="greenland_icesheet_mass"
-                                               fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="greenland_icesheet_timelapsed"
-                                               fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="greenland_icesheet_mass_roc"
-                                               fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="antarctica_icesheet_mass"
-                                               fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="antarctica_icesheet_timelapsed"
-                                               fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="antarctica_icesheet_mass_roc"
-                                               fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="frc_v";          fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="frc_t";          fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="frc_s";          fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="frc_wn_t";       fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="frc_wn_s";       fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="ssh_ini";        fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="e3t_ini";        fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="hc_loc_ini";     fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="sc_loc_ini";     fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="ssh_hc_loc_ini"; fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ssh_sc_loc_ini"; fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="tilde_e3t_b";    fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="tilde_e3t_n";    fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="hdiv_lf";        fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ub2_b";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="vb2_b";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sshbb_e";        fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ubb_e";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="vbb_e";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sshb_e";         fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ub_e";           fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="vb_e";           fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="fwf_isf_b";      fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="isf_sc_b";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="isf_hc_b";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ssh_ibb";        fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="rnf_b";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="rnf_hc_b";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="rnf_sc_b";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="nn_fsbc";        fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="ssu_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ssv_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sst_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="sss_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ssh_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="e3t_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="frq_m";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="avmb";           fields(i)%grid="grid_vector"
-        i = i + 1; fields(i)%vname="avtb";           fields(i)%grid="grid_vector"
-        i = i + 1; fields(i)%vname="ub2_i_b";        fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="vb2_i_b";        fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="ntime";          fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="Dsst";           fields(i)%grid="grid_scalar"
-        i = i + 1; fields(i)%vname="tmask";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="umask";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="vmask";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="smask";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="gdepw_n";        fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="e3t_n";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="e3u_n";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="e3v_n";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="surf_ini";       fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="e3t_b";          fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="hmxl_n";         fields(i)%grid="grid_N_3D"
-        i = i + 1; fields(i)%vname="un_bf";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="vn_bf";          fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="hbl";            fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="hbli";           fields(i)%grid="grid_N"
-        i = i + 1; fields(i)%vname="wn";             fields(i)%grid="grid_N_3D"
-
-        IF( i-1 > max_rst_fields) THEN
-           WRITE(ctmp1,*) 'E R R O R : iom_set_rst_vars SIZE of RST_FIELD array is too small'
-           CALL ctl_stop( 'iom_set_rst_vars:', ctmp1 )
-        ENDIF
-   END SUBROUTINE iom_set_rst_vars
-
-
-   SUBROUTINE iom_set_rstw_active(cdrst_file)
+      IF(lwp) WRITE(numout,*) 'Setting restart filename (for XIOS write) to: ', TRIM(cdrst_file)
+      CALL xios_get_handle("file_definition", filegroup_hdl )
+      CALL xios_add_child(filegroup_hdl, file_hdl, 'wrestart')
+      IF(nxioso.eq.1) THEN 
+         CALL xios_set_file_attr( "wrestart", type="one_file", enabled=.TRUE.,& 
+                                       mode="write", output_freq=xios_timestep) 
+         IF(lwp) write(numout,*) 'OPEN ', trim(cdrst_file), ' in one_file mode' 
+      ELSE  
+         CALL xios_set_file_attr( "wrestart", type="multiple_file", enabled=.TRUE.,& 
+                                            mode="write", output_freq=xios_timestep) 
+         IF(lwp) write(numout,*) 'OPEN ', trim(cdrst_file), ' in multiple_file mode' 
+      ENDIF 
+      CALL xios_set_file_attr( "wrestart", name=trim(cdrst_file))
+#endif
+   END SUBROUTINE iom_set_rstw_file
+
+
+   SUBROUTINE iom_set_rstw_active(sdfield, rd0, rs0, rd1, rs1, rd2, rs2, rd3, rs3)
       !!---------------------------------------------------------------------
       !!                   ***  SUBROUTINE iom_set_rstw_active   ***
@@ -558,51 +479,110 @@
       !!---------------------------------------------------------------------
 !sets enabled = .TRUE. for each field in restart file
-   CHARACTER(len=*) :: cdrst_file
+      CHARACTER(len = *), INTENT(IN)                     :: sdfield
+      REAL(dp), OPTIONAL, INTENT(IN)                     :: rd0
+      REAL(sp), OPTIONAL, INTENT(IN)                     :: rs0
+      REAL(dp), OPTIONAL, INTENT(IN), DIMENSION(:)       :: rd1
+      REAL(sp), OPTIONAL, INTENT(IN), DIMENSION(:)       :: rs1
+      REAL(dp), OPTIONAL, INTENT(IN), DIMENSION(:, :)    :: rd2
+      REAL(sp), OPTIONAL, INTENT(IN), DIMENSION(:, :)    :: rs2
+      REAL(dp), OPTIONAL, INTENT(IN), DIMENSION(:, :, :) :: rd3  
+      REAL(sp), OPTIONAL, INTENT(IN), DIMENSION(:, :, :) :: rs3
 #if defined key_iomput
-   TYPE(xios_field) :: field_hdl
-   TYPE(xios_file) :: file_hdl
-   TYPE(xios_filegroup) :: filegroup_hdl
-   INTEGER :: i
-   CHARACTER(lc)  ::   clpath
-
-!set name of the restart file and enable available fields
-        IF(lwp) WRITE(numout,*) 'Setting restart filename (for XIOS write) to: ',cdrst_file
-        CALL xios_get_handle("file_definition", filegroup_hdl )
-        CALL xios_add_child(filegroup_hdl, file_hdl, 'wrestart')
-        IF(nxioso.eq.1) THEN 
-           CALL xios_set_file_attr( "wrestart", type="one_file", enabled=.TRUE.,& 
-                                    mode="write", output_freq=xios_timestep) 
-           if(lwp) write(numout,*) 'OPEN ', trim(cdrst_file), ' in one_file mode' 
-        ELSE  
-           CALL xios_set_file_attr( "wrestart", type="multiple_file", enabled=.TRUE.,& 
-                                    mode="write", output_freq=xios_timestep) 
-           if(lwp) write(numout,*) 'OPEN ', trim(cdrst_file), ' in multiple_file mode' 
-        ENDIF 
-        CALL xios_set_file_attr( "wrestart", name=trim(cdrst_file))
+      TYPE(xios_field) :: field_hdl
+      TYPE(xios_file) :: file_hdl
+
+      CALL xios_get_handle("wrestart", file_hdl)
 !define fields for restart context
-        DO i = 1, max_rst_fields
-         IF( rst_wfields(i)%active ) THEN
-                CALL xios_add_child(file_hdl, field_hdl, TRIM(rst_wfields(i)%vname))
-                SELECT CASE (TRIM(rst_wfields(i)%grid))
-                 CASE ("grid_N_3D")
-                    CALL xios_set_attr (field_hdl, enabled = .TRUE., name = TRIM(rst_wfields(i)%vname), &
-                        domain_ref="grid_N", axis_ref="nav_lev", prec = 8, operation = "instant")
-                 CASE ("grid_N")
-                    CALL xios_set_attr (field_hdl, enabled = .TRUE., name = TRIM(rst_wfields(i)%vname), &
-                        domain_ref="grid_N", prec = 8, operation = "instant") 
-                 CASE ("grid_vector")
-                    CALL xios_set_attr (field_hdl, enabled = .TRUE., name = TRIM(rst_wfields(i)%vname), &
-                         axis_ref="nav_lev", prec = 8, operation = "instant")
-                 CASE ("grid_scalar")
-                    CALL xios_set_attr (field_hdl, enabled = .TRUE., name = TRIM(rst_wfields(i)%vname), &
-                        scalar_ref = "grid_scalar", prec = 8, operation = "instant")
-                END SELECT
-         ENDIF
-        END DO
+      CALL xios_add_child(file_hdl, field_hdl, sdfield)
+
+      IF(PRESENT(rd3)) THEN
+         CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
+                             domain_ref = "grid_N",                       &
+                             axis_ref = iom_axis(size(rd3, 3)),           &
+                             prec = 8, operation = "instant"              )
+      ELSEIF(PRESENT(rs3)) THEN
+         CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
+                             domain_ref = "grid_N",                       &
+                             axis_ref = iom_axis(size(rd3, 3)),           &
+                             prec = 4, operation = "instant"              )
+      ELSEIF(PRESENT(rd2)) THEN
+         CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
+                             domain_ref = "grid_N", prec = 8,             &
+                             operation = "instant"                        ) 
+      ELSEIF(PRESENT(rs2)) THEN
+         CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
+                             domain_ref = "grid_N", prec = 4,             &
+                             operation = "instant"                        )
+      ELSEIF(PRESENT(rd1)) THEN
+         CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
+                             axis_ref = iom_axis(size(rd1, 1)),           &
+                             prec = 8, operation = "instant"              )
+      ELSEIF(PRESENT(rs1)) THEN
+         CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
+                             axis_ref = iom_axis(size(rd1, 1)),           &
+                             prec = 4, operation = "instant"              )
+      ELSEIF(PRESENT(rd0)) THEN
+         CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
+                             scalar_ref = "grid_scalar", prec = 8,        &
+                             operation = "instant"                        )
+      ELSEIF(PRESENT(rs0)) THEN
+         CALL xios_set_attr (field_hdl, enabled = .TRUE., name = sdfield, &
+                             scalar_ref = "grid_scalar", prec = 4,        &
+                             operation = "instant"                        )
+      ENDIF
 #endif
    END SUBROUTINE iom_set_rstw_active
 
+   FUNCTION iom_axis(idlev) result(axis_ref)
+      !!---------------------------------------------------------------------
+      !!                   ***  FUNCTION  iom_axis  ***
+      !!
+      !! ** Purpose : Used for grid definition when XIOS is used to read/write
+      !!              restart. Returns axis corresponding to the number of levels
+      !!              given as an input variable. Axes are defined in routine 
+      !!              iom_set_rst_context
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(IN) :: idlev
+      CHARACTER(len=lc)   :: axis_ref
+      CHARACTER(len=12)   :: str
+      IF(idlev == jpk) THEN
+         axis_ref="nav_lev"
+#if defined key_si3
+      ELSEIF(idlev == jpl) THEN
+         axis_ref="numcat"
+#endif         
+      ELSE
+         write(str, *) idlev
+         CALL ctl_stop( 'iom_axis', 'Definition for axis with '//TRIM(ADJUSTL(str))//' levels missing')
+      ENDIF
+   END FUNCTION iom_axis
+
+   FUNCTION iom_xios_setid(cdname) result(kid)
+     !!---------------------------------------------------------------------
+      !!                   ***  FUNCTION    ***
+      !!
+      !! ** Purpose : this function returns first available id to keep information about file 
+      !!              sets filename in iom_file structure and sets name
+      !!              of XIOS context depending on cdcomp
+      !!              corresponds to iom_nf90_open
+      !!---------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in   ) :: cdname      ! File name
+      INTEGER                         :: kid      ! identifier of the opened file
+      INTEGER                         :: jl
+
+      kid = 0
+      DO jl = jpmax_files, 1, -1
+         IF( iom_file(jl)%nfid == 0 )   kid = jl
+      ENDDO
+
+      iom_file(kid)%name   = TRIM(cdname)
+      iom_file(kid)%nfid   = 1
+      iom_file(kid)%nvars  = 0
+      iom_file(kid)%irec   = -1
+
+   END FUNCTION iom_xios_setid
+
    SUBROUTINE iom_set_rst_context(ld_rstr) 
-     !!---------------------------------------------------------------------
+      !!---------------------------------------------------------------------
       !!                   ***  SUBROUTINE  iom_set_rst_context  ***
       !!
@@ -611,30 +591,71 @@
       !!               
       !!---------------------------------------------------------------------
-   LOGICAL, INTENT(IN)               :: ld_rstr
-!ld_rstr is true for restart context. There is no need to define grid for 
-!restart read, because it's read from file
+      LOGICAL, INTENT(IN)               :: ld_rstr
+      INTEGER :: ji
 #if defined key_iomput
-   TYPE(xios_domaingroup)            :: domaingroup_hdl 
-   TYPE(xios_domain)                 :: domain_hdl 
-   TYPE(xios_axisgroup)              :: axisgroup_hdl 
-   TYPE(xios_axis)                   :: axis_hdl 
-   TYPE(xios_scalar)                 :: scalar_hdl 
-   TYPE(xios_scalargroup)            :: scalargroup_hdl 
-
-     CALL xios_get_handle("domain_definition",domaingroup_hdl) 
-     CALL xios_add_child(domaingroup_hdl, domain_hdl, "grid_N") 
-     CALL set_grid("N", glamt, gphit, .TRUE., ld_rstr) 
+      TYPE(xios_domaingroup)            :: domaingroup_hdl 
+      TYPE(xios_domain)                 :: domain_hdl 
+      TYPE(xios_axisgroup)              :: axisgroup_hdl 
+      TYPE(xios_axis)                   :: axis_hdl 
+      TYPE(xios_scalar)                 :: scalar_hdl 
+      TYPE(xios_scalargroup)            :: scalargroup_hdl 
+
+      CALL xios_get_handle("domain_definition",domaingroup_hdl) 
+      CALL xios_add_child(domaingroup_hdl, domain_hdl, "grid_N") 
+      CALL set_grid("N", glamt, gphit, .TRUE., ld_rstr) 
  
-     CALL xios_get_handle("axis_definition",axisgroup_hdl) 
-     CALL xios_add_child(axisgroup_hdl, axis_hdl, "nav_lev") 
+      CALL xios_get_handle("axis_definition",axisgroup_hdl) 
+      CALL xios_add_child(axisgroup_hdl, axis_hdl, "nav_lev") 
 !AGRIF fails to compile when unit= is in call to xios_set_axis_attr
-!    CALL xios_set_axis_attr( "nav_lev", long_name="Vertical levels",  unit="m", positive="down") 
-     CALL xios_set_axis_attr( "nav_lev", long_name="Vertical levels in meters", positive="down")
-     CALL iom_set_axis_attr( "nav_lev", paxis = gdept_1d ) 
-
-     CALL xios_get_handle("scalar_definition", scalargroup_hdl) 
-     CALL xios_add_child(scalargroup_hdl, scalar_hdl, "grid_scalar") 
+!     CALL xios_set_axis_attr( "nav_lev", long_name="Vertical levels",  unit="m", positive="down") 
+      CALL xios_set_axis_attr( "nav_lev", long_name = "Vertical levels in meters", positive = "down")
+      CALL iom_set_axis_attr( "nav_lev", paxis = gdept_1d ) 
+#if defined key_si3
+      CALL xios_add_child(axisgroup_hdl, axis_hdl, "numcat")
+      CALL iom_set_axis_attr( "numcat", (/ (REAL(ji,wp), ji=1,jpl) /) )
+#endif
+      CALL xios_get_handle("scalar_definition", scalargroup_hdl) 
+      CALL xios_add_child(scalargroup_hdl, scalar_hdl, "grid_scalar") 
 #endif
    END SUBROUTINE iom_set_rst_context
+
+
+   SUBROUTINE set_xios_context(kdid, cdcont) 
+      !!---------------------------------------------------------------------
+      !!                   ***  SUBROUTINE  iom_set_rst_context  ***
+      !!
+      !! ** Purpose : set correct XIOS context based on kdid
+      !!               
+      !!---------------------------------------------------------------------
+      INTEGER,           INTENT(IN)     :: kdid           ! Identifier of the file
+      CHARACTER(LEN=lc), INTENT(OUT)    :: cdcont         ! name of the context for XIOS read/write
+      
+      cdcont = "NONE"
+
+      IF(lrxios) THEN
+         IF(kdid == numror) THEN
+            cdcont = cr_ocerst_cxt
+         ELSEIF(kdid == numrir) THEN
+            cdcont = cr_icerst_cxt 
+         ELSEIF(kdid == numrtr) THEN
+            cdcont = cr_toprst_cxt
+         ELSEIF(kdid == numrsr) THEN
+            cdcont = cr_sedrst_cxt
+         ENDIF
+      ENDIF
+
+      IF(lwxios) THEN
+         IF(kdid == numrow) THEN
+            cdcont = cw_ocerst_cxt
+         ELSEIF(kdid == numriw) THEN
+            cdcont = cw_icerst_cxt
+         ELSEIF(kdid == numrtw) THEN
+            cdcont = cw_toprst_cxt
+         ELSEIF(kdid == numrsw) THEN
+            cdcont = cw_sedrst_cxt
+         ENDIF
+      ENDIF
+   END SUBROUTINE set_xios_context
+
 
    SUBROUTINE iom_swap( cdname )
@@ -647,5 +668,4 @@
 #if defined key_iomput
       TYPE(xios_context) :: nemo_hdl
-
       IF( TRIM(Agrif_CFixed()) == '0' ) THEN
         CALL xios_get_handle(TRIM(cdname),nemo_hdl)
@@ -897,5 +917,5 @@
    !!                   INTERFACE iom_get
    !!----------------------------------------------------------------------
-   SUBROUTINE iom_g0d_sp( kiomid, cdvar, pvar, ktime, ldxios )
+   SUBROUTINE iom_g0d_sp( kiomid, cdvar, pvar, ktime )
       INTEGER         , INTENT(in   )                 ::   kiomid    ! Identifier of the file
       CHARACTER(len=*), INTENT(in   )                 ::   cdvar     ! Name of the variable
@@ -903,5 +923,4 @@
       REAL(dp)                                        ::   ztmp_pvar ! tmp var to read field
       INTEGER         , INTENT(in   ),     OPTIONAL   ::   ktime     ! record number
-      LOGICAL         , INTENT(in   ),     OPTIONAL   ::   ldxios    ! use xios to read restart
       !
       INTEGER                                         ::   idvar     ! variable id
@@ -911,10 +930,9 @@
       CHARACTER(LEN=100)                              ::   clname    ! file name
       CHARACTER(LEN=1)                                ::   cldmspc   !
-      LOGICAL                                         ::   llxios
-      !
-      llxios = .FALSE.
-      IF( PRESENT(ldxios) ) llxios = ldxios
-
-      IF(.NOT.llxios) THEN  ! read data using default library
+      CHARACTER(LEN=lc)                               ::   context
+      !
+      CALL set_xios_context(kiomid, context)
+
+      IF(context == "NONE") THEN  ! read data using default library
          itime = 1
          IF( PRESENT(ktime) ) itime = ktime
@@ -939,7 +957,7 @@
 #if defined key_iomput
          IF(lwp) WRITE(numout,*) 'XIOS RST READ (0D): ', trim(cdvar)
-         CALL iom_swap( TRIM(crxios_context) )
+         CALL iom_swap(context)
          CALL xios_recv_field( trim(cdvar), pvar)
-         CALL iom_swap( TRIM(cxios_context) )
+         CALL iom_swap(cxios_context)
 #else
          WRITE(ctmp1,*) 'Can not use XIOS in iom_g0d, file: '//trim(clname)//', var:'//trim(cdvar)
@@ -949,10 +967,9 @@
    END SUBROUTINE iom_g0d_sp
 
-   SUBROUTINE iom_g0d_dp( kiomid, cdvar, pvar, ktime, ldxios )
+   SUBROUTINE iom_g0d_dp( kiomid, cdvar, pvar, ktime )
       INTEGER         , INTENT(in   )                 ::   kiomid    ! Identifier of the file
       CHARACTER(len=*), INTENT(in   )                 ::   cdvar     ! Name of the variable
       REAL(dp)        , INTENT(  out)                 ::   pvar      ! read field
       INTEGER         , INTENT(in   ),     OPTIONAL   ::   ktime     ! record number
-      LOGICAL         , INTENT(in   ),     OPTIONAL   ::   ldxios    ! use xios to read restart
       !
       INTEGER                                         ::   idvar     ! variable id
@@ -962,10 +979,9 @@
       CHARACTER(LEN=100)                              ::   clname    ! file name
       CHARACTER(LEN=1)                                ::   cldmspc   !
-      LOGICAL                                         ::   llxios
-      !
-      llxios = .FALSE.
-      IF( PRESENT(ldxios) ) llxios = ldxios
-
-      IF(.NOT.llxios) THEN  ! read data using default library
+      CHARACTER(LEN=lc)                               ::   context
+      !
+      CALL set_xios_context(kiomid, context)
+
+      IF(context == "NONE") THEN  ! read data using default library
          itime = 1
          IF( PRESENT(ktime) ) itime = ktime
@@ -989,7 +1005,7 @@
 #if defined key_iomput
          IF(lwp) WRITE(numout,*) 'XIOS RST READ (0D): ', trim(cdvar)
-         CALL iom_swap( TRIM(crxios_context) )
+         CALL iom_swap(context)
          CALL xios_recv_field( trim(cdvar), pvar)
-         CALL iom_swap( TRIM(cxios_context) )
+         CALL iom_swap(cxios_context)
 #else
          WRITE(ctmp1,*) 'Can not use XIOS in iom_g0d, file: '//trim(clname)//', var:'//trim(cdvar)
@@ -999,5 +1015,5 @@
    END SUBROUTINE iom_g0d_dp
 
-   SUBROUTINE iom_g1d_sp( kiomid, kdom, cdvar, pvar, ktime, kstart, kcount, ldxios )
+   SUBROUTINE iom_g1d_sp( kiomid, kdom, cdvar, pvar, ktime, kstart, kcount )
       INTEGER         , INTENT(in   )                         ::   kiomid    ! Identifier of the file
       INTEGER         , INTENT(in   )                         ::   kdom      ! Type of domain to be read
@@ -1008,5 +1024,4 @@
       INTEGER         , INTENT(in   ), DIMENSION(1), OPTIONAL ::   kstart    ! start axis position of the reading 
       INTEGER         , INTENT(in   ), DIMENSION(1), OPTIONAL ::   kcount    ! number of points in each axis
-      LOGICAL         , INTENT(in   ),               OPTIONAL ::   ldxios    ! read data using XIOS
       !
       IF( kiomid > 0 ) THEN
@@ -1014,6 +1029,5 @@
             ALLOCATE(ztmp_pvar(size(pvar,1)))
             CALL iom_get_123d( kiomid, kdom       , cdvar        , pv_r1d=ztmp_pvar,   &
-              &                                                     ktime=ktime, kstart=kstart, kcount=kcount, &
-              &                                                     ldxios=ldxios )
+              &                                                     ktime=ktime, kstart=kstart, kcount=kcount )
             pvar = ztmp_pvar
             DEALLOCATE(ztmp_pvar)
@@ -1023,5 +1037,5 @@
 
 
-   SUBROUTINE iom_g1d_dp( kiomid, kdom, cdvar, pvar, ktime, kstart, kcount, ldxios )
+   SUBROUTINE iom_g1d_dp( kiomid, kdom, cdvar, pvar, ktime, kstart, kcount )
       INTEGER         , INTENT(in   )                         ::   kiomid    ! Identifier of the file
       INTEGER         , INTENT(in   )                         ::   kdom      ! Type of domain to be read
@@ -1031,14 +1045,12 @@
       INTEGER         , INTENT(in   ), DIMENSION(1), OPTIONAL ::   kstart    ! start axis position of the reading 
       INTEGER         , INTENT(in   ), DIMENSION(1), OPTIONAL ::   kcount    ! number of points in each axis
-      LOGICAL         , INTENT(in   ),               OPTIONAL ::   ldxios    ! read data using XIOS
       !
       IF( kiomid > 0 ) THEN
          IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom       , cdvar        , pv_r1d=pvar,   &
-              &                                                     ktime=ktime, kstart=kstart, kcount=kcount, &
-              &                                                     ldxios=ldxios )
+              &                                                     ktime=ktime, kstart=kstart, kcount=kcount)
       ENDIF
    END SUBROUTINE iom_g1d_dp
 
-   SUBROUTINE iom_g2d_sp( kiomid, kdom, cdvar, pvar, ktime, cd_type, psgn, kfill, kstart, kcount, ldxios)
+   SUBROUTINE iom_g2d_sp( kiomid, kdom, cdvar, pvar, ktime, cd_type, psgn, kfill, kstart, kcount)
       INTEGER         , INTENT(in   )                         ::   kiomid    ! Identifier of the file
       INTEGER         , INTENT(in   )                         ::   kdom      ! Type of domain to be read
@@ -1052,5 +1064,4 @@
       INTEGER         , INTENT(in   ), DIMENSION(2), OPTIONAL ::   kstart    ! start axis position of the reading 
       INTEGER         , INTENT(in   ), DIMENSION(2), OPTIONAL ::   kcount    ! number of points in each axis
-      LOGICAL         , INTENT(in   ),               OPTIONAL ::   ldxios    ! read data using XIOS
       !
       IF( kiomid > 0 ) THEN
@@ -1059,5 +1070,5 @@
             CALL iom_get_123d( kiomid, kdom, cdvar      , pv_r2d = ztmp_pvar  , ktime = ktime,   &
              &                                                      cd_type = cd_type, psgn   = psgn  , kfill = kfill,   &
-             &                                                      kstart  = kstart , kcount = kcount, ldxios=ldxios  )
+             &                                                      kstart  = kstart , kcount = kcount  )
             pvar = ztmp_pvar
             DEALLOCATE(ztmp_pvar)
@@ -1066,5 +1077,5 @@
    END SUBROUTINE iom_g2d_sp
 
-   SUBROUTINE iom_g2d_dp( kiomid, kdom, cdvar, pvar, ktime, cd_type, psgn, kfill, kstart, kcount, ldxios)
+   SUBROUTINE iom_g2d_dp( kiomid, kdom, cdvar, pvar, ktime, cd_type, psgn, kfill, kstart, kcount)
       INTEGER         , INTENT(in   )                         ::   kiomid    ! Identifier of the file
       INTEGER         , INTENT(in   )                         ::   kdom      ! Type of domain to be read
@@ -1077,14 +1088,13 @@
       INTEGER         , INTENT(in   ), DIMENSION(2), OPTIONAL ::   kstart    ! start axis position of the reading 
       INTEGER         , INTENT(in   ), DIMENSION(2), OPTIONAL ::   kcount    ! number of points in each axis
-      LOGICAL         , INTENT(in   ),               OPTIONAL ::   ldxios    ! read data using XIOS
       !
       IF( kiomid > 0 ) THEN
          IF( iom_file(kiomid)%nfid > 0 ) CALL iom_get_123d( kiomid, kdom, cdvar      , pv_r2d = pvar  , ktime = ktime,   &
             &                                                       cd_type = cd_type, psgn   = psgn  , kfill = kfill,   &
-            &                                                       kstart  = kstart , kcount = kcount, ldxios=ldxios  )
+            &                                                       kstart  = kstart , kcount = kcount                )
       ENDIF
    END SUBROUTINE iom_g2d_dp
 
-   SUBROUTINE iom_g3d_sp( kiomid, kdom, cdvar, pvar, ktime, cd_type, psgn, kfill, kstart, kcount, ldxios )
+   SUBROUTINE iom_g3d_sp( kiomid, kdom, cdvar, pvar, ktime, cd_type, psgn, kfill, kstart, kcount )
       INTEGER         , INTENT(in   )                         ::   kiomid    ! Identifier of the file
       INTEGER         , INTENT(in   )                         ::   kdom      ! Type of domain to be read
@@ -1098,5 +1108,4 @@
       INTEGER         , INTENT(in   ), DIMENSION(3), OPTIONAL ::   kstart    ! start axis position of the reading 
       INTEGER         , INTENT(in   ), DIMENSION(3), OPTIONAL ::   kcount    ! number of points in each axis
-      LOGICAL         , INTENT(in   ),               OPTIONAL ::   ldxios    ! read data using XIOS
       !
       IF( kiomid > 0 ) THEN
@@ -1105,5 +1114,5 @@
             CALL iom_get_123d( kiomid, kdom, cdvar      , pv_r3d = ztmp_pvar  , ktime = ktime,   &
             &                                                       cd_type = cd_type, psgn   = psgn  , kfill = kfill,   &
-            &                                                       kstart  = kstart , kcount = kcount, ldxios=ldxios  )
+            &                                                       kstart  = kstart , kcount = kcount                )
             pvar = ztmp_pvar
             DEALLOCATE(ztmp_pvar)
@@ -1112,5 +1121,5 @@
    END SUBROUTINE iom_g3d_sp
 
-   SUBROUTINE iom_g3d_dp( kiomid, kdom, cdvar, pvar, ktime, cd_type, psgn, kfill, kstart, kcount, ldxios )
+   SUBROUTINE iom_g3d_dp( kiomid, kdom, cdvar, pvar, ktime, cd_type, psgn, kfill, kstart, kcount )
       INTEGER         , INTENT(in   )                         ::   kiomid    ! Identifier of the file
       INTEGER         , INTENT(in   )                         ::   kdom      ! Type of domain to be read
@@ -1123,5 +1132,4 @@
       INTEGER         , INTENT(in   ), DIMENSION(3), OPTIONAL ::   kstart    ! start axis position of the reading 
       INTEGER         , INTENT(in   ), DIMENSION(3), OPTIONAL ::   kcount    ! number of points in each axis
-      LOGICAL         , INTENT(in   ),               OPTIONAL ::   ldxios    ! read data using XIOS
       !
       IF( kiomid > 0 ) THEN
@@ -1129,5 +1137,5 @@
             CALL iom_get_123d( kiomid, kdom, cdvar      , pv_r3d = pvar  , ktime = ktime,   &
             &                                                       cd_type = cd_type, psgn   = psgn  , kfill = kfill,   &
-            &                                                       kstart  = kstart , kcount = kcount, ldxios=ldxios  )
+            &                                                       kstart  = kstart , kcount = kcount                )
          END IF
       ENDIF
@@ -1137,5 +1145,5 @@
 
    SUBROUTINE iom_get_123d( kiomid , kdom, cdvar, pv_r1d, pv_r2d, pv_r3d, ktime ,   &
-         &                  cd_type, psgn, kfill, kstart, kcount, ldxios )
+         &                  cd_type, psgn, kfill, kstart, kcount )
       !!-----------------------------------------------------------------------
       !!                  ***  ROUTINE  iom_get_123d  ***
@@ -1157,8 +1165,6 @@
       INTEGER , DIMENSION(:)     , INTENT(in   ), OPTIONAL ::   kstart    ! start position of the reading in each axis 
       INTEGER , DIMENSION(:)     , INTENT(in   ), OPTIONAL ::   kcount    ! number of points to be read in each axis
-      LOGICAL                    , INTENT(in   ), OPTIONAL ::   ldxios    ! use XIOS to read restart
       !
       LOGICAL                        ::   llok        ! true if ok!
-      LOGICAL                        ::   llxios      ! local definition for XIOS read
       INTEGER                        ::   jl          ! loop on number of dimension 
       INTEGER                        ::   idom        ! type of domain
@@ -1187,15 +1193,14 @@
       REAL(dp)                       ::   gma, gmi
       !---------------------------------------------------------------------
-      !
+      CHARACTER(LEN=lc)                               ::   context
+      !
+      CALL set_xios_context(kiomid, context)
       inlev = -1
       IF( PRESENT(pv_r3d) )   inlev = SIZE(pv_r3d, 3)
       !
-      llxios = .FALSE.
-      IF( PRESENT(ldxios) )   llxios = ldxios
-      !
       idom = kdom
       istop = nstop
       !
-      IF(.NOT.llxios) THEN
+      IF(context == "NONE") THEN
          clname = iom_file(kiomid)%name   !   esier to read
          clinfo = '          iom_get_123d, file: '//trim(clname)//', var: '//trim(cdvar)
@@ -1364,18 +1369,28 @@
 #if defined key_iomput
 !would be good to be able to check which context is active and swap only if current is not restart
-         CALL iom_swap( TRIM(crxios_context) ) 
+         idvar = iom_varid( kiomid, cdvar )
+         CALL iom_swap(context)
+         zsgn = 1._wp
+         IF( PRESENT(psgn   ) )   zsgn    = psgn
+         cl_type = 'T'
+         IF( PRESENT(cd_type) )   cl_type = cd_type
+
          IF( PRESENT(pv_r3d) ) THEN
             IF(lwp) WRITE(numout,*) 'XIOS RST READ (3D): ',TRIM(cdvar)
-            CALL xios_recv_field( trim(cdvar), pv_r3d)
-            IF(idom /= jpdom_unknown )   CALL lbc_lnk( 'iom', pv_r3d,'Z', -999., kfillmode = jpfillnothing)
+            CALL xios_recv_field( trim(cdvar), pv_r3d(:, :, :))
+            IF(idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+               CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill)
+            ENDIF
          ELSEIF( PRESENT(pv_r2d) ) THEN
             IF(lwp) WRITE(numout,*) 'XIOS RST READ (2D): ', TRIM(cdvar)
-            CALL xios_recv_field( trim(cdvar), pv_r2d)
-            IF(idom /= jpdom_unknown )   CALL lbc_lnk('iom', pv_r2d,'Z',-999., kfillmode = jpfillnothing)
+            CALL xios_recv_field( trim(cdvar), pv_r2d(:, :))
+            IF(idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+               CALL lbc_lnk('iom', pv_r2d, cl_type, zsgn, kfillmode = kfill)
+            ENDIF
          ELSEIF( PRESENT(pv_r1d) ) THEN
             IF(lwp) WRITE(numout,*) 'XIOS RST READ (1D): ', TRIM(cdvar)
             CALL xios_recv_field( trim(cdvar), pv_r1d)
          ENDIF
-         CALL iom_swap( TRIM(cxios_context) )
+         CALL iom_swap(cxios_context)
 #else
          istop = istop + 1 
@@ -1392,12 +1407,12 @@
       zofs = iom_file(kiomid)%ofs(idvar)      ! offset
       IF(     PRESENT(pv_r1d) ) THEN
-         IF( zscf /= 1. )   pv_r1d(:) = pv_r1d(:) * zscf 
-         IF( zofs /= 0. )   pv_r1d(:) = pv_r1d(:) + zofs
+         IF( zscf /= 1._wp )   pv_r1d(:) = pv_r1d(:) * zscf 
+         IF( zofs /= 0._wp )   pv_r1d(:) = pv_r1d(:) + zofs
       ELSEIF( PRESENT(pv_r2d) ) THEN
-         IF( zscf /= 1.)   pv_r2d(:,:) = pv_r2d(:,:) * zscf
-         IF( zofs /= 0.)   pv_r2d(:,:) = pv_r2d(:,:) + zofs
+         IF( zscf /= 1._wp)   pv_r2d(:,:) = pv_r2d(:,:) * zscf
+         IF( zofs /= 0._wp)   pv_r2d(:,:) = pv_r2d(:,:) + zofs
       ELSEIF( PRESENT(pv_r3d) ) THEN
-         IF( zscf /= 1.)   pv_r3d(:,:,:) = pv_r3d(:,:,:) * zscf
-         IF( zofs /= 0.)   pv_r3d(:,:,:) = pv_r3d(:,:,:) + zofs
+         IF( zscf /= 1._wp)   pv_r3d(:,:,:) = pv_r3d(:,:,:) * zscf
+         IF( zofs /= 0._wp)   pv_r3d(:,:,:) = pv_r3d(:,:,:) + zofs
       ENDIF
       !
@@ -1573,5 +1588,5 @@
    !!                   INTERFACE iom_rstput
    !!----------------------------------------------------------------------
-   SUBROUTINE iom_rp0d_sp( kt, kwrite, kiomid, cdvar, pvar, ktype, ldxios )
+   SUBROUTINE iom_rp0d_sp( kt, kwrite, kiomid, cdvar, pvar, ktype )
       INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
       INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
@@ -1580,16 +1595,26 @@
       REAL(sp)        , INTENT(in)                         ::   pvar     ! written field
       INTEGER         , INTENT(in), OPTIONAL               ::   ktype    ! variable external type
-      LOGICAL, OPTIONAL :: ldxios   ! xios write flag
-      LOGICAL :: llx                ! local xios write flag
-      INTEGER :: ivid   ! variable id
-
-      llx = .FALSE.
-      IF(PRESENT(ldxios)) llx = ldxios
+      !
+      LOGICAL           :: llx                ! local xios write flag
+      INTEGER           :: ivid   ! variable id
+      CHARACTER(LEN=lc) :: context
+      !
+      CALL set_xios_context(kiomid, context)
+
+      llx = .NOT. (context == "NONE")
+
       IF( llx ) THEN
 #ifdef key_iomput
-      IF( kt == kwrite ) THEN
-          IF(lwp) write(numout,*) 'RESTART: write (XIOS 0D) ',trim(cdvar)
-          CALL xios_send_field(trim(cdvar), pvar)
-      ENDIF
+         IF( kt == kwrite ) THEN
+            IF(lwp) write(numout,*) 'RESTART: write (XIOS 0D) ',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_put(trim(cdvar), pvar)
+            CALL iom_swap(cxios_context)
+         ELSE
+            IF(lwp) write(numout,*) 'RESTART: define (XIOS 0D) ',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_set_rstw_active( trim(cdvar), rs0 = pvar ) 
+            CALL iom_swap(cxios_context)
+         ENDIF
 #endif
       ELSE
@@ -1603,5 +1628,5 @@
    END SUBROUTINE iom_rp0d_sp
 
-   SUBROUTINE iom_rp0d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype, ldxios )
+   SUBROUTINE iom_rp0d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype )
       INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
       INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
@@ -1610,16 +1635,26 @@
       REAL(dp)        , INTENT(in)                         ::   pvar     ! written field
       INTEGER         , INTENT(in), OPTIONAL               ::   ktype    ! variable external type
-      LOGICAL, OPTIONAL :: ldxios   ! xios write flag
-      LOGICAL :: llx                ! local xios write flag
-      INTEGER :: ivid   ! variable id
-
-      llx = .FALSE.
-      IF(PRESENT(ldxios)) llx = ldxios
+      !
+      LOGICAL           :: llx                ! local xios write flag
+      INTEGER           :: ivid   ! variable id
+      CHARACTER(LEN=lc) :: context
+      !
+      CALL set_xios_context(kiomid, context)
+
+      llx = .NOT. (context == "NONE")
+
       IF( llx ) THEN
 #ifdef key_iomput
-      IF( kt == kwrite ) THEN
-          IF(lwp) write(numout,*) 'RESTART: write (XIOS 0D) ',trim(cdvar)
-          CALL xios_send_field(trim(cdvar), pvar)
-      ENDIF
+         IF( kt == kwrite ) THEN
+            IF(lwp) write(numout,*) 'RESTART: write (XIOS 0D) ',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_put(trim(cdvar), pvar)
+            CALL iom_swap(cxios_context)
+         ELSE
+            IF(lwp) write(numout,*) 'RESTART: define (XIOS 0D) ',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_set_rstw_active( trim(cdvar), rd0 = pvar ) 
+            CALL iom_swap(cxios_context)
+         ENDIF
 #endif
       ELSE
@@ -1634,5 +1669,5 @@
 
 
-   SUBROUTINE iom_rp1d_sp( kt, kwrite, kiomid, cdvar, pvar, ktype, ldxios )
+   SUBROUTINE iom_rp1d_sp( kt, kwrite, kiomid, cdvar, pvar, ktype )
       INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
       INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
@@ -1641,16 +1676,26 @@
       REAL(sp)        , INTENT(in), DIMENSION(          :) ::   pvar     ! written field
       INTEGER         , INTENT(in), OPTIONAL               ::   ktype    ! variable external type
-      LOGICAL, OPTIONAL                                    ::   ldxios   ! xios write flag
-      LOGICAL :: llx                ! local xios write flag
-      INTEGER :: ivid   ! variable id
-
-      llx = .FALSE.
-      IF(PRESENT(ldxios)) llx = ldxios
+      !
+      LOGICAL           :: llx                ! local xios write flag
+      INTEGER           :: ivid   ! variable id
+      CHARACTER(LEN=lc) :: context
+      !
+      CALL set_xios_context(kiomid, context)
+
+      llx = .NOT. (context == "NONE")
+
       IF( llx ) THEN
 #ifdef key_iomput
-      IF( kt == kwrite ) THEN
-         IF(lwp) write(numout,*) 'RESTART: write (XIOS 1D) ',trim(cdvar)
-         CALL xios_send_field(trim(cdvar), pvar)
-      ENDIF
+         IF( kt == kwrite ) THEN
+            IF(lwp) write(numout,*) 'RESTART: write (XIOS 1D) ',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_put(trim(cdvar), pvar)
+            CALL iom_swap(cxios_context)
+         ELSE
+            IF(lwp) write(numout,*) 'RESTART: define (XIOS 1D)',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_set_rstw_active( trim(cdvar), rs1 = pvar )
+            CALL iom_swap(cxios_context)
+         ENDIF
 #endif
       ELSE
@@ -1664,5 +1709,5 @@
    END SUBROUTINE iom_rp1d_sp
 
-   SUBROUTINE iom_rp1d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype, ldxios )
+   SUBROUTINE iom_rp1d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype )
       INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
       INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
@@ -1671,16 +1716,26 @@
       REAL(dp)        , INTENT(in), DIMENSION(          :) ::   pvar     ! written field
       INTEGER         , INTENT(in), OPTIONAL               ::   ktype    ! variable external type
-      LOGICAL, OPTIONAL                                    ::   ldxios   ! xios write flag
-      LOGICAL :: llx                ! local xios write flag
-      INTEGER :: ivid   ! variable id
-
-      llx = .FALSE.
-      IF(PRESENT(ldxios)) llx = ldxios
+      !
+      LOGICAL           :: llx                ! local xios write flag
+      INTEGER           :: ivid   ! variable id
+      CHARACTER(LEN=lc) :: context
+      !
+      CALL set_xios_context(kiomid, context)
+
+      llx = .NOT. (context == "NONE")
+
       IF( llx ) THEN
 #ifdef key_iomput
-      IF( kt == kwrite ) THEN
-         IF(lwp) write(numout,*) 'RESTART: write (XIOS 1D) ',trim(cdvar)
-         CALL xios_send_field(trim(cdvar), pvar)
-      ENDIF
+         IF( kt == kwrite ) THEN
+            IF(lwp) write(numout,*) 'RESTART: write (XIOS 1D) ',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_put(trim(cdvar), pvar)
+            CALL iom_swap(cxios_context)
+         ELSE
+            IF(lwp) write(numout,*) 'RESTART: define (XIOS 1D)',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_set_rstw_active( trim(cdvar), rd1 = pvar )
+            CALL iom_swap(cxios_context)
+         ENDIF
 #endif
       ELSE
@@ -1695,5 +1750,5 @@
 
 
-   SUBROUTINE iom_rp2d_sp( kt, kwrite, kiomid, cdvar, pvar, ktype, ldxios )
+   SUBROUTINE iom_rp2d_sp( kt, kwrite, kiomid, cdvar, pvar, ktype )
       INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
       INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
@@ -1702,16 +1757,26 @@
       REAL(sp)        , INTENT(in), DIMENSION(:,    :    ) ::   pvar     ! written field
       INTEGER         , INTENT(in), OPTIONAL               ::   ktype    ! variable external type
-      LOGICAL, OPTIONAL :: ldxios   ! xios write flag
-      LOGICAL :: llx
-      INTEGER :: ivid   ! variable id
-
-      llx = .FALSE.
-      IF(PRESENT(ldxios)) llx = ldxios
+      !
+      LOGICAL            :: llx
+      INTEGER            :: ivid   ! variable id
+      CHARACTER(LEN=lc)  :: context
+      !
+      CALL set_xios_context(kiomid, context)
+
+      llx = .NOT. (context == "NONE")
+
       IF( llx ) THEN
 #ifdef key_iomput
-      IF( kt == kwrite ) THEN
-         IF(lwp) write(numout,*) 'RESTART: write (XIOS 2D) ',trim(cdvar)
-         CALL xios_send_field(trim(cdvar), pvar)
-      ENDIF
+         IF( kt == kwrite ) THEN
+            IF(lwp) write(numout,*) 'RESTART: write (XIOS 2D) ',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_put(trim(cdvar), pvar)
+            CALL iom_swap(cxios_context)
+         ELSE
+            IF(lwp) write(numout,*) 'RESTART: define (XIOS 2D)',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_set_rstw_active( trim(cdvar), rs2 = pvar )
+            CALL iom_swap(cxios_context)
+         ENDIF
 #endif
       ELSE
@@ -1725,5 +1790,5 @@
    END SUBROUTINE iom_rp2d_sp
 
-   SUBROUTINE iom_rp2d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype, ldxios )
+   SUBROUTINE iom_rp2d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype )
       INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
       INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
@@ -1732,16 +1797,26 @@
       REAL(dp)        , INTENT(in), DIMENSION(:,    :    ) ::   pvar     ! written field
       INTEGER         , INTENT(in), OPTIONAL               ::   ktype    ! variable external type
-      LOGICAL, OPTIONAL :: ldxios   ! xios write flag
-      LOGICAL :: llx
-      INTEGER :: ivid   ! variable id
-
-      llx = .FALSE.
-      IF(PRESENT(ldxios)) llx = ldxios
+      !
+      LOGICAL           :: llx
+      INTEGER           :: ivid   ! variable id
+      CHARACTER(LEN=lc) :: context
+      !
+      CALL set_xios_context(kiomid, context)
+
+      llx = .NOT. (context == "NONE")
+
       IF( llx ) THEN
 #ifdef key_iomput
-      IF( kt == kwrite ) THEN
-         IF(lwp) write(numout,*) 'RESTART: write (XIOS 2D) ',trim(cdvar)
-         CALL xios_send_field(trim(cdvar), pvar)
-      ENDIF
+         IF( kt == kwrite ) THEN
+            IF(lwp) write(numout,*) 'RESTART: write (XIOS 2D) ',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_put(trim(cdvar), pvar)
+            CALL iom_swap(cxios_context)
+         ELSE
+            IF(lwp) write(numout,*) 'RESTART: define (XIOS 2D)',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_set_rstw_active( trim(cdvar), rd2 = pvar )
+            CALL iom_swap(cxios_context)
+         ENDIF
 #endif
       ELSE
@@ -1756,5 +1831,5 @@
 
 
-   SUBROUTINE iom_rp3d_sp( kt, kwrite, kiomid, cdvar, pvar, ktype, ldxios )
+   SUBROUTINE iom_rp3d_sp( kt, kwrite, kiomid, cdvar, pvar, ktype )
       INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
       INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
@@ -1763,16 +1838,26 @@
       REAL(sp)        , INTENT(in),       DIMENSION(:,:,:) ::   pvar     ! written field
       INTEGER         , INTENT(in), OPTIONAL               ::   ktype    ! variable external type
-      LOGICAL, OPTIONAL :: ldxios   ! xios write flag
-      LOGICAL :: llx                 ! local xios write flag
-      INTEGER :: ivid   ! variable id
-
-      llx = .FALSE.
-      IF(PRESENT(ldxios)) llx = ldxios
+      !
+      LOGICAL           :: llx                 ! local xios write flag
+      INTEGER           :: ivid   ! variable id
+      CHARACTER(LEN=lc) :: context
+      !
+      CALL set_xios_context(kiomid, context)
+
+      llx = .NOT. (context == "NONE")
+
       IF( llx ) THEN
 #ifdef key_iomput
-      IF( kt == kwrite ) THEN
-         IF(lwp) write(numout,*) 'RESTART: write (XIOS 3D) ',trim(cdvar)
-         CALL xios_send_field(trim(cdvar), pvar)
-      ENDIF
+         IF( kt == kwrite ) THEN
+            IF(lwp) write(numout,*) 'RESTART: write (XIOS 3D) ',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_put(trim(cdvar), pvar)
+            CALL iom_swap(cxios_context)
+         ELSE
+            IF(lwp) write(numout,*) 'RESTART: define (XIOS 3D)',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_set_rstw_active( trim(cdvar), rs3 = pvar )
+            CALL iom_swap(cxios_context)
+         ENDIF
 #endif
       ELSE
@@ -1786,5 +1871,5 @@
    END SUBROUTINE iom_rp3d_sp
 
-   SUBROUTINE iom_rp3d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype, ldxios )
+   SUBROUTINE iom_rp3d_dp( kt, kwrite, kiomid, cdvar, pvar, ktype )
       INTEGER         , INTENT(in)                         ::   kt       ! ocean time-step
       INTEGER         , INTENT(in)                         ::   kwrite   ! writing time-step
@@ -1793,16 +1878,26 @@
       REAL(dp)        , INTENT(in),       DIMENSION(:,:,:) ::   pvar     ! written field
       INTEGER         , INTENT(in), OPTIONAL               ::   ktype    ! variable external type
-      LOGICAL, OPTIONAL :: ldxios   ! xios write flag
-      LOGICAL :: llx                 ! local xios write flag
-      INTEGER :: ivid   ! variable id
-
-      llx = .FALSE.
-      IF(PRESENT(ldxios)) llx = ldxios
+      !
+      LOGICAL           :: llx                 ! local xios write flag
+      INTEGER           :: ivid   ! variable id
+      CHARACTER(LEN=lc) :: context
+      !
+      CALL set_xios_context(kiomid, context)
+
+      llx = .NOT. (context == "NONE")
+
       IF( llx ) THEN
 #ifdef key_iomput
-      IF( kt == kwrite ) THEN
-         IF(lwp) write(numout,*) 'RESTART: write (XIOS 3D) ',trim(cdvar)
-         CALL xios_send_field(trim(cdvar), pvar)
-      ENDIF
+         IF( kt == kwrite ) THEN
+            IF(lwp) write(numout,*) 'RESTART: write (XIOS 3D) ',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_put(trim(cdvar), pvar)
+            CALL iom_swap(cxios_context)
+         ELSE
+            IF(lwp) write(numout,*) 'RESTART: define (XIOS 3D)',trim(cdvar)
+            CALL iom_swap(context)
+            CALL iom_set_rstw_active( trim(cdvar), rd3 = pvar )
+            CALL iom_swap(cxios_context)
+         ENDIF
 #endif
       ELSE
@@ -2150,5 +2245,5 @@
       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( cxios_context )   ! return back to nemo context
    END SUBROUTINE iom_setkt
 
@@ -2164,5 +2259,5 @@
          CALL iom_swap( cdname )   ! swap to cdname context
          CALL xios_context_finalize() ! finalize the context
-         IF( cdname /= TRIM(cxios_context) ) CALL iom_swap( TRIM(cxios_context) )   ! return back to nemo context
+         IF( cdname /= cxios_context ) CALL iom_swap( cxios_context )   ! return back to nemo context
       ENDIF
       !
