Changeset 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/ICB/icbutl.F90

Timestamp:

2021-05-05T13:18:04+02:00 (3 years ago)

Author:

mcastril

Message:

[2021/HPC-11_mcastril_HPDAonline_DiagGPU] Update externals

Location:

NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU

Files:

• . (modified) (1 prop)
• src/OCE/ICB/icbutl.F90 (modified) (18 diffs)

NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU

@@ -3 +3 @@
 ^/utils/build/mk@HEAD         mk
 ^/utils/tools@HEAD            tools
-^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF
+^/vendors/AGRIF/dev@HEAD      ext/AGRIF
 ^/vendors/FCM@HEAD            ext/FCM
 ^/vendors/IOIPSL@HEAD         ext/IOIPSL
+^/vendors/PPR@HEAD            ext/PPR
 
 # SETTE
-^/utils/CI/sette@13559        sette
+^/utils/CI/sette@14244        sette

@@ -3 +3 @@
 ^/utils/build/mk@HEAD         mk
 ^/utils/tools@HEAD            tools
-^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF
+^/vendors/AGRIF/dev@HEAD      ext/AGRIF
 ^/vendors/FCM@HEAD            ext/FCM
 ^/vendors/IOIPSL@HEAD         ext/IOIPSL
+^/vendors/PPR@HEAD            ext/PPR
 
 # SETTE
-^/utils/CI/sette@13559        sette
+^/utils/CI/sette@14244        sette

NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/ICB/icbutl.F90

@@ -8 +8 @@
    !!            -    !                            Removal of mapping from another grid
    !!            -    !  2011-04  (Alderson)       Split into separate modules
+   !!           4.2   !  2020-07  (P. Mathiot)     simplification of interpolation routine
+   !!                 !                            and add Nacho Merino work
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   icb_utl_interp   :
-   !!   icb_utl_bilin    :
-   !!   icb_utl_bilin_e  :
+   !!   icb_utl_pos      : compute bottom left corner indice, weight and mask
+   !!   icb_utl_bilin_h  : interpolation field to icb position
+   !!   icb_utl_bilin_e  : interpolation of scale factor to icb position
    !!----------------------------------------------------------------------
    USE par_oce                             ! ocean parameters
+   USE oce,    ONLY: ts, uu, vv
    USE dom_oce                             ! ocean domain
    USE in_out_manager                      ! IO parameters
@@ -31 +35 @@
    PRIVATE
 
+   INTERFACE icb_utl_bilin_h
+      MODULE PROCEDURE icb_utl_bilin_2d_h, icb_utl_bilin_3d_h
+   END INTERFACE
+
    PUBLIC   icb_utl_copy          ! routine called in icbstp module
+   PUBLIC   icb_utl_getkb         ! routine called in icbdyn and icbthm modules
+   PUBLIC   test_icb_utl_getkb    ! routine called in icbdyn and icbthm modules
+   PUBLIC   icb_utl_zavg          ! routine called in icbdyn and icbthm modules
    PUBLIC   icb_utl_interp        ! routine called in icbdyn, icbthm modules
-   PUBLIC   icb_utl_bilin         ! routine called in icbini, icbdyn modules
-   PUBLIC   icb_utl_bilin_x       ! routine called in icbdyn module
+   PUBLIC   icb_utl_bilin_h       ! routine called in icbdyn module
    PUBLIC   icb_utl_add           ! routine called in icbini.F90, icbclv, icblbc and icbrst modules
    PUBLIC   icb_utl_delete        ! routine called in icblbc, icbthm modules
@@ -47 +57 @@
    PUBLIC   icb_utl_heat          ! routine called in icbdia module
 
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -54 +66 @@
 CONTAINS
 
-   SUBROUTINE icb_utl_copy()
+   SUBROUTINE icb_utl_copy( Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_utl_copy  ***
@@ -62 +74 @@
       !! ** Method  : - blah blah
       !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1) :: ztmp
 #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
+      INTEGER :: jk   ! vertical loop index
+      INTEGER :: Kmm  ! ocean time levelindex
+      !
       ! 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(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1)
-      vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1)
+      !
+      ! surface forcing
+      !
+      ssu_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1)
+      ssv_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1)
+      sst_e(1:jpi,1:jpj) = sst_m(:,:)
+      sss_e(1:jpi,1:jpj) = sss_m(:,:)
+      fr_e (1:jpi,1:jpj) = fr_i (:,:)
+      ua_e (1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk
+      va_e (1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk
       ff_e(1:jpi,1:jpj) = ff_f (:,:) 
-      tt_e(1:jpi,1:jpj) = sst_m(:,:)
-      ss_e(1:jpi,1:jpj) = sss_m(:,:)
-      fr_e(1:jpi,1:jpj) = fr_i (:,:)
-      ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk
-      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 )
-      CALL lbc_lnk_icb( 'icbutl', ss_e, 'T', +1._wp, 1, 1 )
+      !
+      CALL lbc_lnk_icb( 'icbutl', ssu_e, 'U', -1._wp, 1, 1 )
+      CALL lbc_lnk_icb( 'icbutl', ssv_e, 'V', -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 )
 #if defined key_si3
       hi_e(1:jpi, 1:jpj) = hm_i (:,:)  
@@ -96 +110 @@
       ssh_e(1:jpi, 1:jpj) = zssh_lead_m(:,:) * tmask(:,:,1)
       !
-      CALL lbc_lnk_icb( 'icbutl', hi_e , '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(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1)
+      ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1)         
 #endif
-      CALL lbc_lnk_icb( 'icbutl', ssh_e, 'T', +1._wp, 1, 1 )
+      !
+      ! (PM) could be improve with a 3d lbclnk gathering both variables
+      ! should be done once extra haloe generalised
+      IF ( ln_M2016 ) THEN
+         DO jk = 1,jpk
+            ! uoce
+            ztmp(1:jpi,1:jpj) = uu(:,:,jk,Kmm)
+            CALL lbc_lnk_icb( 'icbutl', ztmp, 'U', -1._wp, 1, 1 )
+            uoce_e(:,:,jk) = ztmp(:,:)
+            !
+            ! voce
+            ztmp(1:jpi,1:jpj) = vv(:,:,jk,Kmm)
+            CALL lbc_lnk_icb( 'icbutl', ztmp, 'V', -1._wp, 1, 1 )
+            voce_e(:,:,jk) = ztmp(:,:)
+            !
+            e3t_e(1:jpi,1:jpj,jk) = e3t(:,:,jk,Kmm)
+         END DO
+         toce_e(1:jpi,1:jpj,:) = ts(:,:,:,1,Kmm)
+      END IF
       !
    END SUBROUTINE icb_utl_copy
 
 
-   SUBROUTINE icb_utl_interp( pi, pe1, puo, pui, pua, pssh_i,   &
-      &                       pj, pe2, pvo, pvi, pva, pssh_j,   &
-      &                       psst, pcn, phi, pff, psss        )
+   SUBROUTINE icb_utl_interp( pi, pj, pe1 , pssu, pui, pua, pssh_i,         &
+      &                               pe2 , pssv, pvi, pva, pssh_j,         &
+      &                               psst, psss, pcn, phi, pff   ,         &
+      &                               plon, plat, ptoce, puoce, pvoce, pe3t )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_utl_interp  ***
@@ -127 +159 @@
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ) ::   pi , pj                        ! position in (i,j) referential
-      REAL(wp), INTENT(  out) ::   pe1, pe2                       ! i- and j scale factors
-      REAL(wp), INTENT(  out) ::   puo, pvo, pui, pvi, pua, pva   ! ocean, ice and wind speeds
-      REAL(wp), INTENT(  out) ::   pssh_i, pssh_j                 ! ssh i- & j-gradients
-      REAL(wp), INTENT(  out) ::   psst, pcn, phi, pff, psss      ! SST, ice concentration, ice thickness, Coriolis, SSS
-      !
+      REAL(wp), INTENT(  out), OPTIONAL ::   pe1, pe2                       ! i- and j scale factors
+      REAL(wp), INTENT(  out), OPTIONAL ::   pssu, pssv, pui, pvi, pua, pva ! ocean, ice and wind speeds
+      REAL(wp), INTENT(  out), OPTIONAL ::   pssh_i, pssh_j                 ! ssh i- & j-gradients
+      REAL(wp), INTENT(  out), OPTIONAL ::   psst, psss, pcn, phi, pff      ! SST, SSS, ice concentration, ice thickness, Coriolis
+      REAL(wp), INTENT(  out), OPTIONAL ::   plat, plon                     ! position
+      REAL(wp), DIMENSION(jpk), INTENT(  out), OPTIONAL ::   ptoce, puoce, pvoce, pe3t   ! 3D variables
+      !
+      REAL(wp), DIMENSION(4) :: zwT  , zwU  , zwV  , zwF   ! interpolation weight
+      REAL(wp), DIMENSION(4) :: zmskF, zmskU, zmskV, zmskT ! mask
+      REAL(wp), DIMENSION(4) :: zwTp, zmskTp, zwTm, zmskTm
+      REAL(wp), DIMENSION(4,jpk) :: zw1d
+      INTEGER                :: iiT, iiU, iiV, iiF, ijT, ijU, ijV, ijF ! bottom left corner
+      INTEGER                :: iiTp, iiTm, ijTp, ijTm
       REAL(wp) ::   zcd, zmod       ! local scalars
       !!----------------------------------------------------------------------
-
-      pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )      ! scale factors
-      pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
-      !
-      puo  = icb_utl_bilin_h( uo_e, pi, pj, 'U', .false.  )    ! ocean velocities
-      pvo  = icb_utl_bilin_h( vo_e, pi, pj, 'V', .false.  )
-      psst = icb_utl_bilin_h( tt_e, pi, pj, 'T', .true.   )    ! SST
-      psss = icb_utl_bilin_h( ss_e, pi, pj, 'T', .true.   )    ! SSS
-      pcn  = icb_utl_bilin_h( fr_e, pi, pj, 'T', .true.   )    ! ice concentration
-      pff  = icb_utl_bilin_h( ff_e, pi, pj, 'F', .false.  )    ! Coriolis parameter
-      !
-      pua  = icb_utl_bilin_h( ua_e, pi, pj, 'U', .true.   )    ! 10m wind
-      pva  = icb_utl_bilin_h( va_e, pi, pj, 'V', .true.   )    ! 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)  )  )
-      pua  = pua * zmod                                       ! note: stress module=0 necessarly implies ua=va=0
-      pva  = pva * zmod
-
+      !
+      ! get position, weight and mask 
+      CALL icb_utl_pos( pi, pj, 'T', iiT, ijT, zwT, zmskT )
+      CALL icb_utl_pos( pi, pj, 'U', iiU, ijU, zwU, zmskU )
+      CALL icb_utl_pos( pi, pj, 'V', iiV, ijV, zwV, zmskV )
+      CALL icb_utl_pos( pi, pj, 'F', iiF, ijF, zwF, zmskF )
+      !
+      ! metrics and coordinates
+      IF ( PRESENT(pe1 ) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )      ! scale factors
+      IF ( PRESENT(pe2 ) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
+      IF ( PRESENT(plon) ) plon= icb_utl_bilin_h( rlon_e, iiT, ijT, zwT, .true.  )
+      IF ( PRESENT(plat) ) plat= icb_utl_bilin_h( rlat_e, iiT, ijT, zwT, .false. )
+      !
+      IF ( PRESENT(pssu) ) pssu = icb_utl_bilin_h( ssu_e, iiU, ijU, zwU        , .false. ) ! ocean velocities
+      IF ( PRESENT(pssv) ) pssv = icb_utl_bilin_h( ssv_e, iiV, ijV, zwV        , .false. ) !
+      IF ( PRESENT(psst) ) psst = icb_utl_bilin_h( sst_e, iiT, ijT, zwT * zmskT, .false. ) ! sst
+      IF ( PRESENT(psss) ) psss = icb_utl_bilin_h( sss_e, iiT, ijT, zwT * zmskT, .false. ) ! sss
+      IF ( PRESENT(pcn ) ) pcn  = icb_utl_bilin_h( fr_e , iiT, ijT, zwT * zmskT, .false. ) ! ice concentration
+      IF ( PRESENT(pff ) ) pff  = icb_utl_bilin_h( ff_e , iiF, ijF, zwF        , .false. ) ! Coriolis parameter
+      !
+      IF ( PRESENT(pua) .AND. PRESENT(pva) ) THEN
+         pua  = icb_utl_bilin_h( ua_e, iiU, ijU, zwU * zmskU, .false. ) ! 10m wind
+         pva  = icb_utl_bilin_h( va_e, iiV, ijV, zwV * zmskV, .false. ) ! 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)  )  )
+         pua  = pua * zmod                                       ! note: stress module=0 necessarly implies ua=va=0
+         pva  = pva * zmod
+      END IF
+      !
 #if defined key_si3
-      pui = icb_utl_bilin_h( ui_e , pi, pj, 'U', .false. )    ! sea-ice velocities
-      pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V', .false. )
-      phi = icb_utl_bilin_h( hi_e , pi, pj, 'T', .true.  )    ! ice thickness
+      IF ( PRESENT(pui) ) pui = icb_utl_bilin_h( ui_e , iiU, ijU, zwU        , .false. ) ! sea-ice velocities
+      IF ( PRESENT(pvi) ) pvi = icb_utl_bilin_h( vi_e , iiV, ijV, zwV        , .false. )
+      IF ( PRESENT(phi) ) phi = icb_utl_bilin_h( hi_e , iiT, ijT, zwT * zmskT, .false. ) ! ice thickness
 #else
-      pui = 0._wp
-      pvi = 0._wp
-      phi = 0._wp
+      IF ( PRESENT(pui) ) pui = 0._wp
+      IF ( PRESENT(pvi) ) pvi = 0._wp
+      IF ( PRESENT(phi) ) phi = 0._wp
 #endif
-
+      !
       ! Estimate SSH gradient in i- and j-direction (centred evaluation)
-      pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T', .true. ) -   &
-         &       icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T', .true. )  ) / ( 0.2_wp * pe1 )
-      pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T', .true. ) -   &
-         &       icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T', .true. )  ) / ( 0.2_wp * pe2 )
+      IF ( PRESENT(pssh_i) .AND. PRESENT(pssh_j) ) THEN
+         CALL icb_utl_pos( pi+0.1, pj    , 'T', iiTp, ijTp, zwTp, zmskTp )
+         CALL icb_utl_pos( pi-0.1, pj    , 'T', iiTm, ijTm, zwTm, zmskTm )
+         !
+         IF ( .NOT. PRESENT(pe1) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )
+         pssh_i = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) -   &
+            &       icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. )  ) / ( 0.2_wp * pe1 )
+         !
+         CALL icb_utl_pos( pi    , pj+0.1, 'T', iiTp, ijTp, zwTp, zmskTp )
+         CALL icb_utl_pos( pi    , pj-0.1, 'T', iiTm, ijTm, zwTm, zmskTm )
+         !
+         IF ( .NOT. PRESENT(pe2) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
+         pssh_j = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) -   &
+            &       icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. )  ) / ( 0.2_wp * pe2 )
+      END IF
+      !
+      ! 3d interpolation
+      IF ( PRESENT(puoce) .AND. PRESENT(pvoce) ) THEN
+         ! no need to mask as 0 is a valid data for land
+         zw1d(1,:) = zwU(1) ; zw1d(2,:) = zwU(2) ; zw1d(3,:) = zwU(3) ; zw1d(4,:) = zwU(4) ;
+         puoce(:) = icb_utl_bilin_h( uoce_e , iiU, ijU, zw1d )
+
+         zw1d(1,:) = zwV(1) ; zw1d(2,:) = zwV(2) ; zw1d(3,:) = zwV(3) ; zw1d(4,:) = zwV(4) ;
+         pvoce(:) = icb_utl_bilin_h( voce_e , iiV, ijV, zw1d )
+      END IF
+
+      IF ( PRESENT(ptoce) ) THEN
+         ! for temperature we need to mask the weight properly
+         ! no need of extra halo as it is a T point variable
+         zw1d(1,:) = tmask(iiT  ,ijT  ,:) * zwT(1) * zmskT(1)
+         zw1d(2,:) = tmask(iiT+1,ijT  ,:) * zwT(2) * zmskT(2)
+         zw1d(3,:) = tmask(iiT  ,ijT+1,:) * zwT(3) * zmskT(3)
+         zw1d(4,:) = tmask(iiT+1,ijT+1,:) * zwT(4) * zmskT(4)
+         ptoce(:) = icb_utl_bilin_h( toce_e , iiT, ijT, zw1d )
+      END IF
+      !
+      IF ( PRESENT(pe3t)  ) pe3t(:)  = e3t_e(iiT,ijT,:)    ! as in Nacho tarball need to be fix once we are able to reproduce Nacho results
       !
    END SUBROUTINE icb_utl_interp
 
-
-   REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type, plmask )
+   SUBROUTINE icb_utl_pos( pi, pj, cd_type, kii, kij, pw, pmsk )
       !!----------------------------------------------------------------------
       !!                  ***  FUNCTION icb_utl_bilin  ***
@@ -182 +265 @@
       !!
       !!----------------------------------------------------------------------
-      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
-      LOGICAL                             , INTENT(in) ::   plmask    ! special treatment of mask point
-      !
-      INTEGER  ::   ii, ij   ! local integer
-      REAL(wp) ::   zi, zj   ! local real
-      REAL(wp) :: zw1, zw2, zw3, zw4
-      REAL(wp), DIMENSION(4) :: zmask
+      REAL(wp)              , INTENT(IN)  ::   pi, pj    ! targeted coordinates in (i,j) referential
+      CHARACTER(len=1)      , INTENT(IN)  ::   cd_type   ! point type
+      REAL(wp), DIMENSION(4), INTENT(OUT) ::   pw, pmsk  ! weight and mask
+      INTEGER ,               INTENT(OUT) ::   kii, kij  ! bottom left corner position in local domain
+      !
+      REAL(wp) :: zwi, zwj ! distance to bottom left corner
+      INTEGER  :: ierr 
+      !
       !!----------------------------------------------------------------------
       !
@@ -198 +280 @@
          ! since we're looking for four T points containing quadrant we're in of 
          ! current T cell
-         ii = MAX(0, INT( pi     ))
-         ij = MAX(0, INT( pj     ))    ! T-point
-         zi = pi - REAL(ii,wp)
-         zj = pj - REAL(ij,wp)
+         kii = MAX(0, INT( pi        ))
+         kij = MAX(0, INT( pj        ))    ! T-point
+         zwi = pi - REAL(kii,wp)
+         zwj = pj - REAL(kij,wp)
       CASE ( 'U' )
-         ii = MAX(0, INT( pi-0.5_wp ))
-         ij = MAX(0, INT( pj     ))    ! U-point
-         zi = pi - 0.5_wp - REAL(ii,wp)
-         zj = pj - REAL(ij,wp)
+         kii = MAX(0, INT( pi-0.5_wp ))
+         kij = MAX(0, INT( pj        ))    ! U-point
+         zwi = pi - 0.5_wp - REAL(kii,wp)
+         zwj = pj - REAL(kij,wp)
       CASE ( 'V' )
-         ii = MAX(0, INT( pi     ))
-         ij = MAX(0, INT( pj-0.5_wp ))    ! V-point
-         zi = pi - REAL(ii,wp)
-         zj = pj - 0.5_wp - REAL(ij,wp)
+         kii = MAX(0, INT( pi        ))
+         kij = MAX(0, INT( pj-0.5_wp ))    ! V-point
+         zwi = pi - REAL(kii,wp)
+         zwj = pj - 0.5_wp - REAL(kij,wp)
       CASE ( 'F' )
-         ii = MAX(0, INT( pi-0.5_wp ))
-         ij = MAX(0, INT( pj-0.5_wp ))    ! F-point
-         zi = pi - 0.5_wp - REAL(ii,wp)
-         zj = pj - 0.5_wp - REAL(ij,wp)
+         kii = MAX(0, INT( pi-0.5_wp ))
+         kij = MAX(0, INT( pj-0.5_wp ))    ! F-point
+         zwi = pi - 0.5_wp - REAL(kii,wp)
+         zwj = pj - 0.5_wp - REAL(kij,wp)
       END SELECT
+      !
+      ! compute weight
+      pw(1) = (1._wp-zwi) * (1._wp-zwj)
+      pw(2) =        zwi  * (1._wp-zwj)
+      pw(3) = (1._wp-zwi) *        zwj
+      pw(4) =        zwi  *        zwj
+      !
+      ! find position in this processor. Prevent near edge problems (see #1389)
+      !
+      IF (TRIM(cd_type) == 'T' ) THEN
+         ierr = 0
+         IF    ( kii <  mig( 1 ) ) THEN   ;  ierr = ierr + 1
+         ELSEIF( kii >= mig(jpi) ) THEN   ;  ierr = ierr + 1
+         ENDIF
+         !
+         IF    ( kij <  mjg( 1 ) ) THEN   ;   ierr = ierr + 1
+         ELSEIF( kij >= mjg(jpj) ) THEN   ;   ierr = ierr + 1
+         ENDIF
+         !
+         IF ( ierr > 0 ) THEN
+            WRITE(numicb,*) 'bottom left corner T point out of bound'
+            WRITE(numicb,*) pi, kii, mig( 1 ), mig(jpi)
+            WRITE(numicb,*) pj, kij, mjg( 1 ), mjg(jpj)
+            WRITE(numicb,*) pmsk
+            CALL FLUSH(numicb)
+            CALL ctl_stop('STOP','icb_utl_bilin_e: an icebergs coordinates is out of valid range (out of bound error).'       , &
+                 &                                'This can be fixed using rn_speed_limit=0.4 in &namberg.'                   , &
+                 &                                'More details in the corresponding iceberg.stat file (nn_verbose_level > 0).' )
+         END IF
+      END IF
       !
       ! find position in this processor. Prevent near edge problems (see #1389)
       ! (PM) will be useless if extra halo is used in NEMO
       !
-      IF    ( ii <= mig(1)-1 ) THEN   ;   ii = 0
-      ELSEIF( ii  > mig(jpi) ) THEN   ;   ii = jpi
-      ELSE                            ;   ii = mi1(ii)
+      IF    ( kii <= mig(1)-1 ) THEN   ;   kii = 0
+      ELSEIF( kii  > mig(jpi) ) THEN   ;   kii = jpi
+      ELSE                             ;   kii = mi1(kii)
       ENDIF
-      IF    ( ij <= mjg(1)-1 ) THEN   ;   ij = 0
-      ELSEIF( ij  > mjg(jpj) ) THEN   ;   ij = jpj
-      ELSE                            ;   ij = mj1(ij)
+      IF    ( kij <= mjg(1)-1 ) THEN   ;   kij = 0
+      ELSEIF( kij  > mjg(jpj) ) THEN   ;   kij = jpj
+      ELSE                             ;   kij = mj1(kij)
       ENDIF
       !
       ! define mask array 
-      IF (plmask) THEN
-         ! land value is not used in the interpolation
-         SELECT CASE ( cd_type )
-         CASE ( 'T' )
-            zmask = (/tmask_e(ii,ij), tmask_e(ii+1,ij), tmask_e(ii,ij+1), tmask_e(ii+1,ij+1)/)
-         CASE ( 'U' )
-            zmask = (/umask_e(ii,ij), umask_e(ii+1,ij), umask_e(ii,ij+1), umask_e(ii+1,ij+1)/)
-         CASE ( 'V' )
-            zmask = (/vmask_e(ii,ij), vmask_e(ii+1,ij), vmask_e(ii,ij+1), vmask_e(ii+1,ij+1)/)
-         CASE ( 'F' )
-            ! F case only used for coriolis, ff_f is not mask so zmask = 1
-            zmask = 1.
-         END SELECT
-      ELSE
-         ! land value is used during interpolation
-         zmask = 1.
-      END iF
-      !
-      ! compute weight
-      zw1 = zmask(1) * (1._wp-zi) * (1._wp-zj)
-      zw2 = zmask(2) *        zi  * (1._wp-zj)
-      zw3 = zmask(3) * (1._wp-zi) *        zj
-      zw4 = zmask(4) *        zi  *        zj
-      !
-      ! compute interpolated value
-      icb_utl_bilin_h = ( pfld(ii,ij)*zw1 + pfld(ii+1,ij)*zw2 + pfld(ii,ij+1)*zw3 + pfld(ii+1,ij+1)*zw4 ) / MAX(1.e-20, zw1+zw2+zw3+zw4) 
-      !
-   END FUNCTION icb_utl_bilin_h
-
-
-   REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type )
+      ! land value is not used in the interpolation
+      SELECT CASE ( cd_type )
+      CASE ( 'T' )
+         pmsk = (/tmask_e(kii,kij), tmask_e(kii+1,kij), tmask_e(kii,kij+1), tmask_e(kii+1,kij+1)/)
+      CASE ( 'U' )
+         pmsk = (/umask_e(kii,kij), umask_e(kii+1,kij), umask_e(kii,kij+1), umask_e(kii+1,kij+1)/)
+      CASE ( 'V' )
+         pmsk = (/vmask_e(kii,kij), vmask_e(kii+1,kij), vmask_e(kii,kij+1), vmask_e(kii+1,kij+1)/)
+      CASE ( 'F' )
+         ! F case only used for coriolis, ff_f is not mask so zmask = 1
+         pmsk = 1.
+      END SELECT
+   END SUBROUTINE icb_utl_pos
+
+   REAL(wp) FUNCTION icb_utl_bilin_2d_h( pfld, pii, pij, pw, pllon )
       !!----------------------------------------------------------------------
       !!                  ***  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
@@ -272 +369 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), 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)
+      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) ::   pfld      ! field to be interpolated
+      REAL(wp), DIMENSION(4)              , INTENT(in) ::   pw        ! weight
+      LOGICAL                             , INTENT(in) ::   pllon     ! input data is a longitude
+      INTEGER ,                             INTENT(in) ::   pii, pij  ! bottom left corner
+      !
+      REAL(wp), DIMENSION(4) :: zdat ! input data
+      !!----------------------------------------------------------------------
+      !
+      ! data
+      zdat(1) = pfld(pii  ,pij  )
+      zdat(2) = pfld(pii+1,pij  )
+      zdat(3) = pfld(pii  ,pij+1)
+      zdat(4) = pfld(pii+1,pij+1)
+      !
+      IF( pllon .AND. MAXVAL(zdat) - MINVAL(zdat) > 90._wp ) THEN
+         WHERE( zdat < 0._wp ) zdat = zdat + 360._wp
       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 = ( 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
-
-
-   REAL(wp) FUNCTION icb_utl_bilin_x( pfld, pi, pj )
-      !!----------------------------------------------------------------------
-      !!                  ***  FUNCTION icb_utl_bilin_x  ***
+      !
+      ! compute interpolated value
+      icb_utl_bilin_2d_h = ( zdat(1)*pw(1) + zdat(2)*pw(2) + zdat(3)*pw(3) + zdat(4)*pw(4) ) / MAX(1.e-20, pw(1)+pw(2)+pw(3)+pw(4)) 
+      !
+      IF( pllon .AND. icb_utl_bilin_2d_h > 180._wp ) icb_utl_bilin_2d_h = icb_utl_bilin_2d_h - 360._wp
+      !
+   END FUNCTION icb_utl_bilin_2d_h
+
+   FUNCTION icb_utl_bilin_3d_h( pfld, pii, pij, pw )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION icb_utl_bilin  ***
       !!
       !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type
-      !!                Special case for interpolating longitude
+      !!                this version deals with extra halo points
       !!
       !!       !!gm  CAUTION an optional argument should be added to handle
@@ -336 +405 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pfld      ! field to be interpolated
-      REAL(wp)                    , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential
-      !
-      INTEGER                                  ::   ii, ij   ! local integer
-      REAL(wp)                                 ::   zi, zj   ! local real
-      REAL(wp)                                 ::   zret     ! local real
-      REAL(wp), DIMENSION(4)                   ::   z4
-      !!----------------------------------------------------------------------
-      !
-      ! 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)
-      !
-      ! 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
-      !
-      z4(1) = pfld(ii  ,ij  )
-      z4(2) = pfld(ii+1,ij  )
-      z4(3) = pfld(ii  ,ij+1)
-      z4(4) = pfld(ii+1,ij+1)
-      IF( MAXVAL(z4) - MINVAL(z4) > 90._wp ) THEN
-         WHERE( z4 < 0._wp ) z4 = z4 + 360._wp
-      ENDIF
-      !
-      zret = (z4(1) * (1.-zi) + z4(2) * zi) * (1.-zj) + (z4(3) * (1.-zi) + z4(4) * zi) * zj
-      IF( zret > 180._wp ) zret = zret - 360._wp
-      icb_utl_bilin_x = zret
-      !
-   END FUNCTION icb_utl_bilin_x
-
+      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1, jpk), INTENT(in) ::   pfld      ! field to be interpolated
+      REAL(wp), DIMENSION(4,jpk)               , INTENT(in) ::   pw        ! weight
+      INTEGER ,                                  INTENT(in) ::   pii, pij  ! bottom left corner
+      REAL(wp), DIMENSION(jpk) :: icb_utl_bilin_3d_h
+      !
+      REAL(wp), DIMENSION(4,jpk) :: zdat ! input data
+      INTEGER :: jk
+      !!----------------------------------------------------------------------
+      !
+      ! data
+      zdat(1,:) = pfld(pii  ,pij  ,:)
+      zdat(2,:) = pfld(pii+1,pij  ,:)
+      zdat(3,:) = pfld(pii  ,pij+1,:)
+      zdat(4,:) = pfld(pii+1,pij+1,:)
+      !
+      ! compute interpolated value
+      DO jk=1,jpk
+         icb_utl_bilin_3d_h(jk) =   ( zdat(1,jk)*pw(1,jk) + zdat(2,jk)*pw(2,jk) + zdat(3,jk)*pw(3,jk) + zdat(4,jk)*pw(4,jk) ) &
+            &                     /   MAX(1.e-20, pw(1,jk)+pw(2,jk)+pw(3,jk)+pw(4,jk)) 
+      END DO
+      !
+   END FUNCTION icb_utl_bilin_3d_h
 
    REAL(wp) FUNCTION icb_utl_bilin_e( pet, peu, pev, pef, pi, pj )
@@ -390 +437 @@
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:), INTENT(in) ::   pet, peu, pev, pef   ! horizontal scale factor to be interpolated at t-,u-,v- & f-pts
-      REAL(wp)                , INTENT(in) ::   pi, pj               ! targeted coordinates in (i,j) referential
-      !
-      INTEGER  ::   ii, ij, icase, ierr   ! local integer
+      REAL(wp)                , INTENT(IN) ::   pi , pj              ! iceberg position
       !
       ! weights corresponding to corner points of a T cell quadrant
       REAL(wp) ::   zi, zj          ! local real
+      INTEGER  ::   ii, ij          ! bottom left corner coordinate in local domain
       !
       ! values at corner points of a T cell quadrant
@@ -402 +448 @@
       !!----------------------------------------------------------------------
       !
+      ! cannot used iiT because need ii/ij reltaive to global indices not local one
       ii = MAX(1, INT( pi ))   ;   ij = MAX(1, INT( pj ))            ! left bottom T-point (i,j) indices
-
+      ! 
       ! fractional box spacing
       ! 0   <= zi < 0.5, 0   <= zj < 0.5   -->  NW quadrant of current T cell
@@ -413 +460 @@
       zj = pj - REAL(ij,wp)
 
-      ! find position in this processor. Prevent near edge problems (see #1389)
-      !
-      ierr = 0
-      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1       ; ierr = ierr + 1
-      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi     ; ierr = ierr + 1
-      ELSE                           ;   ii = mi1(ii)
-      ENDIF
-      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1       ; ierr = ierr + 1
-      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj     ; ierr = ierr + 1
-      ELSE                           ;   ij  = mj1(ij)
-      ENDIF
-      !
-      IF( ii == jpi ) THEN ; ii = ii-1 ; ierr = ierr + 1 ; END IF     
-      IF( ij == jpj ) THEN ; ij = ij-1 ; ierr = ierr + 1 ; END IF
-      !
-      IF ( ierr > 0 ) CALL ctl_stop('STOP','icb_utl_bilin_e: an icebergs coordinates is out of valid range (out of bound error)')
+      ! conversion to local domain (no need to do a sanity check already done in icbpos)
+      ii = mi1(ii)
+      ij = mj1(ij)
       !
       IF(    0.0_wp <= zi .AND. zi < 0.5_wp   ) THEN
@@ -465 +499 @@
    END FUNCTION icb_utl_bilin_e
 
+   SUBROUTINE icb_utl_getkb( kb, pe3, pD )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE icb_utl_getkb         ***
+      !!
+      !! ** Purpose :   compute the latest level affected by icb
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER,                INTENT(out):: kb
+      REAL(wp), DIMENSION(:), INTENT(in) :: pe3
+      REAL(wp),               INTENT(in) :: pD
+      !!
+      INTEGER  :: jk
+      REAL(wp) :: zdepw
+      !!----------------------------------------------------------------------
+      !!
+      zdepw = pe3(1) ; kb = 2
+      DO WHILE ( zdepw <  pD)
+         zdepw = zdepw + pe3(kb)
+         kb = kb + 1
+      END DO
+      kb = MIN(kb - 1,jpk)
+   END SUBROUTINE
+
+   SUBROUTINE icb_utl_zavg(pzavg, pdat, pe3, pD, kb )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE icb_utl_getkb         ***
+      !!
+      !! ** Purpose :   compute the vertical average of ocean properties affected by icb
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER,                INTENT(in ) :: kb        ! deepest level affected by icb
+      REAL(wp), DIMENSION(:), INTENT(in ) :: pe3, pdat ! vertical profile
+      REAL(wp),               INTENT(in ) :: pD        ! draft
+      REAL(wp),               INTENT(out) :: pzavg     ! z average
+      !!----------------------------------------------------------------------
+      INTEGER  :: jk
+      REAL(wp) :: zdep
+      !!----------------------------------------------------------------------
+      pzavg = 0.0 ; zdep = 0.0
+      DO jk = 1,kb-1
+         pzavg = pzavg + pe3(jk)*pdat(jk)
+         zdep  = zdep  + pe3(jk)
+      END DO
+      ! if kb is limited by mbkt  => bottom value is used between bathy and icb tail
+      ! if kb not limited by mbkt => ocean value over mask is used (ie 0.0 for u, v)
+      pzavg = ( pzavg + (pD - zdep)*pdat(kb)) / pD
+   END SUBROUTINE
 
    SUBROUTINE icb_utl_add( bergvals, ptvals )
@@ -653 +734 @@
       WRITE(numicb, 9200) kt, berg%number(1), &
                    pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel,  &
-                   pt%uo, pt%vo, pt%ua, pt%va, pt%ui, pt%vi
+                   pt%ssu, pt%ssv, pt%ua, pt%va, pt%ui, pt%vi
       CALL flush( numicb )
  9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4))
@@ -679 +760 @@
          WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea
          WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' )   &
-            &         'timestep', 'number', 'xi,yj','lon,lat','u,v','uo,vo','ua,va','ui,vi'
+            &         'timestep', 'number', 'xi,yj','lon,lat','u,v','ssu,ssv','ua,va','ui,vi'
       ENDIF
       DO WHILE( ASSOCIATED(this) )
@@ -823 +904 @@
    END FUNCTION icb_utl_heat
 
+   SUBROUTINE test_icb_utl_getkb
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION test_icb_utl_getkb  ***
+      !!
+      !! ** Purpose : Test routine icb_utl_getkb, icb_utl_zavg
+      !! ** Methode : Call each subroutine with specific input data
+      !!              What should be the output is easy to determined and check 
+      !!              if NEMO return the correct answer.
+      !! ** Comments : not called, if needed a CALL test_icb_utl_getkb need to be added in icb_step
+      !!----------------------------------------------------------------------
+      INTEGER :: ikb
+      REAL(wp) :: zD, zout
+      REAL(wp), DIMENSION(jpk) :: ze3, zin
+      WRITE(numout,*) 'Test icb_utl_getkb : '
+      zD = 0.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      zD = 8000000.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      zD = 80.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      zD = 85.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      zD = 75.0 ; ze3= 20.0
+      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      WRITE(numout,*) 'OUTPUT : kb = ',ikb
+
+      WRITE(numout,*) '=================================='
+      WRITE(numout,*) 'Test icb_utl_zavg'
+      zD = 0.0 ; ze3= 20.0 ; zin=1.0
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
+      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
+      WRITE(numout,*) 'OUTPUT : zout = ',zout
+
+      zD = 50.0 ; ze3= 20.0 ; zin=1.0; zin(3:jpk) = 0.0
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
+      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
+      WRITE(numout,*) 'OUTPUT : zout = ',zout
+      CALL FLUSH(numout)
+
+      zD = 80.0 ; ze3= 20.0 ; zin=1.0; zin(3:jpk) = 0.0
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
+      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
+      WRITE(numout,*) 'OUTPUT : zout = ',zout
+
+      zD = 80 ; ze3= 20.0 ; zin=1.0 ; zin(3:jpk) = 0.0
+      CALL icb_utl_getkb(ikb, ze3, zD)
+      ikb = 2
+      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
+      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
+      WRITE(numout,*) 'OUTPUT : zout = ',zout
+
+      CALL FLUSH(numout)
+
+   END SUBROUTINE test_icb_utl_getkb
+
    !!======================================================================
 END MODULE icbutl

@@ -3 +3 @@
 ^/utils/build/mk@HEAD         mk
 ^/utils/tools@HEAD            tools
-^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF
+^/vendors/AGRIF/dev@HEAD      ext/AGRIF
 ^/vendors/FCM@HEAD            ext/FCM
 ^/vendors/IOIPSL@HEAD         ext/IOIPSL
+^/vendors/PPR@HEAD            ext/PPR
 
 # SETTE
-^/utils/CI/sette@13559        sette
+^/utils/CI/sette@14244        sette