Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/namelist_ref	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/namelist_ref	(revision 11223)
@@ -609,5 +609,4 @@
    ln_vol        = .false.    !  total volume correction (see nn_volctl parameter)
    nn_volctl     =  1         !  = 0, the total water flux across open boundaries is zero
-   nb_jpk_bdy    = -1         ! number of levels in the bdy data (set < 0 if consistent with planned run)
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SPITZ12/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SPITZ12/EXPREF/namelist_cfg	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SPITZ12/EXPREF/namelist_cfg	(revision 11223)
@@ -177,5 +177,4 @@
    nn_rimwidth   = 1          !  width of the relaxation zone
    ln_vol        = .false.    !  total volume correction (see nn_volctl parameter)
-   nb_jpk_bdy    = -1         ! number of levels in the bdy data (set < 0 if consistent with planned run)
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90	(revision 11223)
@@ -73,9 +73,12 @@
    PUBLIC   ice_var_zapneg
    PUBLIC   ice_var_roundoff
-   PUBLIC   ice_var_itd
-   PUBLIC   ice_var_itd2
    PUBLIC   ice_var_bv           
    PUBLIC   ice_var_enthalpy           
    PUBLIC   ice_var_sshdyn
+   PUBLIC   ice_var_itd
+
+   INTERFACE ice_var_itd
+      MODULE PROCEDURE ice_var_itd_1c1c, ice_var_itd_Nc1c, ice_var_itd_1cMc, ice_var_itd_NcMc
+   END INTERFACE
 
    !!----------------------------------------------------------------------
@@ -656,10 +659,154 @@
    END SUBROUTINE ice_var_roundoff
    
+
+   SUBROUTINE ice_var_bv
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE ice_var_bv ***
+      !!
+      !! ** Purpose :   computes mean brine volume (%) in sea ice
+      !!
+      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
+      !!
+      !! References : Vancoppenolle et al., JGR, 2007
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
+      !!-------------------------------------------------------------------
+      !
+!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
+!!   instead of setting everything to zero as just below
+      bv_i (:,:,:) = 0._wp
+      DO jl = 1, jpl
+         DO jk = 1, nlay_i
+            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
+               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
+            END WHERE
+         END DO
+      END DO
+      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
+      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
+      END WHERE
+      !
+   END SUBROUTINE ice_var_bv
+
+
+   SUBROUTINE ice_var_enthalpy
+      !!-------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_var_enthalpy *** 
+      !!                 
+      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
+      !!
+      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ji, jk   ! dummy loop indices
+      REAL(wp) ::   ztmelts  ! local scalar 
+      !!-------------------------------------------------------------------
+      !
+      DO jk = 1, nlay_i             ! Sea ice energy of melting
+         DO ji = 1, npti
+            ztmelts      = - rTmlt  * sz_i_1d(ji,jk)
+            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue
+                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
+            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
+               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
+               &                   - rcp   * ztmelts )
+         END DO
+      END DO
+      DO jk = 1, nlay_s             ! Snow energy of melting
+         DO ji = 1, npti
+            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
+         END DO
+      END DO
+      !
+   END SUBROUTINE ice_var_enthalpy
+
    
-   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
-      !!-------------------------------------------------------------------
-      !!                ***  ROUTINE ice_var_itd   ***
-      !!
-      !! ** Purpose :  converting 1-cat ice to multiple ice categories
+   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_var_sshdyn  ***
+      !!                     
+      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
+      !!
+      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0
+      !!
+      !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
+      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*, 
+      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008
+      !!----------------------------------------------------------------------
+      !
+      ! input
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh            !: ssh [m]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass    !: mass of snow and ice at current  ice time step [Kg/m2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b  !: mass of snow and ice at previous ice time step [Kg/m2]
+      !
+      ! output
+      REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn  ! equivalent ssh in lead [m]
+      !
+      ! temporary
+      REAL(wp) :: zintn, zintb                     ! time interpolation weights []
+      REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m]
+      !
+      ! compute ice load used to define the equivalent ssh in lead
+      IF( ln_ice_embd ) THEN
+         !                                            
+         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
+         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1}
+         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
+         !
+         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
+         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
+         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
+         !
+         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
+         !
+      ELSE
+         zsnwiceload(:,:) = 0.0_wp
+      ENDIF
+      ! compute equivalent ssh in lead
+      ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
+      !
+   END FUNCTION ice_var_sshdyn
+
+   
+   !!-------------------------------------------------------------------
+   !!                ***  INTERFACE ice_var_itd   ***
+   !!
+   !! ** Purpose :  converting N-cat ice to jpl ice categories
+   !!-------------------------------------------------------------------
+   SUBROUTINE ice_var_itd_1c1c( zhti, zhts, zati, zh_i, zh_s, za_i )
+      !!-------------------------------------------------------------------
+      !! ** Purpose :  converting 1-cat ice to 1 ice category
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
+      REAL(wp), DIMENSION(:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
+      !!-------------------------------------------------------------------
+      zh_i(:) = zhti(:)
+      zh_s(:) = zhts(:)
+      za_i(:) = zati(:)
+   END SUBROUTINE ice_var_itd_1c1c
+
+   SUBROUTINE ice_var_itd_Nc1c( zhti, zhts, zati, zh_i, zh_s, za_i )
+      !!-------------------------------------------------------------------
+      !! ** Purpose :  converting N-cat ice to 1 ice category
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
+      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
+      !!-------------------------------------------------------------------
+      !
+      za_i(:) = SUM( zati(:,:), dim=2 )
+      !
+      WHERE( za_i(:) /= 0._wp )
+         zh_i(:) = SUM( zhti(:,:) * zati(:,:), dim=2 ) / za_i(:)
+         zh_s(:) = SUM( zhts(:,:) * zati(:,:), dim=2 ) / za_i(:)
+      ELSEWHERE
+         zh_i(:) = 0._wp
+         zh_s(:) = 0._wp
+      END WHERE
+      !
+   END SUBROUTINE ice_var_itd_Nc1c
+   
+   SUBROUTINE ice_var_itd_1cMc( zhti, zhts, zati, zh_i, zh_s, za_i )
+      !!-------------------------------------------------------------------
+      !!
+      !! ** Purpose :  converting 1-cat ice to jpl ice categories
       !!
       !!                  ice thickness distribution follows a gaussian law
@@ -705,4 +852,9 @@
       zh_s(1:idim,1:jpl) = 0._wp
       za_i(1:idim,1:jpl) = 0._wp
+      !
+      IF( jpl == 1 ) THEN
+         CALL ice_var_itd_1c1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) )
+         RETURN
+      ENDIF
       !
       DO ji = 1, idim
@@ -801,10 +953,8 @@
       END DO
       !
-   END SUBROUTINE ice_var_itd
-
-
-   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
-      !!-------------------------------------------------------------------
-      !!                ***  ROUTINE ice_var_itd2   ***
+   END SUBROUTINE ice_var_itd_1cMc
+
+   SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i )
+      !!-------------------------------------------------------------------
       !!
       !! ** Purpose :  converting N-cat ice to jpl ice categories
@@ -845,197 +995,112 @@
       idim = SIZE( zhti, 1 )
       icat = SIZE( zhti, 2 )
-      !
-      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
-      ALLOCATE( jlmin(idim), jlmax(idim) )
-
-      ! --- initialize output fields to 0 --- !
-      zh_i(1:idim,1:jpl) = 0._wp
-      zh_s(1:idim,1:jpl) = 0._wp
-      za_i(1:idim,1:jpl) = 0._wp
-      !
-      ! --- fill the categories --- !
-      !     find where cat-input = cat-output and fill cat-output fields  
-      jlmax(:) = 0
-      jlmin(:) = 999
-      jlfil(:,:) = 0
-      DO jl1 = 1, jpl
-         DO jl2 = 1, icat
+      !                                 ! ---------------------- !
+      IF( icat == jpl ) THEN            ! input cat = output cat !
+         !                              ! ---------------------- !
+         zh_i(:,:) = zhti(:,:)
+         zh_s(:,:) = zhts(:,:)
+         za_i(:,:) = zati(:,:)
+         !                              ! ---------------------- !
+      ELSEIF( icat == 1 ) THEN          ! specific case if N = 1 !
+         !                              ! ---------------------- !
+         !
+         CALL ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i, zh_s, za_i )
+         !
+         !                              ! ---------------------- !
+      ELSEIF( jpl == 1 ) THEN           ! specific case if M = 1 !
+         !                              ! ---------------------- !
+         !
+         CALL ice_var_itd_Nc1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) )
+         !
+         !                              ! ----------------------- !
+      ELSE                              ! input cat /= output cat !
+         !                              ! ----------------------- !
+         
+         ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
+         ALLOCATE( jlmin(idim), jlmax(idim) )
+
+         ! --- initialize output fields to 0 --- !
+         zh_i(1:idim,1:jpl) = 0._wp
+         zh_s(1:idim,1:jpl) = 0._wp
+         za_i(1:idim,1:jpl) = 0._wp
+         !
+         ! --- fill the categories --- !
+         !     find where cat-input = cat-output and fill cat-output fields  
+         jlmax(:) = 0
+         jlmin(:) = 999
+         jlfil(:,:) = 0
+         DO jl1 = 1, jpl
+            DO jl2 = 1, icat
+               DO ji = 1, idim
+                  IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
+                     ! fill the right category
+                     zh_i(ji,jl1) = zhti(ji,jl2)
+                     zh_s(ji,jl1) = zhts(ji,jl2)
+                     za_i(ji,jl1) = zati(ji,jl2)
+                     ! record categories that are filled
+                     jlmax(ji) = MAX( jlmax(ji), jl1 )
+                     jlmin(ji) = MIN( jlmin(ji), jl1 )
+                     jlfil(ji,jl1) = jl1
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         !
+         ! --- fill the gaps between categories --- !  
+         !     transfer from categories filled at the previous step to the empty ones in between
+         DO ji = 1, idim
+            jl1 = jlmin(ji)
+            jl2 = jlmax(ji)
+            IF( jl1 > 1 ) THEN
+               ! fill the lower cat (jl1-1)
+               za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
+               zh_i(ji,jl1-1) = hi_mean(jl1-1)
+               ! remove from cat jl1
+               za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
+            ENDIF
+            IF( jl2 < jpl ) THEN
+               ! fill the upper cat (jl2+1)
+               za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
+               zh_i(ji,jl2+1) = hi_mean(jl2+1)
+               ! remove from cat jl2
+               za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
+            ENDIF
+         END DO
+         !
+         jlfil2(:,:) = jlfil(:,:) 
+         ! fill categories from low to high
+         DO jl = 2, jpl-1
             DO ji = 1, idim
-               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
-                  ! fill the right category
-                  zh_i(ji,jl1) = zhti(ji,jl2)
-                  zh_s(ji,jl1) = zhts(ji,jl2)
-                  za_i(ji,jl1) = zati(ji,jl2)
-                  ! record categories that are filled
-                  jlmax(ji) = MAX( jlmax(ji), jl1 )
-                  jlmin(ji) = MIN( jlmin(ji), jl1 )
-                  jlfil(ji,jl1) = jl1
+               IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
+                  ! fill high
+                  za_i(ji,jl) = ztrans * za_i(ji,jl-1)
+                  zh_i(ji,jl) = hi_mean(jl)
+                  jlfil(ji,jl) = jl
+                  ! remove low
+                  za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
                ENDIF
             END DO
          END DO
-      END DO
-      !
-      ! --- fill the gaps between categories --- !  
-      !     transfer from categories filled at the previous step to the empty ones in between
-      DO ji = 1, idim
-         jl1 = jlmin(ji)
-         jl2 = jlmax(ji)
-         IF( jl1 > 1 ) THEN
-            ! fill the lower cat (jl1-1)
-            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
-            zh_i(ji,jl1-1) = hi_mean(jl1-1)
-            ! remove from cat jl1
-            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
-         ENDIF
-         IF( jl2 < jpl ) THEN
-            ! fill the upper cat (jl2+1)
-            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
-            zh_i(ji,jl2+1) = hi_mean(jl2+1)
-            ! remove from cat jl2
-            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
-         ENDIF
-      END DO
-      !
-      jlfil2(:,:) = jlfil(:,:) 
-      ! fill categories from low to high
-      DO jl = 2, jpl-1
-         DO ji = 1, idim
-            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
-               ! fill high
-               za_i(ji,jl) = ztrans * za_i(ji,jl-1)
-               zh_i(ji,jl) = hi_mean(jl)
-               jlfil(ji,jl) = jl
-               ! remove low
-               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
-            ENDIF
-         END DO
-      END DO
-      !
-      ! fill categories from high to low
-      DO jl = jpl-1, 2, -1
-         DO ji = 1, idim
-            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
-               ! fill low
-               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
-               zh_i(ji,jl) = hi_mean(jl) 
-               jlfil2(ji,jl) = jl
-               ! remove high
-               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
-            ENDIF
-         END DO
-      END DO
-      !
-      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
-      DEALLOCATE( jlmin, jlmax )
-      !
-   END SUBROUTINE ice_var_itd2
-
-
-   SUBROUTINE ice_var_bv
-      !!-------------------------------------------------------------------
-      !!                ***  ROUTINE ice_var_bv ***
-      !!
-      !! ** Purpose :   computes mean brine volume (%) in sea ice
-      !!
-      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
-      !!
-      !! References : Vancoppenolle et al., JGR, 2007
-      !!-------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      !!-------------------------------------------------------------------
-      !
-!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
-!!   instead of setting everything to zero as just below
-      bv_i (:,:,:) = 0._wp
-      DO jl = 1, jpl
-         DO jk = 1, nlay_i
-            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
-               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
-            END WHERE
-         END DO
-      END DO
-      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
-      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
-      END WHERE
-      !
-   END SUBROUTINE ice_var_bv
-
-
-   SUBROUTINE ice_var_enthalpy
-      !!-------------------------------------------------------------------
-      !!                   ***  ROUTINE ice_var_enthalpy *** 
-      !!                 
-      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
-      !!
-      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
-      !!-------------------------------------------------------------------
-      INTEGER  ::   ji, jk   ! dummy loop indices
-      REAL(wp) ::   ztmelts  ! local scalar 
-      !!-------------------------------------------------------------------
-      !
-      DO jk = 1, nlay_i             ! Sea ice energy of melting
-         DO ji = 1, npti
-            ztmelts      = - rTmlt  * sz_i_1d(ji,jk)
-            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue
-                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
-            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
-               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
-               &                   - rcp   * ztmelts )
-         END DO
-      END DO
-      DO jk = 1, nlay_s             ! Snow energy of melting
-         DO ji = 1, npti
-            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
-         END DO
-      END DO
-      !
-   END SUBROUTINE ice_var_enthalpy
-
-   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
-      !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE ice_var_sshdyn  ***
-      !!                     
-      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
-      !!
-      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0
-      !!
-      !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
-      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*, 
-      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008
-      !!----------------------------------------------------------------------
-      !
-      ! input
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh            !: ssh [m]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass    !: mass of snow and ice at current  ice time step [Kg/m2]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b  !: mass of snow and ice at previous ice time step [Kg/m2]
-      !
-      ! output
-      REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn  ! equivalent ssh in lead [m]
-      !
-      ! temporary
-      REAL(wp) :: zintn, zintb                     ! time interpolation weights []
-      REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m]
-      !
-      ! compute ice load used to define the equivalent ssh in lead
-      IF( ln_ice_embd ) THEN
-         !                                            
-         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
-         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1}
-         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
-         !
-         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
-         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
-         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
-         !
-         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
-         !
-      ELSE
-         zsnwiceload(:,:) = 0.0_wp
+         !
+         ! fill categories from high to low
+         DO jl = jpl-1, 2, -1
+            DO ji = 1, idim
+               IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
+                  ! fill low
+                  za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
+                  zh_i(ji,jl) = hi_mean(jl) 
+                  jlfil2(ji,jl) = jl
+                  ! remove high
+                  za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
+               ENDIF
+            END DO
+         END DO
+         !
+         DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
+         DEALLOCATE( jlmin, jlmax )
+         !
       ENDIF
-      ! compute equivalent ssh in lead
-      ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
-      !
-   END FUNCTION ice_var_sshdyn
-
+      !
+   END SUBROUTINE ice_var_itd_NcMc
 
 #else
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdy_oce.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdy_oce.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdy_oce.F90	(revision 11223)
@@ -42,12 +42,9 @@
    TYPE, PUBLIC ::   OBC_DATA     !: Storage for external data
       INTEGER          , DIMENSION(2)   ::  nread
-      LOGICAL                           ::  ll_ssh
-      LOGICAL                           ::  ll_u2d
-      LOGICAL                           ::  ll_v2d
-      LOGICAL                           ::  ll_u3d
-      LOGICAL                           ::  ll_v3d
-      LOGICAL                           ::  ll_tem
-      LOGICAL                           ::  ll_sal
-      LOGICAL                           ::  ll_fvl
+      LOGICAL                           ::  lneed_ssh
+      LOGICAL                           ::  lneed_dyn2d
+      LOGICAL                           ::  lneed_dyn3d
+      LOGICAL                           ::  lneed_tra
+      LOGICAL                           ::  lneed_ice
       REAL(wp), POINTER, DIMENSION(:)   ::  ssh
       REAL(wp), POINTER, DIMENSION(:)   ::  u2d
@@ -57,12 +54,7 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::  tem
       REAL(wp), POINTER, DIMENSION(:,:) ::  sal
-#if defined key_si3
-      LOGICAL                           ::   ll_a_i
-      LOGICAL                           ::   ll_h_i
-      LOGICAL                           ::   ll_h_s
-      REAL(wp), POINTER, DIMENSION(:,:) ::   a_i    !: now ice leads fraction climatology
-      REAL(wp), POINTER, DIMENSION(:,:) ::   h_i    !: Now ice  thickness climatology
-      REAL(wp), POINTER, DIMENSION(:,:) ::   h_s    !: now snow thickness
-#endif
+      REAL(wp), POINTER, DIMENSION(:,:) ::  a_i    !: now ice leads fraction climatology
+      REAL(wp), POINTER, DIMENSION(:,:) ::  h_i    !: Now ice  thickness climatology
+      REAL(wp), POINTER, DIMENSION(:,:) ::  h_s    !: now snow thickness
 #if defined key_top
       CHARACTER(LEN=20)                   :: cn_obc  !: type of boundary condition to apply
@@ -87,5 +79,4 @@
    !
    INTEGER                    ::   nb_bdy                   !: number of open boundary sets
-   INTEGER, DIMENSION(jp_bdy) ::   nb_jpk_bdy               !: number of levels in the bdy data (set < 0 if consistent with planned run)
    INTEGER, DIMENSION(jp_bdy) ::   nn_rimwidth              !: boundary rim width for Flow Relaxation Scheme
    INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P 
@@ -130,10 +121,4 @@
    INTEGER,  DIMENSION(jp_bdy)                     ::   nn_dta            !: =0 => *all* data is set to initial conditions
                                                                           !: =1 => some data to be read in from data files
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays (unstr.  bdy)
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_z      !: workspace for reading in global depth arrays (unstr.  bdy)
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_dz     !: workspace for reading in global depth arrays (unstr.  bdy)
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy)
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_z     !: workspace for reading in global depth arrays (struct. bdy)
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_dz    !: workspace for reading in global depth arrays (struct. bdy)
 !$AGRIF_DO_NOT_TREAT
    TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process)
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydta.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydta.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydta.F90	(revision 11223)
@@ -43,18 +43,19 @@
    PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90
 
-   INTEGER, ALLOCATABLE, DIMENSION(:)   ::   nb_bdy_fld        ! Number of fields to update for each boundary set.
-   INTEGER                              ::   nb_bdy_fld_sum    ! Total number of fields to update for all boundary sets.
-   LOGICAL,           DIMENSION(jp_bdy) ::   ln_full_vel_array ! =T => full velocities in 3D boundary conditions
-                                                               ! =F => baroclinic velocities in 3D boundary conditions
+   INTEGER , PARAMETER ::   jpbdyfld  = 10    ! maximum number of files to read 
+   INTEGER , PARAMETER ::   jp_bdyssh = 1     ! 
+   INTEGER , PARAMETER ::   jp_bdyu2d = 2     ! 
+   INTEGER , PARAMETER ::   jp_bdyv2d = 3     !
+   INTEGER , PARAMETER ::   jp_bdyu3d = 4     !
+   INTEGER , PARAMETER ::   jp_bdyv3d = 5     !
+   INTEGER , PARAMETER ::   jp_bdytem = 6     ! 
+   INTEGER , PARAMETER ::   jp_bdysal = 7     ! 
+   INTEGER , PARAMETER ::   jp_bdya_i = 8     ! 
+   INTEGER , PARAMETER ::   jp_bdyh_i = 9     ! 
+   INTEGER , PARAMETER ::   jp_bdyh_S = 10    ! 
+                                                              ! =F => baroclinic velocities in 3D boundary conditions
 !$AGRIF_DO_NOT_TREAT
-   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:), TARGET ::   bf        ! structure of input fields (file informations, fields read)
+   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:,:), TARGET ::   bf   ! structure of input fields (file informations, fields read)
 !$AGRIF_END_DO_NOT_TREAT
-   TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr   ! array of pointers to nbmap
-
-#if defined key_si3
-   INTEGER ::   nice_cat                      ! number of categories in the input file
-   INTEGER ::   jfld_hti, jfld_hts, jfld_ai   ! indices of ice thickness, snow thickness and concentration in bf structure
-   INTEGER, DIMENSION(jp_bdy) :: jfld_htit, jfld_htst, jfld_ait
-#endif
 
    !!----------------------------------------------------------------------
@@ -65,5 +66,5 @@
 CONTAINS
 
-      SUBROUTINE bdy_dta( kt, jit, time_offset )
+   SUBROUTINE bdy_dta( kt, kit, kt_offset )
       !!----------------------------------------------------------------------
       !!                   ***  SUBROUTINE bdy_dta  ***
@@ -75,17 +76,18 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in)           ::   kt           ! ocean time-step index 
-      INTEGER, INTENT(in), OPTIONAL ::   jit          ! subcycle time-step index (for timesplitting option)
-      INTEGER, INTENT(in), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit
+      INTEGER, INTENT(in), OPTIONAL ::   kit          ! subcycle time-step index (for timesplitting option)
+      INTEGER, INTENT(in), OPTIONAL ::   kt_offset    ! time offset in units of timesteps. NB. if kit
       !                                               ! is present then units = subcycle timesteps.
-      !                                               ! time_offset = 0 => get data at "now" time level
-      !                                               ! time_offset = -1 => get data at "before" time level
-      !                                               ! time_offset = +1 => get data at "after" time level
+      !                                               ! kt_offset = 0 => get data at "now" time level
+      !                                               ! kt_offset = -1 => get data at "before" time level
+      !                                               ! kt_offset = +1 => get data at "after" time level
       !                                               ! etc.
       !
-      INTEGER ::  jbdy, jfld, jstart, jend, ib, jl  ! dummy loop indices
-      INTEGER ::  ii, ij, ik, igrd                  ! local integers
-      INTEGER,          DIMENSION(jpbgrd) ::   ilen1 
-      INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts
-      TYPE(OBC_DATA), POINTER             ::   dta              ! short cut
+      INTEGER ::  jbdy, jfld, jstart, jend, ib, jl    ! dummy loop indices
+      INTEGER ::  ii, ij, ik, igrd, ipl               ! local integers
+      INTEGER,   DIMENSION(jpbgrd)     ::   ilen1 
+      INTEGER,   DIMENSION(:), POINTER ::   nblen, nblenrim  ! short cuts
+      TYPE(OBC_DATA)         , POINTER ::   dta_alias        ! short cut
+      TYPE(FLD), DIMENSION(:), POINTER ::   bf_alias
       !!---------------------------------------------------------------------------
       !
@@ -94,18 +96,17 @@
       ! Initialise data arrays once for all from initial conditions where required
       !---------------------------------------------------------------------------
-      IF( kt == nit000 .AND. .NOT.PRESENT(jit) ) THEN
+      IF( kt == nit000 .AND. .NOT.PRESENT(kit) ) THEN
 
          ! Calculate depth-mean currents
          !-----------------------------
-         
+
          DO jbdy = 1, nb_bdy
             !
             nblen    => idx_bdy(jbdy)%nblen
             nblenrim => idx_bdy(jbdy)%nblenrim
-            dta      => dta_bdy(jbdy)
             !
             IF( nn_dyn2d_dta(jbdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
-               IF( dta%ll_ssh ) THEN 
+               IF( dta_bdy(jbdy)%lneed_ssh ) THEN 
                   igrd = 1
                   DO ib = 1, ilen1(igrd)
@@ -113,7 +114,7 @@
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
                      dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
-                  END DO 
-               ENDIF
-               IF( dta%ll_u2d ) THEN 
+                  END DO
+               ENDIF
+               IF( dta_bdy(jbdy)%lneed_dyn2d) THEN 
                   igrd = 2
                   DO ib = 1, ilen1(igrd)
@@ -121,7 +122,5 @@
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
                      dta_bdy(jbdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)         
-                  END DO 
-               ENDIF
-               IF( dta%ll_v2d ) THEN 
+                  END DO
                   igrd = 3
                   DO ib = 1, ilen1(igrd)
@@ -129,5 +128,5 @@
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
                      dta_bdy(jbdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)         
-                  END DO 
+                  END DO
                ENDIF
             ENDIF
@@ -135,5 +134,5 @@
             IF( nn_dyn3d_dta(jbdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
-               IF( dta%ll_u3d ) THEN 
+               IF( dta_bdy(jbdy)%lneed_dyn3d ) THEN 
                   igrd = 2 
                   DO ib = 1, ilen1(igrd)
@@ -143,7 +142,5 @@
                         dta_bdy(jbdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)         
                      END DO
-                  END DO 
-               ENDIF
-               IF( dta%ll_v3d ) THEN 
+                  END DO
                   igrd = 3 
                   DO ib = 1, ilen1(igrd)
@@ -152,6 +149,6 @@
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
                         dta_bdy(jbdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)         
-                        END DO
-                  END DO 
+                     END DO
+                  END DO
                ENDIF
             ENDIF
@@ -159,5 +156,5 @@
             IF( nn_tra_dta(jbdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
-               IF( dta%ll_tem ) THEN
+               IF( dta_bdy(jbdy)%lneed_tra ) THEN
                   igrd = 1 
                   DO ib = 1, ilen1(igrd)
@@ -165,17 +162,8 @@
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                        dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)         
+                        dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_bdytem) * tmask(ii,ij,ik)         
+                        dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_bdysal) * tmask(ii,ij,ik)         
                      END DO
-                  END DO 
-               ENDIF
-               IF( dta%ll_sal ) THEN
-                  igrd = 1 
-                  DO ib = 1, ilen1(igrd)
-                     DO ik = 1, jpkm1
-                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
-                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                        dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)         
-                     END DO
-                  END DO 
+                  END DO
                ENDIF
             ENDIF
@@ -184,5 +172,5 @@
             IF( nn_ice_dta(jbdy) == 0 ) THEN    ! set ice to initial values
                ilen1(:) = nblen(:)
-               IF( dta%ll_a_i ) THEN
+               IF( dta_bdy(jbdy)%lneed_ice ) THEN
                   igrd = 1   
                   DO jl = 1, jpl
@@ -191,23 +179,5 @@
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
                         dta_bdy(jbdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1) 
-                     END DO
-                  END DO
-               ENDIF
-               IF( dta%ll_h_i ) THEN
-                  igrd = 1   
-                  DO jl = 1, jpl
-                     DO ib = 1, ilen1(igrd)
-                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
-                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
                         dta_bdy(jbdy)%h_i (ib,jl) =  h_i(ii,ij,jl) * tmask(ii,ij,1) 
-                     END DO
-                  END DO
-               ENDIF
-               IF( dta%ll_h_s ) THEN
-                  igrd = 1   
-                  DO jl = 1, jpl
-                     DO ib = 1, ilen1(igrd)
-                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
-                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
                         dta_bdy(jbdy)%h_s (ib,jl) =  h_s(ii,ij,jl) * tmask(ii,ij,1) 
                      END DO
@@ -222,188 +192,121 @@
       ! update external data from files
       !--------------------------------
-     
-      jstart = 1
-      DO jbdy = 1, nb_bdy   
-         dta => dta_bdy(jbdy)
-         IF( nn_dta(jbdy) == 1 ) THEN ! skip this bit if no external data required
-      
-            IF( PRESENT(jit) ) THEN
-               ! Update barotropic boundary conditions only
-               ! jit is optional argument for fld_read and bdytide_update
-               IF( cn_dyn2d(jbdy) /= 'none' ) THEN
-                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
-                     IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
-                     IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
-                     IF( dta%ll_u3d ) dta%v2d(:) = 0._wp
-                  ENDIF
-                  IF (cn_tra(jbdy) /= 'runoff') THEN
-                     IF( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) THEN
-
-                        jend = jstart + dta%nread(2) - 1
-                        IF( ln_full_vel_array(jbdy) ) THEN
-                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
-                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy(jbdy),   &
-                                     & fvl=ln_full_vel_array(jbdy)  )
-                        ELSE
-                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
-                                     & kit=jit, kt_offset=time_offset  )
-                        ENDIF
-
-                        ! If full velocities in boundary data then extract barotropic velocities from 3D fields
-                        IF( ln_full_vel_array(jbdy) .AND.                                             &
-                          &    ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR.  &
-                          &      nn_dyn3d_dta(jbdy) == 1 ) )THEN
-
-                           igrd = 2                      ! zonal velocity
-                           dta%u2d(:) = 0._wp
-                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
-                              ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-                              ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-                              DO ik = 1, jpkm1
-                                 dta%u2d(ib) = dta%u2d(ib) &
-                       &                          + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
-                              END DO
-                              dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij)
-                           END DO
-                           igrd = 3                      ! meridional velocity
-                           dta%v2d(:) = 0._wp
-                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
-                              ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-                              ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-                              DO ik = 1, jpkm1
-                                 dta%v2d(ib) = dta%v2d(ib) &
-                       &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
-                              END DO
-                              dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij)
-                           END DO
-                        ENDIF                    
-                     ENDIF
-                     IF( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN ! update tidal harmonic forcing
-                        CALL bdytide_update( kt=kt, idx=idx_bdy(jbdy), dta=dta, td=tides(jbdy),   & 
-                          &                 jit=jit, time_offset=time_offset )
-                     ENDIF
-                  ENDIF
-               ENDIF
-            ELSE
-               IF (cn_tra(jbdy) == 'runoff') then      ! runoff condition
-                  jend = nb_bdy_fld(jbdy)
-                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  &
-                               & map=nbmap_ptr(jstart:jend), kt_offset=time_offset )
-                  !
-                  igrd = 2                      ! zonal velocity
-                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
-                     ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-                     ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-                     dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
-                  END DO
-                  !
-                  igrd = 3                      ! meridional velocity
-                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
-                     ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-                     ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-                     dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
-                  END DO
-               ELSE
-                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
-                     IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
-                     IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
-                     IF( dta%ll_v2d ) dta%v2d(:) = 0._wp
-                  ENDIF
-                  IF( dta%nread(1) .gt. 0 ) THEN ! update external data
-                     jend = jstart + dta%nread(1) - 1
-                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), &
-                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy(jbdy),   &
-                                  & fvl=ln_full_vel_array(jbdy) )
-                  ENDIF
-                  ! If full velocities in boundary data then split into barotropic and baroclinic data
-                  IF( ln_full_vel_array(jbdy) .and.                                             &
-                    & ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR. &
-                    &   nn_dyn3d_dta(jbdy) == 1 ) ) THEN
-                     igrd = 2                      ! zonal velocity
-                     dta%u2d(:) = 0._wp
-                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
-                        ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-                        ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-                        DO ik = 1, jpkm1
-                           dta%u2d(ib) = dta%u2d(ib) &
-                 &                       + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
-                        END DO
-                        dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij)
-                        DO ik = 1, jpkm1
-                           dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib)
-                        END DO
-                     END DO
-                     igrd = 3                      ! meridional velocity
-                     dta%v2d(:) = 0._wp
-                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
-                        ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-                        ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-                        DO ik = 1, jpkm1
-                           dta%v2d(ib) = dta%v2d(ib) &
-                 &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
-                        END DO
-                        dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij)
-                        DO ik = 1, jpkm1
-                           dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib)
-                        END DO
-                     END DO
-                  ENDIF
-
-               ENDIF
+
+      DO jbdy = 1, nb_bdy
+
+         dta_alias => dta_bdy(jbdy)
+         bf_alias  => bf(:,jbdy)
+
+         ! read/update all bdy data
+         ! ------------------------
+         CALL fld_read( kt, 1, bf_alias, kit = kit, kt_offset = kt_offset )
+
+         ! apply some corrections in some specific cases...
+         ! --------------------------------------------------
+         !
+         ! if runoff condition: change river flow we read (in m3/s) into barotropic velocity (m/s)
+         IF( cn_tra(jbdy) == 'runoff' .AND. TRIM(bf_alias(jp_bdyu2d)%clrootname) /= 'NOT USED' ) THEN   ! runoff and we read u/v2d
+            !
+            igrd = 2                      ! zonal flow (m3/s) to barotropic zonal velocity (m/s)
+            DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+               dta_alias%u2d(ib) = dta_alias%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
+            END DO
+            igrd = 3                      ! meridional flow (m3/s) to barotropic meridional velocity (m/s)
+            DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+               dta_alias%v2d(ib) = dta_alias%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
+            END DO
+         ENDIF
+
+         ! tidal harmonic forcing ONLY: initialise arrays
+         IF( nn_dyn2d_dta(jbdy) == 2 ) THEN   ! we did not read ssh, u/v2d 
+            IF( dta_alias%lneed_ssh   ) dta_alias%ssh(:) = 0._wp
+            IF( dta_alias%lneed_dyn2d ) dta_alias%u2d(:) = 0._wp
+            IF( dta_alias%lneed_dyn2d ) dta_alias%v2d(:) = 0._wp
+         ENDIF
+
+         ! If full velocities in boundary data, then split it into barotropic and baroclinic component
+         IF( bf_alias(jp_bdyu3d)%ltotvel ) THEN     ! if we read 3D total velocity (can be true only if u3d was read)
+            !
+            igrd = 2                       ! zonal velocity
+            dta_alias%u2d(:) = 0._wp       ! compute barotrope zonal velocity and put it in u2d
+            DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+               DO ik = 1, jpkm1
+                  dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik)
+               END DO
+               dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu_n(ii,ij)
+               DO ik = 1, jpkm1            ! compute barocline zonal velocity and put it in u3d
+                  dta_alias%u3d(ib,ik) = dta_alias%u3d(ib,ik) - dta_alias%u2d(ib)
+               END DO
+            END DO
+            igrd = 3                       ! meridional velocity
+            dta_alias%v2d(:) = 0._wp       ! compute barotrope meridional velocity and put it in v2d
+            DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+               DO ik = 1, jpkm1
+                  dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik)
+               END DO
+               dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv_n(ii,ij)
+               DO ik = 1, jpkm1            ! compute barocline meridional velocity and put it in v3d
+                  dta_alias%v3d(ib,ik) = dta_alias%v3d(ib,ik) - dta_alias%v2d(ib)
+               END DO
+            END DO
+         ENDIF   ! ltotvel
+
+         ! update tidal harmonic forcing
+         IF( PRESENT(kit) .AND. nn_dyn2d_dta(jbdy) .GE. 2 ) THEN
+            CALL bdytide_update( kt = kt, idx = idx_bdy(jbdy), dta = dta_alias, td = tides(jbdy),   & 
+               &                 kit = kit, kt_offset = kt_offset )
+         ENDIF
+
+         !  atm surface pressure : add inverted barometer effect to ssh if it was read
+         IF ( ln_apr_obc .AND. TRIM(bf_alias(jp_bdyssh)%clrootname) /= 'NOT USED' ) THEN
+            igrd = 1
+            DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd)   ! ssh is used only on the rim
+               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+               dta_alias%ssh(ib) = dta_alias%ssh(ib) + ssh_ib(ii,ij)
+            END DO
+         ENDIF
+
 #if defined key_si3
-               ! convert N-cat fields (input) into jpl-cat (output)
-               IF( cn_ice(jbdy) /= 'none' .AND. nn_ice_dta(jbdy) == 1 ) THEN
-                  jfld_hti = jfld_htit(jbdy)
-                  jfld_hts = jfld_htst(jbdy)
-                  jfld_ai  = jfld_ait(jbdy)
-                  IF    ( jpl /= 1 .AND. nice_cat == 1 ) THEN                       ! case input cat = 1
-                     CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
-                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    )
-                  ELSEIF( jpl /= 1 .AND. nice_cat /= 1 .AND. nice_cat /= jpl ) THEN ! case input cat /=1 and /=jpl
-                     CALL ice_var_itd2( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), &
-                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    )
-                  ENDIF
-               ENDIF
+         ! ice: convert N-cat fields (input) into jpl-cat (output)
+         IF( dta_alias%lneed_ice ) THEN
+            ipl = SIZE(bf_alias(jp_bdya_i)%fnow, 3)
+            IF( ipl /= jpl ) THEN      ! ice: convert N-cat fields (input) into jpl-cat (output)
+               CALL ice_var_itd(bf_alias(jp_bdyh_i)%fnow(:,1,:), bf_alias(jp_bdyh_s)%fnow(:,1,:), bf_alias(jp_bdya_i)%fnow(:,1,:), &
+                  &              dta_alias%h_i               , dta_alias%h_s               , dta_alias%a_i                 )
+            ENDIF
+         ENDIF
 #endif
-            ENDIF
-            jstart = jstart + dta%nread(1)
-         ENDIF    ! nn_dta(jbdy) = 1
       END DO  ! jbdy
-
-      IF ( ln_apr_obc ) THEN
-         DO jbdy = 1, nb_bdy
-            IF (cn_tra(jbdy) /= 'runoff')THEN
-               igrd = 1                      ! meridional velocity
-               DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd)
-                  ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-                  ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-                  dta_bdy(jbdy)%ssh(ib) = dta_bdy(jbdy)%ssh(ib) + ssh_ib(ii,ij)
-               END DO
-            ENDIF
-         END DO
-      ENDIF
 
       IF ( ln_tide ) THEN
          IF (ln_dynspg_ts) THEN      ! Fill temporary arrays with slow-varying bdy data                           
-            DO jbdy = 1, nb_bdy    ! Tidal component added in ts loop
-               IF ( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN
+            DO jbdy = 1, nb_bdy      ! Tidal component added in ts loop
+               IF ( nn_dyn2d_dta(jbdy) .GE. 2 ) THEN
                   nblen => idx_bdy(jbdy)%nblen
                   nblenrim => idx_bdy(jbdy)%nblenrim
-                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF 
-                  IF ( dta_bdy(jbdy)%ll_ssh ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1))
-                  IF ( dta_bdy(jbdy)%ll_u2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2))
-                  IF ( dta_bdy(jbdy)%ll_v2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3))
-               ENDIF
-            END DO
-         ELSE ! Add tides if not split-explicit free surface else this is done in ts loop
-            !
-            CALL bdy_dta_tides( kt=kt, time_offset=time_offset )
-         ENDIF
-      ENDIF
-
-      !
-      IF( ln_timing )   CALL timing_stop('bdy_dta')
-      !
-   END SUBROUTINE bdy_dta
+                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN ; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF 
+                     IF ( dta_bdy(jbdy)%lneed_ssh   ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1))
+                     IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2))
+                     IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3))
+                  ENDIF
+               END DO
+            ELSE ! Add tides if not split-explicit free surface else this is done in ts loop
+               !
+               CALL bdy_dta_tides( kt=kt, kt_offset=kt_offset )
+            ENDIF
+         ENDIF
+         !
+         IF( ln_timing )   CALL timing_stop('bdy_dta')
+         !
+      END SUBROUTINE bdy_dta
 
 
@@ -418,32 +321,28 @@
       !!                
       !!----------------------------------------------------------------------
-      INTEGER ::   jbdy, jfld, jstart, jend, ierror, ios     ! Local integers
-      !
+      INTEGER ::   jbdy, jfld    ! Local integers
+      INTEGER ::   ierror, ios     ! 
+      !
+      CHARACTER(len=3)                       ::   cl3           ! 
       CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files
-      CHARACTER(len=100), DIMENSION(nb_bdy)  ::   cn_dir_array  ! Root directory for location of data files
-      CHARACTER(len = 256)::   clname                           ! temporary file name
       LOGICAL                                ::   ln_full_vel   ! =T => full velocities in 3D boundary data
       !                                                         ! =F => baroclinic velocities in 3D boundary data
-      INTEGER                                ::   ilen_global   ! Max length required for global bdy dta arrays
-      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ilen1, ilen3  ! size of 1st and 3rd dimensions of local arrays
-      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ibdy           ! bdy set for a particular jfld
-      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   igrid         ! index for grid type (1,2,3 = T,U,V)
-      INTEGER, POINTER, DIMENSION(:)         ::   nblen, nblenrim  ! short cuts
-      TYPE(OBC_DATA), POINTER                ::   dta           ! short cut
-#if defined key_si3
-      INTEGER               ::   kndims   ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat)
-      INTEGER, DIMENSION(4) ::   kdimsz   ! size   of dimensions
-      INTEGER               ::   inum,id1 ! local integer
-#endif
-      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   blf_i         !  array of namelist information structures
-      TYPE(FLD_N) ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! 
-      TYPE(FLD_N) ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read
-#if defined key_si3
-      TYPE(FLD_N) ::   bn_a_i, bn_h_i, bn_h_s      
-#endif
+      INTEGER                                ::   ipk,ipl           !
+      INTEGER                                ::   idvar         ! variable ID
+      INTEGER                                ::   indims        ! number of dimensions of the variable
+      INTEGER                                ::   iszdim        ! number of dimensions of the variable
+      INTEGER, DIMENSION(4)                  ::   i4dimsz       ! size of variable dimensions 
+      INTEGER                                ::   igrd          ! index for grid type (1,2,3 = T,U,V)
+      LOGICAL                                ::   lluld         ! is the variable using the unlimited dimension
+      LOGICAL                                ::   llneed        !
+      LOGICAL                                ::   llread        !
+      TYPE(FLD_N), DIMENSION(1), TARGET ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill
+      TYPE(FLD_N), DIMENSION(1), TARGET ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read
+      TYPE(FLD_N), DIMENSION(1), TARGET ::   bn_a_i, bn_h_i, bn_h_s      
+      TYPE(FLD_N), DIMENSION(:), POINTER ::   bn_alias                        ! must be an array to be used with fld_fill
+      TYPE(FLD  ), DIMENSION(:), POINTER ::   bf_alias
+      !
       NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 
-#if defined key_si3
       NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s
-#endif
       NAMELIST/nambdy_dta/ ln_full_vel
       !!---------------------------------------------------------------------------
@@ -454,389 +353,205 @@
       IF(lwp) WRITE(numout,*) ''
 
-      ! Set nn_dta
-      DO jbdy = 1, nb_bdy
-         nn_dta(jbdy) = MAX(   nn_dyn2d_dta  (jbdy)    &
-            &                , nn_dyn3d_dta  (jbdy)    &
-            &                , nn_tra_dta    (jbdy)    &
-#if defined key_si3
-            &                , nn_ice_dta    (jbdy)    &
-#endif
-                              )
-         IF(nn_dta(jbdy) > 1)   nn_dta(jbdy) = 1
-      END DO
-
-      ! Work out upper bound of how many fields there are to read in and allocate arrays
-      ! ---------------------------------------------------------------------------
-      ALLOCATE( nb_bdy_fld(nb_bdy) )
-      nb_bdy_fld(:) = 0
-      DO jbdy = 1, nb_bdy         
-         IF( cn_dyn2d(jbdy) /= 'none' .AND. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) THEN
-            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3
-         ENDIF
-         IF( cn_dyn3d(jbdy) /= 'none' .AND. nn_dyn3d_dta(jbdy) == 1 ) THEN
-            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2
-         ENDIF
-         IF( cn_tra(jbdy) /= 'none' .AND. nn_tra_dta(jbdy) == 1  ) THEN
-            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2
-         ENDIF
-#if defined key_si3
-         IF( cn_ice(jbdy) /= 'none' .AND. nn_ice_dta(jbdy) == 1  ) THEN
-            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3
-         ENDIF
-#endif               
-         IF(lwp) WRITE(numout,*) 'Maximum number of files to open =', nb_bdy_fld(jbdy)
-      END DO            
-
-      nb_bdy_fld_sum = SUM( nb_bdy_fld )
-
-      ALLOCATE( bf(nb_bdy_fld_sum), STAT=ierror )
+      ALLOCATE( bf(jpbdyfld,nb_bdy), STAT=ierror )
       IF( ierror > 0 ) THEN   
          CALL ctl_stop( 'bdy_dta: unable to allocate bf structure' )   ;   RETURN  
       ENDIF
-      ALLOCATE( blf_i(nb_bdy_fld_sum), STAT=ierror )
-      IF( ierror > 0 ) THEN   
-         CALL ctl_stop( 'bdy_dta: unable to allocate blf_i structure' )   ;   RETURN  
-      ENDIF
-      ALLOCATE( nbmap_ptr(nb_bdy_fld_sum), STAT=ierror )
-      IF( ierror > 0 ) THEN   
-         CALL ctl_stop( 'bdy_dta: unable to allocate nbmap_ptr structure' )   ;   RETURN  
-      ENDIF
-      ALLOCATE( ilen1(nb_bdy_fld_sum), ilen3(nb_bdy_fld_sum) ) 
-      ALLOCATE( ibdy(nb_bdy_fld_sum) ) 
-      ALLOCATE( igrid(nb_bdy_fld_sum) ) 
-
+      bf(:,:)%clrootname = 'NOT USED'   ! default definition used as a flag in fld_read to do nothing.
+      bf(:,:)%ltotvel = .FALSE.         ! default definition
+ 
       ! Read namelists
       ! --------------
       REWIND(numnam_cfg)
-      jfld = 0 
-      DO jbdy = 1, nb_bdy         
-         IF( nn_dta(jbdy) == 1 ) THEN
-            REWIND(numnam_ref)
-            READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
-901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp )
+      DO jbdy = 1, nb_bdy
+
+         WRITE(ctmp1, '(a,i2)') 'BDY number ', jbdy
+         WRITE(ctmp2, '(a,i2)') 'block nambdy_dta number ', jbdy
+
+         ! There is only one nambdy_dta block in namelist_ref -> use it for each bdy so we do a rewind 
+         REWIND(numnam_ref)
+         READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
+901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp )
+
+         !   by-pass nambdy_dta reading if no input data used in this bdy   
+         IF(       ( dta_bdy(jbdy)%lneed_dyn2d .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1 )   &
+            & .OR. ( dta_bdy(jbdy)%lneed_dyn3d .AND.     nn_dyn3d_dta(jbdy)    == 1 )   &
+            & .OR. ( dta_bdy(jbdy)%lneed_tra   .AND.       nn_tra_dta(jbdy)    == 1 )   &
+            & .OR. ( dta_bdy(jbdy)%lneed_ice   .AND.       nn_ice_dta(jbdy)    == 1 )   )   THEN
+            ! WARNING: we don't do a rewind here, each bdy reads its own nambdy_dta block one after another
             READ  ( numnam_cfg, nambdy_dta, IOSTAT = ios, ERR = 902 )
 902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_dta in configuration namelist', lwp )
-            IF(lwm) WRITE( numond, nambdy_dta )
-
-            cn_dir_array(jbdy) = cn_dir
-            ln_full_vel_array(jbdy) = ln_full_vel
-
-            nblen => idx_bdy(jbdy)%nblen
-            nblenrim => idx_bdy(jbdy)%nblenrim
-            dta => dta_bdy(jbdy)
-            dta%nread(2) = 0
-
-            ! Only read in necessary fields for this set.
-            ! Important that barotropic variables come first.
-            IF( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN 
-
-               IF( dta%ll_ssh ) THEN 
-                  if(lwp) write(numout,*) '++++++ reading in ssh field'
-                  jfld = jfld + 1
-                  blf_i(jfld) = bn_ssh
-                  ibdy(jfld) = jbdy
-                  igrid(jfld) = 1
-                  ilen1(jfld) = nblen(igrid(jfld))
-                  ilen3(jfld) = 1
-                  dta%nread(2) = dta%nread(2) + 1
-               ENDIF
-
-               IF( dta%ll_u2d .and. .not. ln_full_vel_array(jbdy) ) THEN
-                  if(lwp) write(numout,*) '++++++ reading in u2d field'
-                  jfld = jfld + 1
-                  blf_i(jfld) = bn_u2d
-                  ibdy(jfld) = jbdy
-                  igrid(jfld) = 2
-                  ilen1(jfld) = nblen(igrid(jfld))
-                  ilen3(jfld) = 1
-                  dta%nread(2) = dta%nread(2) + 1
-               ENDIF
-
-               IF( dta%ll_v2d .and. .not. ln_full_vel_array(jbdy) ) THEN
-                  if(lwp) write(numout,*) '++++++ reading in v2d field'
-                  jfld = jfld + 1
-                  blf_i(jfld) = bn_v2d
-                  ibdy(jfld) = jbdy
-                  igrid(jfld) = 3
-                  ilen1(jfld) = nblen(igrid(jfld))
-                  ilen3(jfld) = 1
-                  dta%nread(2) = dta%nread(2) + 1
-               ENDIF
-
-            ENDIF
-
-            ! read 3D velocities if baroclinic velocities require OR if
-            ! barotropic velocities required and ln_full_vel set to .true.
-            IF( nn_dyn3d_dta(jbdy) == 1 .OR. &
-           &  ( ln_full_vel_array(jbdy) .AND. ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN
-
-               IF( dta%ll_u3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN
-                  if(lwp) write(numout,*) '++++++ reading in u3d field'
-                  jfld = jfld + 1
-                  blf_i(jfld) = bn_u3d
-                  ibdy(jfld) = jbdy
-                  igrid(jfld) = 2
-                  ilen1(jfld) = nblen(igrid(jfld))
-                  ilen3(jfld) = jpk
-                  IF( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1
-               ENDIF
-
-               IF( dta%ll_v3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN
-                  if(lwp) write(numout,*) '++++++ reading in v3d field'
-                  jfld = jfld + 1
-                  blf_i(jfld) = bn_v3d
-                  ibdy(jfld) = jbdy
-                  igrid(jfld) = 3
-                  ilen1(jfld) = nblen(igrid(jfld))
-                  ilen3(jfld) = jpk
-                  IF( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1
-               ENDIF
-
-            ENDIF
-
-            ! temperature and salinity
-            IF( nn_tra_dta(jbdy) == 1 ) THEN
-
-               IF( dta%ll_tem ) THEN
-                  if(lwp) write(numout,*) '++++++ reading in tem field'
-                  jfld = jfld + 1
-                  blf_i(jfld) = bn_tem
-                  ibdy(jfld) = jbdy
-                  igrid(jfld) = 1
-                  ilen1(jfld) = nblen(igrid(jfld))
-                  ilen3(jfld) = jpk
-               ENDIF
-
-               IF( dta%ll_sal ) THEN
-                  if(lwp) write(numout,*) '++++++ reading in sal field'
-                  jfld = jfld + 1
-                  blf_i(jfld) = bn_sal
-                  ibdy(jfld) = jbdy
-                  igrid(jfld) = 1
-                  ilen1(jfld) = nblen(igrid(jfld))
-                  ilen3(jfld) = jpk
-               ENDIF
-
-            ENDIF
-
-#if defined key_si3
-            ! sea ice
-            IF( nn_ice_dta(jbdy) == 1 ) THEN
-               ! Test for types of ice input (1cat or Xcat) 
-               ! Build file name to find dimensions 
-               clname=TRIM( cn_dir )//TRIM(bn_a_i%clname)
-               IF( .NOT. bn_a_i%ln_clim ) THEN   
-                                                  WRITE(clname, '(a,"_y",i4.4)' ) TRIM( clname ), nyear    ! add year
-                  IF( bn_a_i%cltype /= 'yearly' ) WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname ), nmonth   ! add month
-               ELSE
-                  IF( bn_a_i%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( clname ), nmonth   ! add month
-               ENDIF
-               IF( bn_a_i%cltype == 'daily' .OR. bn_a_i%cltype(1:4) == 'week' ) &
-               &                                  WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), nday     ! add day
+            IF(lwm) WRITE( numond, nambdy_dta )           
+         ENDIF
+
+         ! get the number of ice categories in bdy data file (use a_i information to do this)
+         ipl = jpl   ! default definition
+         IF( dta_bdy(jbdy)%lneed_ice ) THEN    ! if we need ice bdy data
+            IF( nn_ice_dta(jbdy) == 1 ) THEN   ! if we get ice bdy data from netcdf file
+               CALL fld_fill(  bf(jp_bdya_i,jbdy:jbdy), bn_a_i, cn_dir, 'bdy_dta', 'a_i'//' '//ctmp1, ctmp2 )   ! use namelist info
+               CALL fld_clopn( bf(jp_bdya_i,jbdy), nyear, nmonth, nday )   ! not a problem when we call it again after
+               idvar = iom_varid( bf(jp_bdya_i,jbdy)%num, bf(jp_bdya_i,jbdy)%clvar, kndims=indims, kdimsz=i4dimsz, lduld=lluld )
+               IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN   ;   ipl = i4dimsz(3)   ! xylt or xyl
+               ELSE                                                            ;   ipl = 1            ! xy or xyt
+               ENDIF
+            ENDIF
+         ENDIF
+
+         DO jfld = 1, jpbdyfld
+
+            ! =====================
+            !          ssh 
+            ! =====================
+            IF( jfld == jp_bdyssh ) THEN
+               cl3 = 'ssh'
+               igrd = 1                                                    ! T point
+               ipk = 1                                                     ! surface data
+               llneed = dta_bdy(jbdy)%lneed_ssh                            ! dta_bdy(jbdy)%ssh will be needed
+               llread = MOD(nn_dyn2d_dta(jbdy),2) == 1                     ! get data from NetCDF file
+               bf_alias => bf(jp_bdyssh,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
+               bn_alias => bn_ssh                                          ! alias for ssh structure of nambdy_dta 
+            ENDIF
+            ! =====================
+            !         dyn2d
+            ! =====================
+            IF( jfld == jp_bdyu2d ) THEN
+               cl3 = 'u2d'
+               igrd = 2                                                    ! U point
+               ipk = 1                                                     ! surface data
+               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
+               llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get u2d from u3d and read NetCDF file
+               bf_alias => bf(jp_bdyu2d,jbdy:jbdy)                         ! alias for u2d structure of bdy number jbdy
+               bn_alias => bn_u2d                                          ! alias for u2d structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdyv2d ) THEN
+               cl3 = 'v2d'
+               igrd = 3                                                    ! V point
+               ipk = 1                                                     ! surface data
+               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
+               llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get v2d from v3d and read NetCDF file
+               bf_alias => bf(jp_bdyv2d,jbdy:jbdy)                         ! alias for v2d structure of bdy number jbdy
+               bn_alias => bn_v2d                                          ! alias for v2d structure of nambdy_dta 
+            ENDIF
+            ! =====================
+            !         dyn3d
+            ! =====================
+            IF( jfld == jp_bdyu3d ) THEN
+               cl3 = 'u3d'
+               igrd = 2                                                    ! U point
+               ipk = jpk                                                   ! 3d data
+               llneed = dta_bdy(jbdy)%lneed_dyn3d .OR.               &     ! dta_bdy(jbdy)%u3d will be needed
+                  &   ( dta_bdy(jbdy)%lneed_dyn2d .AND. ln_full_vel )      !   u3d needed to compute u2d
+               llread = nn_dyn3d_dta(jbdy) == 1                            ! get data from NetCDF file
+               bf_alias => bf(jp_bdyu3d,jbdy:jbdy)                         ! alias for u3d structure of bdy number jbdy
+               bn_alias => bn_u3d                                          ! alias for u3d structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdyv3d ) THEN
+               cl3 = 'v3d'
+               igrd = 3                                                    ! V point
+               ipk = jpk                                                   ! 3d data
+               llneed = dta_bdy(jbdy)%lneed_dyn3d .OR.               &     ! dta_bdy(jbdy)%v3d will be needed
+                  &   ( dta_bdy(jbdy)%lneed_dyn2d .AND. ln_full_vel )      !   v3d needed to compute v2d
+               llread = nn_dyn3d_dta(jbdy) == 1                            ! get data from NetCDF file
+               bf_alias => bf(jp_bdyv3d,jbdy:jbdy)                         ! alias for v3d structure of bdy number jbdy
+               bn_alias => bn_v3d                                          ! alias for v3d structure of nambdy_dta 
+            ENDIF
+
+            ! =====================
+            !          tra
+            ! =====================
+            IF( jfld == jp_bdytem ) THEN
+               cl3 = 'tem'
+               igrd = 1                                                    ! T point
+               ipk = jpk                                                   ! 3d data
+               llneed = dta_bdy(jbdy)%lneed_tra                            ! dta_bdy(jbdy)%tem will be needed
+               llread = nn_tra_dta(jbdy) == 1                              ! get data from NetCDF file
+               bf_alias => bf(jp_bdytem,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
+               bn_alias => bn_tem                                          ! alias for ssh structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdysal ) THEN
+               cl3 = 'sal'
+               igrd = 1                                                    ! T point
+               ipk = jpk                                                   ! 3d data
+               llneed = dta_bdy(jbdy)%lneed_tra                            ! dta_bdy(jbdy)%sal will be needed
+               llread = nn_tra_dta(jbdy) == 1                              ! get data from NetCDF file
+               bf_alias => bf(jp_bdysal,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
+               bn_alias => bn_sal                                          ! alias for ssh structure of nambdy_dta 
+            ENDIF
+
+            ! =====================
+            !          ice
+            ! =====================
+            IF( jfld == jp_bdya_i ) THEN
+               cl3 = 'a_i'
+               igrd = 1                                                    ! T point
+               ipk = ipl                                                   ! 
+               llneed = dta_bdy(jbdy)%lneed_ice                            ! dta_bdy(jbdy)%a_i will be needed
+               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file
+               bf_alias => bf(jp_bdya_i,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
+               bn_alias => bn_a_i                                          ! alias for ssh structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdyh_i ) THEN
+               cl3 = 'h_i'
+               igrd = 1                                                    ! T point
+               ipk = ipl                                                   ! 
+               llneed = dta_bdy(jbdy)%lneed_ice                            ! dta_bdy(jbdy)%h_i will be needed
+               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file
+               bf_alias => bf(jp_bdyh_i,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
+               bn_alias => bn_h_i                                          ! alias for ssh structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdyh_s ) THEN
+               cl3 = 'h_s'
+               igrd = 1                                                    ! T point
+               ipk = ipl                                                   ! 
+               llneed = dta_bdy(jbdy)%lneed_ice                            ! dta_bdy(jbdy)%h_s will be needed
+               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file
+               bf_alias => bf(jp_bdyh_s,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
+               bn_alias => bn_h_s                                          ! alias for ssh structure of nambdy_dta 
+            ENDIF
+
+
+            IF( llneed ) THEN                                              ! dta_bdy(jbdy)%xxx will be needed
+               !                                                           !   -> must be associated with an allocated target
+               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
+               ALLOCATE( bf_alias(1)%fnow( iszdim, 1, ipk ) )              ! allocate the target
                !
-               CALL iom_open  ( clname, inum )
-               id1 = iom_varid( inum, bn_a_i%clvar, kdimsz=kdimsz, kndims=kndims, ldstop = .FALSE. )
-               CALL iom_close ( inum )
-
-      	       IF ( kndims == 4 ) THEN
-                 nice_cat = kdimsz(4)   ! Xcat input
-               ELSE
-                 nice_cat = 1           ! 1cat input      
-               ENDIF
-               ! End test
-
-               IF( dta%ll_a_i ) THEN
-                  jfld = jfld + 1
-                  blf_i(jfld) = bn_a_i
-                  ibdy(jfld)  = jbdy
-                  igrid(jfld) = 1
-                  ilen1(jfld) = nblen(igrid(jfld))
-                  ilen3(jfld) = nice_cat
-               ENDIF
-
-               IF( dta%ll_h_i ) THEN
-                  jfld = jfld + 1
-                  blf_i(jfld) = bn_h_i
-                  ibdy(jfld)  = jbdy
-                  igrid(jfld) = 1
-                  ilen1(jfld) = nblen(igrid(jfld))
-                  ilen3(jfld) = nice_cat
-               ENDIF
-
-               IF( dta%ll_h_s ) THEN
-                  jfld = jfld + 1
-                  blf_i(jfld) = bn_h_s
-                  ibdy(jfld)  = jbdy
-                  igrid(jfld) = 1
-                  ilen1(jfld) = nblen(igrid(jfld))
-                  ilen3(jfld) = nice_cat
-               ENDIF
-
-            ENDIF
-#endif
-            ! Recalculate field counts
-            !-------------------------
-            IF( jbdy == 1 ) THEN 
-               nb_bdy_fld_sum = 0
-               nb_bdy_fld(jbdy) = jfld
-               nb_bdy_fld_sum     = jfld              
-            ELSE
-               nb_bdy_fld(jbdy) = jfld - nb_bdy_fld_sum
-               nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(jbdy)
-            ENDIF
-
-            dta%nread(1) = nb_bdy_fld(jbdy)
-
-         ENDIF ! nn_dta == 1
-      ENDDO ! jbdy
-
-      DO jfld = 1, nb_bdy_fld_sum
-         ALLOCATE( bf(jfld)%fnow(ilen1(jfld),1,ilen3(jfld)) )
-         IF( blf_i(jfld)%ln_tint ) ALLOCATE( bf(jfld)%fdta(ilen1(jfld),1,ilen3(jfld),2) )
-         nbmap_ptr(jfld)%ptr => idx_bdy(ibdy(jfld))%nbmap(:,igrid(jfld))
-         nbmap_ptr(jfld)%ll_unstruc = ln_coords_file(ibdy(jfld))
-      ENDDO
-
-      ! fill bf with blf_i and control print
-      !-------------------------------------
-      jstart = 1
-      DO jbdy = 1, nb_bdy
-         jend = jstart - 1 + nb_bdy_fld(jbdy) 
-         CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(jbdy), 'bdy_dta',   &
-         &              'open boundary conditions', 'nambdy_dta' )
-         jstart = jend + 1
-      ENDDO
-
-      DO jfld = 1, nb_bdy_fld_sum
-         bf(jfld)%igrd = igrid(jfld)
-         bf(jfld)%ibdy = ibdy(jfld)
-      ENDDO
-
-      ! Initialise local boundary data arrays
-      ! nn_xxx_dta=0 : allocate space - will be filled from initial conditions later
-      ! nn_xxx_dta=1 : point to "fnow" arrays
-      !-------------------------------------
-
-      jfld = 0
-      DO jbdy=1, nb_bdy
-
-         nblen => idx_bdy(jbdy)%nblen
-         dta => dta_bdy(jbdy)
-
-         if(lwp) then
-            write(numout,*) '++++++ dta%ll_ssh = ',dta%ll_ssh
-            write(numout,*) '++++++ dta%ll_u2d = ',dta%ll_u2d
-            write(numout,*) '++++++ dta%ll_v2d = ',dta%ll_v2d
-            write(numout,*) '++++++ dta%ll_u3d = ',dta%ll_u3d
-            write(numout,*) '++++++ dta%ll_v3d = ',dta%ll_v3d
-            write(numout,*) '++++++ dta%ll_tem = ',dta%ll_tem
-            write(numout,*) '++++++ dta%ll_sal = ',dta%ll_sal
-         endif
-
-         IF ( nn_dyn2d_dta(jbdy) == 0 .or. nn_dyn2d_dta(jbdy) == 2 ) THEN
-            if(lwp) write(numout,*) '++++++ dta%ssh/u2d/u3d allocated space'
-            IF( dta%ll_ssh ) ALLOCATE( dta%ssh(nblen(1)) )
-            IF( dta%ll_u2d ) ALLOCATE( dta%u2d(nblen(2)) )
-            IF( dta%ll_v2d ) ALLOCATE( dta%v2d(nblen(3)) )
-         ENDIF
-         IF ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN
-            IF( dta%ll_ssh ) THEN
-               if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow'
-               jfld = jfld + 1
-               dta%ssh => bf(jfld)%fnow(:,1,1)
-            ENDIF
-            IF ( dta%ll_u2d ) THEN
-               IF ( ln_full_vel_array(jbdy) ) THEN
-                  if(lwp) write(numout,*) '++++++ dta%u2d allocated space'
-                  ALLOCATE( dta%u2d(nblen(2)) )
-               ELSE
-                  if(lwp) write(numout,*) '++++++ dta%u2d pointing to fnow'
-                  jfld = jfld + 1
-                  dta%u2d => bf(jfld)%fnow(:,1,1)
-               ENDIF
-            ENDIF
-            IF ( dta%ll_v2d ) THEN
-               IF ( ln_full_vel_array(jbdy) ) THEN
-                  if(lwp) write(numout,*) '++++++ dta%v2d allocated space'
-                  ALLOCATE( dta%v2d(nblen(3)) )
-               ELSE
-                  if(lwp) write(numout,*) '++++++ dta%v2d pointing to fnow'
-                  jfld = jfld + 1
-                  dta%v2d => bf(jfld)%fnow(:,1,1)
-               ENDIF
-            ENDIF
-         ENDIF
-
-         IF ( nn_dyn3d_dta(jbdy) == 0 ) THEN
-            if(lwp) write(numout,*) '++++++ dta%u3d/v3d allocated space'
-            IF( dta%ll_u3d ) ALLOCATE( dta_bdy(jbdy)%u3d(nblen(2),jpk) )
-            IF( dta%ll_v3d ) ALLOCATE( dta_bdy(jbdy)%v3d(nblen(3),jpk) )
-         ENDIF
-         IF ( nn_dyn3d_dta(jbdy) == 1 .or. &
-           &  ( ln_full_vel_array(jbdy) .and. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN
-            IF ( dta%ll_u3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN
-               if(lwp) write(numout,*) '++++++ dta%u3d pointing to fnow'
-               jfld = jfld + 1
-               dta_bdy(jbdy)%u3d => bf(jfld)%fnow(:,1,:)
-            ENDIF
-            IF ( dta%ll_v3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN
-               if(lwp) write(numout,*) '++++++ dta%v3d pointing to fnow'
-               jfld = jfld + 1
-               dta_bdy(jbdy)%v3d => bf(jfld)%fnow(:,1,:)
-            ENDIF
-         ENDIF
-
-         IF( nn_tra_dta(jbdy) == 0 ) THEN
-            if(lwp) write(numout,*) '++++++ dta%tem/sal allocated space'
-            IF( dta%ll_tem ) ALLOCATE( dta_bdy(jbdy)%tem(nblen(1),jpk) )
-            IF( dta%ll_sal ) ALLOCATE( dta_bdy(jbdy)%sal(nblen(1),jpk) )
-         ELSE
-            IF( dta%ll_tem ) THEN
-               if(lwp) write(numout,*) '++++++ dta%tem pointing to fnow'
-               jfld = jfld + 1
-               dta_bdy(jbdy)%tem => bf(jfld)%fnow(:,1,:)
-            ENDIF
-            IF( dta%ll_sal ) THEN 
-               if(lwp) write(numout,*) '++++++ dta%sal pointing to fnow'
-               jfld = jfld + 1
-               dta_bdy(jbdy)%sal => bf(jfld)%fnow(:,1,:)
-            ENDIF
-         ENDIF
-
-#if defined key_si3
-         IF (cn_ice(jbdy) /= 'none') THEN
-            IF( nn_ice_dta(jbdy) == 0 ) THEN
-               ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) )
-               ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) )
-               ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) )
-            ELSE
-               IF ( nice_cat == jpl ) THEN ! case input cat = jpl
-                  jfld = jfld + 1
-                  dta_bdy(jbdy)%a_i => bf(jfld)%fnow(:,1,:)
-                  jfld = jfld + 1
-                  dta_bdy(jbdy)%h_i => bf(jfld)%fnow(:,1,:)
-                  jfld = jfld + 1
-                  dta_bdy(jbdy)%h_s => bf(jfld)%fnow(:,1,:)
-               ELSE                        ! case input cat = 1 OR (/=1 and /=jpl)
-                  jfld_ait(jbdy)  = jfld + 1
-                  jfld_htit(jbdy) = jfld + 2
-                  jfld_htst(jbdy) = jfld + 3
-                  jfld     = jfld + 3
-                  ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) )
-                  ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) )
-                  ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) )
-                  dta_bdy(jbdy)%a_i(:,:) = 0._wp
-                  dta_bdy(jbdy)%h_i(:,:) = 0._wp
-                  dta_bdy(jbdy)%h_s(:,:) = 0._wp
-               ENDIF
-
-            ENDIF
-         ENDIF
-#endif
+               IF( llread ) THEN                                           ! get data from NetCDF file
+                  CALL fld_fill( bf_alias, bn_alias, cn_dir, 'bdy_dta', cl3//' '//ctmp1, ctmp2 )   ! use namelist info
+                  IF( bf_alias(1)%ln_tint ) ALLOCATE( bf_alias(1)%fdta( iszdim, 1, ipk, 2 ) )
+                  bf_alias(1)%imap => idx_bdy(jbdy)%nbmap(1:iszdim,igrd)      ! associate the mapping used for this bdy
+                  bf_alias(1)%igrd = igrd                                     ! used only for vertical integration of 3D arrays
+                  bf_alias(1)%ltotvel = ln_full_vel                           ! T if u3d is full velocity
+               ENDIF
+
+               ! associate the pointer and get rid of the dimensions with a size equal to 1
+               IF( jfld == jp_bdyssh           ) dta_bdy(jbdy)%ssh => bf_alias(1)%fnow(:,1,1)
+               IF( jfld == jp_bdyu2d           ) dta_bdy(jbdy)%u2d => bf_alias(1)%fnow(:,1,1)
+               IF( jfld == jp_bdyv2d           ) dta_bdy(jbdy)%v2d => bf_alias(1)%fnow(:,1,1)
+               IF( jfld == jp_bdyu3d           ) dta_bdy(jbdy)%u3d => bf_alias(1)%fnow(:,1,:)
+               IF( jfld == jp_bdyv3d           ) dta_bdy(jbdy)%v3d => bf_alias(1)%fnow(:,1,:)
+               IF( jfld == jp_bdytem           ) dta_bdy(jbdy)%tem => bf_alias(1)%fnow(:,1,:)
+               IF( jfld == jp_bdysal           ) dta_bdy(jbdy)%sal => bf_alias(1)%fnow(:,1,:)
+               IF( jfld == jp_bdya_i ) THEN
+                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%a_i => bf_alias(1)%fnow(:,1,:)
+                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%a_i(iszdim,jpl) )
+                  ENDIF
+               ENDIF
+               IF( jfld == jp_bdyh_i ) THEN
+                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%h_i => bf_alias(1)%fnow(:,1,:)
+                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%h_i(iszdim,jpl) )
+                  ENDIF
+               ENDIF
+               IF( jfld == jp_bdyh_s ) THEN
+                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%h_s => bf_alias(1)%fnow(:,1,:)
+                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%h_s(iszdim,jpl) )
+                  ENDIF
+               ENDIF
+            ENDIF
+
+         END DO   ! jpbdyfld
          !
       END DO ! jbdy 
       !
    END SUBROUTINE bdy_dta_init
-
+   
    !!==============================================================================
 END MODULE bdydta
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90	(revision 11223)
@@ -37,6 +37,4 @@
 
    INTEGER, PARAMETER ::   jp_nseg = 100   ! 
-   INTEGER, PARAMETER ::   nrimmax =  20   ! maximum rimwidth in structured
-                                               ! open boundary data files
    ! Straight open boundary segment parameters:
    INTEGER  ::   nbdysege, nbdysegw, nbdysegn, nbdysegs 
@@ -70,5 +68,5 @@
          &             cn_ice, nn_ice_dta,                                     &
          &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
-         &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
+         &             ln_vol, nn_volctl, nn_rimwidth
          !
       INTEGER  ::   ios                 ! Local integer output status for namelist read
@@ -81,4 +79,22 @@
       READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp )
+      ! make sur that all elements of the namelist variables have a default definition from namelist_ref
+      ln_coords_file (2:jp_bdy) = ln_coords_file (1)
+      cn_coords_file (2:jp_bdy) = cn_coords_file (1)
+      cn_dyn2d       (2:jp_bdy) = cn_dyn2d       (1)
+      nn_dyn2d_dta   (2:jp_bdy) = nn_dyn2d_dta   (1)
+      cn_dyn3d       (2:jp_bdy) = cn_dyn3d       (1)
+      nn_dyn3d_dta   (2:jp_bdy) = nn_dyn3d_dta   (1)
+      cn_tra         (2:jp_bdy) = cn_tra         (1)
+      nn_tra_dta     (2:jp_bdy) = nn_tra_dta     (1)    
+      ln_tra_dmp     (2:jp_bdy) = ln_tra_dmp     (1)
+      ln_dyn3d_dmp   (2:jp_bdy) = ln_dyn3d_dmp   (1)
+      rn_time_dmp    (2:jp_bdy) = rn_time_dmp    (1)
+      rn_time_dmp_out(2:jp_bdy) = rn_time_dmp_out(1)
+      cn_ice         (2:jp_bdy) = cn_ice         (1)
+      nn_ice_dta     (2:jp_bdy) = nn_ice_dta     (1)
+      rn_ice_tem     (2:jp_bdy) = rn_ice_tem     (1)
+      rn_ice_sal     (2:jp_bdy) = rn_ice_sal     (1)
+      rn_ice_age     (2:jp_bdy) = rn_ice_age     (1)
       REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries
       READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
@@ -87,4 +103,6 @@
 
       IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE.   ! forced for Agrif children
+
+      IF( nb_bdy == 0 ) ln_bdy = .FALSE.
       
       ! -----------------------------------------
@@ -97,5 +115,5 @@
          !
          ! Open boundaries definition (arrays and masks)
-         CALL bdy_segs
+         CALL bdy_def
          !
          ! Open boundaries initialisation of external data arrays
@@ -115,5 +133,5 @@
 
 
-   SUBROUTINE bdy_segs
+   SUBROUTINE bdy_def
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE bdy_init  ***
@@ -130,25 +148,20 @@
       INTEGER  ::   ilen1, ibm1                            !   -       -
       INTEGER  ::   iwe, ies, iso, ino, inum, id_dummy     !   -       -
-      INTEGER  ::   jpbdta, jpbdtau, jpbdtas               !   -       -
+      INTEGER  ::   jpbdta                                 !   -       -
       INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       -
       INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3           !   -       -
       INTEGER  ::   iibe, ijbe, iibi, ijbi                 !   -       -
       INTEGER  ::   flagu, flagv                           ! short cuts
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zz_read      ! work space for 2D global boundary data
       REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask   ! pointer to 2D mask fields
-      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
-      INTEGER, DIMENSION (2)                  ::   kdimsz
+      INTEGER, DIMENSION (4)                  ::   kdimsz
       INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays 
       INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta
       INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
       CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid
-      INTEGER :: jpk_max
       REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                  ! temporary fmask array excluding coastal boundary condition (shlat)
       REAL(wp), DIMENSION(jpi,jpj) ::   ztmask, zumask, zvmask  ! temporary u/v mask array
       !!
-      CHARACTER(LEN=1)                     ::   ctypebdy   !     -        - 
       INTEGER                              ::   nbdyind, nbdybeg, nbdyend
-      !!
-      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
-      INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!----------------------------------------------------------------------
       !
@@ -161,247 +174,197 @@
          &                               ' and general open boundary condition are not compatible' )
 
-      IF( nb_bdy == 0 ) THEN 
-        IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
-      ELSE
-        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy
+      IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy
+
+      DO ib_bdy = 1,nb_bdy
+
+         IF(lwp) THEN
+            WRITE(numout,*) ' ' 
+            WRITE(numout,*) '------ Open boundary data set ',ib_bdy,' ------' 
+            IF( ln_coords_file(ib_bdy) ) THEN
+               WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
+            ELSE
+               WRITE(numout,*) 'Boundary defined in namelist.'
+            ENDIF
+            WRITE(numout,*)
+         ENDIF
+
+         ! barotropic bdy
+         !----------------
+         IF(lwp) THEN
+            WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
+            SELECT CASE( cn_dyn2d(ib_bdy) )                  
+            CASE( 'none' )           ;   WRITE(numout,*) '      no open boundary condition'        
+            CASE( 'frs' )            ;   WRITE(numout,*) '      Flow Relaxation Scheme'
+            CASE( 'flather' )        ;   WRITE(numout,*) '      Flather radiation condition'
+            CASE( 'orlanski' )       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
+            CASE( 'orlanski_npo' )   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
+            CASE DEFAULT             ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' )
+            END SELECT
+         ENDIF
+
+         dta_bdy(ib_bdy)%lneed_ssh   = cn_dyn2d(ib_bdy) == 'flather'
+         dta_bdy(ib_bdy)%lneed_dyn2d = cn_dyn2d(ib_bdy) /= 'none'
+
+         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
+            SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   ! 
+            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'        
+            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file'
+            CASE( 2 )      ;   WRITE(numout,*) '      tidal harmonic forcing taken from file'
+            CASE( 3 )      ;   WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files'
+            CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
+            END SELECT
+         ENDIF
+         IF ( dta_bdy(ib_bdy)%lneed_dyn2d .AND. nn_dyn2d_dta(ib_bdy) .GE. 2  .AND. .NOT.ln_tide ) THEN
+            CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' )
+         ENDIF
+         IF(lwp) WRITE(numout,*)
+
+         ! baroclinic bdy
+         !----------------
+         IF(lwp) THEN
+            WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
+            SELECT CASE( cn_dyn3d(ib_bdy) )                  
+            CASE('none')           ;   WRITE(numout,*) '      no open boundary condition'        
+            CASE('frs')            ;   WRITE(numout,*) '      Flow Relaxation Scheme'
+            CASE('specified')      ;   WRITE(numout,*) '      Specified value'
+            CASE('neumann')        ;   WRITE(numout,*) '      Neumann conditions'
+            CASE('zerograd')       ;   WRITE(numout,*) '      Zero gradient for baroclinic velocities'
+            CASE('zero')           ;   WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
+            CASE('orlanski')       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
+            CASE('orlanski_npo')   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
+            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn3d' )
+            END SELECT
+         ENDIF
+
+         dta_bdy(ib_bdy)%lneed_dyn3d = cn_dyn3d(ib_bdy) == 'frs'      .OR. cn_dyn3d(ib_bdy) == 'specified'   &
+            &                     .OR. cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo'
+
+         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn3d ) THEN
+            SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   ! 
+            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'        
+            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file'
+            CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
+            END SELECT
+         END IF
+
+         IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
+            IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN
+               IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'
+               ln_dyn3d_dmp(ib_bdy) = .false.
+            ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN
+               CALL ctl_stop( 'Use FRS OR relaxation' )
+            ELSE
+               IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone'
+               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
+               IF(rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
+               dta_bdy(ib_bdy)%lneed_dyn3d = .TRUE.
+            ENDIF
+         ELSE
+            IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities'
+         ENDIF
+         IF(lwp) WRITE(numout,*)
+
+         !    tra bdy
+         !----------------
+         IF(lwp) THEN
+            WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
+            SELECT CASE( cn_tra(ib_bdy) )                  
+            CASE('none')           ;   WRITE(numout,*) '      no open boundary condition'        
+            CASE('frs')            ;   WRITE(numout,*) '      Flow Relaxation Scheme'
+            CASE('specified')      ;   WRITE(numout,*) '      Specified value'
+            CASE('neumann')        ;   WRITE(numout,*) '      Neumann conditions'
+            CASE('runoff')         ;   WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity'
+            CASE('orlanski')       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
+            CASE('orlanski_npo')   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
+            CASE DEFAULT           ;   CALL ctl_stop( 'unrecognised value for cn_tra' )
+            END SELECT
+         ENDIF
+
+         dta_bdy(ib_bdy)%lneed_tra = cn_tra(ib_bdy) == 'frs'       .OR. cn_tra(ib_bdy) == 'specified'   &
+            &                   .OR. cn_tra(ib_bdy) == 'orlanski'  .OR. cn_tra(ib_bdy) == 'orlanski_npo' 
+
+         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_tra ) THEN
+            SELECT CASE( nn_tra_dta(ib_bdy) )                   ! 
+            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'        
+            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file'
+            CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
+            END SELECT
+         ENDIF
+
+         IF ( ln_tra_dmp(ib_bdy) ) THEN
+            IF ( cn_tra(ib_bdy) == 'none' ) THEN
+               IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
+               ln_tra_dmp(ib_bdy) = .false.
+            ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN
+               CALL ctl_stop( 'Use FRS OR relaxation' )
+            ELSE
+               IF(lwp) WRITE(numout,*) '      + T/S relaxation zone'
+               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
+               IF(lwp) WRITE(numout,*) '      Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days'
+               IF(lwp.AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
+               dta_bdy(ib_bdy)%lneed_tra = .TRUE.
+            ENDIF
+         ELSE
+            IF(lwp) WRITE(numout,*) '      NO T/S relaxation'
+         ENDIF
+         IF(lwp) WRITE(numout,*)
+
+#if defined key_si3
+         IF(lwp) THEN
+            WRITE(numout,*) 'Boundary conditions for sea ice:  '
+            SELECT CASE( cn_ice(ib_bdy) )                  
+            CASE('none')   ;   WRITE(numout,*) '      no open boundary condition'        
+            CASE('frs')    ;   WRITE(numout,*) '      Flow Relaxation Scheme'
+            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice' )
+            END SELECT
+         ENDIF
+
+         dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none'
+
+         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN 
+            SELECT CASE( nn_ice_dta(ib_bdy) )                   ! 
+            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'        
+            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file'
+            CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_dta must be 0 or 1' )
+            END SELECT
+            WRITE(numout,*)
+            WRITE(numout,*) '      tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)         
+            WRITE(numout,*) '      sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)         
+            WRITE(numout,*) '      age of bdy sea-ice = ', rn_ice_age(ib_bdy)         
+         ENDIF
+#else
+         dta_bdy(ib_bdy)%lneed_ice = .FALSE.
+#endif
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
+         IF(lwp) WRITE(numout,*)
+         !
+      END DO   ! nb_bdy
+
+      IF( lwp ) THEN
+         IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
+            WRITE(numout,*) 'Volume correction applied at open boundaries'
+            WRITE(numout,*)
+            SELECT CASE ( nn_volctl )
+            CASE( 1 )      ;   WRITE(numout,*) '      The total volume will be constant'
+            CASE( 0 )      ;   WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
+            CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
+            END SELECT
+            WRITE(numout,*)
+            !
+            ! sanity check if used with tides        
+            IF( ln_tide ) THEN 
+               WRITE(numout,*) ' The total volume correction is not working with tides. '
+               WRITE(numout,*) ' Set ln_vol to .FALSE. '
+               WRITE(numout,*) ' or '
+               WRITE(numout,*) ' equilibriate your bdy input files '
+               CALL ctl_stop( 'The total volume correction is not working with tides.' )
+            END IF
+         ELSE
+            WRITE(numout,*) 'No volume correction applied at open boundaries'
+            WRITE(numout,*)
+         ENDIF
       ENDIF
-
-      DO ib_bdy = 1,nb_bdy
-        IF(lwp) WRITE(numout,*) ' ' 
-        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------' 
-
-        IF( ln_coords_file(ib_bdy) ) THEN
-           IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
-        ELSE
-           IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'
-        ENDIF
-        IF(lwp) WRITE(numout,*)
-
-        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
-        SELECT CASE( cn_dyn2d(ib_bdy) )                  
-          CASE( 'none' )         
-             IF(lwp) WRITE(numout,*) '      no open boundary condition'        
-             dta_bdy(ib_bdy)%ll_ssh = .false.
-             dta_bdy(ib_bdy)%ll_u2d = .false.
-             dta_bdy(ib_bdy)%ll_v2d = .false.
-          CASE( 'frs' )          
-             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
-             dta_bdy(ib_bdy)%ll_ssh = .false.
-             dta_bdy(ib_bdy)%ll_u2d = .true.
-             dta_bdy(ib_bdy)%ll_v2d = .true.
-          CASE( 'flather' )      
-             IF(lwp) WRITE(numout,*) '      Flather radiation condition'
-             dta_bdy(ib_bdy)%ll_ssh = .true.
-             dta_bdy(ib_bdy)%ll_u2d = .true.
-             dta_bdy(ib_bdy)%ll_v2d = .true.
-          CASE( 'orlanski' )     
-             IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
-             dta_bdy(ib_bdy)%ll_ssh = .false.
-             dta_bdy(ib_bdy)%ll_u2d = .true.
-             dta_bdy(ib_bdy)%ll_v2d = .true.
-          CASE( 'orlanski_npo' ) 
-             IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
-             dta_bdy(ib_bdy)%ll_ssh = .false.
-             dta_bdy(ib_bdy)%ll_u2d = .true.
-             dta_bdy(ib_bdy)%ll_v2d = .true.
-          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' )
-        END SELECT
-        IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN
-           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   ! 
-              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
-              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
-              CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file'
-              CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files'
-              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
-           END SELECT
-           IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.ln_tide)) THEN
-             CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' )
-           ENDIF
-        ENDIF
-        IF(lwp) WRITE(numout,*)
-
-        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
-        SELECT CASE( cn_dyn3d(ib_bdy) )                  
-          CASE('none')
-             IF(lwp) WRITE(numout,*) '      no open boundary condition'        
-             dta_bdy(ib_bdy)%ll_u3d = .false.
-             dta_bdy(ib_bdy)%ll_v3d = .false.
-          CASE('frs')       
-             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
-             dta_bdy(ib_bdy)%ll_u3d = .true.
-             dta_bdy(ib_bdy)%ll_v3d = .true.
-          CASE('specified')
-             IF(lwp) WRITE(numout,*) '      Specified value'
-             dta_bdy(ib_bdy)%ll_u3d = .true.
-             dta_bdy(ib_bdy)%ll_v3d = .true.
-          CASE('neumann')
-             IF(lwp) WRITE(numout,*) '      Neumann conditions'
-             dta_bdy(ib_bdy)%ll_u3d = .false.
-             dta_bdy(ib_bdy)%ll_v3d = .false.
-          CASE('zerograd')
-             IF(lwp) WRITE(numout,*) '      Zero gradient for baroclinic velocities'
-             dta_bdy(ib_bdy)%ll_u3d = .false.
-             dta_bdy(ib_bdy)%ll_v3d = .false.
-          CASE('zero')
-             IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
-             dta_bdy(ib_bdy)%ll_u3d = .false.
-             dta_bdy(ib_bdy)%ll_v3d = .false.
-          CASE('orlanski')
-             IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
-             dta_bdy(ib_bdy)%ll_u3d = .true.
-             dta_bdy(ib_bdy)%ll_v3d = .true.
-          CASE('orlanski_npo')
-             IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
-             dta_bdy(ib_bdy)%ll_u3d = .true.
-             dta_bdy(ib_bdy)%ll_v3d = .true.
-          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn3d' )
-        END SELECT
-        IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
-           SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   ! 
-              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
-              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
-              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
-           END SELECT
-        ENDIF
-
-        IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
-           IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN
-              IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'
-              ln_dyn3d_dmp(ib_bdy)=.false.
-           ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN
-              CALL ctl_stop( 'Use FRS OR relaxation' )
-           ELSE
-              IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone'
-              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
-              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
-              dta_bdy(ib_bdy)%ll_u3d = .true.
-              dta_bdy(ib_bdy)%ll_v3d = .true.
-           ENDIF
-        ELSE
-           IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities'
-        ENDIF
-        IF(lwp) WRITE(numout,*)
-
-        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
-        SELECT CASE( cn_tra(ib_bdy) )                  
-          CASE('none')
-             IF(lwp) WRITE(numout,*) '      no open boundary condition'        
-             dta_bdy(ib_bdy)%ll_tem = .false.
-             dta_bdy(ib_bdy)%ll_sal = .false.
-          CASE('frs')
-             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
-             dta_bdy(ib_bdy)%ll_tem = .true.
-             dta_bdy(ib_bdy)%ll_sal = .true.
-          CASE('specified')
-             IF(lwp) WRITE(numout,*) '      Specified value'
-             dta_bdy(ib_bdy)%ll_tem = .true.
-             dta_bdy(ib_bdy)%ll_sal = .true.
-          CASE('neumann')
-             IF(lwp) WRITE(numout,*) '      Neumann conditions'
-             dta_bdy(ib_bdy)%ll_tem = .false.
-             dta_bdy(ib_bdy)%ll_sal = .false.
-          CASE('runoff')
-             IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity'
-             dta_bdy(ib_bdy)%ll_tem = .false.
-             dta_bdy(ib_bdy)%ll_sal = .false.
-          CASE('orlanski')
-             IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
-             dta_bdy(ib_bdy)%ll_tem = .true.
-             dta_bdy(ib_bdy)%ll_sal = .true.
-          CASE('orlanski_npo')
-             IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
-             dta_bdy(ib_bdy)%ll_tem = .true.
-             dta_bdy(ib_bdy)%ll_sal = .true.
-          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_tra' )
-        END SELECT
-        IF( cn_tra(ib_bdy) /= 'none' ) THEN
-           SELECT CASE( nn_tra_dta(ib_bdy) )                   ! 
-              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
-              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
-              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
-           END SELECT
-        ENDIF
-
-        IF ( ln_tra_dmp(ib_bdy) ) THEN
-           IF ( cn_tra(ib_bdy) == 'none' ) THEN
-              IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
-              ln_tra_dmp(ib_bdy)=.false.
-           ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN
-              CALL ctl_stop( 'Use FRS OR relaxation' )
-           ELSE
-              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone'
-              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
-              IF(lwp) WRITE(numout,*) '      Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days'
-              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
-              dta_bdy(ib_bdy)%ll_tem = .true.
-              dta_bdy(ib_bdy)%ll_sal = .true.
-           ENDIF
-        ELSE
-           IF(lwp) WRITE(numout,*) '      NO T/S relaxation'
-        ENDIF
-        IF(lwp) WRITE(numout,*)
-
-#if defined key_si3
-         IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
-         SELECT CASE( cn_ice(ib_bdy) )                  
-         CASE('none')
-             IF(lwp) WRITE(numout,*) '      no open boundary condition'        
-             dta_bdy(ib_bdy)%ll_a_i = .false.
-             dta_bdy(ib_bdy)%ll_h_i = .false.
-             dta_bdy(ib_bdy)%ll_h_s = .false.
-         CASE('frs')
-             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
-             dta_bdy(ib_bdy)%ll_a_i = .true.
-             dta_bdy(ib_bdy)%ll_h_i = .true.
-             dta_bdy(ib_bdy)%ll_h_s = .true.
-         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice' )
-         END SELECT
-        IF( cn_ice(ib_bdy) /= 'none' ) THEN 
-           SELECT CASE( nn_ice_dta(ib_bdy) )                   ! 
-              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
-              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
-              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_dta must be 0 or 1' )
-           END SELECT
-        ENDIF
-        IF(lwp) WRITE(numout,*)
-        IF(lwp) WRITE(numout,*) '      tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)         
-        IF(lwp) WRITE(numout,*) '      sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)         
-        IF(lwp) WRITE(numout,*) '      age of bdy sea-ice = ', rn_ice_age(ib_bdy)         
-#endif
-
-        IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
-        IF(lwp) WRITE(numout,*)
-         !
-      END DO
-
-     IF( nb_bdy > 0 ) THEN
-        IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
-          IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
-          IF(lwp) WRITE(numout,*)
-          SELECT CASE ( nn_volctl )
-            CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant'
-            CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
-            CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
-          END SELECT
-          IF(lwp) WRITE(numout,*)
-          !
-          ! sanity check if used with tides        
-          IF( ln_tide ) THEN 
-             IF(lwp) WRITE(numout,*) ' The total volume correction is not working with tides. '
-             IF(lwp) WRITE(numout,*) ' Set ln_vol to .FALSE. '
-             IF(lwp) WRITE(numout,*) ' or '
-             IF(lwp) WRITE(numout,*) ' equilibriate your bdy input files '
-             CALL ctl_stop( 'The total volume correction is not working with tides.' )
-          END IF
-        ELSE
-          IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
-          IF(lwp) WRITE(numout,*)
-        ENDIF
-        IF( nb_jpk_bdy(ib_bdy) > 0 ) THEN
-           IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***'
-        ELSE
-           IF(lwp) WRITE(numout,*) '*** open boundary will be read straight onto the native grid without vertical interpolation ***'
-        ENDIF
-     ENDIF
 
       ! -------------------------------------------------
@@ -409,8 +372,5 @@
       ! -------------------------------------------------
 
-      ! Work out global dimensions of boundary data
-      ! ---------------------------------------------
       REWIND( numnam_cfg )     
-
       nblendta(:,:) = 0
       nbdysege = 0
@@ -418,123 +378,32 @@
       nbdysegn = 0
       nbdysegs = 0
-      icount   = 0 ! count user defined segments
-      ! Dimensions below are used to allocate arrays to read external data
-      jpbdtas = 1 ! Maximum size of boundary data (structured case)
-      jpbdtau = 1 ! Maximum size of boundary data (unstructured case)
-
+
+      ! Define all boundaries 
+      ! ---------------------
       DO ib_bdy = 1, nb_bdy
-
-         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
- 
-            icount = icount + 1
-            ! No REWIND here because may need to read more than one nambdy_index namelist.
-            ! Read only namelist_cfg to avoid unseccessfull overwrite 
-            ! keep full control of the configuration namelist
-            READ  ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 )
-904         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp )
-            IF(lwm) WRITE ( numond, nambdy_index )
-
-            SELECT CASE ( TRIM(ctypebdy) )
-              CASE( 'N' )
-                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
-                    nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain.
-                    nbdybeg  = 2
-                    nbdyend  = jpiglo - 1
-                 ENDIF
-                 nbdysegn = nbdysegn + 1
-                 npckgn(nbdysegn) = ib_bdy ! Save bdy package number
-                 jpjnob(nbdysegn) = nbdyind
-                 jpindt(nbdysegn) = nbdybeg
-                 jpinft(nbdysegn) = nbdyend
-                 !
-              CASE( 'S' )
-                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
-                    nbdyind  = 2           ! set boundary to whole side of model domain.
-                    nbdybeg  = 2
-                    nbdyend  = jpiglo - 1
-                 ENDIF
-                 nbdysegs = nbdysegs + 1
-                 npckgs(nbdysegs) = ib_bdy ! Save bdy package number
-                 jpjsob(nbdysegs) = nbdyind
-                 jpisdt(nbdysegs) = nbdybeg
-                 jpisft(nbdysegs) = nbdyend
-                 !
-              CASE( 'E' )
-                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
-                    nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain.
-                    nbdybeg  = 2
-                    nbdyend  = jpjglo - 1
-                 ENDIF
-                 nbdysege = nbdysege + 1 
-                 npckge(nbdysege) = ib_bdy ! Save bdy package number
-                 jpieob(nbdysege) = nbdyind
-                 jpjedt(nbdysege) = nbdybeg
-                 jpjeft(nbdysege) = nbdyend
-                 !
-              CASE( 'W' )
-                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
-                    nbdyind  = 2           ! set boundary to whole side of model domain.
-                    nbdybeg  = 2
-                    nbdyend  = jpjglo - 1
-                 ENDIF
-                 nbdysegw = nbdysegw + 1
-                 npckgw(nbdysegw) = ib_bdy ! Save bdy package number
-                 jpiwob(nbdysegw) = nbdyind
-                 jpjwdt(nbdysegw) = nbdybeg
-                 jpjwft(nbdysegw) = nbdyend
-                 !
-              CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
-            END SELECT
-
-            ! For simplicity we assume that in case of straight bdy, arrays have the same length
-            ! (even if it is true that last tangential velocity points
-            ! are useless). This simplifies a little bit boundary data format (and agrees with format
-            ! used so far in obc package)
-
-            nblendta(1:jpbgrd,ib_bdy) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy)
-            jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1))
-            IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) &
-            & CALL ctl_stop( 'rimwidth must be lower than nrimmax' )
-
-         ELSE            ! Read size of arrays in boundary coordinates file.
+         !
+         IF( .NOT. ln_coords_file(ib_bdy) ) THEN     ! build bdy coordinates with segments defined in namelist
+
+            CALL bdy_read_seg( ib_bdy, nblendta(:,ib_bdy) )
+            ! Now look for crossings in user (namelist) defined open boundary segments:
+            IF( nbdysege > 0 .OR. nbdysegw > 0 .OR. nbdysegn > 0 .OR. nbdysegs > 0)   CALL bdy_ctl_seg
+
+         ELSE                                        ! Read size of arrays in boundary coordinates file.
+            
             CALL iom_open( cn_coords_file(ib_bdy), inum )
             DO igrd = 1, jpbgrd
                id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )  
                nblendta(igrd,ib_bdy) = MAXVAL(kdimsz)
-               jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz))
             END DO
             CALL iom_close( inum )
-            !
-         ENDIF 
+         ENDIF
          !
       END DO ! ib_bdy
-
-      IF (nb_bdy>0) THEN
-         jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
-
-         ! Allocate arrays
-         !---------------
-         ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    &
-            &      nbrdta(jpbdta, jpbgrd, nb_bdy) )
-         
-         jpk_max = MAXVAL(nb_jpk_bdy)
-         jpk_max = MAX(jpk_max, jpk)
-
-         ALLOCATE( dta_global(jpbdtau, 1, jpk_max) )
-         ALLOCATE( dta_global_z(jpbdtau, 1, jpk_max) ) ! needed ?? TODO
-         ALLOCATE( dta_global_dz(jpbdtau, 1, jpk_max) )! needed ?? TODO
-
-         IF ( icount>0 ) THEN
-            ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk_max) )
-            ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk_max) ) ! needed ?? TODO
-            ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk_max) )! needed ?? TODO  
-         ENDIF
-         ! 
-      ENDIF
-
-      ! Now look for crossings in user (namelist) defined open boundary segments:
-      !--------------------------------------------------------------------------
-      IF( icount>0 )   CALL bdy_ctl_seg
-
+      
+      ! Allocate arrays
+      !---------------
+      jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
+      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), nbrdta(jpbdta, jpbgrd, nb_bdy) )
+    
       ! Calculate global boundary index arrays or read in from file
       !------------------------------------------------------------               
@@ -544,17 +413,19 @@
          IF( ln_coords_file(ib_bdy) ) THEN
             !
+            ALLOCATE( zz_read( MAXVAL(nblendta), 1 ) )          
             CALL iom_open( cn_coords_file(ib_bdy), inum )
+            !
             DO igrd = 1, jpbgrd
-               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
+               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
                DO ii = 1,nblendta(igrd,ib_bdy)
-                  nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
+                  nbidta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
                END DO
-               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
+               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
                DO ii = 1,nblendta(igrd,ib_bdy)
-                  nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
+                  nbjdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
                END DO
-               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
+               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
                DO ii = 1,nblendta(igrd,ib_bdy)
-                  nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
+                  nbrdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
                END DO
                !
@@ -564,192 +435,16 @@
                IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
                IF (ibr_max < nn_rimwidth(ib_bdy))   &
-                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
-            END DO
+                  CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
+            END DO
+            !
             CALL iom_close( inum )
+            DEALLOCATE( zz_read )
             !
-         ENDIF 
-         !
-      END DO      
-    
+         ENDIF
+         !
+      END DO
+
       ! 2. Now fill indices corresponding to straight open boundary arrays:
-      ! East
-      !-----
-      DO iseg = 1, nbdysege
-         ib_bdy = npckge(iseg)
-         !
-         ! ------------ T points -------------
-         igrd=1
-         icount=0
-         DO ir = 1, nn_rimwidth(ib_bdy)
-            DO ij = jpjedt(iseg), jpjeft(iseg)
-               icount = icount + 1
-               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
-               nbjdta(icount, igrd, ib_bdy) = ij
-               nbrdta(icount, igrd, ib_bdy) = ir
-            ENDDO
-         ENDDO
-         !
-         ! ------------ U points -------------
-         igrd=2
-         icount=0
-         DO ir = 1, nn_rimwidth(ib_bdy)
-            DO ij = jpjedt(iseg), jpjeft(iseg)
-               icount = icount + 1
-               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir
-               nbjdta(icount, igrd, ib_bdy) = ij
-               nbrdta(icount, igrd, ib_bdy) = ir
-            ENDDO
-         ENDDO
-         !
-         ! ------------ V points -------------
-         igrd=3
-         icount=0
-         DO ir = 1, nn_rimwidth(ib_bdy)
-!            DO ij = jpjedt(iseg), jpjeft(iseg) - 1
-            DO ij = jpjedt(iseg), jpjeft(iseg)
-               icount = icount + 1
-               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
-               nbjdta(icount, igrd, ib_bdy) = ij
-               nbrdta(icount, igrd, ib_bdy) = ir
-            ENDDO
-            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
-            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
-         ENDDO
-      ENDDO
-      !
-      ! West
-      !-----
-      DO iseg = 1, nbdysegw
-         ib_bdy = npckgw(iseg)
-         !
-         ! ------------ T points -------------
-         igrd=1
-         icount=0
-         DO ir = 1, nn_rimwidth(ib_bdy)
-            DO ij = jpjwdt(iseg), jpjwft(iseg)
-               icount = icount + 1
-               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
-               nbjdta(icount, igrd, ib_bdy) = ij
-               nbrdta(icount, igrd, ib_bdy) = ir
-            ENDDO
-         ENDDO
-         !
-         ! ------------ U points -------------
-         igrd=2
-         icount=0
-         DO ir = 1, nn_rimwidth(ib_bdy)
-            DO ij = jpjwdt(iseg), jpjwft(iseg)
-               icount = icount + 1
-               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
-               nbjdta(icount, igrd, ib_bdy) = ij
-               nbrdta(icount, igrd, ib_bdy) = ir
-            ENDDO
-         ENDDO
-         !
-         ! ------------ V points -------------
-         igrd=3
-         icount=0
-         DO ir = 1, nn_rimwidth(ib_bdy)
-!            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
-            DO ij = jpjwdt(iseg), jpjwft(iseg)
-               icount = icount + 1
-               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
-               nbjdta(icount, igrd, ib_bdy) = ij
-               nbrdta(icount, igrd, ib_bdy) = ir
-            ENDDO
-            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
-            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
-         ENDDO
-      ENDDO
-      !
-      ! North
-      !-----
-      DO iseg = 1, nbdysegn
-         ib_bdy = npckgn(iseg)
-         !
-         ! ------------ T points -------------
-         igrd=1
-         icount=0
-         DO ir = 1, nn_rimwidth(ib_bdy)
-            DO ii = jpindt(iseg), jpinft(iseg)
-               icount = icount + 1
-               nbidta(icount, igrd, ib_bdy) = ii
-               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
-               nbrdta(icount, igrd, ib_bdy) = ir
-            ENDDO
-         ENDDO
-         !
-         ! ------------ U points -------------
-         igrd=2
-         icount=0
-         DO ir = 1, nn_rimwidth(ib_bdy)
-!            DO ii = jpindt(iseg), jpinft(iseg) - 1
-            DO ii = jpindt(iseg), jpinft(iseg)
-               icount = icount + 1
-               nbidta(icount, igrd, ib_bdy) = ii
-               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
-               nbrdta(icount, igrd, ib_bdy) = ir
-            ENDDO
-            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
-            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
-         ENDDO
-         !
-         ! ------------ V points -------------
-         igrd=3
-         icount=0
-         DO ir = 1, nn_rimwidth(ib_bdy)
-            DO ii = jpindt(iseg), jpinft(iseg)
-               icount = icount + 1
-               nbidta(icount, igrd, ib_bdy) = ii
-               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir
-               nbrdta(icount, igrd, ib_bdy) = ir
-            ENDDO
-         ENDDO
-      ENDDO
-      !
-      ! South
-      !-----
-      DO iseg = 1, nbdysegs
-         ib_bdy = npckgs(iseg)
-         !
-         ! ------------ T points -------------
-         igrd=1
-         icount=0
-         DO ir = 1, nn_rimwidth(ib_bdy)
-            DO ii = jpisdt(iseg), jpisft(iseg)
-               icount = icount + 1
-               nbidta(icount, igrd, ib_bdy) = ii
-               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
-               nbrdta(icount, igrd, ib_bdy) = ir
-            ENDDO
-         ENDDO
-         !
-         ! ------------ U points -------------
-         igrd=2
-         icount=0
-         DO ir = 1, nn_rimwidth(ib_bdy)
-!            DO ii = jpisdt(iseg), jpisft(iseg) - 1
-            DO ii = jpisdt(iseg), jpisft(iseg)
-               icount = icount + 1
-               nbidta(icount, igrd, ib_bdy) = ii
-               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
-               nbrdta(icount, igrd, ib_bdy) = ir
-            ENDDO
-            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
-            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
-         ENDDO
-         !
-         ! ------------ V points -------------
-         igrd=3
-         icount=0
-         DO ir = 1, nn_rimwidth(ib_bdy)
-            DO ii = jpisdt(iseg), jpisft(iseg)
-               icount = icount + 1
-               nbidta(icount, igrd, ib_bdy) = ii
-               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
-               nbrdta(icount, igrd, ib_bdy) = ir
-            ENDDO
-         ENDDO
-      ENDDO
+      CALL bdy_coords_seg( nbidta, nbjdta, nbrdta )
 
       !  Deal with duplicated points
@@ -765,19 +460,19 @@
                      DO ib2 = 1, nblendta(igrd,ib_bdy2)
                         IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. &
-                        &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
-!                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', & 
-!                                                       &              nbidta(ib1, igrd, ib_bdy1),      & 
-!                                                       &              nbjdta(ib2, igrd, ib_bdy2)
+                           &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
+                           !                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', & 
+                           !                                                       &              nbidta(ib1, igrd, ib_bdy1),      & 
+                           !                                                       &              nbjdta(ib2, igrd, ib_bdy2)
                            ! keep only points with the lowest distance to boundary:
                            IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN
-                             nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2
-                             nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2
+                              nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2
+                              nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2
                            ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN
-                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
-                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
-                           ! Arbitrary choice if distances are the same:
+                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
+                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
+                              ! Arbitrary choice if distances are the same:
                            ELSE
-                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
-                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
+                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
+                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
                            ENDIF
                         END IF
@@ -811,7 +506,7 @@
                   IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
                      CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', &
-                          &        ' in order of distance from edge nbr A utility for re-ordering ', &
-                          &        ' boundary coordinates and data files exists in the TOOLS/OBC directory')
-                  ENDIF    
+                        &        ' in order of distance from edge nbr A utility for re-ordering ', &
+                        &        ' boundary coordinates and data files exists in the TOOLS/OBC directory')
+                  ENDIF
                ENDIF
                ! check if point is in local domain
@@ -911,8 +606,8 @@
                ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same weights
                idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 )      ! tanh formulation
-!               idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.  ! quadratic
-!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy))       ! linear
-            END DO
-         END DO 
+               !               idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.  ! quadratic
+               !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy))       ! linear
+            END DO
+         END DO
 
          ! Compute damping coefficients
@@ -922,9 +617,9 @@
                ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same damping coefficients
                idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 
-               & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
+                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
                idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) & 
-               & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
-            END DO
-         END DO 
+                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
+            END DO
+         END DO
 
       END DO ! ib_bdy
@@ -935,5 +630,5 @@
 
       ! ------------------------------------------
-      ! handel rim0, do as if rim 1 was free ocean
+      ! handle rim0, do as if rim 1 was free ocean
       ! ------------------------------------------
 
@@ -944,6 +639,6 @@
          DO ii = 2, jpim1
             zfmask(ii,ij) = ztmask(ii,ij  ) * ztmask(ii+1,ij  )   &
-           &              * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
-         END DO      
+               &              * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
+         END DO
       END DO
       CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. )
@@ -953,5 +648,5 @@
       ! bdytmask = 1  on the computational domain AND on open boundaries
       !          = 0  elsewhere   
- 
+
       bdytmask(:,:) = ssmask(:,:)
 
@@ -981,5 +676,5 @@
 
       ! ------------------------------------
-      ! handel rim1, do as if rim 0 was land
+      ! handle rim1, do as if rim 0 was land
       ! ------------------------------------
 
@@ -1000,6 +695,6 @@
          DO ii = 2, jpim1
             zfmask(ii,ij) = ztmask(ii,ij  ) * ztmask(ii+1,ij  )   &
-           &              * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
-         END DO      
+               &              * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
+         END DO
       END DO
       CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. )
@@ -1019,5 +714,4 @@
 
       CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. )   ! compute flagu, flagv, ntreat on rim 1
-
 
       !
@@ -1053,6 +747,6 @@
                IF( iibi == 0     .OR. ii1 == 0     .OR. ii2 == 0     .OR. ii3 == 0     )   lrecv_bdyint(ib_bdy,igrd,1,ir) = .true.
                IF( iibi == jpi+1 .OR. ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 )   lrecv_bdyint(ib_bdy,igrd,2,ir) = .true.  
-               IF( iibe == 0     )   lrecv_bdyext(ib_bdy,igrd,1,ir) = .true.
-               IF( iibe == jpi+1 )   lrecv_bdyext(ib_bdy,igrd,2,ir) = .true.  
+               IF( iibe == 0                                                           )   lrecv_bdyext(ib_bdy,igrd,1,ir) = .true.
+               IF( iibe == jpi+1                                                       )   lrecv_bdyext(ib_bdy,igrd,2,ir) = .true.  
                ! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo
                ! :¨¨¨¨¨|¨¨-->    |                                             |    <--¨¨|¨¨¨¨¨: 
@@ -1060,9 +754,9 @@
                ! :.....|_._:_____|   (1) W neighbour         E neighbour (2)   |_____:_._|.....:
                IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. &
-                 & ( iibi == 3     .OR. ii1 == 3     .OR. ii2 == 3     .OR. ii3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,1,ir) = .true.
+                  & ( iibi == 3     .OR. ii1 == 3     .OR. ii2 == 3     .OR. ii3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,1,ir)=.true.
                IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. &
-                 & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) )   lsend_bdyint(ib_bdy,igrd,2,ir) = .true.
-               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. iibe == 3    )   lsend_bdyext(ib_bdy,igrd,1,ir) = .true.
-               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2)   lsend_bdyext(ib_bdy,igrd,2,ir) = .true.
+                  & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) )   lsend_bdyint(ib_bdy,igrd,2,ir)=.true.
+               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. iibe == 3     )   lsend_bdyext(ib_bdy,igrd,1,ir)=.true.
+               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2 )   lsend_bdyext(ib_bdy,igrd,2,ir)=.true.
                !
                ! search neighbour in the north/south direction   
@@ -1073,6 +767,6 @@
                IF( ijbi == 0     .OR. ij1 == 0     .OR. ij2 == 0     .OR. ij3 == 0     )   lrecv_bdyint(ib_bdy,igrd,3,ir) = .true.
                IF( ijbi == jpj+1 .OR. ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 )   lrecv_bdyint(ib_bdy,igrd,4,ir) = .true.
-               IF( ijbe == 0     )   lrecv_bdyext(ib_bdy,igrd,3,ir) = .true.
-               IF( ijbe == jpj+1 )   lrecv_bdyext(ib_bdy,igrd,4,ir) = .true.
+               IF( ijbe == 0                                                           )   lrecv_bdyext(ib_bdy,igrd,3,ir) = .true.
+               IF( ijbe == jpj+1                                                       )   lrecv_bdyext(ib_bdy,igrd,4,ir) = .true.
                ! Check if neighbour has its rim parallel to its mpi subdomain     _________  border and next to its halo
                !   ^  |    o    |                                                :         : 
@@ -1080,9 +774,9 @@
                !      :_________:  (3) S neighbour          N neighbour (4)   v  |    o    |   
                IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. &
-                 & ( ijbi == 3     .OR. ij1 == 3     .OR. ij2 == 3     .OR. ij3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,3,ir) = .true.
+                  & ( ijbi == 3     .OR. ij1 == 3     .OR. ij2 == 3     .OR. ij3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,3,ir)=.true.
                IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. &
-                 & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) )   lsend_bdyint(ib_bdy,igrd,4,ir) = .true.
-               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. ijbe == 3    )   lsend_bdyext(ib_bdy,igrd,3,ir) = .true.
-               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2)   lsend_bdyext(ib_bdy,igrd,4,ir) = .true.
+                  & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) )   lsend_bdyint(ib_bdy,igrd,4,ir)=.true.
+               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. ijbe == 3     )   lsend_bdyext(ib_bdy,igrd,3,ir)=.true.
+               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2 )   lsend_bdyext(ib_bdy,igrd,4,ir)=.true.
             END DO
          END DO
@@ -1091,6 +785,6 @@
       DO ib_bdy = 1,nb_bdy
          IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. &
-           & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. &
-           & cn_tra(ib_bdy)   == 'orlanski' .OR. cn_tra(ib_bdy)   == 'orlanski_npo'      ) THEN
+            & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. &
+            & cn_tra(ib_bdy)   == 'orlanski' .OR. cn_tra(ib_bdy)   == 'orlanski_npo'      ) THEN
             DO igrd = 1, jpbgrd
                DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
@@ -1110,7 +804,7 @@
       ! Tidy up
       !--------
-      IF( nb_bdy>0 )   DEALLOCATE( nbidta, nbjdta, nbrdta )
-      !
-   END SUBROUTINE bdy_segs
+      DEALLOCATE( nbidta, nbjdta, nbrdta )
+      !
+   END SUBROUTINE bdy_def
 
 
@@ -1120,10 +814,10 @@
       !!
       !! ** Purpose :   Initialize structures ( flagu, flagv, ntreat ) indicating how rim points
-      !!                  are to be handeled in the boundary condition treatment
-      !!
-      !! ** Method  :   - to handel rim 0 zmasks must indicate ocean points      (set at one on rim 0 and rim 1 and interior)
+      !!                  are to be handled in the boundary condition treatment
+      !!
+      !! ** Method  :   - to handle rim 0 zmasks must indicate ocean points      (set at one on rim 0 and rim 1 and interior)
       !!                            and bdymasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
       !!                    (as if rim 1 was free ocean)
-      !!                - to handel rim 1 zmasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
+      !!                - to handle rim 1 zmasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
       !!                            and bdymasks must indicate free ocean points (set at one on interior)
       !!                    (as if rim 0 was land)
@@ -1373,7 +1067,94 @@
       CASE( 16 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij
       END SELECT
-    END SUBROUTINE find_neib
+   END SUBROUTINE find_neib
     
 
+   SUBROUTINE bdy_read_seg( kb_bdy, knblendta ) 
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE bdy_coords_seg  ***
+      !!
+      !! ** Purpose :  build bdy coordinates with segments defined in namelist
+      !!
+      !! ** Method  :  read namelist nambdy_index blocks
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                   , INTENT (in   ) ::   kb_bdy           ! bdy number
+      INTEGER, DIMENSION(jpbgrd), INTENT (  out) ::   knblendta        ! length of index arrays 
+      !!
+      INTEGER          ::   ios                 ! Local integer output status for namelist read
+      INTEGER          ::   nbdyind, nbdybeg, nbdyend
+      CHARACTER(LEN=1) ::   ctypebdy   !     -        - 
+      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
+      !!----------------------------------------------------------------------
+
+      ! No REWIND here because may need to read more than one nambdy_index namelist.
+      ! Read only namelist_cfg to avoid unseccessfull overwrite 
+      ! keep full control of the configuration namelist
+      READ  ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 )
+904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nambdy_index )
+      
+      SELECT CASE ( TRIM(ctypebdy) )
+      CASE( 'N' )
+         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
+            nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain.
+            nbdybeg  = 2
+            nbdyend  = jpiglo - 1
+         ENDIF
+         nbdysegn = nbdysegn + 1
+         npckgn(nbdysegn) = kb_bdy ! Save bdy package number
+         jpjnob(nbdysegn) = nbdyind
+         jpindt(nbdysegn) = nbdybeg
+         jpinft(nbdysegn) = nbdyend
+         !
+      CASE( 'S' )
+         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
+            nbdyind  = 2           ! set boundary to whole side of model domain.
+            nbdybeg  = 2
+            nbdyend  = jpiglo - 1
+         ENDIF
+         nbdysegs = nbdysegs + 1
+         npckgs(nbdysegs) = kb_bdy ! Save bdy package number
+         jpjsob(nbdysegs) = nbdyind
+         jpisdt(nbdysegs) = nbdybeg
+         jpisft(nbdysegs) = nbdyend
+         !
+      CASE( 'E' )
+         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
+            nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain.
+            nbdybeg  = 2
+            nbdyend  = jpjglo - 1
+         ENDIF
+         nbdysege = nbdysege + 1 
+         npckge(nbdysege) = kb_bdy ! Save bdy package number
+         jpieob(nbdysege) = nbdyind
+         jpjedt(nbdysege) = nbdybeg
+         jpjeft(nbdysege) = nbdyend
+         !
+      CASE( 'W' )
+         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
+            nbdyind  = 2           ! set boundary to whole side of model domain.
+            nbdybeg  = 2
+            nbdyend  = jpjglo - 1
+         ENDIF
+         nbdysegw = nbdysegw + 1
+         npckgw(nbdysegw) = kb_bdy ! Save bdy package number
+         jpiwob(nbdysegw) = nbdyind
+         jpjwdt(nbdysegw) = nbdybeg
+         jpjwft(nbdysegw) = nbdyend
+         !
+      CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
+      END SELECT
+      
+      ! For simplicity we assume that in case of straight bdy, arrays have the same length
+      ! (even if it is true that last tangential velocity points
+      ! are useless). This simplifies a little bit boundary data format (and agrees with format
+      ! used so far in obc package)
+      
+      knblendta(1:jpbgrd) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy)
+      
+   END SUBROUTINE bdy_read_seg
+
+   
    SUBROUTINE bdy_ctl_seg
       !!----------------------------------------------------------------------
@@ -1719,5 +1500,206 @@
    END SUBROUTINE bdy_ctl_seg
 
-
+   
+   SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta ) 
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE bdy_coords_seg  ***
+      !!
+      !! ** Purpose :  build nbidta, nbidta, nbrdta for bdy built with segments
+      !!
+      !! ** Method  :  
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, DIMENSION(:,:,:), intent(  out)  ::   nbidta, nbjdta, nbrdta   ! Index arrays: i and j indices of bdy dta
+      !!
+      INTEGER  ::   ii, ij, ir, iseg
+      INTEGER  ::   igrd         ! grid type (t=1, u=2, v=3)
+      INTEGER  ::   icount       ! 
+      INTEGER  ::   ib_bdy       ! bdy number
+      !!----------------------------------------------------------------------
+
+      ! East
+      !-----
+      DO iseg = 1, nbdysege
+         ib_bdy = npckge(iseg)
+         !
+         ! ------------ T points -------------
+         igrd=1
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ij = jpjedt(iseg), jpjeft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
+               nbjdta(icount, igrd, ib_bdy) = ij
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+         !
+         ! ------------ U points -------------
+         igrd=2
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ij = jpjedt(iseg), jpjeft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir
+               nbjdta(icount, igrd, ib_bdy) = ij
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+         !
+         ! ------------ V points -------------
+         igrd=3
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            !            DO ij = jpjedt(iseg), jpjeft(iseg) - 1
+            DO ij = jpjedt(iseg), jpjeft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
+               nbjdta(icount, igrd, ib_bdy) = ij
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+         ENDDO
+      ENDDO
+      !
+      ! West
+      !-----
+      DO iseg = 1, nbdysegw
+         ib_bdy = npckgw(iseg)
+         !
+         ! ------------ T points -------------
+         igrd=1
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ij = jpjwdt(iseg), jpjwft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
+               nbjdta(icount, igrd, ib_bdy) = ij
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+         !
+         ! ------------ U points -------------
+         igrd=2
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ij = jpjwdt(iseg), jpjwft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
+               nbjdta(icount, igrd, ib_bdy) = ij
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+         !
+         ! ------------ V points -------------
+         igrd=3
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            !            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
+            DO ij = jpjwdt(iseg), jpjwft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
+               nbjdta(icount, igrd, ib_bdy) = ij
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+         ENDDO
+      ENDDO
+      !
+      ! North
+      !-----
+      DO iseg = 1, nbdysegn
+         ib_bdy = npckgn(iseg)
+         !
+         ! ------------ T points -------------
+         igrd=1
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ii = jpindt(iseg), jpinft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = ii
+               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+         !
+         ! ------------ U points -------------
+         igrd=2
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            !            DO ii = jpindt(iseg), jpinft(iseg) - 1
+            DO ii = jpindt(iseg), jpinft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = ii
+               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+         ENDDO
+         !
+         ! ------------ V points -------------
+         igrd=3
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ii = jpindt(iseg), jpinft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = ii
+               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      ! South
+      !-----
+      DO iseg = 1, nbdysegs
+         ib_bdy = npckgs(iseg)
+         !
+         ! ------------ T points -------------
+         igrd=1
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ii = jpisdt(iseg), jpisft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = ii
+               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+         !
+         ! ------------ U points -------------
+         igrd=2
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            !            DO ii = jpisdt(iseg), jpisft(iseg) - 1
+            DO ii = jpisdt(iseg), jpisft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = ii
+               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+         ENDDO
+         !
+         ! ------------ V points -------------
+         igrd=3
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ii = jpisdt(iseg), jpisft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = ii
+               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+      ENDDO
+
+      
+   END SUBROUTINE bdy_coords_seg
+   
+   
    SUBROUTINE bdy_ctl_corn( ib1, ib2 )
       !!----------------------------------------------------------------------
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdytides.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdytides.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdytides.F90	(revision 11223)
@@ -70,5 +70,4 @@
       INTEGER                                   ::   inum, igrd
       INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays)
-      INTEGER, POINTER, DIMENSION(:)            ::   nblen, nblenrim     ! short cuts
       INTEGER                                   ::   ios                 ! Local integer output status for namelist read
       CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file 
@@ -77,14 +76,11 @@
       !!
       TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
-      TYPE(MAP_POINTER), DIMENSION(jpbgrd)      ::   ibmap_ptr           !: array of pointers to nbmap
       !!
       NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
       !!----------------------------------------------------------------------
       !
-      IF (nb_bdy>0) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
-      ENDIF
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
 
       REWIND(numnam_cfg)
@@ -94,13 +90,12 @@
             !
             td => tides(ib_bdy)
-            nblen => idx_bdy(ib_bdy)%nblen
-            nblenrim => idx_bdy(ib_bdy)%nblenrim
 
             ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
             filtide(:) = ''
 
-            ! Don't REWIND here - may need to read more than one of these namelists. 
+            REWIND( numnam_ref )
             READ  ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
 901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_tide in reference namelist', lwp )
+            ! Don't REWIND here - may need to read more than one of these namelists. 
             READ  ( numnam_cfg, nambdy_tide, IOSTAT = ios, ERR = 902 )
 902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist', lwp )
@@ -125,6 +120,6 @@
             ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
             ! relaxation area      
-            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:)
-            ELSE                                   ;   ilen0(:) = nblenrim(:)
+            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = idx_bdy(ib_bdy)%nblen   (:)
+            ELSE                                   ;   ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:)
             ENDIF
 
@@ -210,9 +205,4 @@
                ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
                !
-               ! Set map structure
-               ibmap_ptr(1)%ptr => idx_bdy(ib_bdy)%nbmap(:,1)   ;   ibmap_ptr(1)%ll_unstruc = ln_coords_file(ib_bdy)
-               ibmap_ptr(2)%ptr => idx_bdy(ib_bdy)%nbmap(:,2)   ;   ibmap_ptr(2)%ll_unstruc = ln_coords_file(ib_bdy)
-               ibmap_ptr(3)%ptr => idx_bdy(ib_bdy)%nbmap(:,3)   ;   ibmap_ptr(3)%ll_unstruc = ln_coords_file(ib_bdy)
-
                ! Open files and read in tidal forcing data
                ! -----------------------------------------
@@ -222,7 +212,7 @@
                   clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
                   CALL iom_open( clfile, inum )
-                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1,  ibmap_ptr(1) )
+                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
                   td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
-                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1,  ibmap_ptr(1) )
+                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
                   td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
                   CALL iom_close( inum )
@@ -230,7 +220,7 @@
                   clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
                   CALL iom_open( clfile, inum )
-                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, ibmap_ptr(2) )
+                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
                   td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
-                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, ibmap_ptr(2) )
+                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
                   td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
                   CALL iom_close( inum )
@@ -238,7 +228,7 @@
                   clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
                   CALL iom_open( clfile, inum )
-                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, ibmap_ptr(3) )
+                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
                   td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
-                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, ibmap_ptr(3) )
+                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
                   td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
                   CALL iom_close( inum )
@@ -272,5 +262,5 @@
 
 
-   SUBROUTINE bdytide_update( kt, idx, dta, td, jit, time_offset )
+   SUBROUTINE bdytide_update( kt, idx, dta, td, kit, kt_offset )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdytide_update  ***
@@ -283,10 +273,10 @@
       TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data
       TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data
-      INTEGER, OPTIONAL, INTENT(in   ) ::   jit         ! Barotropic timestep counter (for timesplitting option)
-      INTEGER, OPTIONAL, INTENT(in   ) ::   time_offset ! time offset in units of timesteps. NB. if jit
+      INTEGER, OPTIONAL, INTENT(in   ) ::   kit         ! Barotropic timestep counter (for timesplitting option)
+      INTEGER, OPTIONAL, INTENT(in   ) ::   kt_offset   ! time offset in units of timesteps. NB. if kit
       !                                                 ! is present then units = subcycle timesteps.
-      !                                                 ! time_offset = 0  => get data at "now"    time level
-      !                                                 ! time_offset = -1 => get data at "before" time level
-      !                                                 ! time_offset = +1 => get data at "after"  time level
+      !                                                 ! kt_offset = 0  => get data at "now"    time level
+      !                                                 ! kt_offset = -1 => get data at "before" time level
+      !                                                 ! kt_offset = +1 => get data at "after"  time level
       !                                                 ! etc.
       !
@@ -303,6 +293,6 @@
 
       zflag=1
-      IF ( PRESENT(jit) ) THEN
-        IF ( jit /= 1 ) zflag=0
+      IF ( PRESENT(kit) ) THEN
+        IF ( kit /= 1 ) zflag=0
       ENDIF
 
@@ -323,10 +313,10 @@
 
       time_add = 0
-      IF( PRESENT(time_offset) ) THEN
-         time_add = time_offset
+      IF( PRESENT(kt_offset) ) THEN
+         time_add = kt_offset
       ENDIF
          
-      IF( PRESENT(jit) ) THEN  
-         z_arg = ((kt-kt_tide) * rdt + (jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) )
+      IF( PRESENT(kit) ) THEN  
+         z_arg = ((kt-kt_tide) * rdt + (kit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) )
       ELSE                              
          z_arg = ((kt-kt_tide)+time_add) * rdt
@@ -361,5 +351,5 @@
 
 
-   SUBROUTINE bdy_dta_tides( kt, kit, time_offset )
+   SUBROUTINE bdy_dta_tides( kt, kit, kt_offset )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_dta_tides  ***
@@ -370,9 +360,9 @@
       INTEGER,           INTENT(in) ::   kt          ! Main timestep counter
       INTEGER, OPTIONAL, INTENT(in) ::   kit         ! Barotropic timestep counter (for timesplitting option)
-      INTEGER, OPTIONAL, INTENT(in) ::   time_offset ! time offset in units of timesteps. NB. if kit
+      INTEGER, OPTIONAL, INTENT(in) ::   kt_offset   ! time offset in units of timesteps. NB. if kit
       !                                              ! is present then units = subcycle timesteps.
-      !                                              ! time_offset = 0  => get data at "now"    time level
-      !                                              ! time_offset = -1 => get data at "before" time level
-      !                                              ! time_offset = +1 => get data at "after"  time level
+      !                                              ! kt_offset = 0  => get data at "now"    time level
+      !                                              ! kt_offset = -1 => get data at "before" time level
+      !                                              ! kt_offset = +1 => get data at "after"  time level
       !                                              ! etc.
       !
@@ -389,6 +379,6 @@
 
       time_add = 0
-      IF( PRESENT(time_offset) ) THEN
-         time_add = time_offset
+      IF( PRESENT(kt_offset) ) THEN
+         time_add = kt_offset
       ENDIF
       
@@ -435,7 +425,7 @@
             ! If time splitting, initialize arrays from slow varying open boundary data:
             IF ( PRESENT(kit) ) THEN           
-               IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
-               IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
-               IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
+               IF ( dta_bdy(ib_bdy)%lneed_ssh   ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
+               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
+               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
             ENDIF
             !
@@ -447,5 +437,5 @@
                z_sist = zramp * SIN( z_sarg )
                !
-               IF ( dta_bdy(ib_bdy)%ll_ssh ) THEN
+               IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN
                   igrd=1                              ! SSH on tracer grid
                   DO ib = 1, ilen0(igrd)
@@ -456,5 +446,5 @@
                ENDIF
                !
-               IF ( dta_bdy(ib_bdy)%ll_u2d ) THEN
+               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
                   igrd=2                              ! U grid
                   DO ib = 1, ilen0(igrd)
@@ -463,7 +453,4 @@
                         &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
                   END DO
-               ENDIF
-               !
-               IF ( dta_bdy(ib_bdy)%ll_v2d ) THEN
                   igrd=3                              ! V grid
                   DO ib = 1, ilen0(igrd) 
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DOM/dommsk.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DOM/dommsk.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DOM/dommsk.F90	(revision 11223)
@@ -101,5 +101,5 @@
          &             cn_ice, nn_ice_dta,                                     &
          &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
-         &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
+         &             ln_vol, nn_volctl, nn_rimwidth
       !!---------------------------------------------------------------------
       !
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90	(revision 11223)
@@ -718,6 +718,6 @@
          !                                                !  ------------------
          ! Update only tidal forcing at open boundaries
-         IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
-         IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset   )
+         IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, kt_offset= noffset+1 )
+         IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, kt_offset= noffset   )
          !
          ! Set extrapolation coefficients for predictor step:
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/IOM/iom.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/IOM/iom.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/IOM/iom.F90	(revision 11223)
@@ -835,5 +835,5 @@
 
 
-   FUNCTION iom_varid ( kiomid, cdvar, kdimsz, kndims, ldstop )  
+   FUNCTION iom_varid ( kiomid, cdvar, kdimsz, kndims, lduld, ldstop )  
       !!-----------------------------------------------------------------------
       !!                  ***  FUNCTION  iom_varid  ***
@@ -844,5 +844,6 @@
       CHARACTER(len=*)     , INTENT(in   )           ::   cdvar    ! name of the variable
       INTEGER, DIMENSION(:), INTENT(  out), OPTIONAL ::   kdimsz   ! size of each dimension
-      INTEGER,               INTENT(  out), OPTIONAL ::   kndims   ! size of the dimensions
+      INTEGER              , INTENT(  out), OPTIONAL ::   kndims   ! number of dimensions
+      LOGICAL              , INTENT(  out), OPTIONAL ::   lduld    ! true if the last dimension is unlimited (time)
       LOGICAL              , INTENT(in   ), OPTIONAL ::   ldstop   ! stop if looking for non-existing variable (default = .TRUE.)
       !
@@ -874,5 +875,5 @@
                iiv = iiv + 1
                IF( iiv <= jpmax_vars ) THEN
-                  iom_varid = iom_nf90_varid( kiomid, cdvar, iiv, kdimsz, kndims )
+                  iom_varid = iom_nf90_varid( kiomid, cdvar, iiv, kdimsz, kndims, lduld )
                ELSE
                   CALL ctl_stop( trim(clinfo), 'Too many variables in the file '//iom_file(kiomid)%name,   &
@@ -892,4 +893,5 @@
                ENDIF
                IF( PRESENT(kndims) )  kndims = iom_file(kiomid)%ndims(iiv)
+               IF( PRESENT( lduld) )  lduld  = iom_file(kiomid)%luld( iiv)
             ENDIF
          ENDIF
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/IOM/iom_nf90.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/IOM/iom_nf90.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/IOM/iom_nf90.F90	(revision 11223)
@@ -187,5 +187,5 @@
 
 
-   FUNCTION iom_nf90_varid ( kiomid, cdvar, kiv, kdimsz, kndims )  
+   FUNCTION iom_nf90_varid ( kiomid, cdvar, kiv, kdimsz, kndims, lduld )  
       !!-----------------------------------------------------------------------
       !!                  ***  FUNCTION  iom_varid  ***
@@ -198,4 +198,5 @@
       INTEGER, DIMENSION(:), INTENT(  out), OPTIONAL ::   kdimsz   ! size of the dimensions
       INTEGER,               INTENT(  out), OPTIONAL ::   kndims   ! size of the dimensions
+      LOGICAL              , INTENT(  out), OPTIONAL ::   lduld    ! true if the last dimension is unlimited (time)
       !
       INTEGER                        ::   iom_nf90_varid   ! iom variable Id
@@ -251,4 +252,5 @@
          ENDIF
          IF( PRESENT(kndims) )  kndims = iom_file(kiomid)%ndims(kiv)
+         IF( PRESENT( lduld) )  lduld  = iom_file(kiomid)%luld(kiv)
       ELSE  
          iom_nf90_varid = -1   !   variable not found, return error code: -1
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/LBC/mppini.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/LBC/mppini.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/LBC/mppini.F90	(revision 11223)
@@ -167,5 +167,5 @@
            &             cn_ice, nn_ice_dta,                                     &
            &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
-           &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
+           &             ln_vol, nn_volctl, nn_rimwidth
       !!----------------------------------------------------------------------
 
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/fldread.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/fldread.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/fldread.F90	(revision 11223)
@@ -80,12 +80,10 @@
       INTEGER                         ::   nreclast     ! last record to be read in the current file
       CHARACTER(len = 256)            ::   lsmname      ! current name of the NetCDF mask file acting as a key
+      !                                                 ! variables related to BDY
       INTEGER                         ::   igrd         ! grid type for bdy data
       INTEGER                         ::   ibdy         ! bdy set id number
+      INTEGER, POINTER, DIMENSION(:)  ::   imap         ! Array of integer pointers to 1D arrays
+      LOGICAL                         ::   ltotvel      ! total velocity or not (T/F)
    END TYPE FLD
-
-   TYPE, PUBLIC ::   MAP_POINTER      !: Map from input data file to local domain
-      INTEGER, POINTER, DIMENSION(:)  ::  ptr           ! Array of integer pointers to 1D arrays
-      LOGICAL                         ::  ll_unstruc    ! Unstructured (T) or structured (F) boundary data file
-   END TYPE MAP_POINTER
 
 !$AGRIF_DO_NOT_TREAT
@@ -129,5 +127,5 @@
 CONTAINS
 
-   SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl )
+   SUBROUTINE fld_read( kt, kn_fsbc, sd, kit, kt_offset )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_read  ***
@@ -144,5 +142,4 @@
       INTEGER  , INTENT(in   )               ::   kn_fsbc   ! sbc computation period (in time step) 
       TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables
-      TYPE(MAP_POINTER),INTENT(in), OPTIONAL, DIMENSION(:) ::   map   ! global-to-local mapping indices
       INTEGER  , INTENT(in   ), OPTIONAL     ::   kit       ! subcycle timestep for timesplitting option
       INTEGER  , INTENT(in   ), OPTIONAL     ::   kt_offset ! provide fields at time other than "now"
@@ -150,6 +147,4 @@
       !                                                     !   kt_offset = +1 => fields at "after"  time level
       !                                                     !   etc.
-      INTEGER  , INTENT(in   ), OPTIONAL     ::   jpk_bdy   ! number of vertical levels in the BDY data
-      LOGICAL  , INTENT(in   ), OPTIONAL     ::   fvl   ! number of vertical levels in the BDY data
       !!
       INTEGER  ::   itmp         ! local variable
@@ -166,5 +161,4 @@
       REAL(wp) ::   ztintb       ! ratio applied to before records when doing time interpolation
       CHARACTER(LEN=1000) ::   clfmt  ! write format
-      TYPE(MAP_POINTER)   ::   imap   ! global-to-local mapping indices
       !!---------------------------------------------------------------------
       ll_firstcall = kt == nit000
@@ -175,6 +169,4 @@
       ENDIF
       IF( PRESENT(kt_offset) )   it_offset = kt_offset
-
-      imap%ptr => NULL()
 
       ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar 
@@ -188,10 +180,6 @@
       IF( ll_firstcall ) THEN                      ! initialization
          DO jf = 1, imf 
-            IF( PRESENT(map) ) imap = map(jf)
-               IF( PRESENT(jpk_bdy) ) THEN
-                  CALL fld_init( kn_fsbc, sd(jf), imap, jpk_bdy, fvl )  ! read each before field (put them in after as they will be swapped)
-               ELSE
-                  CALL fld_init( kn_fsbc, sd(jf), imap )  ! read each before field (put them in after as they will be swapped)
-               ENDIF
+            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE
+            CALL fld_init( kn_fsbc, sd(jf) )       ! read each before field (put them in after as they will be swapped)
          END DO
          IF( lwp ) CALL wgt_print()                ! control print
@@ -202,8 +190,8 @@
          !
          DO jf = 1, imf                            ! ---   loop over field   --- !
-            
+
+            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE
+                      
             IF( isecsbc > sd(jf)%nrec_a(2) .OR. ll_firstcall ) THEN    ! read/update the after data?
-
-               IF( PRESENT(map) )   imap = map(jf)   ! temporary definition of map
 
                sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:)                                  ! swap before record informations
@@ -222,5 +210,5 @@
                   itmp = sd(jf)%nrec_a(1)                       ! temporary storage
                   sd(jf)%nrec_a(1) = sd(jf)%nreclast            ! read the last record of the file currently opened
-                  CALL fld_get( sd(jf), imap )                  ! read after data
+                  CALL fld_get( sd(jf) )                        ! read after data
                   sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! re-swap before record field
                   sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1)           ! update before record informations
@@ -240,5 +228,5 @@
                      &                   .AND. sd(jf)%nrec_b(1) /= sd(jf)%nrec_a(1) - 1 ) THEN   
                      sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - 1       ! move back to before record
-                     CALL fld_get( sd(jf), imap )                  ! read after data
+                     CALL fld_get( sd(jf) )                        ! read after data
                      sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! re-swap before record field
                      sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1)           ! update before record informations
@@ -285,9 +273,6 @@
                   
                ! read after data
-               IF( PRESENT(jpk_bdy) ) THEN
-                  CALL fld_get( sd(jf), imap, jpk_bdy, fvl)
-               ELSE
-                  CALL fld_get( sd(jf), imap )
-               ENDIF
+               CALL fld_get( sd(jf) )
+               
             ENDIF   ! read new data?
          END DO                                    ! --- end loop over field --- !
@@ -296,4 +281,6 @@
 
          DO jf = 1, imf                            ! ---   loop over field   --- !
+            !
+            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE
             !
             IF( sd(jf)%ln_tint ) THEN              ! temporal interpolation
@@ -327,5 +314,5 @@
 
 
-   SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy, fvl)
+   SUBROUTINE fld_init( kn_fsbc, sdjf )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_init  ***
@@ -336,7 +323,4 @@
       INTEGER  , INTENT(in   ) ::   kn_fsbc      ! sbc computation period (in time step) 
       TYPE(FLD), INTENT(inout) ::   sdjf         ! input field related variables
-      TYPE(MAP_POINTER),INTENT(in) ::   map      ! global-to-local mapping indices
-      INTEGER  , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data
-      LOGICAL  , INTENT(in), OPTIONAL :: fvl     ! number of vertical levels in the BDY data
       !!
       LOGICAL :: llprevyr              ! are we reading previous year  file?
@@ -351,4 +335,5 @@
       CHARACTER(LEN=1000) ::   clfmt   ! write format
       !!---------------------------------------------------------------------
+      !
       llprevyr   = .FALSE.
       llprevmth  = .FALSE.
@@ -433,9 +418,5 @@
          !
          ! read before data in after arrays(as we will swap it later)
-         IF( PRESENT(jpk_bdy) ) THEN
-            CALL fld_get( sdjf, map, jpk_bdy, fvl )
-         ELSE
-            CALL fld_get( sdjf, map )
-         ENDIF
+         CALL fld_get( sdjf )
          !
          clfmt = "('   fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')"
@@ -613,5 +594,5 @@
 
 
-   SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl )
+   SUBROUTINE fld_get( sdjf )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_get  ***
@@ -620,7 +601,4 @@
       !!----------------------------------------------------------------------
       TYPE(FLD)        , INTENT(inout) ::   sdjf   ! input field related variables
-      TYPE(MAP_POINTER), INTENT(in   ) ::   map    ! global-to-local mapping indices
-      INTEGER  , INTENT(in), OPTIONAL  ::   jpk_bdy ! number of vertical levels in the bdy data
-      LOGICAL  , INTENT(in), OPTIONAL  ::   fvl     ! number of vertical levels in the bdy data
       !
       INTEGER ::   ipk      ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
@@ -634,16 +612,10 @@
       ipk = SIZE( sdjf%fnow, 3 )
       !
-      IF( ASSOCIATED(map%ptr) ) THEN
-         IF( PRESENT(jpk_bdy) ) THEN
-            IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2),                &
-                                                        sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl )
-            ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ),                &
-                                                        sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl )
-            ENDIF
-         ELSE
-            IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map )
-            ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map )
-            ENDIF
-         ENDIF        
+      IF( ASSOCIATED(sdjf%imap) ) THEN
+         IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2),                &
+            &                                        sdjf%nrec_a(1), sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel )
+         ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ),                &
+            &                                        sdjf%nrec_a(1), sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel )
+         ENDIF
       ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
          CALL wgt_list( sdjf, iw )
@@ -700,5 +672,6 @@
    END SUBROUTINE fld_get
 
-   SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl )
+   
+   SUBROUTINE fld_map( knum, cdvar, pdta, krec, kmap, kgrd, kbdy, ldtotvel )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_map  ***
@@ -707,136 +680,139 @@
       !!                using a general mapping (for open boundaries)
       !!----------------------------------------------------------------------
-
-      USE bdy_oce, ONLY: ln_bdy, idx_bdy, dta_global, dta_global_z, dta_global_dz, dta_global2, dta_global2_z, dta_global2_dz                 ! workspace to read in global data arrays
-
-      INTEGER                   , INTENT(in ) ::   num     ! stream number
-      CHARACTER(LEN=*)          , INTENT(in ) ::   clvar   ! variable name
-      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta     ! output field on model grid (2 dimensional)
-      INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice)
-      TYPE(MAP_POINTER)         , INTENT(in ) ::   map     ! global-to-local mapping indices
-      INTEGER  , INTENT(in), OPTIONAL         ::   igrd, ibdy, jpk_bdy  ! grid type, set number and number of vertical levels in the bdy data
-      LOGICAL  , INTENT(in), OPTIONAL         ::   fvl     ! grid type, set number and number of vertical levels in the bdy data
-      INTEGER                                 ::   jpkm1_bdy! number of vertical levels in the bdy data minus 1
-      !!
-      INTEGER                                 ::   ipi      ! length of boundary data on local process
-      INTEGER                                 ::   ipj      ! length of dummy dimension ( = 1 )
-      INTEGER                                 ::   ipk      ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
-      INTEGER                                 ::   ilendta  ! length of data in file
-      INTEGER                                 ::   idvar    ! variable ID
-      INTEGER                                 ::   ib, ik, ji, jj   ! loop counters
-      INTEGER                                 ::   ierr
-      REAL(wp)                                ::   fv          ! fillvalue 
-      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read    ! work space for global data
-      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read_z  ! work space for global data
-      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read_dz ! work space for global data
-      !!---------------------------------------------------------------------
-      !
-      ipi = SIZE( dta, 1 )
-      ipj = 1
-      ipk = SIZE( dta, 3 )
-      !
-      idvar   = iom_varid( num, clvar )
-      ilendta = iom_file(num)%dimsz(1,idvar)
-
-      IF ( ln_bdy ) THEN
-         ipj = iom_file(num)%dimsz(2,idvar)
-         IF( map%ll_unstruc) THEN   ! unstructured open boundary data file
-            dta_read => dta_global
-            IF( PRESENT(jpk_bdy) ) THEN
-               IF( jpk_bdy>0 ) THEN
-                  dta_read_z => dta_global_z
-                  dta_read_dz => dta_global_dz
-                  jpkm1_bdy = jpk_bdy-1
-               ENDIF
-            ENDIF
-         ELSE                       ! structured open boundary file
-            dta_read => dta_global2
-            IF( PRESENT(jpk_bdy) ) THEN
-               IF( jpk_bdy>0 ) THEN
-                  dta_read_z => dta_global2_z
-                  dta_read_dz => dta_global2_dz
-                  jpkm1_bdy = jpk_bdy-1
-               ENDIF
-            ENDIF
-         ENDIF
-      ENDIF
-
-      IF(lwp) WRITE(numout,*) 'Dim size for ',        TRIM(clvar),' is ', ilendta
-      IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk
-      !
-      SELECT CASE( ipk )
-      CASE(1)        ;   
-      CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec )
-         IF ( map%ll_unstruc) THEN ! unstructured open boundary data file
-            DO ib = 1, ipi
-              DO ik = 1, ipk
-                dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik)
-              END DO
-            END DO
-         ELSE ! we assume that this is a structured open boundary file
-            DO ib = 1, ipi
-               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
-               ji=map%ptr(ib)-(jj-1)*ilendta
-               DO ik = 1, ipk
-                  dta(ib,1,ik) =  dta_read(ji,jj,ik)
-               END DO
-            END DO
-         ENDIF
+      INTEGER                   , INTENT(in   ) ::   knum         ! stream number
+      CHARACTER(LEN=*)          , INTENT(in   ) ::   cdvar        ! variable name
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta         ! bdy output field on model grid
+      INTEGER                   , INTENT(in   ) ::   krec         ! record number to read (ie time slice)
+      INTEGER , DIMENSION(:)    , INTENT(in   ) ::   kmap         ! global-to-local bdy mapping indices
+      ! optional variables used for vertical interpolation:
+      INTEGER, OPTIONAL         , INTENT(in   ) ::   kgrd         ! grid type (t, u, v)
+      INTEGER, OPTIONAL         , INTENT(in   ) ::   kbdy         ! bdy number
+      LOGICAL, OPTIONAL         , INTENT(in   ) ::   ldtotvel     ! true if toal ( = barotrop + barocline) velocity
+      !!
+      INTEGER                                   ::   ipi          ! length of boundary data on local process
+      INTEGER                                   ::   ipj          ! length of dummy dimension ( = 1 )
+      INTEGER                                   ::   ipk          ! number of vertical levels of pdta ( 2D: ipk=1 ; 3D: ipk=jpk )
+      INTEGER                                   ::   ipkb         ! number of vertical levels in boundary data file
+      INTEGER                                   ::   idvar        ! variable ID
+      INTEGER                                   ::   indims       ! number of dimensions of the variable
+      INTEGER, DIMENSION(4)                     ::   idimsz       ! size of variable dimensions 
+      REAL(wp)                                  ::   zfv          ! fillvalue 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zz_read      ! work space for global boundary data
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read    ! work space local data requiring vertical interpolation
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read_z  ! work space local data requiring vertical interpolation
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read_dz ! work space local data requiring vertical interpolation
+      CHARACTER(LEN=1),DIMENSION(3)             ::   clgrid
+      LOGICAL                                   ::   lluld        ! is the variable using the unlimited dimension
+      !!---------------------------------------------------------------------
+      !
+      clgrid = (/'t','u','v'/)
+      !
+      ipi = SIZE( pdta, 1 )
+      ipj = SIZE( pdta, 2 )   ! must be equal to 1
+      ipk = SIZE( pdta, 3 )
+      !
+      idvar = iom_varid( knum, cdvar, kndims = indims, kdimsz = idimsz, lduld = lluld  )
+      IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN   ;   ipkb = idimsz(3)   ! xy(zl)t or xy(zl)
+      ELSE                                                            ;   ipkb = 1           ! xy or xyt
+      ENDIF
+      !
+      ALLOCATE( zz_read( idimsz(1), idimsz(2), ipkb ) )  ! ++++++++ !!! this can be very big...         
+      !
+      IF( ipk == 1 ) THEN
+
+         IF( ipkb /= 1 ) CALL ctl_stop( 'fld_map : we must have ipkb = 1 to read surface data' )
+         CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,1), krec )   ! call iom_get with a 2D file
+         CALL fld_map_core( zz_read, kmap, pdta )
 
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       ! Do we include something here to adjust barotropic velocities !
       ! in case of a depth difference between bdy files and          !
-      ! bathymetry in the case ln_full_vel = .false. and jpk_bdy>0?  !
+      ! bathymetry in the case ln_totvel = .false. and ipkb>0?       !
       ! [as the enveloping and parital cells could change H]         !
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-      CASE DEFAULT   ;
-
-      IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN       ! boundary data not on model grid: vertical interpolation
-         CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar )
-         dta_read(:,:,:) = -ABS(fv)
-         dta_read_z(:,:,:) = 0._wp
-         dta_read_dz(:,:,:) = 0._wp
-         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec )
-         SELECT CASE( igrd )                  
-            CASE(1)
-               CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
-               CALL iom_get ( num, jpdom_unknown, 'e3t',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
-            CASE(2)  
-               CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
-               CALL iom_get ( num, jpdom_unknown, 'e3u',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
-            CASE(3)
-               CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
-               CALL iom_get ( num, jpdom_unknown, 'e3v',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
-         END SELECT
-
-      IF ( ln_bdy ) & 
-         CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta)
-
-      ELSE ! boundary data assumed to be on model grid
-         
-         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )                    
-         IF ( map%ll_unstruc) THEN ! unstructured open boundary data file
-            DO ib = 1, ipi
-              DO ik = 1, ipk
-                dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik)
-              END DO
+      ELSE
+         !
+         CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,:), krec )   ! call iom_get with a 3D file
+         !
+         IF( ipkb /= ipk ) THEN   ! boundary data not on model vertical grid : vertical interpolation
+            !
+            IF( ipk == jpk .AND. iom_varid(knum,'gdep'//clgrid(kgrd)) /= -1 .AND. iom_varid(knum,'e3'//clgrid(kgrd)) /= -1 ) THEN
+               
+               ALLOCATE( zdta_read(ipi,ipj,ipkb), zdta_read_z(ipi,ipj,ipkb), zdta_read_dz(ipi,ipj,ipkb) )
+                
+               CALL fld_map_core( zz_read, kmap, zdta_read )
+               CALL iom_get ( knum, jpdom_unknown, 'gdep'//clgrid(kgrd), zz_read )   ! read only once? Potential temporal evolution?
+               CALL fld_map_core( zz_read, kmap, zdta_read_z )
+               CALL iom_get ( knum, jpdom_unknown,   'e3'//clgrid(kgrd), zz_read )   ! read only once? Potential temporal evolution?
+               CALL fld_map_core( zz_read, kmap, zdta_read_dz )
+               
+               CALL iom_getatt(knum, '_FillValue', zfv, cdvar=cdvar )
+               CALL fld_bdy_interp(zdta_read, zdta_read_z, zdta_read_dz, pdta, kgrd, kbdy, zfv, ldtotvel)
+               DEALLOCATE( zdta_read, zdta_read_z, zdta_read_dz )
+               
+            ELSE
+               IF( ipk /= jpk ) CALL ctl_stop( 'fld_map : this should be an impossible case...' )
+               WRITE(ctmp1,*) 'fld_map : vertical interpolation for bdy variable '//TRIM(cdvar)//' requires ' 
+               IF( iom_varid(knum, 'gdep'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//'gdep'//clgrid(kgrd)//' variable' )
+               IF( iom_varid(knum,   'e3'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//  'e3'//clgrid(kgrd)//' variable' )
+
+            ENDIF
+            !
+         ELSE                            ! bdy data assumed to be the same levels as bdy variables
+            !
+            CALL fld_map_core( zz_read, kmap, pdta )
+            !
+         ENDIF   ! ipkb /= ipk
+      ENDIF   ! ipk == 1
+      
+      DEALLOCATE( zz_read )
+
+   END SUBROUTINE fld_map
+
+     
+   SUBROUTINE fld_map_core( pdta_read, kmap, pdta_bdy )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE fld_map_core  ***
+      !!
+      !! ** Purpose :  inner core of fld_map
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read    ! global boundary data
+      INTEGER,  DIMENSION(:    ), INTENT(in   ) ::   kmap         ! global-to-local bdy mapping indices
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta_bdy     ! bdy output field on model grid
+      !!
+      INTEGER,  DIMENSION(3) ::   idim_read,  idim_bdy            ! arrays dimensions
+      INTEGER                ::   ji, jj, jk, jb                  ! loop counters
+      INTEGER                ::   im1
+      !!---------------------------------------------------------------------
+      !
+      idim_read = SHAPE( pdta_read )
+      idim_bdy  = SHAPE( pdta_bdy  )
+      !
+      ! in all cases: idim_bdy(2) == 1 .AND. idim_read(1) * idim_read(2) == idim_bdy(1)
+      ! structured BDY with rimwidth > 1                     : idim_read(2) == rimwidth /= 1
+      ! structured BDY with rimwidth == 1 or unstructured BDY: idim_read(2) == 1
+      !
+      IF( idim_read(2) > 1 ) THEN    ! structured BDY with rimwidth > 1  
+         DO jk = 1, idim_bdy(3)
+            DO jb = 1, idim_bdy(1)
+               im1 = kmap(jb) - 1
+               jj = im1 / idim_read(1) + 1
+               ji = MOD( im1, idim_read(1) ) + 1
+               pdta_bdy(jb,1,jk) =  pdta_read(ji,jj,jk)
             END DO
-         ELSE ! we assume that this is a structured open boundary file
-            DO ib = 1, ipi
-               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
-               ji=map%ptr(ib)-(jj-1)*ilendta
-               DO ik = 1, ipk
-                  dta(ib,1,ik) =  dta_read(ji,jj,ik)
-               END DO
+         END DO
+      ELSE
+         DO jk = 1, idim_bdy(3)
+            DO jb = 1, idim_bdy(1)   ! horizontal remap of bdy data on the local bdy 
+               pdta_bdy(jb,1,jk) = pdta_read(kmap(jb),1,jk)
             END DO
-         ENDIF
-      ENDIF ! PRESENT(jpk_bdy)
-      END SELECT
-
-   END SUBROUTINE fld_map
+         END DO
+      ENDIF
+      
+   END SUBROUTINE fld_map_core
    
-   SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta)
-
+   
+   SUBROUTINE fld_bdy_interp(pdta_read, pdta_read_z, pdta_read_dz, pdta, kgrd, kbdy, pfv, ldtotvel)
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_bdy_interp  ***
@@ -847,284 +823,159 @@
       USE bdy_oce, ONLY:  idx_bdy         ! indexing for map <-> ij transformation
 
-      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read      ! work space for global data
-      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_z    ! work space for global data
-      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_dz   ! work space for global data
-      REAL(wp) , INTENT(in)                                ::   fv            ! fillvalue and alternative -ABS(fv)
-      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta                        ! output field on model grid (2 dimensional)
-      TYPE(MAP_POINTER)         , INTENT(in ) ::   map                        ! global-to-local mapping indices
-      LOGICAL  , INTENT(in), OPTIONAL         ::   fvl                        ! grid type, set number and number of vertical levels in the bdy data
-      INTEGER  , INTENT(in)                   ::   igrd, ibdy, jpk_bdy        ! number of levels in bdy data
-      INTEGER  , INTENT(in)                   ::   ilendta                    ! length of data in file
-      !!
-      INTEGER                                 ::   ipi                        ! length of boundary data on local process
-      INTEGER                                 ::   ipj                        ! length of dummy dimension ( = 1 )
-      INTEGER                                 ::   ipk                        ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
-      INTEGER                                 ::   jpkm1_bdy                  ! number of levels in bdy data minus 1
-      INTEGER                                 ::   ib, ik, ikk                ! loop counters
-      INTEGER                                 ::   ji, jj, zij, zjj           ! temporary indices
-      REAL(wp)                                ::   zl, zi, zh                 ! tmp variable for current depth and interpolation factor
-      REAL(wp)                                ::   fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv)
-      CHARACTER (LEN=10)                      ::   ibstr
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read       ! data read in bdy file
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_z     ! depth of the data read in bdy file
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_dz    ! thickness of the levels in bdy file
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta            ! output field on model grid (2 dimensional)
+      REAL(wp)                  , INTENT(in   ) ::   pfv             ! fillvalue of the data read in bdy file
+      LOGICAL                   , INTENT(in   ) ::   ldtotvel        ! true if toal ( = barotrop + barocline) velocity
+      INTEGER                   , INTENT(in   ) ::   kgrd            ! grid type (t, u, v)
+      INTEGER                   , INTENT(in   ) ::   kbdy            ! bdy number
+      !!
+      INTEGER                                   ::   ipi             ! length of boundary data on local process
+      INTEGER                                   ::   ipkb            ! number of vertical levels in boundary data file
+      INTEGER                                   ::   jb, ji, jj, jk, jkb   ! loop counters
+      REAL(wp)                                  ::   zcoef
+      REAL(wp)                                  ::   zl, zi, zh      ! tmp variable for current depth and interpolation factor
+      REAL(wp)                                  ::   zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv)
+      REAL(wp), DIMENSION(jpk)                  ::   zdepth, zdhalf  ! level and half-level depth
       !!---------------------------------------------------------------------
      
-
-      ipi       = SIZE( dta, 1 )
-      ipj       = SIZE( dta_read, 2 )
-      ipk       = SIZE( dta, 3 )
-      jpkm1_bdy = jpk_bdy-1
+      ipi  = SIZE( pdta, 1 )
+      ipkb = SIZE( pdta_read, 3 )
       
-      fv_alt = -ABS(fv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later
-      DO ib = 1, ipi
-            zij = idx_bdy(ibdy)%nbi(ib,igrd)
-            zjj = idx_bdy(ibdy)%nbj(ib,igrd)
-         IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj
-      ENDDO
-      !
-      IF ( map%ll_unstruc ) THEN                            ! unstructured open boundary data file
-
-         DO ib = 1, ipi
-            DO ik = 1, jpk_bdy
-               IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN
-                  dta_read_z(map%ptr(ib),1,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data
-                  dta_read_dz(map%ptr(ib),1,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct
+      zfv_alt = -ABS(pfv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later
+      !
+      WHERE( pdta_read == pfv )
+         pdta_read_z  = zfv_alt ! safety: put fillvalue into external depth field so consistent with data
+         pdta_read_dz = 0._wp   ! safety: put 0._wp into external thickness factors to ensure transport is correct
+      ENDWHERE
+      
+      DO jb = 1, ipi
+         ji = idx_bdy(kbdy)%nbi(jb,kgrd)
+         jj = idx_bdy(kbdy)%nbj(jb,kgrd)
+         zh  = SUM(pdta_read_dz(jb,1,:) )
+         !
+         ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary?
+         SELECT CASE( kgrd )                         
+         CASE(1)
+            IF( ABS( (zh - ht_n(ji,jj)) / ht_n(ji,jj)) * tmask(ji,jj,1) > 0.01_wp ) THEN
+               WRITE(ctmp1,"(I10.10)") jb 
+               CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%')
+               !   IF(lwp) WRITE(numout,*) 'DEPTHT', zh, sum(e3t_n(ji,jj,:), mask=tmask(ji,jj,:)==1),  ht_n(ji,jj), jb, jb, ji, jj
+            ENDIF
+         CASE(2)
+            IF( ABS( (zh - hu_n(ji,jj)) * r1_hu_n(ji,jj)) * umask(ji,jj,1) > 0.01_wp ) THEN
+               WRITE(ctmp1,"(I10.10)") jb 
+               CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%')
+               !   IF(lwp) WRITE(numout,*) 'DEPTHU', zh, SUM(e3u_n(ji,jj,:), mask=umask(ji,jj,:)==1),  SUM(umask(ji,jj,:)), &
+               !      &                hu_n(ji,jj), jb, jb, ji, jj, narea-1, pdta_read(jb,1,:)
+            ENDIF
+         CASE(3)
+            IF( ABS( (zh - hv_n(ji,jj)) * r1_hv_n(ji,jj)) * vmask(ji,jj,1) > 0.01_wp ) THEN
+               WRITE(ctmp1,"(I10.10)") jb
+               CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%')
+            ENDIF
+         END SELECT
+         !
+         SELECT CASE( kgrd )                         
+         CASE(1)
+            ! depth of T points:
+            zdepth(:) = gdept_n(ji,jj,:)
+         CASE(2)
+            ! depth of U points: we must not use gdept_n as we don't want to do a communication
+            !   --> copy what is done for gdept_n in domvvl...
+            zdhalf(1) = 0.0_wp
+            zdepth(1) = 0.5_wp * e3uw_n(ji,jj,1)
+            DO jk = 2, jpk                               ! vertical sum
+               !    zcoef = umask - wumask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
+               !                              ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
+               !                              ! 0.5 where jk = mikt     
+               !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
+               zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) )
+               zdhalf(jk) = zdhalf(jk-1) + e3u_n(ji,jj,jk-1)
+               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3uw_n(ji,jj,jk))  &
+                  &         + (1-zcoef) * ( zdepth(jk-1) +       e3uw_n(ji,jj,jk))
+            END DO
+         CASE(3)
+            ! depth of V points: we must not use gdept_n as we don't want to do a communication
+            !   --> copy what is done for gdept_n in domvvl...
+            zdhalf(1) = 0.0_wp
+            zdepth(1) = 0.5_wp * e3vw_n(ji,jj,1)
+            DO jk = 2, jpk                               ! vertical sum
+               !    zcoef = vmask - wvmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
+               !                              ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
+               !                              ! 0.5 where jk = mikt     
+               !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
+               zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) )
+               zdhalf(jk) = zdhalf(jk-1) + e3v_n(ji,jj,jk-1)
+               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3vw_n(ji,jj,jk))  &
+                  &         + (1-zcoef) * ( zdepth(jk-1) +       e3vw_n(ji,jj,jk))
+            END DO
+         END SELECT
+         !
+         DO jk = 1, jpk                      
+            IF(     zdepth(jk) < pdta_read_z(jb,1,          1) ) THEN                ! above the first level of external data
+               pdta(jb,1,jk) =  pdta_read(jb,1,1)
+            ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN                       ! below the last level of external data 
+               pdta(jb,1,jk) =  pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1))
+            ELSE                                                             ! inbetween: vertical interpolation between jkb & jkb+1
+               DO jkb = 1, ipkb-1                                            ! when  gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1)
+                  IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) &
+                     &    .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN   ! linear interpolation between 2 levels
+                     zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) )
+                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read  (jb,1,jkb+1) - pdta_read  (jb,1,jkb) ) * zi
+                  ENDIF
+               END DO
+            ENDIF
+         END DO   ! jpk
+         !
+      END DO   ! ipi
+      
+      IF(kgrd == 2) THEN                                  ! do we need to adjust the transport term?
+         DO jb = 1, ipi
+            ji = idx_bdy(kbdy)%nbi(jb,kgrd)
+            jj = idx_bdy(kbdy)%nbj(jb,kgrd)
+            zh  = SUM(pdta_read_dz(jb,1,:) )
+            ztrans = 0._wp
+            ztrans_new = 0._wp
+            DO jkb = 1, ipkb                              ! calculate transport on input grid
+               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
+            ENDDO
+            DO jk = 1, jpk                                ! calculate transport on model grid
+               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3u_n(ji,jj,jk ) * umask(ji,jj,jk)
+            ENDDO
+            DO jk = 1, jpk                                ! make transport correction
+               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
+                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk)
+               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
+                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hu_n(ji,jj)   * umask(ji,jj,jk)
                ENDIF
             ENDDO
-         ENDDO 
-
-         DO ib = 1, ipi
-            zij = idx_bdy(ibdy)%nbi(ib,igrd)
-            zjj = idx_bdy(ibdy)%nbj(ib,igrd)
-            zh  = SUM(dta_read_dz(map%ptr(ib),1,:) )
-            ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary?
-            SELECT CASE( igrd )                         
-               CASE(1)
-                  IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN
-                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
-                     CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
-                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj
-                  ENDIF
-               CASE(2)
-                  IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN
-                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
-                     CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
-                     IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1),  sum(umask(zij,zjj,:)), &
-                       &                hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1  , &
-                        &                dta_read(map%ptr(ib),1,:)
-                  ENDIF
-               CASE(3)
-                  IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN
-                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
-                     CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
-                  ENDIF
-            END SELECT
-            DO ik = 1, ipk                      
-               SELECT CASE( igrd )                       
-                  CASE(1)
-                     zl =  gdept_n(zij,zjj,ik)                                          ! if using in step could use fsdept instead of gdept_n?
-                  CASE(2)
-                     IF(ln_sco) THEN
-                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
-                     ELSE
-                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 
-                     ENDIF
-                  CASE(3)
-                     IF(ln_sco) THEN
-                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
-                     ELSE
-                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) )
-                     ENDIF
-               END SELECT
-               IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN                                         ! above the first level of external data
-                  dta(ib,1,ik) =  dta_read(map%ptr(ib),1,1)
-               ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN                           ! below the last level of external data 
-                  dta(ib,1,ik) =  dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1))
-               ELSE                                                                                ! inbetween : vertical interpolation between ikk & ikk+1
-                  DO ikk = 1, jpkm1_bdy                                                            ! when  gdept_n(ikk) < zl < gdept_n(ikk+1)
-                     IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) &
-                    &    .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN
-                        zi           = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / &
-                       &               ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) )
-                        dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + &
-                       &               ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi
-                     ENDIF
-                  END DO
-               ENDIF
-            END DO
-         END DO
-
-         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term?
-            DO ib = 1, ipi
-              zij = idx_bdy(ibdy)%nbi(ib,igrd)
-              zjj = idx_bdy(ibdy)%nbj(ib,igrd)
-              zh  = SUM(dta_read_dz(map%ptr(ib),1,:) )
-              ztrans = 0._wp
-              ztrans_new = 0._wp
-              DO ik = 1, jpk_bdy                            ! calculate transport on input grid
-                  ztrans     = ztrans     + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik)
-              ENDDO
-              DO ik = 1, ipk                                ! calculate transport on model grid
-                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik)
-              ENDDO
-              DO ik = 1, ipk                                ! make transport correction
-                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
-                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik)
-                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
-                    IF( ABS(ztrans * r1_hu_n(zij,zjj)) > 0.01_wp ) &
-                   &   CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at')
-                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) * umask(zij,zjj,ik)
-                 ENDIF
-              ENDDO
+         ENDDO
+      ENDIF
+      
+      IF(kgrd == 3) THEN                                  ! do we need to adjust the transport term?
+         DO jb = 1, ipi
+            ji = idx_bdy(kbdy)%nbi(jb,kgrd)
+            jj = idx_bdy(kbdy)%nbj(jb,kgrd)
+            zh  = SUM(pdta_read_dz(jb,1,:) )
+            ztrans = 0._wp
+            ztrans_new = 0._wp
+            DO jkb = 1, ipkb                              ! calculate transport on input grid
+               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
             ENDDO
-         ENDIF
-
-         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term?
-            DO ib = 1, ipi
-              zij = idx_bdy(ibdy)%nbi(ib,igrd)
-              zjj = idx_bdy(ibdy)%nbj(ib,igrd)
-              zh  = SUM(dta_read_dz(map%ptr(ib),1,:) )
-              ztrans = 0._wp
-              ztrans_new = 0._wp
-              DO ik = 1, jpk_bdy                            ! calculate transport on input grid
-                  ztrans     = ztrans     + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik)
-              ENDDO
-              DO ik = 1, ipk                                ! calculate transport on model grid
-                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik)
-              ENDDO
-              DO ik = 1, ipk                                ! make transport correction
-                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
-                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik)
-                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
-                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) * vmask(zij,zjj,ik)
-                 ENDIF
-              ENDDO
+            DO jk = 1, jpk                                ! calculate transport on model grid
+               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3v_n(ji,jj,jk ) * vmask(ji,jj,jk)
             ENDDO
-         ENDIF
-  
-      ELSE ! structured open boundary file
-
-         DO ib = 1, ipi
-            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
-            ji=map%ptr(ib)-(jj-1)*ilendta
-            DO ik = 1, jpk_bdy                      
-               IF( ( dta_read(ji,jj,ik) == fv ) ) THEN
-                  dta_read_z(ji,jj,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data
-                  dta_read_dz(ji,jj,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct
+            DO jk = 1, jpk                                ! make transport correction
+               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
+                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk)
+               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
+                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hv_n(ji,jj)   * vmask(ji,jj,jk)
                ENDIF
             ENDDO
-         ENDDO 
-       
-
-         DO ib = 1, ipi
-            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
-            ji=map%ptr(ib)-(jj-1)*ilendta
-            zij = idx_bdy(ibdy)%nbi(ib,igrd)
-            zjj = idx_bdy(ibdy)%nbj(ib,igrd)
-            zh  = SUM(dta_read_dz(ji,jj,:) )
-            ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary?
-            SELECT CASE( igrd )                         
-               CASE(1)
-                  IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN
-                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
-                     CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
-                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj
-                  ENDIF
-               CASE(2)
-                  IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN
-                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
-                     CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
-                  ENDIF
-               CASE(3)
-                  IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN
-                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
-                     CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
-                  ENDIF
-            END SELECT
-            DO ik = 1, ipk                      
-               SELECT CASE( igrd )                          ! coded for sco - need zco and zps option using min
-                  CASE(1)
-                     zl =  gdept_n(zij,zjj,ik)              ! if using in step could use fsdept instead of gdept_n?
-                  CASE(2)
-                     IF(ln_sco) THEN
-                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
-                     ELSE
-                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) )
-                     ENDIF
-                  CASE(3)
-                     IF(ln_sco) THEN
-                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
-                     ELSE
-                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) )
-                     ENDIF
-               END SELECT
-               IF( zl < dta_read_z(ji,jj,1) ) THEN                                      ! above the first level of external data
-                  dta(ib,1,ik) =  dta_read(ji,jj,1)
-               ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN                        ! below the last level of external data 
-                  dta(ib,1,ik) =  dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1))
-               ELSE                                                                     ! inbetween : vertical interpolation between ikk & ikk+1
-                  DO ikk = 1, jpkm1_bdy                                                 ! when  gdept_n(ikk) < zl < gdept_n(ikk+1)
-                     IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) &
-                    &    .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN
-                        zi           = ( zl - dta_read_z(ji,jj,ikk) ) / &
-                       &               ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) )
-                        dta(ib,1,ik) = dta_read(ji,jj,ikk) + &
-                       &               ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi
-                     ENDIF
-                  END DO
-               ENDIF
-            END DO
-         END DO
-
-         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term?
-            DO ib = 1, ipi
-               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
-               ji=map%ptr(ib)-(jj-1)*ilendta
-               zij = idx_bdy(ibdy)%nbi(ib,igrd)
-               zjj = idx_bdy(ibdy)%nbj(ib,igrd)
-               zh = SUM(dta_read_dz(ji,jj,:) )
-               ztrans = 0._wp
-               ztrans_new = 0._wp
-               DO ik = 1, jpk_bdy                            ! calculate transport on input grid
-                  ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik)
-               ENDDO
-               DO ik = 1, ipk                                ! calculate transport on model grid
-                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik)
-               ENDDO
-               DO ik = 1, ipk                                ! make transport correction
-                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
-                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik)
-                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
-                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term?
-            DO ib = 1, ipi
-               jj  = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
-               ji  = map%ptr(ib)-(jj-1)*ilendta
-               zij = idx_bdy(ibdy)%nbi(ib,igrd)
-               zjj = idx_bdy(ibdy)%nbj(ib,igrd)
-               zh  = SUM(dta_read_dz(ji,jj,:) )
-               ztrans = 0._wp
-               ztrans_new = 0._wp
-               DO ik = 1, jpk_bdy                            ! calculate transport on input grid
-                  ztrans     = ztrans     + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik)
-               ENDDO
-               DO ik = 1, ipk                                ! calculate transport on model grid
-                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik)
-               ENDDO
-               DO ik = 1, ipk                                ! make transport correction
-                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
-                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik)
-                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
-                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik)
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-      ENDIF ! endif unstructured or structured
-
+         ENDDO
+      ENDIF
+      
    END SUBROUTINE fld_bdy_interp
 
@@ -1151,4 +1002,5 @@
       imf = SIZE( sd )
       DO ju = 1, imf
+         IF( TRIM(sd(ju)%clrootname) == 'NOT USED' )   CYCLE
          ill = LEN_TRIM( sd(ju)%vcomp )
          DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2
@@ -1159,4 +1011,5 @@
                   iv = -1
                   DO jv = 1, imf
+                     IF( TRIM(sd(jv)%clrootname) == 'NOT USED' )   CYCLE
                      IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) )   iv = jv
                   END DO
@@ -1299,5 +1152,6 @@
       !
       DO jf = 1, SIZE(sdf)
-         sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
+         sdf(jf)%clrootname = sdf_n(jf)%clname
+         IF( TRIM(sdf_n(jf)%clname) /= 'NOT USED' )   sdf(jf)%clrootname = TRIM( cdir )//sdf(jf)%clrootname
          sdf(jf)%clname     = "not yet defined"
          sdf(jf)%nfreqh     = sdf_n(jf)%nfreqh
@@ -1308,7 +1162,7 @@
          sdf(jf)%num        = -1
          sdf(jf)%wgtname    = " "
-         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
+         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//sdf_n(jf)%wname
          sdf(jf)%lsmname = " "
-         IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//TRIM( sdf_n(jf)%lname )
+         IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//sdf_n(jf)%lname
          sdf(jf)%vcomp      = sdf_n(jf)%vcomp
          sdf(jf)%rotn(:)    = .TRUE.   ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
@@ -1317,5 +1171,9 @@
          IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim )   &
             &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
-         sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
+         sdf(jf)%nreclast   = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
+         sdf(jf)%igrd       = 0
+         sdf(jf)%ibdy       = 0
+         sdf(jf)%imap       => NULL()
+         sdf(jf)%ltotvel    = .FALSE.
       END DO
       !
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/updtide.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/updtide.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/updtide.F90	(revision 11223)
@@ -27,5 +27,5 @@
 CONTAINS
 
-   SUBROUTINE upd_tide( kt, kit, time_offset )
+   SUBROUTINE upd_tide( kt, kit, kt_offset )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE upd_tide  ***
@@ -39,9 +39,9 @@
       INTEGER, INTENT(in)           ::   kt      ! ocean time-step index
       INTEGER, INTENT(in), OPTIONAL ::   kit     ! external mode sub-time-step index (lk_dynspg_ts=T)
-      INTEGER, INTENT(in), OPTIONAL ::   time_offset ! time offset in number 
+      INTEGER, INTENT(in), OPTIONAL ::   kt_offset ! time offset in number 
                                                      ! of internal steps             (lk_dynspg_ts=F)
                                                      ! of external steps             (lk_dynspg_ts=T)
       !
-      INTEGER  ::   joffset      ! local integer
+      INTEGER  ::   ioffset      ! local integer
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zt, zramp    ! local scalar
@@ -52,11 +52,11 @@
       zt = ( kt - kt_tide ) * rdt
       !
-      joffset = 0
-      IF( PRESENT( time_offset ) )   joffset = time_offset
+      ioffset = 0
+      IF( PRESENT( kt_offset ) )   ioffset = kt_offset
       !
       IF( PRESENT( kit ) )   THEN
-         zt = zt + ( kit +  joffset - 1 ) * rdt / REAL( nn_baro, wp )
+         zt = zt + ( kit +  ioffset - 1 ) * rdt / REAL( nn_baro, wp )
       ELSE
-         zt = zt + joffset * rdt
+         zt = zt + ioffset * rdt
       ENDIF
       !
@@ -70,5 +70,5 @@
       IF( ln_tide_ramp ) THEN         ! linear increase if asked
          zt = ( kt - nit000 ) * rdt
-         IF( PRESENT( kit ) )   zt = zt + ( kit + joffset -1) * rdt / REAL( nn_baro, wp )
+         IF( PRESENT( kit ) )   zt = zt + ( kit + ioffset -1) * rdt / REAL( nn_baro, wp )
          zramp = MIN(  MAX( zt / (rdttideramp*rday) , 0._wp ) , 1._wp  )
          pot_astro(:,:) = zramp * pot_astro(:,:)
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/step.F90	(revision 11222)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/step.F90	(revision 11223)
@@ -112,5 +112,5 @@
       IF( ln_tide    )   CALL sbc_tide( kstp )                   ! update tide potential
       IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                   ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 
-      IF( ln_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
+      IF( ln_bdy     )   CALL bdy_dta ( kstp, kt_offset = +1 )   ! update dynamic & tracer data at open boundaries
                          CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice)
 
