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 8586)
+++ 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 )
