Index: NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icb_oce.F90
===================================================================
--- NEMO/trunk/src/OCE/ICB/icb_oce.F90	(revision 10570)
+++ NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icb_oce.F90	(revision 10679)
@@ -63,4 +63,6 @@
       REAL(wp) ::   uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi    ! properties of iceberg environment 
       REAL(wp) ::   mass_of_bits, heat_density
+      REAL(wp), DIMENSION(4) :: xRK1, xRK2, xRK3, xRK4
+      REAL(wp), DIMENSION(4) :: yRK1, yRK2, yRK3, yRK4
    END TYPE point
 
@@ -85,11 +87,5 @@
    ! Extra arrays with bigger halo, needed when interpolating forcing onto iceberg position
    ! particularly for MPP when iceberg can lie inside T grid but outside U, V, or f grid
-   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   uo_e, vo_e
-   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ff_e, tt_e, fr_e, hicth
-   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ua_e, va_e
-   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_e
-#if defined key_si3 || defined key_cice
-   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ui_e, vi_e
-#endif
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hicth
 
    !!gm almost all those PARAM ARE defined in NEMO
@@ -168,15 +164,5 @@
       icb_alloc = icb_alloc + ill
       !
-      ! expanded arrays for bilinear interpolation
-      ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,   &
-         &      vo_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,   &
-#if defined key_si3 || defined key_cice
-         &      ui_e(0:jpi+1,0:jpj+1) ,                            &
-         &      vi_e(0:jpi+1,0:jpj+1) ,                            &
-#endif
-         &      ff_e(0:jpi+1,0:jpj+1) , fr_e(0:jpi+1,0:jpj+1)  ,   &
-         &      tt_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,   &
-         &      hicth(0:jpi+1,0:jpj+1),                            &
-         &      first_width(nclasses) , first_length(nclasses) ,   &
+      ALLOCATE( first_width(nclasses) , first_length(nclasses) ,   &
          &      src_calving (jpi,jpj) ,                            &
          &      src_calving_hflx(jpi,jpj) , STAT=ill)
Index: NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbdyn.F90
===================================================================
--- NEMO/trunk/src/OCE/ICB/icbdyn.F90	(revision 10570)
+++ NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbdyn.F90	(revision 10679)
@@ -18,4 +18,8 @@
    USE icbutl         ! iceberg utility routines
    USE icbdia         ! iceberg budget routines
+   USE icblbc         ! iceberg: lateral boundary routines (including mpp)
+   !
+   USE in_out_manager
+   USE lib_mpp        ! massively parallel library 
 
    IMPLICIT NONE
@@ -68,99 +72,198 @@
       zdt_6 = zdt / 6._wp
 
-      berg => first_berg                    ! start from the first berg
-      !
-      DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
+      IF( ASSOCIATED(first_berg) ) THEN
+         berg => first_berg                    ! start from the first berg
          !
-         pt => berg%current_point
-
-         ll_bounced = .FALSE.
-
-
-         ! STEP 1 !
-         ! ====== !
-         zxi1 = pt%xi   ;   zuvel1 = pt%uvel     !**   X1 in (i,j)  ;  V1 in m/s
-         zyj1 = pt%yj   ;   zvvel1 = pt%vvel
-
-
-         !                                         !**   A1 = A(X1,V1)
-         CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     &
-            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 )
-         !
-         zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt
-         zv1 = zvvel1 / ze2
-
-         ! STEP 2 !
-         ! ====== !
-         !                                         !**   X2 = X1+dt/2*V1   ;   V2 = V1+dt/2*A1
-         ! position using di/dt & djdt   !   V2  in m/s
-         zxi2 = zxi1 + zdt_2 * zu1          ;   zuvel2 = zuvel1 + zdt_2 * zax1
-         zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1
-         !
-         CALL icb_ground( zxi2, zxi1, zu1,   &
-            &             zyj2, zyj1, zv1, ll_bounced )
-
-         !                                         !**   A2 = A(X2,V2)
-         CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
-            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )
-         !
-         zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt
-         zv2 = zvvel2 / ze2
-         !
-         ! STEP 3 !
-         ! ====== !
-         !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3)
-         zxi3  = zxi1  + zdt_2 * zu2   ;   zuvel3 = zuvel1 + zdt_2 * zax2
-         zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2
-         !
-         CALL icb_ground( zxi3, zxi1, zu3,   &
-            &                zyj3, zyj1, zv3, ll_bounced )
-
-         !                                         !**   A3 = A(X3,V3)
-         CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
-            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt )
-         !
-         zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt
-         zv3 = zvvel3 / ze2
-
-         ! STEP 4 !
-         ! ====== !
-         !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3
-         zxi4 = zxi1 + zdt * zu3   ;   zuvel4 = zuvel1 + zdt * zax3
-         zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3
-
-         CALL icb_ground( zxi4, zxi1, zu4,   &
-            &             zyj4, zyj1, zv4, ll_bounced )
-
-         !                                         !**   A4 = A(X4,V4)
-         CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    &
-            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt )
-
-         zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt
-         zv4 = zvvel4 / ze2
-
-         ! FINAL STEP !
-         ! ========== !
-         !                                         !**   Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
-         !                                         !**   Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
-         zxi_n   = pt%xi   + zdt_6 * (  zu1  + 2.*(zu2  + zu3 ) + zu4  )
-         zyj_n   = pt%yj   + zdt_6 * (  zv1  + 2.*(zv2  + zv3 ) + zv4  )
-         zuvel_n = pt%uvel + zdt_6 * (  zax1 + 2.*(zax2 + zax3) + zax4 )
-         zvvel_n = pt%vvel + zdt_6 * (  zay1 + 2.*(zay2 + zay3) + zay4 )
-
-         CALL icb_ground( zxi_n, zxi1, zuvel_n,   &
-            &             zyj_n, zyj1, zvvel_n, ll_bounced )
-
-         pt%uvel = zuvel_n                        !** save in berg structure
-         pt%vvel = zvvel_n
-         pt%xi   = zxi_n
-         pt%yj   = zyj_n
-
-         ! update actual position
-         pt%lon  = icb_utl_bilin_x(glamt, pt%xi, pt%yj )
-         pt%lat  = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' )
-
-         berg => berg%next                         ! switch to the next berg
-         !
-      END DO                                  !==  end loop over all bergs  ==!
+         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
+            !
+            pt => berg%current_point
+            !
+            ! STEP 1 !
+            ! ====== !
+            zxi1 = pt%xi   ;   zuvel1 = pt%uvel     !**   X1 in (i,j)  ;  V1 in m/s
+            zyj1 = pt%yj   ;   zvvel1 = pt%vvel
+            !
+            !                                         !**   A1 = A(X1,V1)
+            CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     &
+               &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 )
+            !
+            zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt
+            zv1 = zvvel1 / ze2
+            !
+            ! STEP 2 !
+            ! ====== !
+            !                                         !**   X2 = X1+dt/2*V1   ;   V2 = V1+dt/2*A1
+            ! position using di/dt & djdt   !   V2  in m/s
+            zxi2 = zxi1 + zdt_2 * zu1          ;   zuvel2 = zuvel1 + zdt_2 * zax1
+            zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1
+            !
+            ! check ground iceberg
+            CALL icb_ground( zxi2, zxi1, zu1,   &
+               &             zyj2, zyj1, zv1, ll_bounced )
+            !
+            ! SAVE intermediaite properties
+            pt%xRK1 = (/ zxi1, zuvel1, zu1, zax1 /); pt%yRK1 = (/ zyj1, zvvel1, zv1, zay1 /); 
+            pt%xRK2 = (/ zxi2, zuvel2, 0.0, 0.0  /); pt%yRK2 = (/ zyj2, zvvel2, 0.0, 0.0  /); 
+            pt%xi = zxi2 ; pt%yj = zyj2 ;
+            ! 
+            berg => berg%next                         ! switch to the next berg
+            !
+         END DO
+      END IF
+      !
+      ! EXCHANGE BERG (in case first step fall out of the inner domain)
+      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs
+      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case
+      ENDIF
+      !
+      IF( ASSOCIATED(first_berg) ) THEN
+         berg => first_berg                    ! start from the first berg
+         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
+            !
+            pt => berg%current_point
+            ! 
+            ! LOAD intermediaite properties
+            zxi1 = pt%xRK1(1) ; zuvel1 = pt%xRK1(2) ; zu1 = pt%xRK1(3) ; zax1 = pt%xRK1(4) 
+            zxi2 = pt%xRK2(1) ; zuvel2 = pt%xRK2(2) ; zu2 = pt%xRK2(3) ; zax2 = pt%xRK2(4) 
+            zyj1 = pt%yRK1(1) ; zvvel1 = pt%yRK1(2) ; zv1 = pt%yRK1(3) ; zay1 = pt%yRK1(4) 
+            zyj2 = pt%yRK2(1) ; zvvel2 = pt%yRK2(2) ; zv2 = pt%yRK2(3) ; zay2 = pt%yRK2(4) 
+            !
+            !                                         !**   A2 = A(X2,V2)
+            CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    &
+               &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )
+            !
+            zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt
+            zv2 = zvvel2 / ze2
+            !
+            ! STEP 3 !
+            ! ====== !
+            !                                         !**  X3 = X1+dt/2*V2  ;   V3 = V1+dt/2*A2; A3=A(X3)
+            zxi3  = zxi1  + zdt_2 * zu2   ;   zuvel3 = zuvel1 + zdt_2 * zax2
+            zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2
+            !
+            CALL icb_ground( zxi3, zxi1, zu3,   &
+               &             zyj3, zyj1, zv3, ll_bounced )
+            !
+            ! SAVE intermediaite properties
+            pt%xRK2(3:4) =               (/ zu2, zax2 /); pt%yRK2(3:4) =               (/ zv2, zay2 /);
+            pt%xRK3      = (/ zxi3, zuvel3, 0.0, 0.0  /); pt%yRK3      = (/ zyj3, zvvel3, 0.0, 0.0  /); 
+            pt%xi = zxi3 ; pt%yj = zyj3 ;
+            ! 
+            berg => berg%next                         ! switch to the next berg
+            !
+         END DO
+      END IF
+      !
+      ! EXCHANGE BERG (in case first step fall out of the inner domain)
+      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs
+      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case
+      ENDIF
+      !
+      IF( ASSOCIATED(first_berg) ) THEN
+         berg => first_berg                    ! start from the first berg
+         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
+            !
+            pt => berg%current_point
+            ! 
+            ! LOAD intermediaite properties
+            zxi1 = pt%xRK1(1) ; zuvel1 = pt%xRK1(2) ; zu1 = pt%xRK1(3) ; zax1 = pt%xRK1(4) 
+            zxi3 = pt%xRK3(1) ; zuvel3 = pt%xRK3(2) ; zu3 = pt%xRK3(3) ; zax3 = pt%xRK3(4) 
+            zyj1 = pt%yRK1(1) ; zvvel1 = pt%yRK1(2) ; zv1 = pt%yRK1(3) ; zay1 = pt%yRK1(4) 
+            zyj3 = pt%yRK3(1) ; zvvel3 = pt%yRK3(2) ; zv3 = pt%yRK3(3) ; zay3 = pt%yRK3(4) 
+            !
+            !                                         !**   A3 = A(X3,V3)
+            CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    &
+               &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt )
+            !
+            zu3 = zuvel3 / ze1                        !**   V3 in d(i,j)/dt
+            zv3 = zvvel3 / ze2
+            !
+            ! STEP 4 !
+            ! ====== !
+            !                                         !**   X4 = X1+dt*V3   ;   V4 = V1+dt*A3
+            zxi4 = zxi1 + zdt * zu3   ;   zuvel4 = zuvel1 + zdt * zax3
+            zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3
+            !
+            CALL icb_ground( zxi4, zxi1, zu4,   &
+               &             zyj4, zyj1, zv4, ll_bounced )
+            !
+            ! SAVE intermediaite properties
+            pt%xRK3(3:4) =               (/ zu3, zax3 /); pt%yRK3(3:4) =               (/ zv3, zay3 /); 
+            pt%xRK4      = (/ zxi4, zuvel4, 0.0, 0.0  /); pt%yRK4      = (/ zyj4, zvvel4, 0.0, 0.0  /); 
+            pt%xi = zxi4 ; pt%yj = zyj4 ;
+            !
+            berg => berg%next                         ! switch to the next berg
+            !
+         END DO
+      END IF
+      !
+      ! EXCHANGE BERG (in case first step fall out of the inner domain)
+      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs
+      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case
+      ENDIF
+      !
+      IF( ASSOCIATED(first_berg) ) THEN
+         berg => first_berg                    ! start from the first berg
+         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
+            !
+            pt => berg%current_point
+            ! 
+            ! LOAD intermediaite properties
+            zxi1 = pt%xRK1(1) ; zuvel1 = pt%xRK1(2) ; zu1 = pt%xRK1(3) ; zax1 = pt%xRK1(4) 
+            zxi2 = pt%xRK2(1) ; zuvel2 = pt%xRK2(2) ; zu2 = pt%xRK2(3) ; zax2 = pt%xRK2(4) 
+            zxi3 = pt%xRK3(1) ; zuvel3 = pt%xRK3(2) ; zu3 = pt%xRK3(3) ; zax3 = pt%xRK3(4) 
+            zxi4 = pt%xRK4(1) ; zuvel4 = pt%xRK4(2) ; zu4 = pt%xRK4(3) ; zax4 = pt%xRK4(4) 
+            zyj1 = pt%yRK1(1) ; zvvel1 = pt%yRK1(2) ; zv1 = pt%yRK1(3) ; zay1 = pt%yRK1(4) 
+            zyj2 = pt%yRK2(1) ; zvvel2 = pt%yRK2(2) ; zv2 = pt%yRK2(3) ; zay2 = pt%yRK2(4) 
+            zyj3 = pt%yRK3(1) ; zvvel3 = pt%yRK3(2) ; zv3 = pt%yRK3(3) ; zay3 = pt%yRK3(4) 
+            zyj4 = pt%yRK4(1) ; zvvel4 = pt%yRK4(2) ; zv4 = pt%yRK4(3) ; zay4 = pt%yRK4(4) 
+            !
+            !                                         !**   A4 = A(X4,V4)
+            CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    &
+               &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt )
+            !
+            zu4 = zuvel4 / ze1                        !**   V4 in d(i,j)/dt
+            zv4 = zvvel4 / ze2
+            !
+            ! FINAL STEP !
+            ! ========== !
+            !                                         !**   Xn = X1+dt*(V1+2*V2+2*V3+V4)/6
+            !                                         !**   Vn = V1+dt*(A1+2*A2+2*A3+A4)/6
+            zxi_n   = zxi1   + zdt_6 * (  zu1  + 2._wp*(zu2  + zu3 ) + zu4  )
+            zyj_n   = zyj1   + zdt_6 * (  zv1  + 2._wp*(zv2  + zv3 ) + zv4  )
+            zuvel_n = zuvel1 + zdt_6 * (  zax1 + 2._wp*(zax2 + zax3) + zax4 )
+            zvvel_n = zvvel1 + zdt_6 * (  zay1 + 2._wp*(zay2 + zay3) + zay4 )
+            !
+            CALL icb_ground( zxi_n, zxi1, zuvel_n,   &
+               &             zyj_n, zyj1, zvvel_n, ll_bounced )
+            !
+            !                                        !** save in berg structure
+            pt%xi   = zxi_n   ; pt%yj   = zyj_n ;
+            pt%uvel = zuvel_n ; pt%vvel = zvvel_n
+            !
+            berg => berg%next                        ! switch to the next berg
+            !
+         END DO
+      END IF
+      !
+      ! EXCHANGE BERG (in case first step fall out of the inner domain)
+      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs
+      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case
+      ENDIF
+      !
+      IF( ASSOCIATED(first_berg) ) THEN
+         berg => first_berg                    ! start from the first berg
+         DO WHILE ( ASSOCIATED(berg) )          !==  loop over all bergs  ==!
+            !
+            pt => berg%current_point
+            !
+            ! update actual position
+            pt%lon  = icb_utl_bilin_x(glamt, pt%xi, pt%yj )
+            pt%lat  = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' )
+            !
+            berg => berg%next                         ! switch to the next berg
+            !
+         END DO                                  !==  end loop over all bergs  ==!
+      END IF
       !
    END SUBROUTINE icb_dyn
@@ -182,21 +285,21 @@
       LOGICAL , INTENT(  out) ::   ld_bounced   ! bounced indicator
       !
-      INTEGER  ::   ii, ii0
-      INTEGER  ::   ij, ij0
+      INTEGER  ::   ii, iig, iig0
+      INTEGER  ::   ij, ijg, ijg0
       INTEGER  ::   ibounce_method
       !!----------------------------------------------------------------------
       !
+      ! (PM) is this variable used somewhere ?
       ld_bounced = .FALSE.
       !
-      ii0 = INT( pi0+0.5 )   ;   ij0 = INT( pj0+0.5 )       ! initial gridpoint position (T-cell)
-      ii  = INT( pi +0.5 )   ;   ij  = INT( pj +0.5 )       ! current     -         -
-      !
-      IF( ii == ii0  .AND.  ij == ij0  )   RETURN           ! berg remains in the same cell
+      ! work with global position in case pi and pi0 are on different domains
+      iig0 = INT( pi0 + 0.5_wp )   ;   ijg0 = INT( pj0 + 0.5_wp )       ! initial gridpoint position (T-cell)
+      iig  = INT( pi  + 0.5_wp )   ;   ijg  = INT( pj  + 0.5_wp )       ! current     -         -
+      !
+      IF( iig == iig0  .AND.  ijg == ijg0  )   RETURN           ! berg remains in the same cell
       !
       ! map into current processor
-      ii0 = mi1( ii0 )
-      ij0 = mj1( ij0 )
-      ii  = mi1( ii  )
-      ij  = mj1( ij  )
+      ii  = mi1( iig  )
+      ij  = mj1( ijg  )
       !
       IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one
@@ -221,9 +324,9 @@
          pv = 0._wp
       CASE ( 2 )
-         IF( ii0 /= ii ) THEN
+         IF( iig0 /= iig ) THEN
             pi = pi0                   ! return back to the initial position
             pu = 0._wp                 ! zeroing of velocity in the direction of the grounding
          ENDIF
-         IF( ij0 /= ij ) THEN
+         IF( ijg0 /= ijg ) THEN
             pj = pj0                   ! return back to the initial position
             pv = 0._wp                 ! zeroing of velocity in the direction of the grounding
Index: NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icblbc.F90
===================================================================
--- NEMO/trunk/src/OCE/ICB/icblbc.F90	(revision 10570)
+++ NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icblbc.F90	(revision 10679)
@@ -56,5 +56,5 @@
 
    INTEGER, PARAMETER, PRIVATE ::   jp_delta_buf = 25             ! Size by which to increment buffers
-   INTEGER, PARAMETER, PRIVATE ::   jp_buffer_width = 15+nkounts  ! items to store for each berg
+   INTEGER, PARAMETER, PRIVATE ::   jp_buffer_width = 47+nkounts  ! items to store for each berg
 
 #endif
@@ -90,5 +90,5 @@
          DO WHILE( ASSOCIATED(this) )
             pt => this%current_point
-            iine = INT( pt%xi + 0.5 )
+            iine = INT( pt%xi + 0.5_wp )
             IF( iine > mig(nicbei) ) THEN
                pt%xi = ricb_right + MOD(pt%xi, 1._wp ) - 1._wp
@@ -125,8 +125,8 @@
       DO WHILE( ASSOCIATED(this) )
          pt => this%current_point
-         ijne = INT( pt%yj + 0.5 )
+         ijne = INT( pt%yj + 0.5_wp )
          IF( ijne .GT. mjg(nicbej) ) THEN
             !
-            iine = INT( pt%xi + 0.5 )
+            iine = INT( pt%xi + 0.5_wp )
             ipts  = nicbfldpts (mi1(iine))
             !
@@ -232,5 +232,5 @@
          DO WHILE (ASSOCIATED(this))
             pt => this%current_point
-            iine = INT( pt%xi + 0.5 )
+            iine = INT( pt%xi + 0.5_wp )
             IF( ipe_E >= 0 .AND. iine > mig(nicbei) ) THEN
                tmpberg => this
@@ -242,5 +242,9 @@
                ENDIF
                ! deal with periodic case
-               tmpberg%current_point%xi = ricb_right + MOD(tmpberg%current_point%xi, 1._wp ) - 1._wp
+               tmpberg%current_point%xi      = ricb_right + MOD(tmpberg%current_point%xi, 1._wp ) - 1._wp
+               tmpberg%current_point%xRK1(1) = ricb_right + MOD(tmpberg%current_point%xRK1(1), 1._wp ) - 1._wp
+               tmpberg%current_point%xRK2(1) = ricb_right + MOD(tmpberg%current_point%xRK2(1), 1._wp ) - 1._wp
+               tmpberg%current_point%xRK3(1) = ricb_right + MOD(tmpberg%current_point%xRK3(1), 1._wp ) - 1._wp
+               tmpberg%current_point%xRK4(1) = ricb_right + MOD(tmpberg%current_point%xRK4(1), 1._wp ) - 1._wp
                ! now pack it into buffer and delete from list
                CALL icb_pack_into_buffer( tmpberg, obuffer_e, ibergs_to_send_e)
@@ -255,5 +259,9 @@
                ENDIF
                ! deal with periodic case
-               tmpberg%current_point%xi = ricb_left + MOD(tmpberg%current_point%xi, 1._wp )
+               tmpberg%current_point%xi      = ricb_left + MOD(tmpberg%current_point%xi, 1._wp )
+               tmpberg%current_point%xRK1(1) = ricb_left + MOD(tmpberg%current_point%xRK1(1), 1._wp )
+               tmpberg%current_point%xRK2(1) = ricb_left + MOD(tmpberg%current_point%xRK2(1), 1._wp )
+               tmpberg%current_point%xRK3(1) = ricb_left + MOD(tmpberg%current_point%xRK3(1), 1._wp )
+               tmpberg%current_point%xRK4(1) = ricb_left + MOD(tmpberg%current_point%xRK4(1), 1._wp )
                ! now pack it into buffer and delete from list
                CALL icb_pack_into_buffer( tmpberg, obuffer_w, ibergs_to_send_w)
@@ -370,5 +378,5 @@
          DO WHILE (ASSOCIATED(this))
             pt => this%current_point
-            ijne = INT( pt%yj + 0.5 )
+            ijne = INT( pt%yj + 0.5_wp )
             IF( ipe_N >= 0 .AND. ijne .GT. mjg(nicbej) ) THEN
                tmpberg => this
@@ -537,6 +545,6 @@
          DO WHILE (ASSOCIATED(this))
             pt => this%current_point
-            iine = INT( pt%xi + 0.5 )
-            ijne = INT( pt%yj + 0.5 )
+            iine = INT( pt%xi + 0.5_wp )
+            ijne = INT( pt%yj + 0.5_wp )
             IF( iine .LT. mig(nicbdi) .OR. &
                 iine .GT. mig(nicbei) .OR. &
@@ -544,7 +552,8 @@
                 ijne .GT. mjg(nicbej)) THEN
                i = i + 1
-               WRITE(numicb,*) 'berg lost in halo: ', this%number(:),iine,ijne
-               WRITE(numicb,*) '                   ', nimpp, njmpp
-               WRITE(numicb,*) '                   ', nicbdi, nicbei, nicbdj, nicbej
+               WRITE(numicb,*) 'berg lost in halo:  ', this%number(:),iine,ijne, pt%xi, pt%yj
+               WRITE(numicb,*) 'nimpp, njmpp        ', nimpp, njmpp
+               WRITE(numicb,*) 'nicb di, ei, dj, ej ', nicbdi     , nicbei     , nicbdj     , nicbej
+               WRITE(numicb,*) 'mijg di, ei, dj, ej ', mig(nicbdi), mig(nicbei), mjg(nicbdj), mjg(nicbej)
                CALL flush( numicb )
             ENDIF
@@ -611,6 +620,6 @@
                DO WHILE (ASSOCIATED(this))
                   pt => this%current_point
-                  iine = INT( pt%xi + 0.5 )
-                  ijne = INT( pt%yj + 0.5 )
+                  iine = INT( pt%xi + 0.5_wp )
+                  ijne = INT( pt%yj + 0.5_wp )
                   iproc = nicbflddest(mi1(iine))
                   IF( ijne .GT. mjg(nicbej) ) THEN
@@ -691,6 +700,6 @@
                DO WHILE (ASSOCIATED(this))
                   pt => this%current_point
-                  iine = INT( pt%xi + 0.5 )
-                  ijne = INT( pt%yj + 0.5 )
+                  iine = INT( pt%xi + 0.5_wp )
+                  ijne = INT( pt%yj + 0.5_wp )
                   ipts  = nicbfldpts (mi1(iine))
                   iproc = nicbflddest(mi1(iine))
@@ -791,4 +800,5 @@
       IF( .NOT. ASSOCIATED(pbuff) ) CALL icb_increase_buffer( pbuff, jp_delta_buf )
       IF( kb .GT. pbuff%size ) CALL icb_increase_buffer( pbuff, jp_delta_buf )
+      IF( kb .GT. pbuff%size ) PRINT *, 'SHOULD NOT SEE THIS'
 
       !! pack points into buffer
@@ -809,7 +819,17 @@
       pbuff%data(14,kb) = berg%current_point%heat_density
 
-      pbuff%data(15,kb) = berg%mass_scaling
+      pbuff%data(15:18,kb) = berg%current_point%xRK1
+      pbuff%data(19:22,kb) = berg%current_point%xRK2
+      pbuff%data(23:26,kb) = berg%current_point%xRK3
+      pbuff%data(27:30,kb) = berg%current_point%xRK4
+
+      pbuff%data(31:34,kb) = berg%current_point%yRK1
+      pbuff%data(35:38,kb) = berg%current_point%yRK2
+      pbuff%data(39:42,kb) = berg%current_point%yRK3
+      pbuff%data(43:46,kb) = berg%current_point%yRK4
+
+      pbuff%data(47,kb) = berg%mass_scaling
       DO k=1,nkounts
-         pbuff%data(15+k,kb) = REAL( berg%number(k), wp )
+         pbuff%data(47+k,kb) = REAL( berg%number(k), wp )
       END DO
       !
@@ -844,7 +864,17 @@
       pt%heat_density   =      pbuff%data(14,kb)
 
-      currentberg%mass_scaling =      pbuff%data(15,kb)
+      pt%xRK1 = pbuff%data(15:18,kb)
+      pt%xRK2 = pbuff%data(19:22,kb) 
+      pt%xRK3 = pbuff%data(23:26,kb)
+      pt%xRK4 = pbuff%data(27:30,kb)
+
+      pt%yRK1 = pbuff%data(31:34,kb)
+      pt%yRK2 = pbuff%data(35:38,kb)
+      pt%yRK3 = pbuff%data(39:42,kb)
+      pt%yRK4 = pbuff%data(43:46,kb)
+
+      currentberg%mass_scaling =      pbuff%data(47,kb)
       DO ik = 1, nkounts
-         currentberg%number(ik) = INT( pbuff%data(15+ik,kb) )
+         currentberg%number(ik) = INT( pbuff%data(47+ik,kb) )
       END DO
       !
Index: NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbstp.F90
===================================================================
--- NEMO/trunk/src/OCE/ICB/icbstp.F90	(revision 10570)
+++ NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbstp.F90	(revision 10679)
@@ -29,5 +29,4 @@
    USE icbclv         ! iceberg: calving routines
    USE icbthm         ! iceberg: thermodynamics routines
-   USE icblbc         ! iceberg: lateral boundary routines (including mpp)
    USE icbtrj         ! iceberg: trajectory I/O routines
    USE icbdia         ! iceberg: budget
@@ -91,8 +90,4 @@
  9100 FORMAT('kt= ',i8, ' day= ',i8,' secs=',i8)
       !
-      !                                   !* copy nemo forcing arrays into iceberg versions with extra halo
-      CALL icb_utl_copy()                 ! only necessary for variables not on T points
-      !
-      !
       !                       !==  process icebergs  ==!
       !                              !
@@ -104,9 +99,5 @@
       !                       !==  For each berg, evolve  ==!
       !
-      IF( ASSOCIATED(first_berg) )   CALL icb_dyn( kt )       ! ice berg dynamics
-
-      IF( lk_mpp ) THEN   ;          CALL icb_lbc_mpp()       ! Send bergs to other PEs
-      ELSE                ;          CALL icb_lbc()           ! Deal with any cyclic boundaries in non-mpp case
-      ENDIF
+                                     CALL icb_dyn( kt )       ! ice berg dynamics
 
       IF( ASSOCIATED(first_berg) )   CALL icb_thm( kt )       ! Ice berg thermodynamics (melting) + rolling
Index: NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbutl.F90
===================================================================
--- NEMO/trunk/src/OCE/ICB/icbutl.F90	(revision 10570)
+++ NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbutl.F90	(revision 10679)
@@ -31,5 +31,4 @@
    PRIVATE
 
-   PUBLIC   icb_utl_copy          ! routine called in icbstp module
    PUBLIC   icb_utl_interp        ! routine called in icbdyn, icbthm modules
    PUBLIC   icb_utl_bilin         ! routine called in icbini, icbdyn modules
@@ -54,66 +53,4 @@
 CONTAINS
 
-   SUBROUTINE icb_utl_copy()
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE icb_utl_copy  ***
-      !!
-      !! ** Purpose :   iceberg initialization.
-      !!
-      !! ** Method  : - blah blah
-      !!----------------------------------------------------------------------
-#if defined key_si3
-      REAL(wp), DIMENSION(jpi,jpj) :: zssh_lead_m    !    ocean surface (ssh_m) if ice is not embedded
-      !                                              !    ocean surface in leads if ice is embedded   
-#endif
-      ! copy nemo forcing arrays into iceberg versions with extra halo
-      ! only necessary for variables not on T points
-      ! and ssh which is used to calculate gradients
-
-      uo_e(:,:) = 0._wp   ;   uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1)
-      vo_e(:,:) = 0._wp   ;   vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1)
-      ff_e(:,:) = 0._wp   ;   ff_e(1:jpi,1:jpj) = ff_f (:,:) 
-      tt_e(:,:) = 0._wp   ;   tt_e(1:jpi,1:jpj) = sst_m(:,:)
-      fr_e(:,:) = 0._wp   ;   fr_e(1:jpi,1:jpj) = fr_i (:,:)
-      ua_e(:,:) = 0._wp   ;   ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk
-      va_e(:,:) = 0._wp   ;   va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk
-      !
-      CALL lbc_lnk_icb( 'icbutl', uo_e, 'U', -1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbutl', vo_e, 'V', -1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbutl', ff_e, 'F', +1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbutl', ua_e, 'U', -1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbutl', va_e, 'V', -1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbutl', fr_e, 'T', +1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbutl', tt_e, 'T', +1._wp, 1, 1 )
-#if defined key_si3
-      hicth(:,:) = 0._wp ;  hicth(1:jpi,1:jpj) = hm_i (:,:)  
-      ui_e(:,:) = 0._wp ;   ui_e(1:jpi, 1:jpj) = u_ice(:,:)
-      vi_e(:,:) = 0._wp ;   vi_e(1:jpi, 1:jpj) = v_ice(:,:)
-      !      
-      ! compute ssh slope using ssh_lead if embedded
-      zssh_lead_m(:,:) = ice_var_sshdyn(ssh_m, snwice_mass, snwice_mass_b)
-      ssh_e(:,:) = 0._wp ;  ssh_e(1:jpi, 1:jpj) = zssh_lead_m(:,:) * tmask(:,:,1)
-      !
-      CALL lbc_lnk_icb( 'icbutl', hicth, 'T', +1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbutl', ui_e , 'U', -1._wp, 1, 1 )
-      CALL lbc_lnk_icb( 'icbutl', vi_e , 'V', -1._wp, 1, 1 )
-#else
-      ssh_e(:,:) = 0._wp ;  ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1)
-#endif
-
-      !! special for ssh which is used to calculate slope
-      !! so fudge some numbers all the way around the boundary
-      ssh_e(0    ,    :) = ssh_e(1  ,  :)
-      ssh_e(jpi+1,    :) = ssh_e(jpi,  :)
-      ssh_e(:    ,    0) = ssh_e(:  ,  1)
-      ssh_e(:    ,jpj+1) = ssh_e(:  ,jpj)
-      ssh_e(0,0)         = ssh_e(1,1)
-      ssh_e(jpi+1,0)     = ssh_e(jpi,1)
-      ssh_e(0,jpj+1)     = ssh_e(1,jpj)
-      ssh_e(jpi+1,jpj+1) = ssh_e(jpi,jpj)
-      CALL lbc_lnk_icb( 'icbutl', ssh_e, 'T', +1._wp, 1, 1 )
-      !
-   END SUBROUTINE icb_utl_copy
-
-
    SUBROUTINE icb_utl_interp( pi, pe1, puo, pui, pua, pssh_i,   &
       &                       pj, pe2, pvo, pvi, pva, pssh_j,   &
@@ -140,4 +77,5 @@
       !
       REAL(wp) ::   zcd, zmod       ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj) :: zssh
       !!----------------------------------------------------------------------
 
@@ -145,12 +83,12 @@
       pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
       !
-      puo  = icb_utl_bilin_h( uo_e, pi, pj, 'U' )             ! ocean velocities
-      pvo  = icb_utl_bilin_h( vo_e, pi, pj, 'V' )
-      psst = icb_utl_bilin_h( tt_e, pi, pj, 'T' )             ! SST
-      pcn  = icb_utl_bilin_h( fr_e , pi, pj, 'T' )            ! ice concentration
-      pff  = icb_utl_bilin_h( ff_e , pi, pj, 'F' )            ! Coriolis parameter
-      !
-      pua  = icb_utl_bilin_h( ua_e , pi, pj, 'U' )            ! 10m wind
-      pva  = icb_utl_bilin_h( va_e , pi, pj, 'V' )            ! here (ua,va) are stress => rough conversion from stress to speed
+      puo  = icb_utl_bilin( ssu_m(:,:) * umask(:,:,1), pi, pj, 'U' )             ! ocean velocities
+      pvo  = icb_utl_bilin( ssv_m(:,:) * vmask(:,:,1), pi, pj, 'V' )
+      psst = icb_utl_bilin( sst_m(:,:), pi, pj, 'T' )             ! SST
+      pcn  = icb_utl_bilin( fr_i(:,:) , pi, pj, 'T' )            ! ice concentration
+      pff  = icb_utl_bilin( ff_f(:,:) , pi, pj, 'F' )            ! Coriolis parameter
+      !
+      pua  = icb_utl_bilin( utau(:,:) * umask(:,:,1), pi, pj, 'U' )            ! 10m wind
+      pva  = icb_utl_bilin( vtau(:,:) * vmask(:,:,1), pi, pj, 'V' )            ! here (ua,va) are stress => rough conversion from stress to speed
       zcd  = 1.22_wp * 1.5e-3_wp                              ! air density * drag coefficient
       zmod = 1._wp / MAX(  1.e-20, SQRT(  zcd * SQRT( pua*pua + pva*pva)  )  )
@@ -159,86 +97,25 @@
 
 #if defined key_si3
-      pui = icb_utl_bilin_h( ui_e , pi, pj, 'U' )              ! sea-ice velocities
-      pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V' )
-      phi = icb_utl_bilin_h( hicth, pi, pj, 'T' )              ! ice thickness
+      pui = icb_utl_bilin( u_ice , pi, pj, 'U' )              ! sea-ice velocities
+      pvi = icb_utl_bilin( v_ice , pi, pj, 'V' )
+      phi = icb_utl_bilin( hm_i  , pi, pj, 'T' )              ! ice thickness
+      ! Estimate SSH gradient in i- and j-direction (centred evaluation)
+      ! compute ssh slope using ssh_lead if embedded
+      zssh(:,:) = ice_var_sshdyn(ssh_m, snwice_mass, snwice_mass_b)
 #else
       pui = 0._wp
       pvi = 0._wp
       phi = 0._wp
+      zssh(:,:) = ssh_m(:,:)
 #endif
+      zssh(:,:) = zssh(:,:) * tmask(:,:,1)
 
       ! Estimate SSH gradient in i- and j-direction (centred evaluation)
-      pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T' ) -   &
-         &       icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' )  ) / ( 0.2_wp * pe1 )
-      pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T' ) -   &
-         &       icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' )  ) / ( 0.2_wp * pe2 )
+      pssh_i = ( icb_utl_bilin( zssh(:,:), pi+0.1_wp, pj, 'T' ) -   &
+         &       icb_utl_bilin( zssh(:,:), pi-0.1_wp, pj, 'T' )  ) / ( 0.2_wp * pe1 )
+      pssh_j = ( icb_utl_bilin( zssh(:,:), pi, pj+0.1_wp, 'T' ) -   &
+         &       icb_utl_bilin( zssh(:,:), pi, pj-0.1_wp, 'T' )  ) / ( 0.2_wp * pe2 )
       !
    END SUBROUTINE icb_utl_interp
-
-
-   REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type )
-      !!----------------------------------------------------------------------
-      !!                  ***  FUNCTION icb_utl_bilin  ***
-      !!
-      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type
-      !!                this version deals with extra halo points
-      !!
-      !!       !!gm  CAUTION an optional argument should be added to handle
-      !!             the slip/no-slip conditions  ==>>> to be done later
-      !!
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) ::   pfld      ! field to be interpolated
-      REAL(wp)                            , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential
-      CHARACTER(len=1)                    , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points
-      !
-      INTEGER  ::   ii, ij   ! local integer
-      REAL(wp) ::   zi, zj   ! local real
-      !!----------------------------------------------------------------------
-      !
-      SELECT CASE ( cd_type )
-      CASE ( 'T' )
-         ! note that here there is no +0.5 added
-         ! since we're looking for four T points containing quadrant we're in of 
-         ! current T cell
-         ii = MAX(1, INT( pi     ))
-         ij = MAX(1, INT( pj     ))    ! T-point
-         zi = pi - REAL(ii,wp)
-         zj = pj - REAL(ij,wp)
-      CASE ( 'U' )
-         ii = MAX(1, INT( pi-0.5 ))
-         ij = MAX(1, INT( pj     ))    ! U-point
-         zi = pi - 0.5 - REAL(ii,wp)
-         zj = pj - REAL(ij,wp)
-      CASE ( 'V' )
-         ii = MAX(1, INT( pi     ))
-         ij = MAX(1, INT( pj-0.5 ))    ! V-point
-         zi = pi - REAL(ii,wp)
-         zj = pj - 0.5 - REAL(ij,wp)
-      CASE ( 'F' )
-         ii = MAX(1, INT( pi-0.5 ))
-         ij = MAX(1, INT( pj-0.5 ))    ! F-point
-         zi = pi - 0.5 - REAL(ii,wp)
-         zj = pj - 0.5 - REAL(ij,wp)
-      END SELECT
-      !
-      ! find position in this processor. Prevent near edge problems (see #1389)
-      !
-      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1
-      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi
-      ELSE                           ;   ii = mi1(ii)
-      ENDIF
-      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1
-      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj
-      ELSE                           ;   ij  = mj1(ij)
-      ENDIF
-      !
-      IF( ii == jpi )   ii = ii-1      
-      IF( ij == jpj )   ij = ij-1
-      !
-      icb_utl_bilin_h = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   &
-         &            + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) *     zj
-      !
-   END FUNCTION icb_utl_bilin_h
-
 
    REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type )
@@ -270,32 +147,22 @@
             zj = pj - REAL(ij,wp)
          CASE ( 'U' )
-            ii = MAX(1, INT( pi-0.5 ))
-            ij = MAX(1, INT( pj     ))    ! U-point
-            zi = pi - 0.5 - REAL(ii,wp)
+            ii = MAX(1, INT( pi-0.5_wp ))
+            ij = MAX(1, INT( pj        ))    ! U-point
+            zi = pi - 0.5_wp - REAL(ii,wp)
             zj = pj - REAL(ij,wp)
          CASE ( 'V' )
-            ii = MAX(1, INT( pi     ))
-            ij = MAX(1, INT( pj-0.5 ))    ! V-point
+            ii = MAX(1, INT( pi        ))
+            ij = MAX(1, INT( pj-0.5_wp ))    ! V-point
             zi = pi - REAL(ii,wp)
-            zj = pj - 0.5 - REAL(ij,wp)
+            zj = pj - 0.5_wp - REAL(ij,wp)
          CASE ( 'F' )
-            ii = MAX(1, INT( pi-0.5 ))
-            ij = MAX(1, INT( pj-0.5 ))    ! F-point
-            zi = pi - 0.5 - REAL(ii,wp)
-            zj = pj - 0.5 - REAL(ij,wp)
+            ii = MAX(1, INT( pi-0.5_wp ))
+            ij = MAX(1, INT( pj-0.5_wp ))    ! F-point
+            zi = pi - 0.5_wp - REAL(ii,wp)
+            zj = pj - 0.5_wp - REAL(ij,wp)
       END SELECT
       !
-      ! find position in this processor. Prevent near edge problems (see #1389)
-      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1
-      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi
-      ELSE                           ;   ii = mi1(ii)
-      ENDIF
-      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1
-      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj
-      ELSE                           ;   ij  = mj1(ij)
-      ENDIF
-      !
-      IF( ii == jpi )   ii = ii-1      
-      IF( ij == jpj )   ij = ij-1
+      ! find position in this processor.
+      ii = mi1(ii) ; ij  = mj1(ij)
       !
       icb_utl_bilin = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   &
@@ -333,16 +200,6 @@
       zj = pj - REAL(ij,wp)
       !
-      ! find position in this processor. Prevent near edge problems (see #1389)
-      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1
-      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi
-      ELSE                           ;   ii = mi1(ii)
-      ENDIF
-      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1
-      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj
-      ELSE                           ;   ij  = mj1(ij)
-      ENDIF
-      !
-      IF( ii == jpi )   ii = ii-1      
-      IF( ij == jpj )   ij = ij-1
+      ! find position in this processor.
+      ii = mi1(ii) ; ij  = mj1(ij)
       !
       z4(1) = pfld(ii  ,ij  )
@@ -392,17 +249,7 @@
       zi = pi - REAL(ii,wp)          !!gm use here mig, mjg arrays
       zj = pj - REAL(ij,wp)
-
-      ! find position in this processor. Prevent near edge problems (see #1389)
-      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1
-      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi
-      ELSE                           ;   ii = mi1(ii)
-      ENDIF
-      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1
-      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj
-      ELSE                           ;   ij  = mj1(ij)
-      ENDIF
-      !
-      IF( ii == jpi )   ii = ii-1      
-      IF( ij == jpj )   ij = ij-1
+      !
+      ! find position in this processor.
+      ii = mi1(ii) ; ij  = mj1(ij)
       !
       IF(    0.0_wp <= zi .AND. zi < 0.5_wp   ) THEN
