Index: NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdy_oce.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdy_oce.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdy_oce.F90	(revision 13540)
@@ -63,4 +63,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::  aip    !: now ice  pond concentration
       REAL(wp), POINTER, DIMENSION(:,:) ::  hip    !: now ice  pond depth
+      REAL(wp), POINTER, DIMENSION(:,:) ::  hil    !: now ice  pond lid depth
 #if defined key_top
       CHARACTER(LEN=20)                   :: cn_obc  !: type of boundary condition to apply
@@ -115,4 +116,5 @@
    REAL(wp), DIMENSION(jp_bdy) ::   rice_apnd               !: pond conc.  of incoming sea ice
    REAL(wp), DIMENSION(jp_bdy) ::   rice_hpnd               !: pond thick. of incoming sea ice
+   REAL(wp), DIMENSION(jp_bdy) ::   rice_hlid               !: pond lid thick. of incoming sea ice
    !
    !!----------------------------------------------------------------------
Index: NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdydta.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdydta.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdydta.F90	(revision 13540)
@@ -43,5 +43,5 @@
    PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90
 
-   INTEGER , PARAMETER ::   jpbdyfld  = 16    ! maximum number of files to read 
+   INTEGER , PARAMETER ::   jpbdyfld  = 17    ! maximum number of files to read 
    INTEGER , PARAMETER ::   jp_bdyssh = 1     ! 
    INTEGER , PARAMETER ::   jp_bdyu2d = 2     ! 
@@ -60,4 +60,5 @@
    INTEGER , PARAMETER ::   jp_bdyaip = 15    ! 
    INTEGER , PARAMETER ::   jp_bdyhip = 16    ! 
+   INTEGER , PARAMETER ::   jp_bdyhil = 17    ! 
 #if ! defined key_si3
    INTEGER , PARAMETER ::   jpl = 1
@@ -70,4 +71,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -91,6 +93,4 @@
       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
@@ -108,12 +108,8 @@
          DO jbdy = 1, nb_bdy
             !
-            nblen    => idx_bdy(jbdy)%nblen
-            nblenrim => idx_bdy(jbdy)%nblenrim
-            !
             IF( nn_dyn2d_dta(jbdy) == 0 ) THEN 
-               ilen1(:) = nblen(:)
                IF( dta_bdy(jbdy)%lneed_ssh ) THEN 
                   igrd = 1
-                  DO ib = 1, ilen1(igrd)
+                  DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd)   ! ssh is allocated and used only on the rim
                      ii = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
@@ -121,13 +117,15 @@
                   END DO
                ENDIF
-               IF( dta_bdy(jbdy)%lneed_dyn2d) THEN 
+               IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain
                   igrd = 2
-                  DO ib = 1, ilen1(igrd)
+                  DO ib = 1, SIZE(dta_bdy(jbdy)%u2d)      ! u2d is used either over the whole bdy or only on the rim
                      ii = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
                      dta_bdy(jbdy)%u2d(ib) = uu_b(ii,ij,Kmm) * umask(ii,ij,1)         
                   END DO
+               ENDIF
+               IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain
                   igrd = 3
-                  DO ib = 1, ilen1(igrd)
+                  DO ib = 1, SIZE(dta_bdy(jbdy)%v2d)      ! v2d is used either over the whole bdy or only on the rim
                      ii = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
@@ -138,8 +136,7 @@
             !
             IF( nn_dyn3d_dta(jbdy) == 0 ) THEN 
-               ilen1(:) = nblen(:)
                IF( dta_bdy(jbdy)%lneed_dyn3d ) THEN 
                   igrd = 2 
-                  DO ib = 1, ilen1(igrd)
+                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
                      DO ik = 1, jpkm1
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
@@ -149,5 +146,5 @@
                   END DO
                   igrd = 3 
-                  DO ib = 1, ilen1(igrd)
+                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
                      DO ik = 1, jpkm1
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
@@ -160,8 +157,7 @@
 
             IF( nn_tra_dta(jbdy) == 0 ) THEN 
-               ilen1(:) = nblen(:)
                IF( dta_bdy(jbdy)%lneed_tra ) THEN
                   igrd = 1 
-                  DO ib = 1, ilen1(igrd)
+                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
                      DO ik = 1, jpkm1
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
@@ -176,9 +172,8 @@
 #if defined key_si3
             IF( nn_ice_dta(jbdy) == 0 ) THEN    ! set ice to initial values
-               ilen1(:) = nblen(:)
                IF( dta_bdy(jbdy)%lneed_ice ) THEN
                   igrd = 1   
                   DO jl = 1, jpl
-                     DO ib = 1, ilen1(igrd)
+                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
@@ -193,4 +188,5 @@
                         dta_bdy(jbdy)%aip(ib,jl) =  a_ip(ii,ij,jl) * tmask(ii,ij,1) 
                         dta_bdy(jbdy)%hip(ib,jl) =  h_ip(ii,ij,jl) * tmask(ii,ij,1) 
+                        dta_bdy(jbdy)%hil(ib,jl) =  h_il(ii,ij,jl) * tmask(ii,ij,1) 
                      END DO
                   END DO
@@ -218,25 +214,29 @@
          !
          ! 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
+         IF( cn_tra(jbdy) == 'runoff' ) THEN   ! runoff
             !
-            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
+            IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain
+               igrd = 2                         ! zonal flow (m3/s) to barotropic zonal velocity (m/s)
+               DO ib = 1, SIZE(dta_alias%u2d)   ! u2d is used either over the whole bdy or only on the rim
+                  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
+            ENDIF
+            IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain
+               igrd = 3                         ! meridional flow (m3/s) to barotropic meridional velocity (m/s)
+               DO ib = 1, SIZE(dta_alias%v2d)   ! v2d is used either over the whole bdy or only on the rim
+                  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
          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
+            IF( ASSOCIATED(dta_alias%ssh) ) dta_alias%ssh(:) = 0._wp
+            IF( ASSOCIATED(dta_alias%u2d) ) dta_alias%u2d(:) = 0._wp
+            IF( ASSOCIATED(dta_alias%v2d) ) dta_alias%v2d(:) = 0._wp
          ENDIF
 
@@ -245,10 +245,11 @@
             !
             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)
+               dta_alias%u2d(ib) = 0._wp   ! compute barotrope zonal velocity and put it in u2d
                DO ik = 1, jpkm1
-                  dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik)
+                  dta_alias%u2d(ib) = dta_alias%u2d(ib)   &
+                     &              + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik)
                END DO
                dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu(ii,ij,Kmm)
@@ -258,10 +259,11 @@
             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)
+               dta_alias%v2d(ib) = 0._wp   ! compute barotrope meridional velocity and put it in v2d
                DO ik = 1, jpkm1
-                  dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik)
+                  dta_alias%v2d(ib) = dta_alias%v2d(ib)   &
+                     &              + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik)
                END DO
                dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv(ii,ij,Kmm)
@@ -283,5 +285,5 @@
 
 #if defined key_si3
-         IF( dta_alias%lneed_ice ) THEN
+         IF( dta_alias%lneed_ice .AND. idx_bdy(jbdy)%nblen(1) > 0 ) THEN
             ! fill temperature and salinity arrays
             IF( TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyt_i)%fnow(:,1,:) = rice_tem (jbdy)
@@ -289,8 +291,9 @@
             IF( TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' )   bf_alias(jp_bdytsu)%fnow(:,1,:) = rice_tem (jbdy)
             IF( TRIM(bf_alias(jp_bdys_i)%clrootname) == 'NOT USED' )   bf_alias(jp_bdys_i)%fnow(:,1,:) = rice_sal (jbdy)
-            IF( TRIM(bf_alias(jp_bdyaip)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyaip)%fnow(:,1,:) = rice_apnd(jbdy) * & ! rice_apnd is the pond fraction
-               &                                                                         bf_alias(jp_bdya_i)%fnow(:,1,:)     !   ( a_ip = rice_apnd * a_i )
+            IF( TRIM(bf_alias(jp_bdyaip)%clrootname) == 'NOT USED' )   &               ! rice_apnd is the pond fraction
+               &   bf_alias(jp_bdyaip)%fnow(:,1,:) = rice_apnd(jbdy) * bf_alias(jp_bdya_i)%fnow(:,1,:)   ! ( a_ip = rice_apnd*a_i )
             IF( TRIM(bf_alias(jp_bdyhip)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyhip)%fnow(:,1,:) = rice_hpnd(jbdy)
-            
+            IF( TRIM(bf_alias(jp_bdyhil)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyhil)%fnow(:,1,:) = rice_hlid(jbdy)
+
             ! if T_i is read and not T_su, set T_su = T_i
             IF( TRIM(bf_alias(jp_bdyt_i)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' ) &
@@ -316,4 +319,8 @@
                bf_alias(jp_bdyaip)%fnow(:,1,:) = 0._wp
                bf_alias(jp_bdyhip)%fnow(:,1,:) = 0._wp
+               bf_alias(jp_bdyhil)%fnow(:,1,:) = 0._wp
+            ENDIF
+            IF ( .NOT.ln_pnd_lids ) THEN
+               bf_alias(jp_bdyhil)%fnow(:,1,:) = 0._wp
             ENDIF
             
@@ -321,12 +328,12 @@
             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                  , &
-                  &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), &
-                  &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), &
-                  &              bf_alias(jp_bdyaip)%fnow(:,1,:), bf_alias(jp_bdyhip)%fnow(:,1,:), &
-                  &              dta_alias%t_i                  , dta_alias%t_s                  , &
-                  &              dta_alias%tsu                  , dta_alias%s_i                  , &
-                  &              dta_alias%aip                  , dta_alias%hip )
+               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,:), & ! in
+                  &              dta_alias%h_i                  , dta_alias%h_s                  , dta_alias%a_i                  , & ! out
+                  &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), &                                  ! in (optional)
+                  &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), &                                  ! in     -
+                  &              bf_alias(jp_bdyaip)%fnow(:,1,:), bf_alias(jp_bdyhip)%fnow(:,1,:), bf_alias(jp_bdyhil)%fnow(:,1,:), & ! in     -
+                  &              dta_alias%t_i                  , dta_alias%t_s                  , &                                  ! out    -
+                  &              dta_alias%tsu                  , dta_alias%s_i                  , &                                  ! out    -
+                  &              dta_alias%aip                  , dta_alias%hip                  , dta_alias%hil )                    ! out    -
             ENDIF
          ENDIF
@@ -338,17 +345,11 @@
             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)%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))
+                  IF( ASSOCIATED(dta_bdy(jbdy)%ssh) ) dta_bdy_s(jbdy)%ssh(:) = dta_bdy(jbdy)%ssh(:)
+                  IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) dta_bdy_s(jbdy)%u2d(:) = dta_bdy(jbdy)%u2d(:)
+                  IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) dta_bdy_s(jbdy)%v2d(:) = dta_bdy(jbdy)%v2d(:)
                ENDIF
             END DO
          ELSE ! Add tides if not split-explicit free surface else this is done in ts loop
             !
-            ! BDY: use pt_offset=1.0 as applied at the end of the step and bdy_dta_tides is referenced at the middle of the step
             CALL bdy_dta_tides( kt=kt, pt_offset = 1._wp )
          ENDIF
@@ -358,5 +359,5 @@
       !
    END SUBROUTINE bdy_dta
-
+   
 
    SUBROUTINE bdy_dta_init
@@ -380,5 +381,5 @@
       !                                                         ! =F => baroclinic velocities in 3D boundary data
       LOGICAL                                ::   ln_zinterp    ! =T => requires a vertical interpolation of the bdydta
-      REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd 
+      REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd, rn_ice_hlid
       INTEGER                                ::   ipk,ipl       !
       INTEGER                                ::   idvar         ! variable ID
@@ -390,14 +391,15 @@
       LOGICAL                                ::   llneed        !
       LOGICAL                                ::   llread        !
+      LOGICAL                                ::   llfullbdy     !
       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, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip       
+      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip, bn_hil       
       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 
-      NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip
-      NAMELIST/nambdy_dta/ rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd
-      NAMELIST/nambdy_dta/ ln_full_vel, ln_zinterp
+      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d,                 &
+                         & bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip, bn_hil, &
+                         & rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd, rn_ice_hlid,      &
+                         & ln_full_vel, ln_zinterp
       !!---------------------------------------------------------------------------
       !
@@ -469,6 +471,9 @@
 #if defined key_si3
          IF( .NOT.ln_pnd ) THEN
-            rn_ice_apnd = 0. ; rn_ice_hpnd = 0.
-            CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 when no ponds' )
+            rn_ice_apnd = 0. ; rn_ice_hpnd = 0. ; rn_ice_hlid = 0.
+            CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 & rn_ice_hlid = 0 when no ponds' )
+         ENDIF
+         IF( .NOT.ln_pnd_lids ) THEN
+            rn_ice_hlid = 0.
          ENDIF
 #endif
@@ -480,5 +485,6 @@
          rice_apnd(jbdy) = rn_ice_apnd
          rice_hpnd(jbdy) = rn_ice_hpnd
-         
+         rice_hlid(jbdy) = rn_ice_hlid
+
          
          DO jfld = 1, jpbdyfld
@@ -504,10 +510,11 @@
                igrd = 2                                                    ! U point
                ipk = 1                                                     ! surface data
-               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
+               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%u2d 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
-               IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from u3d -> need on the full bdy
-               ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim
+               llfullbdy = ln_full_vel .OR. cn_dyn2d(jbdy) == 'frs'        ! need u2d over the whole bdy or only over the rim?
+               IF( llfullbdy ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)
+               ELSE                  ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)
                ENDIF
             ENDIF
@@ -516,10 +523,11 @@
                igrd = 3                                                    ! V point
                ipk = 1                                                     ! surface data
-               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
+               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%v2d 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 
-               IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from v3d -> need on the full bdy
-               ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim
+               llfullbdy = ln_full_vel .OR. cn_dyn2d(jbdy) == 'frs'        ! need v2d over the whole bdy or only over the rim?
+               IF( llfullbdy ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)
+               ELSE                  ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)
                ENDIF
             ENDIF
@@ -579,5 +587,5 @@
             IF(  jfld == jp_bdya_i .OR. jfld == jp_bdyh_i .OR. jfld == jp_bdyh_s .OR. &
                & jfld == jp_bdyt_i .OR. jfld == jp_bdyt_s .OR. jfld == jp_bdytsu .OR. &
-               & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip      ) THEN
+               & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip .OR. jfld == jp_bdyhil ) THEN
                igrd = 1                                                    ! T point
                ipk = ipl                                                   ! jpl-cat data
@@ -630,4 +638,9 @@
                bf_alias => bf(jp_bdyhip,jbdy:jbdy)                         ! alias for hip structure of bdy number jbdy
                bn_alias => bn_hip                                          ! alias for hip structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdyhil ) THEN
+               cl3 = 'hil'
+               bf_alias => bf(jp_bdyhil,jbdy:jbdy)                         ! alias for hil structure of bdy number jbdy
+               bn_alias => bn_hil                                          ! alias for hil structure of nambdy_dta 
             ENDIF
 
@@ -699,4 +712,9 @@
                   ENDIF
                ENDIF
+               IF( jfld == jp_bdyhil ) THEN
+                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%hil => bf_alias(1)%fnow(:,1,:)
+                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%hil(iszdim,jpl) )
+                  ENDIF
+               ENDIF
             ENDIF
 
Index: NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdydyn.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdydyn.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdydyn.F90	(revision 13540)
@@ -30,4 +30,6 @@
    PUBLIC   bdy_dyn    ! routine called in dyn_nxt
 
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Index: NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdydyn2d.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdydyn2d.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdydyn2d.F90	(revision 13540)
@@ -102,8 +102,8 @@
          END DO
          IF( ANY(llsend2) .OR. ANY(llrecv2) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdydyn2d', pua2d, 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
+            CALL lbc_lnk( 'bdydyn2d', pua2d, 'U', -1.0_wp, kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
          END IF
          IF( ANY(llsend3) .OR. ANY(llrecv3) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdydyn2d', pva2d, 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
+            CALL lbc_lnk( 'bdydyn2d', pva2d, 'V', -1.0_wp, kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
          END IF
          !
@@ -324,5 +324,5 @@
          IF( nn_hls > 1 .AND. ir == 1 ) CYCLE   ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
          IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdydyn2d', zssh(:,:,1), 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
+            CALL lbc_lnk( 'bdydyn2d', zssh(:,:,1), 'T',  1.0_wp, kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
          END IF
       END DO
Index: NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdydyn3d.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdydyn3d.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdydyn3d.F90	(revision 13540)
@@ -99,8 +99,8 @@
          !
          IF( ANY(llsend2) .OR. ANY(llrecv2) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdydyn2d', puu(:,:,:,Kaa), 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
+            CALL lbc_lnk( 'bdydyn2d', puu(:,:,:,Kaa), 'U', -1.0_wp, kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
          END IF
          IF( ANY(llsend3) .OR. ANY(llrecv3) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdydyn2d', pvv(:,:,:,Kaa), 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
+            CALL lbc_lnk( 'bdydyn2d', pvv(:,:,:,Kaa), 'V', -1.0_wp, kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
          END IF
       END DO   ! ir
Index: NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdyice.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdyice.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdyice.F90	(revision 13540)
@@ -94,11 +94,11 @@
          IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
             ! exchange 3d arrays
-            CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1. &
-                 &                      , a_ip, 'T', 1., v_ip, 'T', 1., s_i , 'T', 1., t_su, 'T', 1. &
-                 &                      , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1.                &
-                 &                      , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1      )
+            CALL lbc_lnk_multi('bdyice', a_i , 'T', 1._wp, h_i , 'T', 1._wp, h_s , 'T', 1._wp, oa_i, 'T', 1._wp                   &
+               &                       , s_i , 'T', 1._wp, t_su, 'T', 1._wp, v_i , 'T', 1._wp, v_s , 'T', 1._wp, sv_i, 'T', 1._wp &
+               &                       , a_ip, 'T', 1._wp, v_ip, 'T', 1._wp, v_il, 'T', 1._wp                                     &
+               &                       , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
             ! exchange 4d arrays :   third dimension = 1   and then   third dimension = jpk
-            CALL lbc_lnk_multi( 'bdyice', t_s , 'T', 1., e_s , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
-            CALL lbc_lnk_multi( 'bdyice', t_i , 'T', 1., e_i , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
+            CALL lbc_lnk_multi('bdyice', t_s , 'T', 1._wp, e_s , 'T', 1._wp, kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
+            CALL lbc_lnk_multi('bdyice', t_i , 'T', 1._wp, e_i , 'T', 1._wp, kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
          END IF
       END DO   ! ir
@@ -163,4 +163,5 @@
             a_ip(ji,jj,  jl) = ( a_ip(ji,jj,  jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond concentration
             h_ip(ji,jj,  jl) = ( h_ip(ji,jj,  jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond depth
+            h_il(ji,jj,  jl) = ( h_il(ji,jj,  jl) * zwgt1 + dta%hil(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond lid depth
             !
             sz_i(ji,jj,:,jl) = s_i(ji,jj,jl)
@@ -170,4 +171,9 @@
                a_ip(ji,jj,jl) = 0._wp
                h_ip(ji,jj,jl) = 0._wp
+               h_il(ji,jj,jl) = 0._wp
+            ENDIF
+
+            IF( .NOT.ln_pnd_lids ) THEN
+               h_il(ji,jj,jl) = 0._wp
             ENDIF
             !
@@ -231,4 +237,5 @@
                a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl)
                h_ip(ji,jj,  jl) = h_ip(ib,jb,  jl)
+               h_il(ji,jj,  jl) = h_il(ib,jb,  jl)
                !
                sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl)
@@ -265,10 +272,6 @@
                !
                ! melt ponds
-               IF( a_i(ji,jj,jl) > epsi10 ) THEN
-                  a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i (ji,jj,jl)
-               ELSE
-                  a_ip_frac(ji,jj,jl) = 0._wp
-               ENDIF
                v_ip(ji,jj,jl) = h_ip(ji,jj,jl) * a_ip(ji,jj,jl)
+               v_il(ji,jj,jl) = h_il(ji,jj,jl) * a_ip(ji,jj,jl)
                !
             ELSE   ! no ice at the boundary
@@ -278,14 +281,11 @@
                h_s (ji,jj,  jl) = 0._wp
                oa_i(ji,jj,  jl) = 0._wp
-               a_ip(ji,jj,  jl) = 0._wp
-               v_ip(ji,jj,  jl) = 0._wp
                t_su(ji,jj,  jl) = rt0
                t_s (ji,jj,:,jl) = rt0
                t_i (ji,jj,:,jl) = rt0 
 
-               a_ip_frac(ji,jj,jl) = 0._wp
-               h_ip     (ji,jj,jl) = 0._wp
-               a_ip     (ji,jj,jl) = 0._wp
-               v_ip     (ji,jj,jl) = 0._wp
+               a_ip(ji,jj,jl) = 0._wp
+               h_ip(ji,jj,jl) = 0._wp
+               h_il(ji,jj,jl) = 0._wp
                
                IF( nn_icesal == 1 ) THEN     ! if constant salinity
@@ -303,4 +303,6 @@
                e_s (ji,jj,:,jl) = 0._wp
                e_i (ji,jj,:,jl) = 0._wp
+               v_ip(ji,jj,  jl) = 0._wp
+               v_il(ji,jj,  jl) = 0._wp
 
             ENDIF
@@ -436,5 +438,5 @@
             END DO
             IF( ANY(llsend2) .OR. ANY(llrecv2) ) THEN   ! if need to send/recv in at least one direction
-               CALL lbc_lnk( 'bdyice', u_ice, 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
+               CALL lbc_lnk( 'bdyice', u_ice, 'U', -1.0_wp, kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
             END IF
          CASE ( 'V' )
@@ -450,5 +452,5 @@
             END DO
             IF( ANY(llsend3) .OR. ANY(llrecv3) ) THEN   ! if need to send/recv in at least one direction
-               CALL lbc_lnk( 'bdyice', v_ice, 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
+               CALL lbc_lnk( 'bdyice', v_ice, 'V', -1.0_wp, kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
             END IF
          END SELECT
Index: NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdyini.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdyini.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdyini.F90	(revision 13540)
@@ -19,9 +19,10 @@
    USE oce            ! ocean dynamics and tracers variables
    USE dom_oce        ! ocean space and time domain
+   USE sbc_oce , ONLY: nn_ice
    USE bdy_oce        ! unstructured open boundary conditions
    USE bdydta         ! open boundary cond. setting   (bdy_dta_init routine)
    USE bdytides       ! open boundary cond. setting   (bdytide_init routine)
    USE tide_mod, ONLY: ln_tide ! tidal forcing
-   USE phycst   , ONLY: rday
+   USE phycst  , ONLY: rday
    !
    USE in_out_manager ! I/O units
@@ -316,4 +317,9 @@
          dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none'
 
+         IF( dta_bdy(ib_bdy)%lneed_ice .AND. nn_ice /= 2 ) THEN
+            WRITE(ctmp1,*) 'bdy number ', ib_bdy,', needs ice model but nn_ice = ', nn_ice
+            CALL ctl_stop( ctmp1 )
+         ENDIF
+
          IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN 
             SELECT CASE( nn_ice_dta(ib_bdy) )                   ! 
@@ -410,9 +416,9 @@
                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) = NINT( zz_read(ii,1) )
+                  nbidta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) + nn_hls
                END DO
                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) = NINT( zz_read(ii,1) )
+                  nbjdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) + nn_hls
                END DO
                CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
@@ -632,5 +638,5 @@
          END DO
       END DO
-      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. )
+      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1.0_wp )
 
       ! Read global 2D mask at T-points: bdytmask
@@ -648,5 +654,5 @@
          END DO
       END DO
-      CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1., bdyvmask, 'V', 1. )   ! Lateral boundary cond.
+      CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1.0_wp , bdyvmask, 'V', 1.0_wp )   ! Lateral boundary cond. 
 
       ! bdy masks are now set to zero on rim 0 points:
@@ -689,5 +695,5 @@
          END DO
       END DO
-      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. )
+      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1.0_wp )
 
       ! bdy masks are now set to zero on rim1 points:
@@ -865,7 +871,7 @@
             ENDIF 
             SELECT CASE( igrd )
-               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
-               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 
-               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 
+               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
+               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
+               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
             END SELECT 
             DO ib = ibeg, iend
@@ -913,7 +919,7 @@
             ENDIF
             SELECT CASE( igrd )
-               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
-               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 
-               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 
+               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
+               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
+               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
             END SELECT 
             DO ib = ibeg, iend
@@ -1001,7 +1007,7 @@
             END DO
             SELECT CASE( igrd )
-               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
-               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 
-               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 
+               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
+               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
+               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
             END SELECT 
             DO ib = ibeg, iend
@@ -1105,13 +1111,13 @@
       CASE( 'N' )
          IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
-            nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain.
+            nbdyind  = Nj0glo - 2  ! set boundary to whole side of model domain.
             nbdybeg  = 2
-            nbdyend  = jpiglo - 1
+            nbdyend  = Ni0glo -1
          ENDIF
          nbdysegn = nbdysegn + 1
          npckgn(nbdysegn) = kb_bdy ! Save bdy package number
-         jpjnob(nbdysegn) = nbdyind
-         jpindt(nbdysegn) = nbdybeg
-         jpinft(nbdysegn) = nbdyend
+         jpjnob(nbdysegn) = nbdyind + nn_hls
+         jpindt(nbdysegn) = nbdybeg + nn_hls
+         jpinft(nbdysegn) = nbdyend + nn_hls
          !
       CASE( 'S' )
@@ -1119,23 +1125,23 @@
             nbdyind  = 2           ! set boundary to whole side of model domain.
             nbdybeg  = 2
-            nbdyend  = jpiglo - 1
+            nbdyend  = Ni0glo - 1
          ENDIF
          nbdysegs = nbdysegs + 1
          npckgs(nbdysegs) = kb_bdy ! Save bdy package number
-         jpjsob(nbdysegs) = nbdyind
-         jpisdt(nbdysegs) = nbdybeg
-         jpisft(nbdysegs) = nbdyend
+         jpjsob(nbdysegs) = nbdyind + nn_hls
+         jpisdt(nbdysegs) = nbdybeg + nn_hls
+         jpisft(nbdysegs) = nbdyend + nn_hls
          !
       CASE( 'E' )
          IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
-            nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain.
+            nbdyind  = Ni0glo - 2  ! set boundary to whole side of model domain.
             nbdybeg  = 2
-            nbdyend  = jpjglo - 1
+            nbdyend  = Nj0glo - 1
          ENDIF
          nbdysege = nbdysege + 1 
          npckge(nbdysege) = kb_bdy ! Save bdy package number
-         jpieob(nbdysege) = nbdyind
-         jpjedt(nbdysege) = nbdybeg
-         jpjeft(nbdysege) = nbdyend
+         jpieob(nbdysege) = nbdyind + nn_hls
+         jpjedt(nbdysege) = nbdybeg + nn_hls
+         jpjeft(nbdysege) = nbdyend + nn_hls
          !
       CASE( 'W' )
@@ -1143,11 +1149,11 @@
             nbdyind  = 2           ! set boundary to whole side of model domain.
             nbdybeg  = 2
-            nbdyend  = jpjglo - 1
+            nbdyend  = Nj0glo - 1
          ENDIF
          nbdysegw = nbdysegw + 1
          npckgw(nbdysegw) = kb_bdy ! Save bdy package number
-         jpiwob(nbdysegw) = nbdyind
-         jpjwdt(nbdysegw) = nbdybeg
-         jpjwft(nbdysegw) = nbdyend
+         jpiwob(nbdysegw) = nbdyind + nn_hls
+         jpjwdt(nbdysegw) = nbdybeg + nn_hls
+         jpjwft(nbdysegw) = nbdyend + nn_hls
          !
       CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
@@ -1186,4 +1192,5 @@
       IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn
       IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs
+      !
       ! 1. Check bounds
       !----------------
@@ -1223,5 +1230,4 @@
          IF (jpjwft(ib).gt.jpjglo)     CALL ctl_stop( 'End index out of domain' )
       ENDDO
-      !
       !      
       ! 2. Look for segment crossings
@@ -1372,8 +1378,6 @@
          DO ji = 1, jpi
             DO jj = 1, jpj             
-              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
-               &  ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
-              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
-               &  ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1)  
+              IF( mig(ji) == jpiwob(ib) .AND. mjg(jj) == jpjwdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
+              IF( mig(ji) == jpiwob(ib) .AND. mjg(jj) == jpjwft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
             END DO
          END DO
@@ -1410,8 +1414,6 @@
          DO ji = 1, jpi
             DO jj = 1, jpj             
-              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
-               &  ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
-              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
-               &  ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1)  
+              IF( mig(ji) == jpieob(ib)+1 .AND. mjg(jj) == jpjedt(ib) )   ztestmask(1) = tmask(ji,jj,1)
+              IF( mig(ji) == jpieob(ib)+1 .AND. mjg(jj) == jpjeft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
             END DO
          END DO
@@ -1448,8 +1450,6 @@
          DO ji = 1, jpi
             DO jj = 1, jpj             
-              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
-               &  ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
-              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
-               &  ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1)  
+              IF( mjg(jj) == jpjsob(ib) .AND. mig(ji) == jpisdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
+              IF( mjg(jj) == jpjsob(ib) .AND. mig(ji) == jpisft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
             END DO
          END DO
@@ -1472,8 +1472,6 @@
          DO ji = 1, jpi
             DO jj = 1, jpj             
-              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
-               &  ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
-              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
-               &  ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1)  
+               IF( mjg(jj) == jpjnob(ib)+1 .AND. mig(ji) == jpindt(ib) )   ztestmask(1) = tmask(ji,jj,1)
+               IF( mjg(jj) == jpjnob(ib)+1 .AND. mig(ji) == jpinft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
             END DO
          END DO
Index: NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdylib.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdylib.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdylib.F90	(revision 13540)
@@ -44,5 +44,5 @@
       !!----------------------------------------------------------------------
       TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
-      REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
+      REAL(wp), DIMENSION(:,:), POINTER,   INTENT(in) ::   dta  ! OBC external data
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend
       !!
@@ -73,5 +73,5 @@
       !!----------------------------------------------------------------------
       TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
-      REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
+      REAL(wp), DIMENSION(:,:), POINTER,   INTENT(in) ::   dta  ! OBC external data
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend
       !!
@@ -100,10 +100,10 @@
       !!
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
-      REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phib  ! before tracer field
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend
-      LOGICAL                 , OPTIONAL,  INTENT(in) ::   lrim0   ! indicate if rim 0 is treated
-      LOGICAL,                             INTENT(in) ::   ll_npo  ! switch for NPO version
+      TYPE(OBC_INDEX),                   INTENT(in   ) ::   idx  ! OBC indices
+      REAL(wp), DIMENSION(:,:), POINTER, INTENT(in   ) ::   dta  ! OBC external data
+      REAL(wp), DIMENSION(jpi,jpj,jpk),  INTENT(inout) ::   phib  ! before tracer field
+      REAL(wp), DIMENSION(jpi,jpj,jpk),  INTENT(inout) ::   phia  ! tracer trend
+      LOGICAL ,                          INTENT(in   ) ::   lrim0   ! indicate if rim 0 is treated
+      LOGICAL ,                          INTENT(in   ) ::   ll_npo  ! switch for NPO version
       !!
       INTEGER  ::   igrd                                    ! grid index
@@ -128,11 +128,11 @@
       !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)    
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX),          INTENT(in   ) ::   idx      ! BDY indices
-      INTEGER ,                 INTENT(in   ) ::   igrd     ! grid index
-      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   phib     ! model before 2D field
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   phia     ! model after 2D field (to be updated)
-      REAL(wp), DIMENSION(:)  , INTENT(in   ) ::   phi_ext  ! external forcing data
-      LOGICAL, OPTIONAL,        INTENT(in   ) ::   lrim0    ! indicate if rim 0 is treated
-      LOGICAL ,                 INTENT(in   ) ::   ll_npo   ! switch for NPO version
+      TYPE(OBC_INDEX),                   INTENT(in   ) ::   idx      ! BDY indices
+      INTEGER ,                          INTENT(in   ) ::   igrd     ! grid index
+      REAL(wp), DIMENSION(:,:),          INTENT(in   ) ::   phib     ! model before 2D field
+      REAL(wp), DIMENSION(:,:),          INTENT(inout) ::   phia     ! model after 2D field (to be updated)
+      REAL(wp), DIMENSION(:  ), POINTER, INTENT(in   ) ::   phi_ext  ! external forcing data
+      LOGICAL ,                          INTENT(in   ) ::   lrim0    ! indicate if rim 0 is treated
+      LOGICAL ,                          INTENT(in   ) ::   ll_npo   ! switch for NPO version
       !
       INTEGER  ::   jb                                     ! dummy loop indices
@@ -188,10 +188,7 @@
       END SELECT
       !
-      IF( PRESENT(lrim0) ) THEN
-         IF( lrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)   ! rim 0
-         ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)    ! rim 1
-         END IF
-      ELSE                  ;   ibeg = 1                       ;   iend = idx%nblenrim(igrd)    ! both
-      END IF
+      IF( lrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)   ! rim 0
+      ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)    ! rim 1
+      ENDIF
       !
       DO jb = ibeg, iend
@@ -249,5 +246,5 @@
 !!$         zdy_centred = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1)
          ! upstream differencing for tangential derivatives
-         zsign_ups = sign( 1., zdt * zdy_centred )
+         zsign_ups = sign( 1.0_wp, zdt * zdy_centred )
          zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
          zdy = zsign_ups * zdy_1 + (1. - zsign_ups) * zdy_2
@@ -257,5 +254,5 @@
          zrx = zdt * zdx / ( zex1 * znor2 ) 
 !!$         zrx = min(zrx,2.0_wp)
-         zout = sign( 1., zrx )
+         zout = sign( 1.0_wp, zrx )
          zout = 0.5*( zout + abs(zout) )
          zwgt = 2.*rn_Dt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) )
@@ -266,5 +263,5 @@
            &                            + zwgt * ( phi_ext(jb) - phib(ii,ij) ) ) / ( 1. + zrx ) 
          else                  !! full oblique radiation !!
-            zsign_ups = sign( 1., zdt * zdy )
+            zsign_ups = sign( 1.0_wp, zdt * zdy )
             zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
             zey = zsign_ups * zey1 + (1.-zsign_ups) * zey2 
@@ -275,5 +272,5 @@
            &                    - (1.-zsign_ups) * zry * ( phib(iijp1,ijjp1) - phib(ii   ,ij    ) ) &
            &                    + zwgt * ( phi_ext(jb) - phib(ii,ij) ) ) / ( 1. + zrx ) 
-         end if
+         endif
          phia(ii,ij) = phia(ii,ij) * zmask(ii,ij)
       END DO
@@ -293,11 +290,11 @@
       !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)    
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX),            INTENT(in   ) ::   idx      ! BDY indices
-      INTEGER ,                   INTENT(in   ) ::   igrd     ! grid index
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phib     ! model before 3D field
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phia     ! model after 3D field (to be updated)
-      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   phi_ext  ! external forcing data
-      LOGICAL, OPTIONAL,          INTENT(in   ) ::   lrim0    ! indicate if rim 0 is treated
-      LOGICAL ,                   INTENT(in   ) ::   ll_npo   ! switch for NPO version
+      TYPE(OBC_INDEX),                     INTENT(in   ) ::   idx      ! BDY indices
+      INTEGER ,                            INTENT(in   ) ::   igrd     ! grid index
+      REAL(wp), DIMENSION(:,:,:),          INTENT(in   ) ::   phib     ! model before 3D field
+      REAL(wp), DIMENSION(:,:,:),          INTENT(inout) ::   phia     ! model after 3D field (to be updated)
+      REAL(wp), DIMENSION(:,:  ), POINTER, INTENT(in   ) ::   phi_ext  ! external forcing data
+      LOGICAL ,                            INTENT(in   ) ::   lrim0    ! indicate if rim 0 is treated
+      LOGICAL ,                            INTENT(in   ) ::   ll_npo   ! switch for NPO version
       !
       INTEGER  ::   jb, jk                                 ! dummy loop indices
@@ -353,10 +350,7 @@
       END SELECT
       !
-      IF( PRESENT(lrim0) ) THEN
-         IF( lrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)   ! rim 0
-         ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)    ! rim 1
-         END IF
-      ELSE                  ;   ibeg = 1                       ;   iend = idx%nblenrim(igrd)    ! both
-      END IF
+      IF( lrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)   ! rim 0
+      ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)    ! rim 1
+      ENDIF
       !
       DO jk = 1, jpk
@@ -414,5 +408,5 @@
 !!$            zdy_centred = phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1jm1,ijbm1jm1,jk)
             ! upstream differencing for tangential derivatives
-            zsign_ups = sign( 1., zdt * zdy_centred )
+            zsign_ups = sign( 1.0_wp, zdt * zdy_centred )
             zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
             zdy = zsign_ups * zdy_1 + (1. - zsign_ups) * zdy_2
@@ -423,5 +417,5 @@
             zrx = zdt * zdx / ( zex1 * znor2 )
 !!$            zrx = min(zrx,2.0_wp)
-            zout = sign( 1., zrx )
+            zout = sign( 1.0_wp, zrx )
             zout = 0.5*( zout + abs(zout) )
             zwgt = 2.*rn_Dt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) )
@@ -432,5 +426,5 @@
               &                            + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) ) ) / ( 1. + zrx ) 
             else                  !! full oblique radiation !!
-               zsign_ups = sign( 1., zdt * zdy )
+               zsign_ups = sign( 1.0_wp, zdt * zdy )
                zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
                zey = zsign_ups * zey1 + (1.-zsign_ups) * zey2 
@@ -441,5 +435,5 @@
               &                       - (1.-zsign_ups) * zry * ( phib(iijp1,ijjp1,jk) - phib(ii   ,ij   ,jk) ) &
               &                       + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) ) ) / ( 1. + zrx ) 
-            end if
+            endif
             phia(ii,ij,jk) = phia(ii,ij,jk) * zmask(ii,ij,jk)
          END DO
@@ -466,5 +460,5 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated), must be masked
       TYPE(OBC_INDEX),            INTENT(in   )  ::   idx      ! OBC indices
-      LOGICAL, OPTIONAL,          INTENT(in   )  ::   lrim0    ! indicate if rim 0 is treated
+      LOGICAL ,                   INTENT(in   )  ::   lrim0    ! indicate if rim 0 is treated
       !! 
       REAL(wp) ::   zweight
@@ -486,10 +480,7 @@
       END SELECT
       !
-      IF( PRESENT(lrim0) ) THEN
-         IF( lrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)   ! rim 0
-         ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)    ! rim 1
-         END IF
-      ELSE                  ;   ibeg = 1                       ;   iend = idx%nblenrim(igrd)    ! both
-      END IF
+      IF( lrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)   ! rim 0
+      ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)    ! rim 1
+      ENDIF
       !
       DO ib = ibeg, iend
Index: NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdytides.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdytides.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdytides.F90	(revision 13540)
@@ -65,19 +65,20 @@
       !! namelist variables
       !!-------------------
-      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
-      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data
+      CHARACTER(len=80)                         ::   filtide             ! Filename root for tidal input files
+      LOGICAL                                   ::   ln_bdytide_2ddta    ! If true, read 2d harmonic data
       !!
-      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices
-      INTEGER                                   ::   ii, ij              !: dummy loop indices
+      INTEGER                                   ::   ib_bdy, itide, ib   ! dummy loop indices
+      INTEGER                                   ::   ii, ij              ! dummy loop indices
       INTEGER                                   ::   inum, igrd
-      INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays)
+      INTEGER                                   ::   isz                 ! bdy data size
       INTEGER                                   ::   ios                 ! Local integer output status for namelist read
       INTEGER                                   ::   nbdy_rdstart, nbdy_loc
-      CHARACTER(LEN=50)                         ::   cerrmsg             !: error string
-      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file 
-      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data
-      REAL(wp),ALLOCATABLE, DIMENSION(:,:)      ::   ztr, zti            !:  "     "    "   "   "   "        "      " 
+      CHARACTER(LEN=50)                         ::   cerrmsg             ! error string
+      CHARACTER(len=80)                         ::   clfile              ! full file name for tidal input file 
+      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            ! work space to read in tidal harmonics data
+      REAL(wp),ALLOCATABLE, DIMENSION(:,:)      ::   ztr, zti            !  "     "    "   "   "   "        "      " 
       !!
-      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
+      TYPE(TIDES_DATA), POINTER                 ::   td                  ! local short cut   
+      TYPE(  OBC_DATA), POINTER                 ::   dta                 ! local short cut
       !!
       NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta
@@ -93,6 +94,7 @@
          IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
             !
-            td => tides(ib_bdy)
-
+            td  => tides(ib_bdy)
+            dta => dta_bdy(ib_bdy)
+         
             ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
             filtide(:) = ''
@@ -130,30 +132,28 @@
             IF(lwp) WRITE(numout,*) ' '
 
-            ! Allocate space for tidal harmonics data - get size from OBC data arrays
+            ! Allocate space for tidal harmonics data - get size from BDY data arrays
+            ! Allocate also slow varying data in the case of time splitting:
+            ! Do it anyway because at this stage knowledge of free surface scheme is unknown
             ! -----------------------------------------------------------------------
-
-            ! 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(:) = idx_bdy(ib_bdy)%nblen   (:)
-            ELSE                                   ;   ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:)
-            ENDIF
-
-            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
-            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
-
-            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
-            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
-
-            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
-            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
-
-            td%ssh0(:,:,:) = 0._wp
-            td%ssh (:,:,:) = 0._wp
-            td%u0  (:,:,:) = 0._wp
-            td%u   (:,:,:) = 0._wp
-            td%v0  (:,:,:) = 0._wp
-            td%v   (:,:,:) = 0._wp
-
+            IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
+               isz = SIZE(dta%ssh)
+               ALLOCATE( td%ssh0( isz, nb_harmo, 2 ), td%ssh( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%ssh( isz ) )
+               dta_bdy_s(ib_bdy)%ssh(:) = 0._wp   ! needed?
+            ENDIF
+            IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
+               isz = SIZE(dta%u2d)
+               ALLOCATE( td%u0  ( isz, nb_harmo, 2 ), td%u  ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%u2d( isz ) )
+               dta_bdy_s(ib_bdy)%u2d(:) = 0._wp   ! needed?
+            ENDIF
+            IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
+               isz = SIZE(dta%v2d)
+               ALLOCATE( td%v0  ( isz, nb_harmo, 2 ), td%v  ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%v2d( isz ) )
+               dta_bdy_s(ib_bdy)%v2d(:) = 0._wp   ! needed?
+            ENDIF
+
+            ! fill td%ssh0, td%u0, td%v0
+            ! -----------------------------------------------------------------------
             IF( ln_bdytide_2ddta ) THEN
+               !
                ! It is assumed that each data file contains all complex harmonic amplitudes
                ! given on the global domain (ie global, jpiglo x jpjglo)
@@ -162,53 +162,56 @@
                !
                ! SSH fields
-               clfile = TRIM(filtide)//'_grid_T.nc'
-               CALL iom_open( clfile , inum ) 
-               igrd = 1                       ! Everything is at T-points here
-               DO itide = 1, nb_harmo
-                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) )
-                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) ) 
-                  DO ib = 1, ilen0(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
-                     td%ssh0(ib,itide,1) = ztr(ii,ij)
-                     td%ssh0(ib,itide,2) = zti(ii,ij)
-                  END DO
-               END DO 
-               CALL iom_close( inum )
+               IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
+                  clfile = TRIM(filtide)//'_grid_T.nc'
+                  CALL iom_open( clfile , inum ) 
+                  igrd = 1                       ! Everything is at T-points here
+                  DO itide = 1, nb_harmo
+                     CALL iom_get( inum, jpdom_auto, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) )
+                     CALL iom_get( inum, jpdom_auto, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) ) 
+                     DO ib = 1, SIZE(dta%ssh)
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        td%ssh0(ib,itide,1) = ztr(ii,ij)
+                        td%ssh0(ib,itide,2) = zti(ii,ij)
+                     END DO
+                  END DO
+                  CALL iom_close( inum )
+               ENDIF
                !
                ! U fields
-               clfile = TRIM(filtide)//'_grid_U.nc'
-               CALL iom_open( clfile , inum ) 
-               igrd = 2                       ! Everything is at U-points here
-               DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) )
-                  DO ib = 1, ilen0(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
-                     td%u0(ib,itide,1) = ztr(ii,ij)
-                     td%u0(ib,itide,2) = zti(ii,ij)
-                  END DO
-               END DO
-               CALL iom_close( inum )
+               IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
+                  clfile = TRIM(filtide)//'_grid_U.nc'
+                  CALL iom_open( clfile , inum ) 
+                  igrd = 2                       ! Everything is at U-points here
+                  DO itide = 1, nb_harmo
+                     CALL iom_get(inum, jpdom_auto, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:),cd_type='U',psgn=-1._wp)
+                     CALL iom_get(inum, jpdom_auto, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:),cd_type='U',psgn=-1._wp)
+                     DO ib = 1, SIZE(dta%u2d)
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        td%u0(ib,itide,1) = ztr(ii,ij)
+                        td%u0(ib,itide,2) = zti(ii,ij)
+                     END DO
+                  END DO
+                  CALL iom_close( inum )
+               ENDIF
                !
                ! V fields
-               clfile = TRIM(filtide)//'_grid_V.nc'
-               CALL iom_open( clfile , inum ) 
-               igrd = 3                       ! Everything is at V-points here
-               DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) )
-                  DO ib = 1, ilen0(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
-                     td%v0(ib,itide,1) = ztr(ii,ij)
-                     td%v0(ib,itide,2) = zti(ii,ij)
-                  END DO
-               END DO  
-               CALL iom_close( inum )
+               IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
+                  clfile = TRIM(filtide)//'_grid_V.nc'
+                  CALL iom_open( clfile , inum ) 
+                  igrd = 3                       ! Everything is at V-points here
+                  DO itide = 1, nb_harmo
+                     CALL iom_get(inum, jpdom_auto, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:),cd_type='V',psgn=-1._wp)
+                     CALL iom_get(inum, jpdom_auto, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:),cd_type='V',psgn=-1._wp)
+                     DO ib = 1, SIZE(dta%v2d)
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        td%v0(ib,itide,1) = ztr(ii,ij)
+                        td%v0(ib,itide,2) = zti(ii,ij)
+                     END DO
+                  END DO
+                  CALL iom_close( inum )
+               ENDIF
                !
                DEALLOCATE( ztr, zti ) 
@@ -218,5 +221,5 @@
                ! Read tidal data only on bdy segments
                ! 
-               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
+               ALLOCATE( dta_read( MAXVAL( idx_bdy(ib_bdy)%nblen(:) ), 1, 1 ) )
                !
                ! Open files and read in tidal forcing data
@@ -225,27 +228,36 @@
                DO itide = 1, nb_harmo
                   !                                                              ! SSH fields
-                  clfile = TRIM(filtide)//TRIM(tide_harmonics(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, 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, idx_bdy(ib_bdy)%nbmap(:,1) )
-                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
-                  CALL iom_close( inum )
+                  IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
+                     isz = SIZE(dta%ssh)
+                     clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc'
+                     CALL iom_open( clfile, inum )
+                     CALL fld_map( inum, 'z1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
+                     td%ssh0(:,itide,1) = dta_read(1:isz,1,1)
+                     CALL fld_map( inum, 'z2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
+                     td%ssh0(:,itide,2) = dta_read(1:isz,1,1)
+                     CALL iom_close( inum )
+                  ENDIF
                   !                                                              ! U fields
-                  clfile = TRIM(filtide)//TRIM(tide_harmonics(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, 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, idx_bdy(ib_bdy)%nbmap(:,2) )
-                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
-                  CALL iom_close( inum )
+                  IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
+                     isz = SIZE(dta%u2d)
+                     clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc'
+                     CALL iom_open( clfile, inum )
+                     CALL fld_map( inum, 'u1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
+                     td%u0(:,itide,1) = dta_read(1:isz,1,1)
+                     CALL fld_map( inum, 'u2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
+                     td%u0(:,itide,2) = dta_read(1:isz,1,1)
+                     CALL iom_close( inum )
+                  ENDIF
                   !                                                              ! V fields
-                  clfile = TRIM(filtide)//TRIM(tide_harmonics(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, 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, idx_bdy(ib_bdy)%nbmap(:,3) )
-                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
-                  CALL iom_close( inum )
+                  IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
+                     isz = SIZE(dta%v2d)
+                     clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc'
+                     CALL iom_open( clfile, inum )
+                     CALL fld_map( inum, 'v1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
+                     td%v0(:,itide,1) = dta_read(1:isz,1,1)
+                     CALL fld_map( inum, 'v2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
+                     td%v0(:,itide,2) = dta_read(1:isz,1,1)
+                     CALL iom_close( inum )
+                  ENDIF
                   !
                END DO ! end loop on tidal components
@@ -254,13 +266,4 @@
                !
             ENDIF ! ln_bdytide_2ddta=.true.
-            !
-            ! Allocate slow varying data in the case of time splitting:
-            ! Do it anyway because at this stage knowledge of free surface scheme is unknown
-            ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
-            ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
-            ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
-            dta_bdy_s(ib_bdy)%ssh(:) = 0._wp
-            dta_bdy_s(ib_bdy)%u2d(:) = 0._wp
-            dta_bdy_s(ib_bdy)%v2d(:) = 0._wp
             !
          ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2
@@ -283,7 +286,5 @@
       !
       LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step
-      INTEGER  ::   itide, ib_bdy, ib, igrd   ! loop indices
-      INTEGER, DIMENSION(jpbgrd)   ::   ilen0 
-      INTEGER, DIMENSION(1:jpbgrd) ::   nblen, nblenrim  ! short cuts
+      INTEGER  ::   itide, ib_bdy, ib         ! loop indices
       REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist, zt_offset   
       !!----------------------------------------------------------------------
@@ -310,11 +311,4 @@
          IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
             !
-            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
-            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
-            !
-            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:)
-            ELSE                                   ;   ilen0(:) = nblenrim(:)
-            ENDIF     
-            !
             ! We refresh nodal factors every day below
             ! This should be done somewhere else
@@ -337,7 +331,7 @@
             ! If time splitting, initialize arrays from slow varying open boundary data:
             IF ( PRESENT(kit) ) THEN           
-               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))
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) dta_bdy(ib_bdy)%ssh(:) = dta_bdy_s(ib_bdy)%ssh(:)
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) dta_bdy(ib_bdy)%u2d(:) = dta_bdy_s(ib_bdy)%u2d(:)
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) dta_bdy(ib_bdy)%v2d(:) = dta_bdy_s(ib_bdy)%v2d(:)
             ENDIF
             !
@@ -349,7 +343,6 @@
                z_sist = zramp * SIN( z_sarg )
                !
-               IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN
-                  igrd=1                              ! SSH on tracer grid
-                  DO ib = 1, ilen0(igrd)
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) THEN   ! SSH on tracer grid
+                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%ssh)
                      dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
                         &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
@@ -358,13 +351,14 @@
                ENDIF
                !
-               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
-                  igrd=2                              ! U grid
-                  DO ib = 1, ilen0(igrd)
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) THEN  ! U grid
+                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%u2d)
                      dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
                         &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
                         &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
                   END DO
-                  igrd=3                              ! V grid
-                  DO ib = 1, ilen0(igrd) 
+               ENDIF
+               !
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) THEN   ! V grid
+                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%v2d)
                      dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
                         &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
@@ -372,6 +366,7 @@
                   END DO
                ENDIF
+               !
             END DO             
-         END IF
+         ENDIF
       END DO
       !
@@ -386,31 +381,31 @@
       TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
       !
-      INTEGER ::   itide, igrd, ib       ! dummy loop indices
-      INTEGER, DIMENSION(1) ::   ilen0   ! length of boundary data (from OBC arrays)
+      INTEGER ::   itide, isz, ib       ! dummy loop indices
       REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
       !!----------------------------------------------------------------------
       !
-      igrd=1   
-                              ! SSH on tracer grid.
-      ilen0(1) =  SIZE(td%ssh0(:,1,1))
-      !
-      ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) )
-      !
-      DO itide = 1, nb_harmo
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
-            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
+      IF( ASSOCIATED(td%ssh0) ) THEN   ! SSH on tracer grid.
+         !
+         isz = SIZE( td%ssh0, dim = 1 )
+         ALLOCATE( mod_tide(isz), phi_tide(isz) )
+         !
+         DO itide = 1, nb_harmo
+            DO ib = 1, isz
+               mod_tide(ib)=SQRT( td%ssh0(ib,itide,1)*td%ssh0(ib,itide,1) + td%ssh0(ib,itide,2)*td%ssh0(ib,itide,2) )
+               phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
+            END DO
+            DO ib = 1, isz
+               mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+               phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u
+            END DO
+            DO ib = 1, isz
+               td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
+               td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
+            END DO
          END DO
-         DO ib = 1 , ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
-            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u
-         ENDDO
-         DO ib = 1 , ilen0(igrd)
-            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
-            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
-         ENDDO
-      END DO
-      !
-      DEALLOCATE( mod_tide, phi_tide )
+         !
+         DEALLOCATE( mod_tide, phi_tide )
+         !
+      ENDIF
       !
    END SUBROUTINE tide_init_elevation
@@ -424,55 +419,57 @@
       TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
       !
-      INTEGER ::   itide, igrd, ib       ! dummy loop indices
-      INTEGER, DIMENSION(3) ::   ilen0   ! length of boundary data (from OBC arrays)
+      INTEGER ::   itide, isz, ib        ! dummy loop indices
       REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
       !!----------------------------------------------------------------------
       !
-      ilen0(2) =  SIZE(td%u0(:,1,1))
-      ilen0(3) =  SIZE(td%v0(:,1,1))
-      !
-      igrd=2                                 ! U grid.
-      !
-      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
-      !
-      DO itide = 1, nb_harmo
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
-            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
+      IF( ASSOCIATED(td%u0) ) THEN   ! U grid. we use bdy u2d on this mpi subdomain
+         !
+         isz = SIZE( td%u0, dim = 1 )
+         ALLOCATE( mod_tide(isz), phi_tide(isz) )
+         !
+         DO itide = 1, nb_harmo
+            DO ib = 1, isz
+               mod_tide(ib)=SQRT( td%u0(ib,itide,1)*td%u0(ib,itide,1) + td%u0(ib,itide,2)*td%u0(ib,itide,2) )
+               phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
+            END DO
+            DO ib = 1, isz
+               mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+               phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
+            END DO
+            DO ib = 1, isz
+               td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
+               td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
+            END DO
          END DO
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
-            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
-         ENDDO
-         DO ib = 1, ilen0(igrd)
-            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
-            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
-         ENDDO
-      END DO
-      !
-      DEALLOCATE( mod_tide , phi_tide )
-      !
-      igrd=3                                 ! V grid.
-      !
-      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
-
-      DO itide = 1, nb_harmo
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
-            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
+         !
+         DEALLOCATE( mod_tide, phi_tide )
+         !
+      ENDIF
+      !
+      IF( ASSOCIATED(td%v0) ) THEN   ! V grid. we use bdy u2d on this mpi subdomain
+         !
+         isz = SIZE( td%v0, dim = 1 )
+         ALLOCATE( mod_tide(isz), phi_tide(isz) )
+         !
+         DO itide = 1, nb_harmo
+            DO ib = 1, isz
+               mod_tide(ib)=SQRT( td%v0(ib,itide,1)*td%v0(ib,itide,1) + td%v0(ib,itide,2)*td%v0(ib,itide,2) )
+               phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
+            END DO
+            DO ib = 1, isz
+               mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+               phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
+            END DO
+            DO ib = 1, isz
+               td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
+               td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
+            END DO
          END DO
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
-            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
-         ENDDO
-         DO ib = 1, ilen0(igrd)
-            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
-            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
-         ENDDO
-      END DO
-      !
-      DEALLOCATE( mod_tide, phi_tide )
-      !
-  END SUBROUTINE tide_init_velocities
+         !
+         DEALLOCATE( mod_tide, phi_tide )
+         !
+      ENDIF
+      !
+   END SUBROUTINE tide_init_velocities
 
    !!======================================================================
Index: NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdytra.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdytra.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/OCE/BDY/bdytra.F90	(revision 13540)
@@ -61,5 +61,5 @@
          IF( ir == 0 ) THEN   ;   llrim0 = .TRUE.
          ELSE                 ;   llrim0 = .FALSE.
-         END IF
+         ENDIF
          DO ib_bdy=1, nb_bdy
             !
@@ -69,16 +69,16 @@
             DO jn = 1, jpts
                !
-               SELECT CASE( TRIM(cn_tra(ib_bdy)) )
+               SELECT CASE( cn_tra(ib_bdy) )
                CASE('none'        )   ;   CYCLE
                CASE('frs'         )   ! treat the whole boundary at once
-                  IF( ir == 0 ) CALL bdy_frs ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), zdta(jn)%tra )
+                  IF( ir == 0 )           CALL bdy_frs ( idx_bdy(ib_bdy),                    pts(:,:,:,jn,Kaa), zdta(jn)%tra )
                CASE('specified'   )   ! treat the whole rim      at once
-                  IF( ir == 0 ) CALL bdy_spe ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), zdta(jn)%tra )
-               CASE('neumann'     )   ;   CALL bdy_nmn ( idx_bdy(ib_bdy), igrd         , pts(:,:,:,jn,Kaa), llrim0 )   ! tsa masked
-               CASE('orlanski'    )   ;   CALL bdy_orl ( idx_bdy(ib_bdy), pts(:,:,:,jn,Kbb), pts(:,:,:,jn,Kaa), &
-                    & zdta(jn)%tra, llrim0, ll_npo=.false. )
-               CASE('orlanski_npo')   ;   CALL bdy_orl ( idx_bdy(ib_bdy), pts(:,:,:,jn,Kbb), pts(:,:,:,jn,Kaa), &
-                    & zdta(jn)%tra, llrim0, ll_npo=.true.  )
-               CASE('runoff'      )   ;   CALL bdy_rnf ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), jn, llrim0 )
+                  IF( ir == 0 )           CALL bdy_spe ( idx_bdy(ib_bdy),                    pts(:,:,:,jn,Kaa), zdta(jn)%tra )
+               CASE('neumann'     )   ;   CALL bdy_nmn ( idx_bdy(ib_bdy), igrd             , pts(:,:,:,jn,Kaa), llrim0 )   ! tsa masked
+               CASE('orlanski'    )   ;   CALL bdy_orl ( idx_bdy(ib_bdy), pts(:,:,:,jn,Kbb), pts(:,:,:,jn,Kaa), zdta(jn)%tra,   &
+                  &                                      llrim0, ll_npo=.FALSE. )
+               CASE('orlanski_npo')   ;   CALL bdy_orl ( idx_bdy(ib_bdy), pts(:,:,:,jn,Kbb), pts(:,:,:,jn,Kaa), zdta(jn)%tra,   &
+                  &                                      llrim0, ll_npo=.TRUE.  )
+               CASE('runoff'      )   ;   CALL bdy_rnf ( idx_bdy(ib_bdy),                    pts(:,:,:,jn,Kaa), jn, llrim0 )
                CASE DEFAULT           ;   CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' )
                END SELECT
@@ -88,7 +88,7 @@
          !
          IF( nn_hls > 1 .AND. ir == 1 ) CYCLE   ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
-         IF( nn_hls == 1 ) THEN   ;   llsend1(:) = .false.   ;   llrecv1(:) = .false.   ;   END IF
+         IF( nn_hls == 1 ) THEN   ;   llsend1(:) = .false.   ;   llrecv1(:) = .false.   ;   ENDIF
          DO ib_bdy=1, nb_bdy
-            SELECT CASE( TRIM(cn_tra(ib_bdy)) )
+            SELECT CASE( cn_tra(ib_bdy) )
             CASE('neumann','runoff')
                llsend1(:) = llsend1(:) .OR. lsend_bdyint(ib_bdy,1,:,ir)   ! possibly every direction, T points
@@ -100,6 +100,6 @@
          END DO
          IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdytra', pts(:,:,:,jn,Kaa), 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
-         END IF
+            CALL lbc_lnk( 'bdytra', pts(:,:,:,jn,Kaa), 'T',  1.0_wp, kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
+         ENDIF
          !
       END DO   ! ir
@@ -135,5 +135,5 @@
             pt(ii,ij,1:jpkm1) = 0.1 * tmask(ii,ij,1:jpkm1)
          END DO
-      END IF
+      ENDIF
       !
    END SUBROUTINE bdy_rnf
